Modeling the Electrocardiogram as Oscillatory Almost-Cyclostationary Process

A new model for the electrocardiogram (ECG) signal is proposed. Specifically, the ECG signal is modeled as an amplitude-modulated and time-warped version of a cyclostationary process which is the sum of a periodic signal and a zero-mean cyclostationary term. For the proposed model, the second-order characterization is derived in both time and frequency domains. The autocorrelation function is shown to be the superposition of amplitude- and angle-modulated sine waves and the Loève bifrequency spectrum a spread version of that of the underlying cyclostationary process. The signal model belongs the recently introduced class of the oscillatory almost-cyclostationary processes. A procedure for estimating the second-order statistical functions in both time and frequency domains is outlined. The effectiveness of the proposed model and of the estimation procedure is corroborated by measurements on real ECG signals which are in full agreement with the theoretical analytical expressions. The proposed model is shown to be effective with observation intervals much larger than those adopted up to now with the classical cyclostationary model and is suitable to be exploited for arrhythmia modeling and characterization, and for diagnosis and biometric purposes.


I. INTRODUCTION
T HE analysis of the electrocardiogram (ECG) signal is at the basis of the diagnosis of cardiovascular diseases. A first step in the analysis consists in providing an accurate signal model for the ECG.
Several models have been proposed in the literature. A pole-zero model is presented in [30]. Hidden Markov models are adopted to describe accurately segmented ECG signals in [25]. In [27], a dynamical model based on three coupled ordinary differential equations is introduced which is capable of generating realistic ECG signals. In [39], a model based on wavelet neural networks is proposed. Sequential Bayesian methods are considered in [6]. The ECG signal is modeled as generated by three nonlinear oscillators in [40]. A simple mathematical model based on the sum of two Gaussian functions is proposed in [1]. Surveys of ECG signal analysis techniques are made in [22] and [44]. The usefulness of the discriminative characteristics of the ECG signal in the biometric field is addressed in [42].
Cyclostationary processes are an appropriate model for signals arising from the interaction of periodic and random phenomena. In such a case, signals are not periodic functions but have statistical functions such as distribution, mean, autocorrelation, moments, and cumulants that are periodic functions of time [9, Part II], [10], [35,Chaps. 1,2]. More generally, if more incommensurate periodicities are present in the generating mechanisms, the resulting statistical functions are poly-periodic or almost-periodic.
The structure of the electrocardiogram (ECG) signal is consequence of a periodic generating mechanism, the heart pulsation, and other less predictable phenomena like the propagation of the electrical wave throughout the heart and the body and possibly artifacts. This combination of periodic and random phenomena suggests that the cyclostationary model could be exploited for ECG signal modeling [18], [19], [29], [43], and [46]. The cyclostationary model has been exploited in [17] for the fetal PQRST extraction, in [15] for arrhythmia classification, in [16] for denoising the ECG signal, in [23] and [24] for heart and respiration rate monitoring, in [12] for removing ballistocardiogram artifacts, in [13] for heart and lung sound separation, and in [26] for heart sound selection. In all the works cited above, the length of the data collection time does not exceed 10 seconds. In fact, only within such a time interval, the heart rate can be considered approximately constant, a necessary assumption for the validity of the cyclostationary model. The chirp-cyclostationary model considered in [28] allows one to consider a data-record length of 30 seconds.
A generalization of the cyclostationary processes that accounts for an irregular pace of a not exactly periodic generating mechanism are the amplitude-modulated (AM) timewarped (TW) almost-cyclostionary (ACS) processes [11], [34], [36], [35,Sec. 14.3]. Preliminary results have shown that such a model allows one to consider for the ECG signal data-record lengths larger than those for which the cyclostationary model is valid [34]. Other approaches to model signals with time-varying cyclostationary characteristics are considered in [2], [38].
In this paper, the ECG signals are modeled as AM-TW ACS processes, which are a special case of the recently introduced class of the oscillatory almost-cyclostationary (OACS) processes [32,Sec. 6], [35,Chap. 14]. For an OACS process, the autocorrelation function is the superposition of amplitude-and angle-modulated sine waves whose frequencies are the cycle frequencies of an underlying ACS process. In particular, the ACS process can be cyclostationary. For the ECG signal, the underlying cyclostationary process is modeled as the sum of a periodic signal with period equal to the reciprocal of the average heart rate and a zero-mean cyclostationary process exhibiting the same cyclostationarity period. The second-order characterization of the proposed model is derived in the time domain in terms of time-varying autocorrelation function and in the frequency domain in terms of a smoothed version of the Loève bifrequency spectrum which is the correlation between two spectral components of the process. For the underlying cyclostationary process, the analytical expressions of cyclic spectra and cyclic autocorrelation functions are derived. Due to the presence of the periodic signal in the ECG signal model, the cyclic autocorrelation functions contain finite-strength additive sine waves and the cyclic spectra contain Dirac impulses. The cyclic covariance and the second-order cyclic polyspectrum of the underlying cyclostationary process are the cyclic autocorrelation and the cyclic spectrum, respectively, of the zeromean cyclostationary component.
The effectiveness of the proposed model is corroborated by measurements made on several ECG data available in the PhysioNet database https://www.physionet.org [14]. The amplitude-modulation and time-warping functions are estimated by amplitude and angle demodulation of one of the modulated sine waves present in the second-order lag product. Then, their effects are compensated in order to reconstruct the underlying cyclostationary signal and its cyclic spectral analysis is carried out. The proposed model for the underlying cyclostationary signal is confirmed by the presence of sinusoidal terms in the estimated cyclic autocorrelations and of impulses in the estimated cyclic spectra.
The analysis made on the ECG signals of several patients confirms the general validity of the proposed model. In addition, the analysis shows that the shapes of cyclic autocorrelations, spectra, autocovariances, and 2nd-order polyspectra are different for each patient. Therefore, biometric information and information on possible diseases can be deducted from the Fourier coefficients of the periodic component and from the cyclic statistics of the zero-mean cyclostationary component. Transitory effects and arrhythmia are linked to the behavior of the time-warping function.
The paper is organized as follows. In Section II, the proposed signal model for the ECG signal is presented (Sec. II-A) and its 2nd-order characterization in both time (Sec. II-B) and frequency (Sec. II-C) domains is derived. The statistical function estimation is addressed in Section III. Numerical results are presented in Section IV and conclusions are drawn in Section V. Proofs are reported in Appendix A.

A. SIGNAL MODEL
In this section, the proposed mathematical model for the ECG signal is presented.
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 proposed model for the electrocardiogram is the realvalued signal y(t) = a(t) x(ψ(t)) . (2.1) In (2.1), 1) x(t) is a real-valued cyclostationary process, with cyclostationarity period T 0 , that can be decomposed into the sum of a periodic component and a zero-mean term. That is, real-valued (deterministic) periodic signal with period T 0 and Fourier coefficients x h/T0 (such that x −h/T0 = x * h/T0 ), and x r (t) real-valued zero-mean cyclostationary process with cyclostationarity period T 0 . 2) ψ(t) is a deterministic time-warping function with (t) slowly varying, that is, such that where . (t) is the first-order derivative of (t). Condition (2.5) means that time variations of (t) are small if compared with 1 which is the derivative of t.
3) a(t) is a deterministic time-varying amplitude.
The function (t) describes the heart-rate variability with time. Hearth-rate variability can be due to physical activity, emotion, or arrhythmia. Both functions (t) and a(t) can be affected by several factors including variation in the propagation of the electrical wave throughout the heart and possibly artifacts. Undesired artifacts are caused by patient movement, loose or defective electrodes, improper grounding, faulty ECG apparatus.

B. SECOND-ORDER TIME-DOMAIN CHARACTERIZATION
Since in model (2.2) E{x r (t)} = 0, we have be the autocorrelation function of the cyclostationary zeromean term x r (t). In (2.7), denotes the cyclic autocorrelation at cycle frequency k/T 0 of x r (t). It is the Fourier coefficient at frequency k/T 0 of the , it follows that the cyclic autocorrelation function of x(t) is given by The signal x(t) is cyclostationary with cycle frequencies α = k/T 0 , k ∈ Z. The presence of the additive periodic term in the expression of x(t) reflects into the presence of sinusoidal terms in the cyclic autocorrelation (2.9b).
The cyclic autocovariance of x(t) is given by that is, it is coincident with the cyclic autocorrelation of the zero-mean term x r (t).
The signal in (2.1) is an AM-TW ACS process. It is a special case of oscillatory almost-cyclostationary (OACS) process whose autocorrelation function can be expressed as (Appendix A) where the functions ρ α y (t, τ ) are referred to as evolutionary cyclic autocorrelation functions and are given by (Appendix A) · e j2π(h/T0)(τ + (t+τ )− (t)) (2.14) The autocorrelation (2.11) is the superposition of amplitude-and angle-modulated sine waves with frequencies that are integer multiples of the average heart rate α 0 = 1/T 0 . The amplitude-and angle-modulation functions are the magnitude and phase, respectively, of the complex valued evolutionary cyclic autocorrelation functions ρ α y (t, τ ), α = k/T 0 . The spectral analysis of y(t) by the Loève bifrequency spectrum [35,Sec. A.3], that is the double Fourier transform with respect to t 1 , t 2 of the autocorrelation E {y(t 1 ) y(t 2 )} is very challenging. The same holds for the rotated Loève bifrequency spectrum, that is, the double Fourier transform with respect to t, τ of the autocorrelation E {y(t + τ ) y(t)}.
In the following, under mild conditions, an approximate expression of the autocorrelation E {y(t + τ ) y(t)} is derived for sufficiently small values of |τ |. From this expression, a smoothed version of the Loève bifrequency spectrum and its rotated version are obtained in Section II-C.
In Appendix A, it is shown that, for slowly varying (t) (condition (2.5)), the time-varying terms in the argument of R k/T0 xr (·) in (2.13) can be neglected for |τ | sufficiently small. Specifically, the following approximate expression holds where B r is the bandwidth of x r (t).
Let h max > 0 be the order of the maximum significant harmonic of the periodic signal x p (t). That is, h max is such that |x h/T0 | is negligible for |h| > h max . In Appendix A, it is shown that with τ max such that Such a condition, accounting for (2.5), reduces to Under the assumption that a(t) is slowly varying with respect to x p (t) and x r (t), we have a(t + τ ) a(t) for |τ | < min(τ max , 1/B r ). Thus, in (2.12) it results Therefore, under conditions (2.5) and (2.17), and accounting for (2.15), (2.16), and (2.19), with R α x (τ ) cyclic autocorrelation of x(t) given in (2.9b), and for the oscillatory cyclic autocorrelation ρ α y (t, τ ) in (2.11) the following approximate expression hold By replacing (2.21) and (2.22) into (2.11), one obtains the following approximate expression for the autocorrelation of the ECG signal y(t) This expression can be interpreted as the autocorrelation of a cyclostationary process with "time-varying cyclic autocorrelations" and "time-varying instantaneous cycle frequencies". Thus, the instantaneous cycle frequency for k = 1 in (2.23) provides the time-varying heart-rate of the proposed ECG signal model

C. SECOND-ORDER FREQUENCY-DOMAIN CHARACTERIZATION
The cyclic spectrum of the underlying cyclostationary process x(t) is the Fourier transform of the cyclic autocorrelation (2.9b): where S k/T0 xr (f ) is the cyclic spectrum of the zero-mean term x r (t), that is, the Fourier transform of R k/T0 xr (τ ), and the presence of Dirac deltas is consequence of the additive periodic term in (2.2).
For a finite-memory or practically finite-memory process x r (t), the cyclic autocorrelation function is summable. Hence, the second-order cyclic polyspectrum of x(t), which is the Fourier transform, defined in ordinary sense, of the cyclic autocovariance (2.10a) of x(t) is coincident with the cyclic spectrum of x r (t) and does not contain Dirac impulses [35,Sec. 1.4].
The Loève bifrequency spectrum of x(t) and its rotated version are given by [35, is the Loève bifrequency spectrum of y(t). In (2.30), Y (f ) formally denotes the Fourier transform of y(t).
The calculation of the double Fourier transform (2.29a) of the autocorrelation function (2.11) with (2.12)-(2.14) substituted into is a formidable problem. However, a smoothed version of the rotated Loève bifrequency spectrum can be easily obtained starting from (2.11) and substituting for (2.21). It should be noted that the approximate expression (2.21) holds for |τ | < min(τ max , 1/B r ).
By doubly Fourier transforming with respect to t and τ the windowed autocorrelation one obtains the following smoothed version of the rotated Loève bifrequency spectrum (Appendix A) A smoothed version of the Loève bifrequency spectrum of y(t) is given by (2.31), it follows that the smoothed rotated Loève bifrequency spectrum (2.31) of the ECG signal y(t) is given by the sum of the cyclic spectra of the underlying cyclostationary signal x(t) smoothed in f and spread in the δf dimension around the cycle frequencies k/T 0 . The spreading is due to the amplitude-and anglemodulation of the sine waves in the autocorrelation function. For this reason, Dirac deltas in (2.28b) are replaced by the functions M α y (·) in (2.31) that have nonzero bandwidth and depend on α. The smoothing along the f dimension is due to the convolution of the cyclic spectra with the spectral window W (f ) as a consequence of the windowing of the autocorrelation.

III. STATISTICAL FUNCTION ESTIMATION
In this section, the problem of estimating the second-order statistical functions of the ECG signal modeled as AM-TW ACS process is addressed.
The autocorrelation function and the Loève bifrequency spectrum are not directly estimated. Rather, the functions that appear in their analytical description are estimated separately. Then, they are used to build estimates of the autocorrelation function and the Loève bifrequency spectrum.
Specifically, the functions involved in the time-and frequency-domain characterizations in Sections II-B and II-C are derived as follows: (a) The amplitude-and angle-modulation functions a(t) and (t) (or ψ(t) = t + (t)) are estimated (Sec. III-A); (b) Using the estimates of a(t) and (t), amplitude modulation is compensated and de-warping is performed in order to reconstruct an estimate x(t) of the underlying cyclostationary signal x(t) (Sec. III-B); (c) Cyclic spectral analysis is carried out on the signal x(t) obtaining estimates of the cyclic autocorrela- 25a), and 2ndorder cyclic polyspectrum P α x (f ) in (2.26a) of the underlying cyclostationary process x(t) (Sec. III-C). Starting from the estimates obtained in steps (a)-(c), the estimate of the autocorrelation function (2.11) of y(t) is obtained as follows: (i) The estimates of a(t) and (t) provide the estimates of amplitude and phase (divided by 2πα), respectively, of the function m α y (τ ) in (2.22); (ii) The estimate of m α y (τ ), jointly with the estimate of R α x (τ ), provides an estimate of the oscillatory cyclic autocorrelation ρ α y (t, τ ) in (2.21); (iii) The estimates of the functions ρ α y (t, τ ), considered for a sufficient number of cycle frequencies α = k/T 0 , provide an estimate of the autocorrelation function (2.11) of y(t).
Starting from the estimates obtained in steps (a)-(c), the estimate of the smoothed Loève bifrequency spectrum (2.33b) of y(t) is obtained as follows: (iv) The estimate of the Fourier transform M α y (f ) in (2.32a) is obtained by Fourier transforming the estimate of m α y (τ ) (see (i)); (v) The estimate of the cyclic spectrum S α x (f ) is obtained at step (c); (vi) The estimates of M α y (f ) and S α x (f ) considered for a sufficient number of cycle frequencies α = k/T 0 , provide an estimate of the the smoothed Loève bifrequency spectrum (2.33b) of y(t). The above outlined procedures allow the estimation of the autocorrelation function and the Loève bifrequency spectrum of y(t). However, it is worthwhile to underline that for the objectives of the numerical analysis conducted in Section IV, only the estimates obtained in steps (a)-(c) are exploited.

A. AMPLITUDE-AND ANGLE-DEMODULATION
This function is the superposition of amplitude-and anglemodulated sine waves centered around the frequencies k/T 0 , k integer. Let us assume that the bandwidths of these modulated terms are such that their power spectra do not significantly VOLUME 4, 2016 overlap the power spectrum of the modulated sine wave with k = 1. Such an assumption has been verified in all the observed ECG signals (see Section IV).
Let α 0 be an estimate of the average heart rate α 0 = 1/T 0 obtained by counting the number n(T ) of pulses in the observation interval of width T and consider the signal is the impulse-response function of a lowpass filter whose bandwidth W is such that the term in (3.3b) passes unaltered and all the other terms in y 2 (t)e −j2π α0t are filtered out.
Under the mild assumption that the constant term in the square brackets in (3.3b) is non zero, the time-varying amplitude a(t) can be estimated as but for a nonzero constant multiplicative factor which is the magnitude of the (unknown) term in the square brackets in (3.3b).
Since a 2 (t) > 0, from (3.3b) it follows that where q denotes the phase of the (unknown) term in the square brackets in (3.3b). The affine term mt + q, with m = 2π(α 0 − α 0 ), accounts for the error in the cycle frequency estimate α 0 .
The function (t) can be estimated but for a constant additive term as where arg uw [·] denotes the unwrapped phase. In (3.6), m and q are estimates of m and q obtained by least-squares linear regression on the available data arg z ( α0,W ) (nT s ) , n = 0, 1, . . . , N − 1, where T s is the sampling period and T = N T s is the data-record length. By letting m 2π( α 0 − α 0 ), the following refined cycle frequency estimate is obtained By replacing α 0 with α 0 into (3.3a)-(3.6), refined estimates of a(t) and (t) are obtained [35, Sec. 14.3.3].

B. AMPLITUDE-MODULATION COMPENSATION AND DE-WARPING
Let us consider the signal model (2.1). If time warping never reverses time, then the warping-function ψ(t) is invertible [11]. Denoting ψ −1 (t) as its inverse, x(t) can be obtained by compensating the amplitude modulation in y(t) and then de-warping the result. Specifically, from (2.1) it follows that for all values of t such that the denominator is non zero. An estimate x(t) of the underlying ACS signal x(t) is obtained by replacing into (3.9) a(t) and ψ −1 (t) with their estimates. Let us consider ψ(t) = t + (t) with (t) slowly varying (see (2.5)). Let a(t) and (t) be the estimates of a(t) and (t) obtained by the procedure described in Section III-A. It can be shown that, for the purpose of de-warping, an estimate of ψ −1 (t) is given by t − (t) [35, Sec. 14.8]. Thus, accounting for (3.9), an estimate of the signal x(t) is obtained by y(t) as provided that a(t)/ a(t) 1, that is, the estimate a(t) is sufficiently accurate, and the estimation error on (t) is sufficiently small in the sense that The samples of the time-warped versions of y(·) and a(·) in (3.10) are obtained from those of y(·), (·), and a(·) by interpolation [31,Sec. 7.7.4]. For example, the samples with sampling period T s of the numerator in (3.10) are computed as where N is the number of available samples of y(·).

C. STATISTICAL FUNCTION ESTIMATION OF THE UNDERLYING CYCLOSTATIONARY PROCESS
Estimates of statistical functions of the underlying ACS signal x(t) are obtained through measurements of those of the de-warped signal x(t) obtained by (3.10). Specifically (see [8], [9,Chap. 13]): 1) The cyclic autocorrelation of x(t) is estimated by the cyclic correlogram of x(t) 2) The cyclic spectrum of x(t) is estimated by the the frequency-smoothed cyclic periodogram of x(t) where ∆f is the width of the frequency-smoothing window and is coincident with the spectral frequency resolution, and is the finite-time Fourier transform of x(t) on the datarecord of length T . Reliable estimates are obtained for T ∆f 1 [8], [9,Chap. 13], [35,Chap. 5]. The additive periodic component x p (t) in the signal x(t) (see (2.2) and (2.3)) gives rise to first-order cyclostationarity and, at second-order, reflects into the presence of finitestrength additive sine-wave components in the cyclic autocorrelation (see (2.9b)) and Dirac deltas in the cyclic spectrum (see (2.25b)). Consequently, the frequency-smoothed cyclic periodogram of x(t) presents spikes whose width is equal to the frequency-smoothing window-width ∆f .
Estimates of the cyclic autocovariance (2.10a) and the second-order cyclic polyspectrum (2.26a) are obtained by removing the effects of the additive periodic component. This is realized by removing the spikes in the frequency-smoothed cyclic periodogram of x(t) by a sliding median-filtering window [33,Sec. 3], [35,Sec. 5.2.5]. This method allows one to remove the effects of an almost-periodic additive component with unknown frequencies and was originally proposed in [37] in the context of cyclic polyspectrum estimation.
Specifically, the 2nd-order cyclic polyspectrum of x(t), that is, the cyclic spectrum of the zero-mean component x r (t) (see (2.26b)), is estimated by the median-filtered frequencysmoothed cyclic periodogram of x(t) [37], [ The inverse Fourier transform of the median-filtered frequency-smoothed cyclic periodogram of x(t) is adopted as estimator of the cyclic autocovariance of x(t), that is, the cyclic autocorrelation function of the zero-mean component x r (t) (see (2.10b)). Reliable estimates are obtained for T ∆f 1 and Df > ∆f (the value Df = 6 ∆f is adopted in the numerical experiments of Section IV).
Alternative techniques for removing the effects of an additive periodic component are proposed and analyzed in [20], [21], [45].
The spectral density of the rotated Loève bifrequency spectrum is obtained by considering cyclic spectrum estimates over a grid of values of α (see (2.31)).
The knowledge of the estimates a(t), (t), and of the cyclic correlograms of x(t), provides an estimate of the autocorrelation function of the amplitude-modulated timewarped ACS signal y(t).

D. EFFECTS OF NOISE AND OTHER DISTURBANCES ON THE MODEL AND THE ESTIMATION PROCEDURE
In this section, the effects of noise and other disturbances on the proposed model and the estimation procedure are briefly discussed.

1) Additive Noise
Assume that the disturbance can be modeled as additive stationary noise independent of the ECG signal. Thus, as a consequence of the well-known signal selectivity property of cyclostationarity-based signal processing algorithms [35,Sec. 9.2], the noise, also after the dewarping procedure of Section III-B, influences the estimated cyclic statistical functions only at α = 0. In particular, it does not influence these functions at cycle frequency equal to the heart rate. The smallest signal-to-noise ratio (SNR) that can be tolerated in practice depends on the considered data-record length [35, Sec. 9.2].
The same result holds under mild assumptions if the additive noise is nonstationary. In fact, the only case in which the noise has cyclic features that overlap those of the underlying cyclostationary signal x(t) is when the noise second-order lag product contains finite-strength additive sine-wave components with the same cycle frequencies of x(t) that have experienced the same time-warping of x(t).

2) Disturbance on the Amplitude-and Angle-Modulation Functions
Baseline wander and electrode movements are expected to influence the amplitude modulation function a(t) and the angle modulation function (t) in model (2.1), mainly if the data-record length is large.
The recovery of the underlying cyclostationary process x(t) can be carried out, independently of the causes of the amplitude and angle modulations, by the procedures outlined in Sections III-A and III-B, provided that the made assumptions are valid. In particular, the angle modulation introduced by (t) must be not too deep in order to have that condition (2.5) holds. In such a case, the statistical functions of x(t) can be still estimated as explained in Section III-C.
The estimated functions a(t) and (t), however, are not representative only of physiological factors but also of arti- VOLUME 4, 2016 facts. Thus, they cannot be exploited for arrhythmia or other disease analysis and diagnosis.

IV. NUMERICAL RESULTS
In this section, experiments are conducted aimed at corroborating the effectiveness of the mathematical model of ECG signal proposed in Section II. Several ECG data available in the PhysioNet database https://www.physionet.org [14] have been analyzed.
The analysis is carried out on observation intervals of several minutes. In such a case, the data-record length is such that the ACS model is not appropriate and the time-warping must be accounted for in the model. Subsampling analysis of cyclic statistic estimates [3], [4], [5] is made in both cases of non compensated and compensated time warping.

A. EXPERIMENT 1
The ECG signal is extracted from the data file m001.dat of the CEBS database (combined measurement of ECG, breathing and seismocardiogram) available at the URL https://physionet.org/physiobank/database /cebsdb/ [7]. Signals in m001.dat are recorded from Subject 001, a 30 years Caucasian male, non-smoker, with sedentary lifestyle, presumed healthy, after recent coffee intake, while listening music.
Data was acquired by a Biopac MP36 data acquisition system (Santa Barbara, CA, USA). The considered signal is channel 1 (lead I) which is devoted to measure the conventional ECG with a bandwidth between 0.05 Hz and 150 Hz. For the ECG measurement, monitoring electrodes were used with foam tape and sticky gel (3M Red Dot 2560). More details on the database description are available at the above mentioned URL and in [7].
The temporal mean value y(t) t is removed, that is, the signal y(t) − y(t) t is analyzed. Data in m001.dat are obtained with a sampling frequency of 5kHz. These data are decimated by a factor 25. Thus, the considered digitalized signal y(t)| t=nTs corresponds to a sampling frequency f s = 1/T s = 200 Hz (sampling period T s = 5 ms). The data is converted into Matlab/Octave format by the wfdb toolbox [41].
In the experiment, N = 160 · 10 3 samples are taken. This corresponds to a data-record length T = N T s = 800 s. The first 3000 samples are drawn in Fig. 1. The considered datarecord length is significantly larger than that adopted in [18], [29], [43], and [46] (less than 10 s), where a cyclostationary model is assumed for the ECG signal y(t).
A coarse estimate of the average heart rate α 0 1/T 0 is obtained by counting the spikes within the observation interval (see (3.2)). In this case, one obtains α 0 = 0.0055f s . Thus, the estimated average cardiac cycle is T 0 = 1/ α 0 182 T s 0.91 s. This corresponds to 60/ T 0 66 beats per minute.
The discrete-time counterparts of the cyclic correlogram R (α, f ) as a function of α/f s and f /f s are computed. The width of the frequency-smoothing window is ∆f = f s /256. The magnitude of the estimates are reported in Fig. 2. According to the theoretical results for the autocorrelation (2.11) and rotated Loève bifrequency spectrum (2.31) derived for the proposed model, the ECG signal y(t) (which is not ACS) exhibits cyclic features which are spread around the cycle frequencies k/T 0 of the underlying cyclostationary process x(t). The estimate of the time-averaged autocorrelation and power spectral density (PSD) of y(t) are the slices for α = 0 in Fig. 2 left and right, respectively. The time-averaged autocorrelation of y(t) is not periodic in τ since y(t) does not contain an additive periodic component but, rather, contains modulated sine waves.
In Fig. 3, the energy of the cyclic correlogram as a function of α/f s (left) and the energy of the frequency-smoothed cyclic periodogram as a function of α/f s (right) are reported, where T is a set in which the cyclic correlogram is significantly non zero and B = (−f s /2, f s /2). If the signal y(t) were cyclostationary, one should obtain functions λ  Fig. 3 is significantly larger than 1/T and is compatible with irregular cyclicity, that is, with the presence of additive amplitude-and/or angle modulated sine waves in the autocorrelation function (see (2.11)). Since a(t) 1 (see Fig. 6 right), for α = 0 we have M 0 y (f ) δ(f ) and there is no spread of the PSD around the line α = 0.
A further confirmation of the presence of irregular cyclicity in the autocorrelation function is given by the subsampling analysis [3], [4], [5] of the cyclic correlogram. The block size is b = 100 T 0 (with T 0 estimated cardiac cycle) and the block-overlap factor is 1/4. In Fig. 4, (left) real part and (right) imaginary part of the cyclic correlogram R (T ) y (α, τ ) as a function of τ /T s at α = α 0 (estimated average heart rate) are reported (thick line). The shaded area is the region within 95% confidence interval. The high variability of the estimate is in accordance with the fact that the (regular) cyclostationary model is not appropriate for the ECG signal y(t) if the data-record length exceeds few seconds.
In order to estimate (t), a(t) and then to recover an estimate of x(t) by de-warping, the procedure described in Section III is carried out. The low-pass filter h W (t) has monolateral bandwidth W = 0.003 f s whose value is obtained by visual inspection of the PSD of y 2 (t) e −j2π α0t (Fig. 5).  The condition that the amplitude-and angle-modulated sine waves with k = 1 in (3.1) have spectra that practically do not overlap the power spectrum of the modulated sine wave with k = 1 (so that (3.3b) holds) is verified.
The estimates (t) and a(t) obtained by (3.6) and (3.4) are reported in Fig. 6 left and right, respectively. From this figure, we see that sup t | (t)| 500 T s which is greater than T 0 . That is, the time-warping is relevant and cannot be neglected.
The refined cycle frequency estimate obtained at the end of the angle-demodulation procedure (see (3.7)) is α 0 = 1/ T 0 = 0.00558018f s .
The estimate x(t) of the underlying cyclostationary signal x(t) is obtained by amplitude-modulation compensation and de-warping from y(t) by (3.10) (Section III-B).
In order to get estimates of the cyclic autocorrelation and the cyclic spectrum of x(t), the discrete-time counterparts of the cyclic correlogram R  Fig. 7. Unlike the ECG signal y(t) (Fig. 2), x(t) has cyclic features well concentrated around the cycle frequencies k/ T 0 by confirming the conjectured model (2.1) with underlying cyclostationary signal x(t).
In Fig. 8, the energy of the cyclic correlogram of x as a function of α/f s (left) and the energy of the frequencysmoothed cyclic periodogram of x as a function of α/f s (right) defined according to (4.1) and (4.2), respectively, with y replaced by x, are reported. As further confirmation of the proposed cyclostationary model for x(t), unlike the ECG signal y(t) (Fig. 3), x(t) has cyclic features well concentrated around the cycle frequencies k/ T 0 .
The slices at α = α 0 of the cyclic correlogram and of the frequency-smoothed cyclic periodogram are reported in Fig. 9. The presence of an additive periodic component in R x (α, τ ) ( Fig. 9 (left)) is in agreement with (2.9b). The spikes in S (T,∆f ) x (α, f ) ( Fig. 9 (right)), whose shape is the frequency-smoothing window with width ∆f , correspond to the Dirac pulses in (2.25b).
The estimates of the cyclic autocorrelation and the cyclic spectrum of the zero-mean cyclostationary component x r (t) are obtained by the median-filtering procedure described in Sec. III-C and are reported in Fig. 10. They are obtained as the estimate C (T,∆f,Df ) x (α, τ ) of the cyclic autocovariance and the estimate P (T,∆f,Df ) x (α, f ) of the 2nd-order cyclic polyspectrum of the underlying cyclostationary signal x(t).
As a further confirmation of the effectiveness of the proposed model for the ECG signal y(t) with an underlying cyclostationary signal x(t), a subsampling analysis of the cyclic correlogram of x(t) is carried out. The block size is b = 100 T 0 and the block-overlap factor is 1/4. In Fig. 11, (left) real part and (right) imaginary part of the cyclic correlogram R (T ) x (α, τ ) as a function of τ /T s at α = α 0 (estimated average heart rate) are reported. The shaded area is the region within 95% confidence interval. The very low variability of the estimate confirms the (regular) cyclostationary model for x(t).

B. EXPERIMENT 2
In order to corroborate the effectiveness of the proposed model for the ECG signal, the analysis made in Section IV-A is repeated for all the subjects whose ECG measurements are available in the CEBS database. For the sake of brevity, only results of three subjects are reported here. Results for the other subjects are similar.
The analysis parameters are the same as those in the experiment of Section IV-A. The only difference is in the value of the estimated average heart rate α 0 since it depends on the subject.
In Figs. 12-14, results are reported for the subjects 002, 004, and 010. For all the subjects the results are in accordance with the proposed theoretical model described in Section II.
In the left-top sub-figures of Figs. 12-14, the thick line represents the PSD of y 2 (t) e −j2π α0t as a functions of f /f s VOLUME 4, 2016 and the dotted line delimits the portion of spectrum within the band of the low-pass filter h W (t). For all the subjects, the spectral region corresponding to the modulated first harmonic in the autocorrelation is well separated from the spectral regions of the neighboring harmonics and can be easily extracted by the filtering procedure (3.3a), (3.3b). In the leftbottom sub-figures, the estimated (t) as a function of t/T s is reported. Spectral analysis of the estimated (t) is suitable to be exploited for arrhythmia diagnosis and characterization.
In the middle-top sub-figures of Figs. 12-14, the energy λ (T ) y (α) (see (4.1)) of the cyclic correlogram of the ECG signal y(t) as a function of α/f s is reported. According to (2.31), the cyclic features are spread around the cycle frequencies of the underlying cyclostationary signal x(t). In the middle-bottom sub-figures, the energy λ (T ) x (α) of the cyclic correlogram of the estimated underlying cyclostationary signal x(t) as a function of α/f s is reported. For all the subjects, the cyclic features of x(t) are well concentrated around the integer multiples of the estimated average heart rate α 0 confirming the validity of the proposed model for the ECG signal.
In the right-top sub-figures of Figs. 12-14, the magnitude of the estimate (3.13) of the cyclic autocorrelation of x(t) at cycle frequency α = α 0 as a function of τ /T s is reported. In all plots, that is, for all the subjects, an approximately periodic pattern can be recognized, with superimposed a function which is nonzero around τ = 0. This behavior is in agreement with the analytical model (2.9a)-(2.10b). In fact, the cyclic autocorrelation R α x (τ ) is the sum (see (2.9b)) of a periodic function and the autocorrelation R α xr (τ ) of the residual term which is coincident with the autocovariance C α x (τ ) (see (2.10b)). The function C α x (τ ) is summable and, hence, approaches zero for large values of |τ |.
By comparing the right-top subfigures of Figs. 12-14 (that correspond to different subjects), it can be observed that different subjects are characterized by different shapes of the cyclic autocoavriance and different periodic patterns. Thus, these features can be exploited for biometric purposes. In the right-bottom sub-figures, the magnitude of the estimate (3.16) of the second-order cyclic polyspectrum of x(t) at cycle frequency α = α 0 as a function of f /f s is reported. Also this feature is different for different subjects and is suitable to be exploited for biometric purposes.      (α, f ) as a function of α/fs and f /fs. Unlike the ECG signal y(t) (Fig. 2), x(t) has cyclic features concentrated around the cycle frequencies k/T0. Napolitano: Modeling the Electrocardiogram as Oscillatory Almost-Cyclostationary Process   A new model for the electrocardiogram (ECG) signal is proposed. The ECG signal is modeled as an amplitudemodulated and time-warped cyclostationary process. The underlying cyclostationary process is decomposed into the sum of a deterministic periodic signal and a zero-mean cyclostationary term. The period of the periodic signal and the cyclostationarity period of the zero-mean term are both equal to the reciprocal of the average heart rate. The time-warping function describes the heart rate variability due to different factors like physical activity, emotions and arrhythmia. Both amplitude-modulating and time-warping functions are also consequence of variations in the propagation of the electrical wave throughout the heart. For the proposed model, the second-order characterization is derived in both time and frequency domains. The autocorrelation function is found to be the superposition of amplitude-and angle-modulated sine waves whose frequecies are the cycle frequencies of the periodic autocorrelation function of the underlying cyclostationary signal. Thus, the proposed ECG signal model belongs to the recently introduced class of the oscillatory almost-cyclostationary processes. Since the derivation of the Loève bifrequency spectrum is a formidable problem, a smoothed version is derived. The presence of the additive periodic term in the signal model reflects into the presence of finite-strength additive sine waves in the cyclic autocorrelations and Dirac impulses in the cyclic spectra of the underlying cyclostationary signal.
Estimators of the amplitude-and angle-modulation functions are derived and a de-warping and amplitude-modulation compensation procedure for estimating the underlying cyclostationary signal is outlined. Thus, for such a signal, cyclic autocorrelation, autocovariance, spectrum, and 2ndorder polyspectrum are estimated.
Real ECG signals taken from the CEBS database of the PhysioNet database are analyzed. The assumptions made to develop the de-warping procedure are found to be satisfied by all the analyzed ECG signals. Thus, the de-warping procedure is successfully adopted to restore the underlying cyclostationary signal and its cyclic statistical functions are estimated. The measurement results are in full agreement with the proposed theoretical model. In particular, the estimates of the cyclic autocorrelation functions contain additive sinusoidal terms and, accordingly, the estimates of the cyclic spectra contain spikes. The estimates of the cyclic autocovariance and of the cyclic second-order polyspectrum, that are coincident with the cyclic autocorrelation and cyclic spectrum of the zero-mean residual term, respectively, are also obtained.
The considered data-record lengths are 800 seconds. They are much larger than those of few seconds adopted up to now with the classical (regular) cyclostationary model. The larger data-record length, which is compatible with the proposed ECG signal model, allows the analysis of diseases like the arrhythmia.
The subsampling analysis of the cyclic autocorrelation estimates provides a further confirmation of the effectiveness of the proposed oscillatory almost-cyclostationary model for data-record lengths exceeding few seconds. In fact, the estimate of the cyclic autocorrelation made on the original ECG signal exhibits a wide 95% confidence interval. It is consequence of the large variability of the estimate along the blocks in accordance with the fact that the classical cyclostationary model is compatible only with a constant heart rate within the observation interval. In contrast, if time-warping is compensated and the estimate is made on the reconstructed underlying cyclostationary signal, a narrow 95% confidence interval for the cyclic autocorrelation estimate is found. Measurements made on the ECG signals of different subjects give rise to underlying cyclostationary signals whose additive periodic components have different Fourier coefficients and whose zero-mean terms have different cyclic statistical functions. Thus, these cyclic features are suitable to be exploited for biometric purposes.
Further numerical results not presented here have confirmed that the cyclic features of the underlying cyclostationary process remain the same for the same subject when measured in different time intervals. In contrast, the amplitudemodulation and time-warping functions, which are consequence of physical activity, emotions, and arrhythmia, depend on the time interval in which they are measured.
Future work will explore the spectral analysis of the timewarping function for the arrhythmia diagnosis and characterization and the exploitation of the cyclostationary features for biometric applications. .

DERIVATION OF (2.31)
By using the autocorrelation expression (2.11) with (2.21) replaced into, one has