Cyclostationary Analysis for Heart Rate Variability

Cyclostationary Analysis for Heart Rate Variability

The Open Bioinformatics Journal 31 Jul 2018 LETTER DOI: 10.2174/1875036201811010164



During the last years, cyclostationarity has emerged as a new approach for the analysis of a certain type of non-stationary signals. This theoretical tool allows us to identify periodicity in signals which cannot be identified easily but also to separate useful signals for other interfering contributions that overlap in the spectral support.


The aim of this work is the exploitation of cyclostationary theory to enhance standard methodologies for the study of heart rate variability. In this framework, a preliminary analysis on healthy patients is proposed to be extended further on pathological patients with the perspective to improve (hopefully) the diagnostic power of some cardiac dysfunctions due to the more complete set of information provided by this analysis.


The proposed approach involves an initial band-pass filtering step in the range 0.5 - 40 Hz of the recorded ECG signal, followed by a first-order derivative filter to reduce the effects of P and T waves and to emphasise the QRS contribution. After that, the auto-correlation function is evaluated and the Cyclic Power Spectrum (CPS) is computed. From this two-dimensional information, a one-dimensional plot is derived via the evaluation of a folded-projected CPS to be compared with standard Lomb-Scargle spectrum.


The proposed analysis has been tested on both numerical simulations as well as for the processing of real data which are available online in the Physionet database.


The proposed cyclostationary analysis has shown a good agreement with the results provided by the classical Lomb-Scargle spectrum in the processing of real data, underlining some contributions in the high-frequency bandwidth which are not visible by means of standard processing.

Keywords: Electrocardiogram (ECG), Heart Rate Variability (HRV), Cyclostationarity signal analysis, Cyclic power spectrum, Spectral support, Lomb-Scargle spectrum.


Heart Rate Variability (HRV) is the variation in the time interval between consecutive heartbeats. It is a physiological phenomenon related to the regulation of the cardiac activity produced by the Autonomic Nervous System (ANS) [1]. More in detail, the accelerations and decelerations of the heard activity are due to the competing activities of the sympathetic and the parasympathetic nervous system branches [2]. The analysis of the HRV has the advantage of being noninvasive, easy to be performed and with good reproducibility, and is useful for the determination of the ANS status [3] and the cardiac activity [4]. In particular, it has shown to have a beneficial role in the diagnosis and analysis of several pathologies related to blood pressure [5], myocardial infarction [6, 7], brain damage [8], depression [9], cardiac arrhythmia [10], diabetes [11] and renal failure [12]. HRV analysis has also shown correlations with sleep [13], drugs [14] or alcohol [15] assumption and smoking [16].

A first important step in order to provide HRV measures consists of the R peak detection in the QRS wave. Regarding this aim, many robust algorithms are known in the literature, among which the well-known Pan-Tompkins algorithm [17] represents a gold standard for the scientific literature.

Variation in heart rate can be evaluated by using four main classes of methods. The first one is based on the evaluation of some synthetic parameters in the time domain. In a continuous electrocardiographic (ECG) record, each QRS complex is detected, and the so-called normal-to-normal (NN) or RR intervals (that is all intervals between adjacent QRS complexes resulting from sinus node depolarizations), or the instantaneous heart rate is determined. Simple time-domain variables that can be calculated include the mean RR interval, the mean heart rate, the difference between the longest and shortest RR interval, and many others [18, 19]. Other time-domain measurements analyse the variations in instantaneous heart rate secondary to respiration, tilt, Valsalva maneuver, or secondary to phenylephrine infusion. These differences can be described as either differences in heart rate or cycle length.

A second important class is represented by frequency domain methods, which are mainly based on the estimation of the Power Spectral Density (PSD) [20]. The analysis of this function provides the basic information of how signal power distributes as a function of frequency. Both parametric and non-parametric methods are available for the PSD estimation, each one with its advantages and limitations. Non-parametric methods are simple and fast, since they are mostly based on the use of the fast Fourier transform, while parametric methods provide smoother spectral components and a more accurate estimation of the PSD even on a small number of samples on which the signal is supposed to maintain the stationarity hypothesis. Among these techniques, the Lomb-Scargle (LS) periodogram is perhaps one of the most well-known techniques employed to compute the periodicity of unequally-spaced data, and it provides a good estimate of the PSD of an ECG signal [21-25]. LS method avoids the major problem of classical approaches related to the low-pass effect due to re-sampling operation. Therefore, the Lomb method is more suitable than fast Fourier transform or autoregressive estimate with linear or cubic interpolation for PSD estimation of unevenly sampled signals [26]. However, in extreme situations (low heart rate or high-frequency components), the Lomb estimate still introduces high-frequency contamination that suggests further studies on superior performance interpolators.

Both time-domain and spectral methods share some limitations imposed by the irregularity of the RR series, since they assume the same trends of increasing or decreasing in the cycle length, which is not always realistic [27, 28]. In practise, this reflects on the amplitude of the peaks at fundamental frequencies in the spectral analysis and enlarge their basis. In order to overcome these limitations, some alternative techniques aiming at analysing the rhythm pattern via considering blocks of RR intervals without considering the internal variability have been proposed. The well-known interval spectrum and spectrum of counts methods are well suited to investigate the relationship between HRV and the variability of other physiological measures, like blood pressure, respiration and arrhythmia events [29].

Last class of approaches for HRV analysis is based on non-linear methods [30-34]. The motivation for using this kind of methodologies is based on the non-linear phenomena involved in the genesis of the HRV signal, which is the result of complex interactions among haemodynamic, electrophysiological and humoral variables. However, although in principle, these techniques have been shown to be powerful tools for characterization of various complex systems, their application to HRV (and biomedical data more in general) still needs to be validated. An exhaustive overview of the main publications for non-linear HRV analysis for the past 25 years has been proposed in reference [35]. In their work, Sassi et al. present a critical review of the state of the art and new methodologies tested in sufficiently sized populations, with particular attention paid to the long-range correlation and fractal analysis, entropy and regularity, non-linear dynamical systems and chaotic behaviour.

Conversely from the approaches described previously in the text, this short communication aims at testing and analysing the performance of cyclostationary (CS) analysis for HRV signals. In most cases, the stationarity hypothesis is often an assumption of convenience rather than a realistic one, and for many biomedical signals (e.g., ECG, EMG, etc.) the use of cyclostationarity, i.e. the cyclic variation of the statistical properties, is more suitable rather than conventional stationary assumption [36-38]. Compared to standard analysis, it provides some extra information due to the hidden periodicities in the signal. Taking advantage of these properties obviously leads to more powerful processing than possible with the stationary approach. Moreover, this extra information often counter-balances the complication that it may involve as compared to other standard approaches.

Methodologies exploiting CS theory for heart monitoring have been proposed by several authors. In reference [39], a CS approach for heart and respiration rates monitoring exploiting a 2.4 GHz Doppler radar is presented. CS theory has also been applied to ECG signals for several applications, such as foetal ECG extraction [40] or for the study of non-linearity in the HRV signal [41]. In this framework, we propose to apply the CS methodology to the ECG signal in order to analyse the HRV. The main advantage of the CS analysis for HRV signals consists in the possibility of studying simultaneously the standard spectral coverage of the ECG signal in the classical frequency domain, but also the cyclic spectral components related to the physiological behaviour of the HRV. Another important advantage is that the analysis is performed directly in the time domain after the ECG recording and does not require any RR interval extraction, which makes this tool more robust.

The remainder of the paper is organised as follows: Section 2 provides a mathematical overview of CS theory, Section 3 proposes the methodology involved in this paper and Section 4 presents some numerical results and real data processing. Finally, Section 5 closes the paper and draws some conclusions.


CS extends the class of stationary signals to those signals whose statistical properties change periodically with time. In this theoretical framework, the minimum “period” of a CS signal is called cycle [38].

Conversely from standard non-stationary signals, CS is a well-defined property and can exploit a powerful spectral analysis in a wide sense, employing the same tools that have been developed historically for stationary signals [42]. The instantaneous auto-correlation function of a random signal , with Δ denoting the sampling period, can be defined as:


where n and τ are two time variables, is the statistical mean operator and β The β. The parameter in Eq. (1) allows a general formulation of various equivalent definitions found in the literature (for example typical values are β = 1/2 for the symmetric instantaneous autocorrelation function, and β = 1 or β = 0 for the asymmetric case) [38].

By definition, the instantaneous auto-correlation function of a (quasi-) CS signal is (quasi-) periodic, and therefore it can be decomposed in a Fourier series:


over the spectrum A = {αi} of the cyclic frequencies associated to non-zero Fourier coefficients:


Coefficients in Eq. (3), which are function of the time τ and of the cyclic frequency ai, are called the cyclic auto-correlation function of the random signal X[n]. Finally, by exploiting Eq. (3), the Cyclic Power Spectrum (CPS) can be defined as:


As well known in the spectral analysis, there is no consistent estimator of the CPS, i.e. an estimator whose variance tends to zero as the length of observation time increases [42, 43] Although not consistent, the averaged cyclic periodogram is one of the most adopted spectral estimators. Given a sequence x{n} of the random CS signal and a positive and smooth Nw-long window W[n], and Wk [n] = W[n - KR] its shifted version of multiple of R samples, the averaged cyclic periodogram can be computed via:




is the short Discrete-Time Fourier Transform (DTFT) of the K-th weighted sequence (where the operator. rounds the argument to rounds the argument to the next smaller integer).

In this manuscript, we want to evaluate the CPS of an ECG signal, therefore, we apply Eq. (5) for the estimation of CPS in order to perform a more complete spectral analysis.


Let us consider a recorded ECG signal:


where z (t) is the noise free signal and n (t) is the additive noise term which can be modelled as uncorrelated and Gaussian distributed.

As first step after the acquisition, a bandpass filter is applied in order to remove the frequency components below 0.5 Hz, which are mainly related to breathing, and the components above 40 Hz, which are due to external interferences. Next, a derivative filter is exploited to isolate and strengthen the QRS information. The motivation behind this choice lies in the attempt to reduce the effects of the P and T waves compared to the QRS complex since the main information regarding HRV is related to the QRS peaks.

Subsequently, the CS spectrum is computed. As described in Section 2, the procedure requires the computation of the auto-correlation function followed by the estimation of the cyclic power spectrum . In other words, Eqs. (1, 5 and 6) are computed.

Conversely from classical HRV spectra, the information arising from the cyclostationary analysis has two dimensions, i.e. the spectral frequencies and the cyclic frequencies f, and spread on both positive and negative frequency a values. Therefore, in order to isolate the HRV information and improve its readability, a further spectrum, namely folded-projected cyclic power spectrum (FPCPS), is computed. The first step consists in integrating the CPS over the spectral frequencies:


The signal , with respect to the CPS, has lost the information about the frequency components of the ECG signal y (t), keeping only the dependency on the cyclic frequencies α which contains information about the HRV.

Another transformation is required in order to refer to the shift of the cyclic frequency values α from the central cyclic frequency αav, i.e. the mean heart rhythm, instead of considering the absolute values. Moreover, since we are interested in the amplitude of the cyclic frequency shifts from the average value, both positive and negative shifts are worth for this evaluation. Thus, a “folding-and-sum” operation around the zero cyclic frequency is performed, obtaining the FPCPS:


in which u(.) represents the Heaviside step function.


In order to show the potentiality of cyclostationary analysis for the study of Heart Rate Variability (HRV), results related to both simulated and real case study are presented.

4.1. Simulated Data: Gaussian Template

A synthetic signal z (t) was generated by replicating M times a template function over time. The template function simulated the signal related to one heart beat, while the delays of the replicas take into account the heart rhythm and the HRV. Analytically, the simulated signal can be expressed as:




and RR [i] is the time distance between the i-th and the (i + 1)- th heartbeat. For this simulation, we considered as template p(t) a Gaussian function and the following model for the RR intervals:


in which RRav is the average RR interval (constant), and α1 represents the amplitude of the deviations from the average RR value. In other words, we considered a sinusoidal behaviour of the HRV. In this example, we assumed RRav = 1 s and a1 = 0.07 , corresponding to a variation of approximately compared to the RRav f1 refers to the frequency involved in the numerical model and the parameter m identifies the m-th RR interval. This choice seems to be realistic, as proved by averaged oscillations in healthy patients. We assumed the sampling frequency fs = 200 Hz, a number of heartbeats M = 100 and an HRV frequency f1 = 0.1 Hz. Fig. (1). reports the Gaussian template function employed for the synthetic signal generation and the RR-interval amplitudes described by the model in Eq. (12).

Fig. (1). Gaussian template function (a) and RR-interval amplitudes (b) of model in Eq. (14).

The computed CPS is reported as surface and as an image in Fig. (2a and b), respectively. On the spectral frequency axis (f), information about the components of the signal template can be appreciated, while on the cyclic frequency axis (α), the HRV can be evaluated. For the sake of clarity, the figures show the cyclic frequency axis (α) centred at αav = 1/RRav, which corresponds to the zero cyclic frequency. In can be noticed that the peaks at are in correspondence with the considered HRV frequency f1 Hz. Of course, harmonics at 0.2 Hz and 0.3 Hz are present.

Fig. (2). Numerical simulation considering a Gaussian template function and one cosine RR model (at 0.1 Hz): 3D (a) and 2D (b) representation of the CPS, (c) FPCPS and (d) LS spectrum.

In Fig. (2c), the FPCPS computed according to Eq. (9) is reported. It can be observed that both the contributions at f = 0 (corresponding to RRav) and f1 = 0.1 Hz are present in the spectrum, even though also some other harmonics at multiple of f1 can be detected in the spectrum. In order to provide a reference, the LS spectrum is reported in Fig. (2d). Such a result confirms the presence of the HRV frequency at 0.1 Hz. It should be noted that the DC component does not appear in the LS spectrum due to the fact that the Matlab function used for this comparison does not evaluate the component centred at f = 0 Hz.

A Monte Carlo (MC) simulation has been implemented in order to evaluate the robustness with respect to noise. First, a white Gaussian noise has been added to the ECG signal Z (t) of Eq. (10), i.e.


where with the noise standard deviation σ fixed. One hundred signal realizations have been considered, and for each, the CS spectrum has been computed and the position of the contribution close to f = 0.1 Hz has been measured. The analysis has been repeated for σecg varying between 1 (SNR of about 50 dB) and 50 (SNR of about 18 dB), and the results are reported in Table 1. It can be appreciated that the mean value of the peak position slightly changes (in case of σecg = 50 it increases of about 11%). Of course, in case of high noise level, the stability of the peak position greatly decreases, with a standard deviation of 0.04 in case of σecg = 30 and higher. Nevertheless, it can be noted that, in case of noise σecg = 10 (SNR of about 30dB) or below, the standard deviation of the estimation is in the order of magnitude of 1.10-3 or lower.

Table 1.
Mean and STD of the frequency peak position in case of noise corrupting the ECG signal.
Noise Level Peak Position (Mean) Peak Position (STD)
σecg = 1 0.0992 0.0000
σecg = 2 0.0992 0.0001
σecg = 3 0.0993 0.0006
σecg = 5 0.0995 0.0006
σecg = 10 0.0996 0.0013
σecg = 20 0.0994 0.0250
σecg = 30 0.1053 0.0407
σecg = 50 0.1132 0.0466

A second MC simulation has been implemented in order to evaluate the effect of noise corrupting the RR intervals. Samples from a white Gaussian random variable have been added to the time instants Tm defined in Eq. (11). As in the previous simulation, 100 MC iterations have been considered and for each, the frequency peak close to f = 0.1 Hz has been identified. The performances were evaluated in case noise standard deviation σt ranging from 5.10-3 to 2.10-1, and results are reported in Table 2. Compared to the previous MC simulation, a different behaviour can be appreciated: When the noise standard deviation is greater than, performances rapidly deteriorate, while remaining good in case of lower values of σt.

Table 2.
Mean and STD of the frequency peak position in case of noise corrupting the RR interval values.
Noise Level Peak Position (Mean) Peak Position (STD)
σt = 0.005 0.0995 0.0012
σt = 0.001 0.0997 0.0022
σt = 0.002 0.1001 0.0048
σt = 0.003 0.0995 0.0810
σt = 0.004 0.1002 0.0105
σt = 0.005 0.0978 0.0180
σt = 0.01 0.0903 0.0419
σt = 0.02 0.1258 0.0591

As a second numerical example, a second cosine component was added to the RR model:


with α2 = 0.07 and f2 = 0.18 Hz and the other parameters equal to the previous case. As in the previous case, the FPCPS and the LS periodogram are reported in Fig. (3).

Again, there is a good agreement between the LS spectrum and the FPCPS, as shown in Fig. (3c, d). Conversely from the previous case study, the presence of more than one frequency is responsible for the intermodulation effects, which results in more lobes than what was expected from the case of a single frequency. Intermodulation effects which are visible in the LS case are still visible in the FPCPS case.

Fig. (3). Numerical simulation considering a Gaussian template function and two cosines RR model (at 0.1 Hz and 0.18 Hz): 3D (a) and 2D (b) representation of the CPS, (c) FPCPS and (d) LS spectrum.

4.2. Simulated data: QRS template

A more realistic numerical example in which the signal template p(t) is a real QRS complex is considered. Again, we considered the RR model containing two sinusoidal components reported in Eq. (14) with α1 = α2 = 0.07, f1 = 0.1 and f2 = 0.18 Hz.

The results are shown in Fig. (4). The HRV information carried by the cyclic frequencies α was close to the one reported in the previous examples, but, as expected, in this case, the CPS was richer in the spectral components. More in detail, the adopted QRS signal template p(t) was characterized by a much wider and complex spectrum in comparison with the Gaussian-shaped template considered in Section 4.1, and therefore the CPS shows several components up to f = 30Hz.

Fig. (4). Numerical simulation considering an ECG template function and two cosines RR model (at 0.1 Hz and 0.18 Hz): 3D (a) and 2D (b) representation of the CPS, (c) FPCPS and (d) LS spectrum.

4.3. Real data

To validate the proposed numerical analysis, some real datasets from the Physionet repository have been considered [44]. Three different patients with physiological ECG were analysed (ID: 16272, 16483 and 16539) from the MIT-BIH normal sinus rhythm database [45]. It includes 18 long-term ECG recordings of subjects which were found to have had no significant arrhythmias. The people involved in this study include people aged from 20 to 50. Signals have been acquired with a sampling frequency equal to 128 Hz for a duration of 300 seconds.

As previously pointed out in Section 3, the ECG signals have been pre-processed by means of a band-pass filter and a derivative filter in order to emphasize the effect of QRS complex compared to the other waves. Then, the CPS of these signals was evaluated via Eq. (5) and reported in Fig, (5).

Fig. (5). CPS of the ECG signals filtered as shown in Section 3 for three different healthy patients. (a) ID: 16272, (b) ID: 16483 and (c) ID: 16539.

These kind of images provide both spectral information as well as cyclic spectral representation, but it is not easy to be interpreted as it is. A way to ease the analysis of this information lies in the use of Eq. (9), which is reported in Fig. (6) and compared with standard LS spectrum.

Fig. (6). Comparison between FPCPS and LS spectrum for three different clinical datasets (downloaded from the physionet database [47]). Patients ID: 16272 (a-b), 16483 (c-d), and 16539 (e-f). The dotted-dash line refers to the envelope of the signals which can help the understanding of the figure.

In the resting physiological subject, three main spectral components for the HRV analysis can be distinguished in short-term ECG recordings: Very Low Frequency (VLF), Low Frequency (LF) and High Frequency (HF). The distribution of the power and the central frequency of LF and HF are not fixed and may vary due to changes in the modulation by the autonomic nervous system [46], whereas the HF components, synchronous with the respiration, occurs at 0.25 Hz approximately. The study of VLF ( Hz) phenomena, which might contain clinically relevant information, requires long-period uninterrupted data; thus, the DC component and the whole range of the VLF have not been addressed in this manuscript, and therefore the frequency axis will be limited in the useful range (0.02-0-32) Hz. Both total power of HRV spectrum as well as the LF-to-HF ratio has proved to be selective indices of cardiac parasympathetic activity [44], and this motivates the interest in the study of these components.

In order to ease the analysis of the two-dimensional spectra proposed in Fig. (5), an easier and more straight comparison has been carried out via the evaluation of Eq. (9) and comparing these results with the standard LS spectrum. The proposed comparison has been reported in Fig. (6), in which the FPCPS per each patient is compared with the corresponding LS spectrum. From a first analysis, it can be observed that for all the considered cases, most of the signal power is located in the VLF and LF bandwidths. All the cases reported in Fig. (6) are characterised by a main lobe in the VLF range and by another strong contribution of around 0.1 Hz, both in the FPCPS as well as LS spectra. Conversely from the standard LS analysis, the FPCPS also underlines some non-negligible contributions in the HF bandwidth, but still fits with the clinical considerations arising from standard analysis of HRV LS spectrum. Hopefully, this apparent higher sensitivity to HF bandwidth could be beneficial in the early diagnosis of some pathological dysfunctions. However, further investigation is required to test the proposed formulation for pathological patients prior to clinical use.


In this paper, the performance of CS for HRV signal analysis has been evaluated. More in detail, the performance of a CPS estimator has been analysed and tested for synthetic numerical experiments with Gaussian function replica and ECG-like waveforms, proving a good agreement with the LS spectrum. Moreover, CS is also able to evaluate signal frequency components and does not require the RR interval extraction from the ECG signal.

Finally, CS analysis has also been conducted on real data, confirming the interesting performance observed numerically.


Not applicable.


The authors declare no conflict of interest, financial or otherwise.


Declared none.


Rajendra Acharya U, Paul Joseph K, Kannathal N, Lim CM, Suri JS. Heart rate variability: A review. Med Biol Eng Comput 2006; 44(12): 1031-51.
Saul JP. Beat-to-beat variations of heart rate reflect modulation of cardiac autonomic outflow. Physiology (Bethesda) 1990; 5(1): 32-7.
Esler M. The autonomic nervous system and cardiac arrhythmias. Clin Auton Res 1992; 2(2): 133-5.
Kamath MV, Fallen EL. Correction of the heart rate variability signal for ectopics and missing beats Heart rate variability 1995; 75-85.
Westerhof BE, Gisolf J, Stok WJ, Wesseling KH, Karemaker JM. Time-domain cross-correlation baroreflex sensitivity: Performance on the EUROBAVAR data set. J Hypertens 2004; 22(7): 1371-80.
Duru F, Candinas R, Dziekan G, Goebbels U, Myers J, Dubach P. Effect of exercise training on heart rate variability in patients with new-onset left ventricular dysfunction after myocardial infarction. Am Heart J 2000; 140(1): 157-61.
Carney RM, Blumenthal JA, Freedland KE, et al. Low heart rate variability and the effect of depression on post-myocardial infarction mortality. Arch Intern Med 2005; 165(13): 1486-91.
Lowensohn RI, Weiss M, Hon EH. Heart-rate variability in brain-damaged adults. Lancet 1977; 1(8012): 626-8.
Carney RM, Blumenthal JA, Stein PK, et al. Depression, heart rate variability, and acute myocardial infarction. Circulation 2001; 104(17): 2024-8.
Ge D, Srinivasan N, Krishnan SM. Cardiac arrhythmia classification using autoregressive modeling. Biomed Eng Online 2002; 1(1): 5.
Pfeifer MA, Cook D, Brodsky J, et al. Quantitative evaluation of cardiac parasympathetic activity in normal and diabetic man. Diabetes 1982; 31(4 Pt 1): 339-45.
Axelrod S, Lishner M, Oz O, Bernheim J, Ravid M. Spectral analysis of fluctuations in heart rate: an objective evaluation of autonomic nervous control in chronic renal failure. Nephron 1987; 45(3): 202-6.
Scherz WD, Fritz D, Velicu OR, Seepold R, Madrid NM. Heart rate spectrum analysis for sleep quality detection. EURASIP J Embed Syst 2017; 2017(1): 26. [].
Pater C, Compagnone D, Luszick J, Verboom C-N. Effect of Omacor on HRV parameters in patients with recent uncomplicated myocardial infarction - A randomized, parallel group, double-blind, placebo-controlled trial: study design [ISRCTN75358739]. Curr Control Trials Cardiovasc Med 2003; 4(1): 2.
Malpas S C, Whiteside E A, Maling T J. Heart rate variability and cardiac autonomic function in men with chronic alcohol dependence. 1991; 65(2)(84): 8.
Zeskind PS, Gingras JL. Maternal cigarette-smoking during pregnancy disrupts rhythms in fetal heart rate. J Pediatr Psychol 2006; 31(1): 5-14.
Pan J, Tompkins WJ. A real-time QRS detection algorithm. IEEE Trans Biomed Eng 1985; 32(3): 230-6.
Saul JP, Albrecht P, Berger RD, Cohen RJ. Analysis of long term heart rate variability: methods, 1/f scaling and implications. Comput Cardiol 1988; 14: 419-22. [PMID: 11542156].
Malik M, Farrell T, Cripps T, Camm AJ. Heart rate variability in relation to prognosis after myocardial infarction: selection of optimal processing techniques. Eur Heart J 1989; 10(12): 1060-74.
Kay SM, Marple SL. Spectrum analysis—a modern perspective. Proc IEEE 1981; 69(11): 1380-419.
Scargle JD. Studies in astronomical time series analysis. ii-statistical aspects of spectral analysis of unevenly spaced data. Astrophys J 1982; 263: 835-53.
Lomb NR. Least-squares frequency analysis of unequally spaced data. Astrophys Space Sci 1976; 39(2): 447-62.
Mateo J, Laguna P. Analysis of heart rate variability in the presence of ectopic beats using the heart timing signal. IEEE Trans Biomed Eng 2003; 50(3): 334-43.
Li L, Li K, Liu C, Liu C-Y. Comparison of detrending methods in spectral analysis of heart rate variability. Res J Appl Sci Eng Technol 2011; 3(9): 1014-21.
Mateo J, Laguna P. Improved heart rate variability signal analysis from the beat occurrence times according to the IPFM model. IEEE Trans Biomed Eng 2000; 47(8): 985-96.
Laguna P, Moody GB, Mark RG. Power spectral density of unevenly sampled data by least-square analysis: performance and application to heart rate signals. IEEE Trans Biomed Eng 1998; 45(6): 698-715.
Katona PG, Jih F. Respiratory sinus arrhythmia: noninvasive measure of parasympathetic cardiac control. J Appl Physiol 1975; 39(5): 801-5.
Eckberg DL. Human sinus arrhythmia as an index of vagal cardiac outflow. J Appl Physiol 1983; 54(4): 961-6.
Berger RD, Akselrod S, Gordon D, Cohen RJ. An efficient algorithm for spectral analysis of heart rate variability. IEEE Trans Biomed Eng 1986; 33(9): 900-4.
Kobayashi M, Musha T. 1/f fluctuation of heartbeat period. IEEE Trans Biomed Eng 1982; 29(6): 456-7.
Yamamoto Y, Hughson RL. Coarse-graining spectral analysis: new method for studying heart rate variability. J Appl Physiol 1991; 71(3): 1143-50.
Perkiömäki JS, Bloch Thomsen PE, Kiviniemi AM, Messier MD, Huikuri HV. Risk factors of self-terminating and perpetuating ventricular tachyarrhythmias in post-infarction patients with moderately depressed left ventricular function, a CARISMA sub-analysis. Europace 2011; 13(11): 1604-11.
Gang UJ, Jøns C, Jørgensen RM, et al. Risk markers of late high-degree atrioventricular block in patients with left ventricular dysfunction after an acute myocardial infarction: a CARISMA substudy. Europace 2011; 13(10): 1471-7.
Cygankiewicz I, Zareba W, Vazquez R, et al. Relation of heart rate turbulence to severity of heart failure. Am J Cardiol 2006; 98(12): 1635-40.
Sassi R, Cerutti S, Lombardi F, et al. Advances in heart rate variability signal analysis: joint position statement by the e-Cardiology ESC Working Group and the European Heart Rhythm Association co-endorsed by the Asia Pacific Heart Rhythm Society. Europace 2015; 17(9): 1341-53.
Napolitano A. Generalizations of cyclostationarity: A new paradigm for signal processing for mobile communications, radar, and sonar. IEEE Signal Process Mag 2013; 30(6): 53-63.
Napolitano A. Cyclic statistic estimators with uncertain cycle frequencies. IEEE Trans Inf Theory 2017; 63(1): 649-75.
Antoni J. Cyclic spectral analysis in practice. Mech Syst Signal Process 2007; 21(2): 597-630.
Kazemi S, Ghorbani A, Amindavar H, Li C. Cyclostationary approach for heart and respiration rates monitoring with body movement cancellation using radar doppler system. 2013. Available from: http://arXiv preprint arXiv:1310.2293
Haritopoulos M, Capdessus C, Nandi AK. Foetal pqrst extraction from ECG recordings using cyclostationarity-based source separation method. Engineering in Medicine and Biology Society (EMBC) 2010; 1910-3.
Seydnejad S. Detection of nonlinearity in cardiovascular variability signals using cyclostationary analysis. Ann Biomed Eng 2007; 35(5): 744-54.
Gardner W. Measurement of spectral correlation. IEEE Trans Acoust Speech Signal Process 1986; 34(5): 1111-23.
Hurd HL. Nonparametric time series analysis for periodically correlated processes. IEEE Trans Inf Theory 1989; 35(2): 350-9.
Goldberger AL, Amaral LAN, Glass L, et al. PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals. Circulation 2000; 101(23): E215-20.
PhysioBank - clinical database Available from:
Malliani A, Pagani M, Lombardi F, Cerutti S. Cardiovascular neural regulation explored in the frequency domain. Circulation 1991; 84(2): 482-92.
Kleiger RE, Miller JP, Bigger JT Jr, Moss AJ. Decreased heart rate variability and its association with increased mortality after acute myocardial infarction. Am J Cardiol 1987; 59(4): 256-62. []. [PMID: 3812275].