On characterization and application of oscillatory almost-cyclostationary processes

The class of the second-order oscillatory almost-cyclostationary processes is characterized. These processes have autocorrelation function which is the superposition of amplitude- and angle-modulated complex sinewaves, where the modulating functions, referred to as evolutionary cyclic autocorrelation functions, depend on both time and lag parameter. This class of processes includes that of the almost-cyclostationary processes. The problem of statistical function measurements is addressed for the special case of amplitude-modulated time-warped almost-cyclostationary processes. These processes are shown to be a suitable model for the electrocardiogram.


I. INTRODUCTION
The cyclostationary and almost-cyclostationary (ACS) models have been suitably exploited for describing many natural and man-made processes [1], [3], [7], [11], [17]. For secondorder wide-sense ACS processes, the autocorrelation is an almost-periodic function of time. Its (generalized) Fourier series expansion has frequencies, called cycle frequencies, belonging to a countable set of possibly incommensurate values and coefficients, called cyclic autocorrelation functions, depending only on the lag parameter. Spectral components of ACS processes are correlated when their separation is equal to one of the cycle frequencies. This corresponds to a Loève bifrequency spectrum [8] whose support is constituted by lines with unit slope in the bifrequency plane.
By generalizing Priestley's evolutionary spectral analysis [16], the more general model of the oscillatory almostcyclostationary (OACS) processes has been recently proposed in [10,Sec. 6]. For OACS processes, the autocorrelation function is given by the superposition of amplitude-and angle-modulated sinewaves whose frequencies are called cycle frequencies. The modulating functions, referred to as evolutionary cyclic autocorrelation functions, depend on both time and lag parameter. In the special case where the evolutionary cyclic autocorrelation functions do not depend on time, there is no amplitude and/or angle modulation and the OACS process reduce to an ACS one. In such a case, the evolutionary cyclic autocorrelation functions are coincident with the cyclic autocorrelation functions.
OACS processes turn out to be a suitable model for describing phenomena where the timing is irregular. Examples are several biological signals and, in communications, the received signal when there is relative motion between transmitter and receiver with generic motion law.
In this paper, an OACS process is characterized by an underlying ACS process and a linear time-variant (LTV) system that transforms this ACS process into the OACS one. Then, an OACS process is characterized in the frequency domain by its Loève bifrequency spectrum. It is shown that it can be seen as obtained by that of the underlying ACS process by spreading the spectral correlation densities around the support lines.
As special case of OACS processes, the amplitudemodulated (AM) time-warped (TW) ACS processes are considered and their characterization is provided. Results for timewarped ACS processes [12] with no amplitude modulation are extended here to estimate the time-warping and amplitudemodulation functions of AM-TW ACS processes and for recovering the underlying ACS process.
As an example of application, the AM-TW ACS model is proposed for the electrocardiogram (ECG). In the model, the underlying ACS signal is given by the superposition of a periodic signal and a zero-mean ACS signal. The timewarping and amplitude-modulation functions are linked to the variability with time of the heart rate and the variation in propagation of the electric wave throughout the heart. Numerical results are presented to corroborate the proposed model.
The paper is organized as follows. In Section II, widesense stationary and oscillatory processes are briefly reviewed to introduce notation. The class of the OACS processes is characterized is Section III. The special case of AM-TW ACS processes is treated in Section IV and the recovery of the underlying ACS process addressed in Section V. Numerical results on the ECG signal are reported in Section VI and conclusions are drawn in Section VII.

II. WIDE-SENSE STATIONARY AND OSCILLATORY PROCESSES
Let x(t) be a continuous-time second-order harmonizable wide-sense stationary (WSS) process. It admits the Cramér representation [8] where the integrated complex spectrum Z(f ) is an orthogonalincrement process, that is, In (2), δ(·) is Dirac delta, ( * ) represents an optional complex conjugation, (−) is an optional minus sign linked to ( * ), and µ x (f 1 ) is a bounded measure on R. Subscript x denotes xx ( * ) . WSS processes have (conjugate) autocorrelation depending only on the lag parameter τ and constant average power In [16], the class of the oscillatory processes is introduced aimed at describing processes with (slowly) time-varying average power. They admit the representation where Z(f ) is an orthogonal-increment process (see (2)) and {A t (f )} is a family (labeled by f ) of low-pass functions of t (i.e., slowly varying with respect to t). Oscillatory processes have time-varying (conjugate) autocorrelation function and time-varying average power The function dσ xx * ,t (f ) |A t (f )| 2 dµ xx * (f ) is referred to as evolutionary spectrum.

III. OSCILLATORY ALMOST-CYCLOSTATIONARY PROCESSES
A process y(t) is said to be oscillatory ACS if it admits representation (5) with Z(f ) integrated complex spectrum of an ACS process x(t) [10,Sec. 6]. That is, x(t) has Cramér representation (1) with is a family of complex measures and A is the countable set of possibly incommensurate (conjugate) cycle frequencies of the almost-periodic (conjugate) autocorrelation function where R α x (τ ) are the (conjugate) cyclic autocorrelation functions of x(t).
If µ α x (f 1 ) does not contain the singular component, then cyclic spectrum of x(t) that possibly contains Dirac deltas in correspondence of the jumps of µ α x (f 1 ). The (conjugate) cyclic spectrum S α x (f 1 ) is the Fourier transform of the (conjugate) cyclic autocorrelation function R α x (τ ). Starting from (5) and (8), one obtains the (conjugate) autocorrelation function of the OACS processes where the functions are referred to as evolutionary (conjugate) cyclic autocorrelation functions. Subscript y denotes yy ( * ) . For ACS processes, the (conjugate) autocorrelation function is the superposition of complex sinewaves whose frequencies are the (conjugate) cycle frequencies and whose complex amplitudes are the (conjugate) cyclic autocorrelation functions (that depend only on the lag parameter τ ) (see (9)). For OACS processes, the (conjugate) autocorrelation function is the superposition of amplitude-and angle-modulated complex sinewaves (see (10)) whose frequencies are the (conjugate) cycle frequencies of the underlying ACS process (see (8)) and the modulating functions are the evolutionary (conjugate) cyclic autocorrelation functions.
Note that the OACS model is useful if the evolutionary (conjugate) cyclic autocorrelation functions ρ α y (t, τ ) defined in (11) are slowly varying with respect to e j2παt . That is, if where B(α, τ ) is the bandwidth of ρ α y (t, τ ) as function of t (with α and τ fixed).
The functions are named evolutionary (conjugate) cyclic spectra and the time-varying averaged power of y(t) is given by The Fourier transform of y(t) in (5) Thus, from (15) and (8), one obtains the Loève bifrequency spectrum for the OACS process y(t) The ACS processes are a special case of OACS processes and the Loève bifrequency spectrum of y(t) is coincident with that of x(t) (see (8)) Therefore, the Loève bifrequency spectrum (17) can be seen as obtained by the Loève bifrequency spectrum of the underlying ACS process x(t) by spreading around the support lines Every OACS process can be obtained by LTV filtering of an ACS process x(t) with Cramér representation (1). In fact, replacing (1) into where h(t, u) is the impulse-response function of the LTV system, one has that y(t) admits the representation (5) with If x(t) is ACS (that is, its integrated complex spectrum Z(f ) satisfies (8)), then y(t) is OACS.
The process x(t) is named the underlying ACS process of the OACS process y(t). The LTV system h(t, u) is named the evolution-inducing LTV system of the OACS process y(t). In fact, when h(t, u) is not almost-periodically time variant, it modifies the (conjugate) cyclic spectra of the underlying ACS process x(t) into the evolutionary (conjugate) cyclic spectra of the OACS process y(t).
Statistical function estimation for the general case of OACS processes is still an open problem. In the following sections, the special class of AM-TW ACS processes is characterized and a procedure for their statistical function measurements outlined.

IV. AMPLITUDE-MODULATED TIME-WARPED ACS PROCESSES
A LTV system that modifies time scale and amplitude of the input signal x(t) is refereed to as time-warping and amplitude-modulation transformation. In (21), h(t, u) = a(t) δ(u − ψ(t)), a(t) is the possibly complex-valued amplitude-modulation function, and ψ(t) is the real-valued time-warping function. Both a(t) and ψ(t) are deterministic. ψ(t) is assumed asymptotically unbounded and nondecreasing. Note that ψ(t) should not be confused with the timing jitter which, in general, is described as a discrete-time random process. If x(t) is ACS with (conjugate) autocorrelation (9), y(t) is referred to as AM-TW ACS. Its (conjugate) autocorrelation is The process y(t) can be modeled as OACS with modulating functions A evolutionary (conjugate) cyclic spectra and evolutionary (conjugate) cyclic autocorrelation functions (25) Let us assume that with ǫ(t) differentiable and slowly varying with respect to t (whose derivative is 1), that is, where . ǫ (t) is the first-order derivative of ǫ(t). In addition, let us assume that a(t + τ ) ≃ a(t) for all τ such that R α x (τ ) is significantly non zero.
In such a case, the OACS model characterized by (23), (24), and (25) reduces to a modulated cyclical model [10, Sec. 6.2.2], [15]. Expression (22) of the (conjugate) autocorrelation simplifies into The proof, not reported here for lack of space, is a generalization of that for Theorem 3.6 in [12].

V. RECOVERY OF THE UNDERLYING ACS PROCESS
Let us consider the AM-TW ACS process defined in (21) under the assumptions such that expression (28) for the autocorrelation holds and with a(t) real and positive. In this section, a procedure for the recovery of the underlying ACS process x(t) is outlined. For this purpose, first functions a(t) and ǫ(t) are estimated. Then, by de-warping and amplitudemodulation compensation, an estimate of the underlying ACS signal x(t) is obtained. This procedure generalizes to the case of a(t) not constant the procedure proposed in [12] for the case a(t) = 1.
Let α 0 be a (conjugate) cycle frequency of x(t) and let us assume that only its coarse estimate, say α 0 is available, for example, by locating the peaks of the power spectral density (PSD) of y(t + τ ) y ( * ) (t). In fact, from (28) it follows that spectral components of E{y(t + τ ) y ( * ) (t)} are concentrated around the cycle frequencies α ∈ A. The value of τ such that |R α x (τ )| peaks should be chosen. In most cases τ = 0 is the right choice for all α.
Let h W (t) be the impulse response function of an ideal lowpass filter with monolateral bandwidth W and let us define It results that where B α is the monolateral bandwidth of a 2 (t) e j2παǫ(t) . In (31), the first inequality assures that a 2 (t) e j2πα0ǫ(t) e j2π(α0− α0)t passes unaltered through h W (t) and the second inequality assures that the terms a 2 (t) e j2παǫ(t) e j2π(α− α0)t (see (28)) with α = α 0 are filtered out. The value of W is chosen as the width of the peak around the zero frequency of the PSD of the function t → y(t + τ ) y ( * ) (t) e −j2π α0t .
From (30), it follows that a(t) can be estimated, but for a constant (with respect to t) multiplicative factor, as In addition, from (30) we have where arg[·] denotes the argument of the complex number in the brackets and mod 2π is the modulo 2π operation. Therefore, ǫ(t) can be estimated (to within an unknown constant representing a fixed time delay) as where arg uw denotes the unwrapped phase and m and q are estimates of the coefficients of the affine term in (33). They are obtained, for example, by least-squares linear regression on the available data arg uw [z α0 (t, τ )].
Under the assumption of ǫ(t) small and slowly varying (see (27)), it can be shown that the inverse function of ψ(t) = t + ǫ(t) can be closely approximated as ψ −1 (t) = t − ǫ(t). Thus, accounting for (21), an estimate of x(t) can be recovered from the available data y(t) by de-warping and compensating the amplitude modulation for all t such that a(t − ǫ(t)) = 0.
Starting from x(t), estimates of (conjugate) cyclic autocorrelation functions and (conjugate) cyclic spectra can be built. From estimates of the (conjugate) cyclic autocorrelations of x(t) and estimates of a(t) and ǫ(t), an estimate of the (conjugate) autocorrelation of the OACS signal y(t) can be obtained (see (28)). In addition, once x(t) is available, all the well known signal processing algorithms for ACS signals (e.g., FRESH filtering to remove undesired additive signals) [3], [11] can be suitably exploited.
An alternative procedure for recovering the underlying ACS signal from a time-warped ACS signal is proposed in [2]. The comparison of the techniques proposed in [2] and [12] is made in [13].

VI. NUMERICAL RESULTS
In this section, the ECG signal is modeled as AM-TW ACS and an experiment is conducted aimed at estimating the time-warping and amplitude-modulation functions and at recovering the underlying ACS signal by exploiting the estimation procedure described in Section V.
Let T 0 be the average cardiac cycle, that is, the reciprocal of the average heart rate α 0 . In practice, the average is made within the observation interval. The ECG signal is modeled here as the AM-TW ACS signal y(t) in (21), where x(t) is a real-valued cyclostationary signal with period T 0 that can be decomposed into a periodic function with period T 0 and a zero-mean cyclostationary signal with period of cyclostationarity T 0 . Thus, the cycle frequencies of x(t) are k/T 0 , k integer. The time-warping function ψ(t) = t + ǫ(t), with ǫ(t) satisfying (27), describes the variability with time of the heart rate and possibly accounts for the presence of artifacts. The time-varying amplitude a(t) is due to variation in the propagation of the electrical wave throughout the heart and possibly to artifacts. The validity of the assumptions made in Sec. IV to obtain the autocorrelation (28) are verified a posteriori. The evolutionary cyclic autocorrelation functions of the conjectured signal model for the ECG signal are given by the terms of the sum over α in (28), where the cyclic autocorrelation functions are non summable due to the presence of the periodic term in x(t).
The signal stored in the data file m001.dat taken form the CEBS database available at https://physionet.org [5] is considered. It corresponds to a healthy 30 year male. The signal y(t)− y(t) t is analyzed, where y(t) t is the temporal mean of y(t). The digitalized signal y(t)| t=nTs is obtained with sampling frequency f s = 1/T s = 250 Hz. The data is converted into Matlab/Octave format by the wfdb toolbox [18]. In the experiment, N = 2 15 samples are taken. This corresponds to a data-record length T = N T s ≃ 131 s.
In Fig. 1 (top), the magnitude of the frequency-smoothed cyclic periodogram [1,Chap. 13] of the ECG signal y(t) is reported. Accordingly with (17), cyclic features of the underlying cyclostationary signal x(t) are spread around the cycle frequencies k/T 0 .
A coarse estimate of the average heart rate α 0 1/T 0 starting from the observation interval [0, T ] is derived as α 0 = n(T )/T , where n(T ) is the number of spikes in y(t) within the observation interval. Then, the time-warping and amplitude modulation functions are estimated by the procedure outlined in Section V and an estimate of the underlying cyclostationary signal x(t) is obtained by (35).
In Fig. 1 (middle), the magnitude of the frequencysmoothed cyclic periodogram of the de-warped and amplitudemodulation compensated signal x(t) is reported. According with the conjectured model, x(t) is cyclostationary since cyclic features are evident in correspondence of cycle frequencies k/T 0 . Moreover, spikes whose shape is that of the frequencysmoothing window are present. They correspond to the Dirac deltas in the cyclic spectrum as a consequence of the additive periodic component in x(t) [11,Sec. 2]. These spikes are removed by median filtering the cyclic periodogram. In Fig. 1  (bottom), the magnitude of the median-frequency-smoothed cyclic periodogram [14] of x(t) is reported. It is an estimate of the 2nd-order cyclic polyspectrum, that is, the inverse Fourier transform of the cyclic covariance which is the cyclic autocorrelation function of the zero-mean ACS component of x(t).
The measurement results reported in Fig. 1, corroborate the conjectured model for the ECG signal. The same behavior for the estimated cyclic statistics is found using ECGs in the database measured on different subjects.

VII. CONCLUSION
The recently introduced class of oscillatory almostcyclostationary processes is characterized in the frequency domain by the Loève bifrequency spectrum. It is shown that it can be seen as obtained by that of an underlying ACS signal by spreading the spectral correlation densities around the support lines. As special case of OACS processes, the AM-TW ACS processes are analyzed and a procedure is proposed to estimate the time-warping and amplitude-modulation functions and for recovering the underlying ACS process. Then, the ECG signal is suitably modeled as AM-TW ACS process. Numerical results on real data confirm the conjectured model that the underlying ACS signal is given by the superposition of a periodic signal and a zero-mean ACS signal.