Spectral analysis (SA) has been extensively applied to the assessment of heart rate variability. Traditional methods require the transformation of the original non-uniformly spaced electrocardiogram RR interval series into regularly spaced ones using interpolation or other approaches. The Lomb-Scargle (L-S) method uses the raw original RR series, avoiding different artifacts introduced by traditional SA methods, but it has been scarcely used in clinical settings. An RR series was recorded from 120 healthy participants (17–25 years) of both genders during a resting condition using four SA methods, including the Classic modified periodogram, the Welch procedure, the autoregressive Burg method and the L-S method. The efficient implementation of the L-S algorithm with the added possibility of the application of a set of options for the RR series pre-processing developed by Eleuteri et al., and also the results obtained in this study, show that the L-S method can be a good choice for future clinical studies. The L-S method seems particularly useful when the heart rates of studied participants will be below 60 or over 120 beats per minute. Nevertheless, it is important to the development of a smoothing procedure for the L-S spectra to avoid the picky behavior of the L-S power spectrum. The implementation of the L-S algorithm used in this study has been recently published by other authors included in our references, and brings some particular filtering features. The results obtained, comparing the four spectral methods, show that this implementation seems particularly useful when the heart rates of studied participants will be below 60 or over 120 beats per minute. Nevertheless, it is important to recommend for all existing L-S software implementations, the development of a smoothing procedure to avoid the picky behavior of the L-S power spectrum.
Spectral analysis (SA) of heart rate variability (HRV) has been widely applied for the assessment of cardiovascular control by the autonomic nervous system (ANS). A query to the Pubmed NLM database (07-13-2013) using the EndNote XV software of Thomson Reuters and the MeSH “Heart Rate Variability” returned a total of 15,892 items. From these, the SA method was included in 3288 studies (20.69%).
The introduction of SA in the assessment of HRV was acknowledged to Solange Akselrod et al.  by a task force  specifically created to establish standards of measurement, physiological interpretation and clinical use of HRV studies.
To be fair, since the early 1960s, Soviet scientists were using and publishing on the use of SA as a tool for HRV assessment [3–6], and in 1981 Baevskii et al.  developed a cybernetic model about the ANS regulation of the cardiovascular system. Later, Baevskii et al. emphasized this fact . Another Soviet scientist, I.G. Nidekker [6, 9] was, to our knowledge, the first to demonstrate that the use of raw RR intervals of the electrocardiogram (ECG) for the direct application of SA was inconsequential because the raw RR series were, in fact, ordinal sequences and not temporal ones. Hence, this author recommended the use of sliding time-windows and measuring the number of cardiac cycles or fractions of it for every one of these windows for the whole RR series. Using this procedure, it was possible to obtain a real temporal RR series with units expressed as cardiac cycles per unit of time that could then be studied by applying the fast Fourier transform (FFT).
The first author of the present paper worked in the Medico-Biological Spatial Intercosmos Program of the then USSR from 1978 to 1992 and was a personal witness to those studies [10–16]. We feel that the language barrier played an important role in the unawareness of Western countries about the pioneering contributions that those scientists played in the study of HRV. Notwithstanding, there is no doubt that the publication of S. Akselrod et al.  became a milestone in this field, as well as another paper published by R.D. Berger et al. .
It is not possible in such short space to mention the most relevant scientists in this field that made possible the improvement of the particular methods and the interpretation of the physiological meaning of the HRV spectral bands and indices, but perhaps it is enough to say that if in only 15 years after the publication of Akselrod, the Task Force  could establish standards and recommendations, the authors included in the references of that publication must be considered in that selected group. At that time (1996), almost 420 publications had reported the use of SA of the HRV.
In general, to estimate the power spectral density (PSD) of the HRV, the most used SA methods have been those using a non-parametric approach, including the classic windowed periodogram [18, 19] and a smoothed variant of this periodogram developed by Welch . Some other non-parametric spectral methods have been developed but not frequently used . Parametric methods, including the Burg autoregressive method and others [21, 22] have been also often used as an alternative.
A common practice for the previously mentioned spectral methods is the pre-processing of the interval series, including some of the following actions: the inspection and edition of the correct detection of the peaks of the consecutive R waves of the ECG; the subtraction of the mean value of the series to avoid the effect of the high energy values of the zero-frequency or DC frequency of the spectrum; the application of some detrending method to suppress drifts or very low frequency processes that can distort the RR time-series; the limitation of the frequencies with recognized physiological meaning with high pass, low pass or band pass filters; the interpolation of the original RR interval series; or the use of other particular procedures  to transform the non-uniformly spaced RR series into regularly spaced ones.
In recent years, interest in the use of methods proposing the usage of raw unevenly spaced RR series has been progressively growing [23–37]. Some of these methods have their sources in the “Cosinor” method  and others derive from the concepts of the Point Process Theory that were initially developed in the field of Astronomy and Astrophysics [39, 40].
Recently, A. Eleuteri et al.  published a method for detrending and bandwidth-limiting the HRV series without resampling and also developed and published, as open source and free downloadable, an implementation of their method and of a particularly fast implementation of the L-S algorithm, both written in MATLAB code. The L-S method for SA in human neonatal HRV  has been shown to be more useful than other traditional methods.
The aims of this work are to develop a comparative analysis between the application of the L-S spectral method using the approach developed by Eleuteri et al.  and those obtained using two non-parametric methods and a parametric autoregressive method in a group of 120 young healthy volunteers during a resting state and breathing spontaneously, and on the basis of the results to consider the convenience of the introduction of the L-S method as an alternative to the more conventional methods of spectral analysis for the RR series of the ECG.
The study included a group of 120 participants including students and young medical doctors or graduated nurses of one of the medical faculties of Havana who gave their written informed consent to participate in the study. Participants with history of cardiorespiratory or neurological disorders, diabetes mellitus or medications with known action on the ANS were excluded. All the volunteers were non-smokers, and their weight was no more or <10% from the accepted standards for each gender. The group was integrated with 60 female and 60 male participants paired by age with ages from 17 to 25 years. The protocol fulfilled the principles of the Declaration of Helsinki. The Ethical Committee of the Institute of Neurology and Neurosurgery approved the study.
All participants were studied in the period from 08:00 to 12:00 h. Volunteers were instructed to abstain of physical efforts during the day before the study, to avoid drinking beverages containing caffeine, to sleep for at least 7 h the night before, to have their usual breakfast and to drink a glass of fruit juice at least 1 h before the study. For 30 min, they had to rest sitting on an easy chair in a semi-recumbent position while the ECG electrodes were placed and the equipment setting and the quality of records were tested. The Laboratory room temperature was controlled to be between 23°C and 25°C with a noise level attenuated to comfortable levels. Only a specialized technician and one of us were present during the studies.
Electrocardiogram (ECG) recordings
Discoid silver electrodes were situated at the positions CM2 and V5 after rubbing and cleaning the skin under the electrodes with an alcoholic solution and later applying an electrolytic paste. An earth electrode was fixed in the left wrist. After 25 min, blood pressure was measured using the sphygmomanometric method. After 30 min at a resting condition, an ECG was continually recorded using a DC amplifier of a polygraphic device (Neurofax, Nihon-Kohden). The ECG signal was digitally monitored and stored in dedicated files on a PC after A/D conversion using a commercial card (12-bits, Advantech) and special tailored software developed with Matlab R2011a (The MathWorks, Inc.). The signal was recorded for 6 min.
RR measurements and editing
The digitally recorded ECG of every participant was submitted to visual inspection to ascertain the quality of the records using the software MultiTools, version 3.2.2, developed by us in Borland Delphi-7 code (Davihmed, Havana, Cuba). Non-physiological artifacts were marked and discarded, and participants with isolated ventricular ectopic beats or supraventricular events were not included in this study. We rejected seven participants, three for isolated ventricular ectopic beats and four for excessive movement artifacts affecting the quality of the ECG record. So, the study included 127 volunteer participants. Automatic detections of the R peaks were visually inspected to avoid false detections. Once the physiologist was satisfied with the R detections, the values of the sequential RR intervals were stored in ASCII files for further analysis.
Pre-processing of the RR series
HRV series to be processed with the classic, Welch and Burg periodogram methods were first submitted to a subtraction of their mean values and later were linearly detrended, computing a least-squares fit of a straight line sequence to the data and subtracting the resulting function from the RR series. Posteriorly, a Butterworth high pass filter of order 6 with a cutoff frequency of 0.02 Hz was applied to the series using a special zero-phase-shift digital filtering by processing the input data in both the forward and reverse directions. A low pass Butterworth filter of order 6 and cutoff frequency at 0.4 Hz was also applied to the series using the same process (zero-phase-shifting), and finally an interpolation process using the piece-wise cubic Hermite function was carried out. As the time duration to analyze the RR series was limited to 5 min (300 seconds) and for the SA 1024 samples, the sampling period for the interpolation was calculated as
The sampling frequency of the new RR interpolated series was
and the Nyquist frequency was
The series to be processed using the L-S method was pre-processed using the method proposed by Eleuteri et al. .
One series of RR data was reserved without any type of pre-processing to develop comparisons between the effects of the pre-processing actions when using the Welch or the L-S methods.
Calculation of the periodograms and spectral indices
As the L-S method uses different specific input parameters, differing from the other methods studied, we defined some criterion to achieve the best homology for the comparison of the four spectral methods. We decided then to use the real number “1.1” for the input value of the IOFAC (oversampling factor) for the L-S method. Other input values for the L-S algorithm were the default ones defined by the implementation of Eleuteri et al. . In Table 1, the number of discrete frequencies included for the L-S method and for the classic, Welch and Burg methods are shown. The spectral resolution of the L-S method was 0.0030 Hz and 0.0033 Hz for the other three methods. This was the best approximation we were able to achieve.
|Spectral bands||Spectral range, Hz||Lomb-scargle method (IOFAC=1.1)||Traditional methods||Lomb-scargle method (IOFAC=4)|
We limited the spectral range studied from 0.02 to 0.40 Hz, and the spectral band limits and the acronyms used were those recommended by the Task Force , excepting the corresponding to the mid-frequency sub-band (MF). Nevertheless, the complete range of the LF band was also considered, but using the acronym LFAll (Table 1). For the calculations of the classic and Welch periodograms, we applied a Hann window to the input RR series. The Burg method does not use data windowing, but it is necessary to assign the order of the autoregressive model to apply. There are different criteria about the best order to use, but recently Dantas and colleagues  have shown that it can be used any value between 9 and 25. As other authors have found that the use of higher values for the autoregressive methods are better to define the different peaks of the spectra , we selected order 25 for the Burg method.
To assess the effect of using a value of “1.1” for the IOFAC, we also carried out the same process for each participant, but using the value “4” for this input value. The number of spectral discrete frequencies with an IOFAC of “4” is shown in Table 1.
The estimation of the PSD using the periodogram of the RR series was one of the standards used by Matlab, and the mathematical expression is presented below:
where f, corresponding discrete spectral frequency; Fs, sampling frequency; N, number of samples in the series; and n, the number of order of the sample in the RR series.
The normalized values (%) of the spectral bands were calculated using the following expression:
where TotalPSD corresponded to the value of the PSD in the spectral range of 0.02 to 0.40 Hz, and AbsVal (Band) was the PSD calculated for the different bands.
For the calculation of the mean frequency for every one of the four bands, the following expression was used.
where m, ordinal number of the discrete spectral frequency at the inferior limit of the band; n, the ordinal number at the superior limit of the band; fk, value of the corresponding discrete spectral frequency expressed in Hz; Pk, PSD for that k discrete spectral frequency (DFi); and TotalPSD (Band), the value of the absolute PSD for that band. The results were expressed in Hz. The LF/HF ratio was calculated as the ratio between the PSD in the frequency range from 0.04 to 0.15 by the PSD in the range from 0.15 to 0.40 Hz.
Grand averages of the 120 spectra were calculated for each method , and the standard error of the mean was obtained for every DFi in the spectral range from 0.02 Hz to 0.40 Hz. The expression for these calculations can be expressed as follows:
where M is the total number of spectra (120), Pk is the k-th individual spectrum, and f represents all the DFi corresponding to each k-th spectrum in the studied spectral range. Graphics for every grand average for each method were obtained for visual comparative analysis.
Analysis of indices of the HRV in the time domain
HRV indices in the time domain were calculated as detailed elsewhere [44, 45]. Standard recommended HRV indices were considered: the mean RRI value (M), the SD and the square root of the mean squared differences of successive RRI intervals RMSSD. The mode (Mo), the median (Me), the amplitude of the mode (AMo), the delta range, and the stress index were calculated by the expression:
with the triangular index calculated with a bin of 5 milliseconds and not as originally described with a bin of 7.8125 [2, 46] and, finally, the index pMean2% calculated by the proportion of the RR intervals larger than the precedent by 2% of the mean RR value for the whole RR sequence.
All the calculations in this study were carried out using routines written in Matlab.
For each quantitative index in this study, the Kolmogorov-Smirnov test was applied to verify the normality of its distribution and accepted if the results were p>0.05. When necessary, data were transformed until normality was achieved and properly showed in tables and figures.
Comparisons between the indices in the time domain obtained for males and females were compared using Student’s t-test for independent samples. To compare the values of the spectral indices obtained with the four methods, we used criteria based on the general linear model assumption. An ANOVA test design comprising four repeated measures (within factor “Method”) and two between factors (“Age” and “Gender”) was used for all the spectral indices. Three subgroups of age were considered for the between factor age: 17 to 19.99; 20 to 21.99; and 22 to 25. As the participants were paired by age and gender, the three subgroups were identical in number. To avoid assumptions of sphericity and compound symmetry, we used the criteria of the multivariate tests of Wilk’s Lambda, Pillai-Bartlett Trace, Hotelling-Lawley Trace and Roy’s larger root [47, 48]. We checked the values obtained for the between factors “Age”, “Gender” and their interactions, and in case of finding significant values for the results of the F test, and at least two of the before mentioned multivariate tests, with a statistical power over 0.9 for an alpha error of p<0.05, the results were accepted as significant. If these conditions were achieved, we considered justified the use of the post hoc procedure of Scheffé .
For comparisons between the effects of pre-processing (or not), the RR series between the results using the L-S and the Welch methods, the RR series were not modified. For the estimations of the PSD using the L-S method, the pre-processing was the one recommended by Eleuteri et al. , and for the Welch method, all those described before in the “Pre-Processing of the RR Series” section of this paper. Comparisons of the spectral indices were carried out using Student’s t-test for dependent samples.
The results obtained using different values for the IOFAC factor (1.1 or 4) with the L-S method were assessed statistically using also Student’s t-test for dependent samples.
Data obtained from the analysis were tabulated, and graphics were created using the statistical package STATISTICA (data analysis software system, version 8.0. StatSoft, Inc.).
Male and female volunteers did not show significant differences for age, systolic and diastolic blood pressure or any indices of HRV calculated in the time domain as shown, but the body mass index (BMI) was significantly higher in male participants (Table 2).
|#||Indices||Values all group n=120||Females n=60||Males n=60||pa|
|1||Decimal age, years||21.05±2.23||20.82±2.2||21.29±2.3||0.265|
|4||Body mass index, Kg/m2||–||20.4±1.3||21.8±1.5||0.001|
|8||Delta range, ms||344.20±99.8||341.32±111.2||347.10±87.7||0.752|
aProbability values associated with comparisons between male and female participants using Student’s t-test for independent samples. bBlood pressure in mmHg.
The ANOVA test results for the values of the absolute PSD are shown in Table 3. Significant differences were not found between the different indices for the factor “gender”. The factor “age” showed significant differences for the index MF (F(2, 114)=3.4272, p=0.03586) with higher values for the group 20–22.99 years of age, and HF (F(2, 114)=3.2354, p=0.04298), with lower values for the group 23–25 years of age vs. the other two subgroups, but the multivariate analysis did not find any significant value for the two between factors, and therefore those values may not be considered statistically significant. The differences obtained from the post hoc Scheffé test are indicated with symbols in Table 3. In general, the PSD values for the L-S and the Burg methods resulted in values lower than those of the classic and Welch methods. There were only significant differences between the L-S method and the Burg method for the VLF absolute values. There were not found significant differences between the classic and the Welch methods.
|Indices||L-S method||Welch method||Classic method||Burg method|
|Ln(LFAll), ms2||11.733±0.42||12.609±0.73a||12.586±0.77a||11.664±0.76b, c|
ap<0.001 [comparisons vs. Lomb-Scargle (L-S) method]. bp<0.001 (comparisons of classic and Burg methods vs. Welch). cp<0.001 (comparisons between Burg and the classic method).
The morphology of the grand averages obtained for the classic method and the L-S method showed a great resemblance (Figure 1), although the values of the PSD were significantly higher for the classic method (Table 3). The Welch grand average showed a great similitude to the L-S and the classic grand averages, but with the smoothing effect that is inherent to this method.
However, the morphology of the grand average obtained for the Burg method (Figure 1) did not show the peak for the VLF band and emphasized the amplitude of the peak in the frequencies between 0.17 and 0.25 Hz.
The ANOVA test showed only significant differences for the within factor “method” (repeated measures). The between factors “gender”, “age” or their interactions did not show statistical significance, and their values of statistical power did not reach the limit of 0.9. In Figure 2, the results of the differences between the four methods are graphically showed.
The normalized values for the VLF band were higher than those obtained with the conventional methods, whereas there were not differences between the Welch and the classic method. The Burg method showed the lower values compared with the other three methods.
The MF band expressed in normalized values did not show significant differences between any of the methods. However, the LF and LFAll bands showed a similar pattern in their statistical behavior, with lower significant values for the classic and Welch methods vs. the L-S method, and higher values vs. the other three methods.
The normalized values of the HF band were significantly lower than those of the three traditional methods, while no differences were observed between them. The LF/HF ratio showed lower values for the Welch and classic methods, but there were no differences between the L-S and the Burg methods.
The ANOVA test for the indices of the mean frequencies of the considered bands also did not show statistical differences for the between factors “age” or “gender” or their interactions. The detailed results are shown in Figure 3.
The mean frequency of the VLF band (mf VLF) did not differ between the L-S and the Welch method and was significantly higher with the classic method vs. all the other methods, and the Burg method showed the lowest values of the four methods.
The mean frequency of the LF band was significantly lower for the Welch and classic methods vs. the L-S method, but the values for the Burg method were higher than any other method. The mean frequency of the MF band was not different between the L-S, the Welch or the classic method, but the values for the Burg method showed higher significant values than the other three.
For the mean frequency of the HF band, the values obtained for the three traditional methods were higher than those of the L-S method, and particularly higher were the values for the Burg method, which differed from those of the L-S and the classic methods.
In Figure 4, the differences observed for particularly important spectral indices that can reflect the effect of the pre-processing that was used for the L-S method as developed by Eleuteri et al.  and the actions that were described for the Welch method are shown.
For the normalized values (nu) of the VLF band, the values were significantly lower when the procedures of pre-processing (PP) were carried out for both the L-S and the Welch methods. The same results were obtained for the nu of the HF band. Nevertheless, the magnitude of the reduction of the high frequency band was higher for the use of the PP with the L-S method.
The mean frequencies of the VLF band were significantly incremented when using PP with the two methods. On the contrary, the mean frequencies of the HF band were reduced significantly when using PP.
The PP procedure produced a significant increment of the LF/HF ratio for both methods, but the magnitude was clearly higher for the L-S method. The comparisons of the effects of using a value of 4 or 1.1 for the IOFAC for the L-S method are shown in Table 4.
|Indices||Lomb-scargle (IOFAC=1.1)||Lomb-scargle (IOFAC=4)|
ap<0.01 probability associated with Student’s t-test for dependent samples. bp<0.05 idem.
All the absolute PSD values for the studied bands were significantly higher with IOFAC=4, but the normalized values (%) were lower for the VLF band and the LF band. The MF band was higher. Considering the spectral range from 0.04 to 0.15 (VLFAll), there were not found significant differences nor for the for the HF band and for the mean frequencies in the four studied bands.
In this study, we did not find any statistically significant differences between male and female volunteers for HRV indices calculated in the time domain, although indices of mean, median and mode were slightly higher in females (lesser values of heart rate). This can be considered as a physiological tendency  related to the gender factor. Differences in BMI were also related to gender.
Grand averages of the spectral methods
Morphology and variability of the grand averages of spectra obtained for the 120 participants showed one possible disadvantage of the L-S method that is also attributed to the classic modified method, its “picky” appearance. The simple visual analysis showed a great similarity of the grand averages obtained for the classic method and the L-S method.
The general variability of the four methods for every discrete spectral frequency of the four grand averages was not visually very different, as shown by the standard error of the means. The Welch method, although in a smoothed way, clearly depicted the main peaks detected by the classic and the L-S method.
The Burg autoregressive method failed in the detection of the VLF peak, which is also a limitation for this method, but despite this limitation the two peaks corresponding to the LF band, as considered in our study from 0.04 to 0.085 Hz and for the MF sub-band (0.085 to 0.15 Hz), were clearly defined by this method.
The approach of averaging spectra of different individuals to produce a grand average for a group of participants was described by Piskorski et al. . We have used this approach in a previous paper  to visually differentiate the spectral characteristics of comatose patients with different Glasgow coma scale scores, and now in this paper we show its utility in detecting the morphology and the main frequencies present in each of the studied spectral bands or sub-bands of the HRV in the range from 0.02 to 0.04 Hz.
It was also shown that the range of averaged amplitudes were different for the methods, with higher mean values for the classic and Welch methods and lesser and approximately similar mean averaged range of amplitudes for the Burg and the L-S methods.
It is important to emphasize that this approach is not the same as the averaging process used for the calculations of the Welch procedure  or intents of producing the same effect for the L-S method . So, this analytic approach must be limited to visual analysis because the statistical properties of the HRV spectra grand averages have not been fully developed until this moment, although the procedure seems very promising for comparisons of groups or subgroups of participants or between physiological or testing conditions.
ANOVA test results
The results of the application of the ANOVA test for the absolute PSD values corresponds very well with the visual analytic approach discussed previously but to our criteria they do not provide much value for our objectives. The statistical differences, shown in Table 3 can only be the result of the inherent mathematics of each method for estimating the PSD.
Nevertheless, the standardizations used for the periodograms of the three traditional methods were those recommended by Matlab , and in the L-S method we adjusted our results to the recommendations of Eleuteri et al. . Of course, these results only have a formal value in the analysis, but they show differences that must be considered, and hence, we show them and their statistical significance.
No significant differences were found, as we had expected, between the values of the normalized PSD estimations obtained with the Welch and classic methods (Figure 2), although the nuVLF, nuLF and nuLFAll estimations with the L-S method were higher than the other three methods, and the nuHF was statistically lower. In the spectral range from 0.085 to 0.15 Hz. (nuMF), no statistical differences were found between the four methods. The LF/HF ratio was found higher for the L-S method.
As differences between the classic and Welch methods with the Burg method are known and have been recognized by other authors [53–56], we are not going to focus on this.
Although statistically significant, the previously described differences were not so striking as to lead us to consider that they could modify the main knowledge accumulated using the traditional methods. For example, the mean increment with the L-S method for the nuVLF was from 2%–3% higher than the classic and Welch methods, for the nuLF was from 3%–4% and for the nuLFAll was from 2%–3% incrementally.
Nevertheless, the importance of a fair estimation of the VLF band is very important in many specific conditions in neurological patients with deep disturbances of consciousness [51, 57–60], and we consider it better to overestimate the VLF presence than the contrary for those patients and conditions. This fact may be considered an advantage for the L-S method.
From another point of view, the lesser values of the nuHF band of the L-S method, using the same arguments, could be a disadvantage because it is precisely the HF band that is the first to suffer the effects of different noxious processes, particularly those associated with conscience disturbances.
The LF/HF ratio was higher with the L-S method, as should had been expected after the previous results, but the classic, Welch and Burg methods also showed an increment of the HRV low frequencies. The Burg method did not show significant differences from the L-S method. The basic finding for healthy participants in resting conditions, consisting of a slight increment of the low frequencies over the high respiratory and vagal frequencies, was correctly detected by the four methods (Figure 2).
The mean frequency indices of the spectral bands (Figure 3) are not very frequently used in research in HRV analysis, but we thought that they could help us in the whole analysis. It called our attention to the fact that for the mean frequency of the VLF band, the results obtained with the L-S and Welch methods did not differ. So, the peak of the band was similarly detected quantitatively.
Now, focusing on the other extreme of the spectrum of HRV, the mean frequency of the HF band was found at lesser value displacements to the left of the energy of the spectral band with the L-S method and at higher values for the traditional methods. The number of spectral discrete frequencies used for the calculations were not very different between the four methods, as shown in Table 1 (83 for the L-S and 74 for the other methods). The possible origin for those differences should be related to the methods themselves and the pre-processing actions to the RR series.
We will now analyze the effect of overestimation of the high frequencies of the HRV spectrum produced by the effect of oversampling the original RR series used by the traditional procedures . The oversampling to obtain uniformly spaced RR intervals behaves itself as a low-pass filter and produces the undesired phenomenon of aliasing in the high frequencies.
In the method proposed by Berger et al.  that we have used for many years, there is, in part, a way to avoid that effect. However, when using the interpolation procedures, the only ways to attenuate this inconvenience are the use of a procedure for resampling (as we consider the one used in our study a piece-wise cubic Hermite interpolation), the application of a low-pass filtering to the RR series before applying the process for estimation the PSD and the use as resampling frequency of a value at least double the mean sampling frequency.
We used a resampling frequency for the traditional methods near 3.4 Hz for a sampling period of approximately 0.293 seconds. Some authors have stated  that for heart rate signals, the mean sampling frequency can be fixed around 2 Hz. In our study, the resampling frequency is not the double of this value (3.41333/2=1.71), but is close to it. Nevertheless, we did not limit the higher limit of the spectral range studied to 0.5, but to 0.4 Hz, and then we could easily calculate that this value is more than four times below the Nyquist frequency (See the “Pre-Processing of the RR Series” section) of our resampling process, as shown in the next calculation.
The L-S method does not use resampling but also introduces some high-frequency contamination, as shown by the same authors . Schaffer et al.  found that the L-S method had limitations for resolving the frequencies in the HF band.
It is then very difficult, or even impossible, to decide which method performs better for detecting the HF band because as stated by Laguna et al. , real heart rate signals are not adequate to estimate experimentally the performance of the L-S PSD estimate because we do not know the real spectra that must be obtained. Therefore, the discussion about whether the observed higher values of the HF band obtained with the traditional methods are a consequence of artifacts or not is futile.
Analysis of the pre-processing actions
To accomplish our aims, it was convenient to establish a comparative analysis of the effects of the PP actions applied for the traditional method of Welch and the L-S method. Those results are presented in Figure 4 and are helpful in the analysis of the problem in question.
For nuVLF values, the effects of the PP actions produced for both methods significant proportional effects: a reduction of the leakage and aliasing dependent on different known factors [18, 19, 21]. So, the set of PP actions can be considered appropriate for both methods.
For nuHF values, although when the PP actions were applied a reduction of the possible leakage in the high frequencies for both methods occurred, the magnitude of the reduction was more evident for the L-S than with the Welch method.
The mean frequencies in both extremes of the studied HRV spectra were found significant and proportional significant effects were seen with both methods (L-S and Welch), consisting of higher values of the mfVLF on account of the effect of the PP actions, and the high frequencies of the nuHF values were significantly reduced and the magnitude of the reductions was higher in both cases. When observing the results of the effects of the PP actions on the LF/HF ratio, although for both methods the expected increment was obtained, the magnitude was obviously more intense for the L-S method.
So, it seems that the PP actions developed by Eleuteri et al.  are obviously effective, and their action in the low frequencies is of the same magnitude as those used in this study for the traditional Welch method. The results are not so clear as to adopt a conclusion with regard to the high frequencies, but there is no doubt that the effects of the PP actions for the L-S method were clearly more intense.
The default value recommended by Eleuteri et al.  for the IOFAC input factor for the algorithm of the L-S method is 4, but in this study, to homologate the conditions for the comparisons of the different methods, we used 1.1 instead. The comparisons of the results (Table 4) confirmed that the value introduced did not produce significant bias to the general results analyzed in this study.
The use of the L-S method in neonates  results is almost mandatory because heart rates are commonly very high, frequently over 120 beats per minute (bpm), and the respiratory frequency is also very high compared with that of adults. Resampling such RR series requires the use of sampling frequencies that introduce significant artifacts into the results of the estimations of PSD. When applying the traditional methods to participants with heart rates with values under 60 bpm, the same problem is present, but the use of the L-S method can be used without evidently introducing artifacts, as showed by Laguna et al.  in their paper.
Recently the L-S method has been promoted as an alternative for the study of very short RR series of 30 beats [36, 37]. This is a questionable but very interesting fact, because it is the situation we encounter during physiological maneuvers such as respiratory gasps, yawning and the Valsalva test in the first 15 seconds after the forced expiration. We will not focus on these aspects, but they are items that represent benefits of the utility of the L-S method.
To our knowledge, this is the first study applying the L-S algorithm using the PP procedure developed by Eleuteri et al.  in a cohort of 120 healthy young participants of both genders in a resting state and breathing spontaneously.
The use of the L-S method in studies of SA of HRV is still very low (approximately 0.6% of papers applying SA in the Pubmed NLM database). The efficient implementation of the algorithm by Eleuteri et al. , with the possibility of the application of a set of options for the RR series pre-processing and also the results obtained in the assessment of HRV indices of the frequency domain, carried out in this study shows that the L-S method can be a good choice for future studies. The L-S method seems particularly useful in those conditions where the heart rates of studied participants will be below 60 or over 120 beats per minute. Nevertheless, the development of a smoothing procedure would be important, such as that used with the Welch method, to avoid the picky behavior of the L-S power spectrum.
The authors would like to thank A. Eleuteri et al.  for the fast implementation of the L-S algorithm and for the pre-processing procedure developed by these authors, developed in Matlab code and generously presented as an open-source, free download from http://clinengnhs.liv.ac.uk/links.htm. G.L. Is supported by the Kamea Dor-Bet program of the State of Israel and by the Children’s Autism Hope Project. None of the authors have any conflicts of interests to declare.
Conflict of interest statement
Author contributions: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.
Research funding: This project was funded in part by the Government of Israel under the Kamea-Dor-Bet program and by the Children’s Autism Hope Project.
Employment or leadership: None declared.
Honorarium: None declared.
Competing interests: The funding organization(s) played no role in the study design; in the collection, analysis, and interpretation of data; in the writing of the report; or in the decision to submit the report for publication.
1. Akselrod S, Gordon D, Ubel FA, Shannon DC, Berger AC, Cohen RJ. Power spectrum analysis of heart rate fluctuation: a quantitative probe of beat-to-beat cardiovascular control. Science 1981;213:220–2.10.1126/science.6166045Search in Google Scholar
2. Malik M, Bigger JT, Cam AJ, Kleiger RE, Malliani A, Moss AJ, et al. Task force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Heart rate variability. Standards of measurement, physiological interpretation, and clinical use. Eur Heart J 1996;17:354–81.10.1093/oxfordjournals.eurheartj.a014868Search in Google Scholar
3. Parin VV, Baevskii RM, Gazenko OG. Achievements and problems of modern space cardiology. Kardiologiia 1965;22:3–11.Search in Google Scholar
4. Parin VV, Baevskii RM, Gazenko OG. Heart and circulation under space conditions. Cor Vasa 1965;32:165–84.Search in Google Scholar
5. Zhemaitite D. Statistical analysis of the sinusal activity in normal and pathological states. In: Methods of mathematical analysis of the heart rhythm. Nauka: Moscow, Russia, 1967.Search in Google Scholar
6. Baeovskii RM, Zamotaev IP, Nidekker IG. Mathematical analysis of sinus automatism for prognosis of rhythm disorders. Kardiologiia 1971;11:65–8.Search in Google Scholar
7. Baevskii RM, Barsukova ZhV, Tazetdinov IG. Cybernetic analysis of the heart rhythm in the measured physical loading test with crew members of the Saliut-6 orbital station. Kardiologiia 1981;21:100–4.Search in Google Scholar
8. Baevskii RM. Analysis of heart rate variability in space medicine. Fiziol Cheloveka 2002;28:70–82.Search in Google Scholar
9. Nidekker IG. Method of spectral analysis for long-term recordings of physiological curves. Kosm Biol Aviakosm Med 1981;15:78–82.Search in Google Scholar
10. Estévez Báez M, Casanova-Sotolongo P, Peñalver JC, Zayas-López F. El electroencefalograma en la evaluación del estado funcional del piloto de caza después de los vuelos. Revista de Medicina Militar 1983;2:95–104.Search in Google Scholar
11. Estévez Báez M, Peñalver JC, Trápaga Ortega M, Cintas González E. Estudio de la excitabilidad cortical del cerebro del hombre en condiciones de hipodinamia y antiortostasis. Boletín Academia de Ciencias de Cuba Suppl Especial 1983;141–3.Search in Google Scholar
12. Estévez Báez M, Asyamolova NM, Brodyetskaya L, Rodríguez F, Peñalver JC, Grachev VA, et al. Investigación de la actividad bioeléctrica cerebral de los cosmonautas en estado de impesantez. Experimento ‘Córtex’. Orbita Extraordinario, 1985;54–76.Search in Google Scholar
13. Estévez Báez M, Matveev AD, Kornilova LN, Tarasov IK, Borunov NL, Peñalver JC, et al. Estudio comparativo de la enfermedad del movimiento y la estabilidad vestíbulo-vegetativa en condiciones terrestres y del vuelo cósmico. Experimento ‘Encuesta,’ órbita Extraordinario 1985;82–91.Search in Google Scholar
14. Gazenko OG, Papenfuss W, Hideg J, Estévez Báez M. Results of the medical experiments during the spatial flights of the international crews. Rev Med Lotnicza 1988;1:6–14.Search in Google Scholar
15. Estévez Báez M. A medicina espacial cubana e o sistema neuromega. Saude para Todos 1993;20–6.Search in Google Scholar
16. Estévez-Báez M, Iglesias-Alfonso J, Villar-Olivera C, Cabana-González J, Fernández-Pérez L, Pujol-García J. El sistema neuromega en la evaluación de la influencia del estrés. Rev Cubana de Med Militar 1994;23:42–55.Search in Google Scholar
17. Berger RD, Akselrod S, Gordon D, Cohen RJ. An efficient algorithm for spectral analysis of heart rate variability. IEEE Trans Biomed Eng 1986;33:900–4.10.1109/TBME.1986.325789Search in Google Scholar
18. Oppenheim AV, Schafer RW. Discrete-time signal processing. Upper Saddle River, New Jersey: Prentice-Hall, 1999.Search in Google Scholar
19. Marple SL. Digital spectral analysis. Engelwood Cliffs, NJ: Prentice Hall, 1987.Search in Google Scholar
20. Welch PD. The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Trans Audio Electroacoustics 1967;AU-15:70–3.10.1109/TAU.1967.1161901Search in Google Scholar
21. The MathWorks. MATLAB signal processing toolbox – help. Availble at: http://www.mathworks.com/help/releases/R2011a/pdf_doc/allpdf.html#signal. Accessed 4 February, 2014.Search in Google Scholar
22. Semmlow JL. Biosignal and biomedical image processing: MATLAB-based applications. New York: Marcel Dekker, Inc., 2014.Search in Google Scholar
23. 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:698–715.10.1109/10.678605Search in Google Scholar
24. Rudiger H, Klinghammer L, Scheuch K. The trigonometric regressive spectral analysis--a method for mapping of beat-to-beat recorded cardiovascular parameters on to frequency domain in comparison with Fourier transformation. Comp Meth Prog Biomed 1999;58:1–15.10.1016/S0169-2607(98)00070-4Search in Google Scholar
25. Clifford GD, Tarassenko L. Segmenting cardiac-related data using sleep stages increases separation between normal subjects and apnoeic patients. Physiol Meas 2004;25:N27–35.10.1088/0967-3334/25/6/N03Search in Google Scholar
26. Laude D, Elghozi JL, Girard A, Bellard E, Bouhaddi M, Castiglioni P, et al. Comparison of various techniques used to estimate spontaneous baroreflex sensitivity (the EuroBaVar study). Am J Physiol Regul Integr Comp Physiol 2004;286:R226–31.10.1152/ajpregu.00709.2002Search in Google Scholar PubMed
27. Rudiger H, Seibt R, Scheuch K, Krause M, Alam S. Sympathetic and parasympathetic activation in heart rate variability in male hypertensive patients under mental stress. J Hum Hypertens 2004;18:307–15.10.1038/sj.jhh.1001671Search in Google Scholar PubMed
28. Barbieri R, Matten EC, Abdul Rasheed AA, Brown EM. A point-process model of human heartbeat intervals: new definitions of heart rate and heart rate variability. Am J Physiol Heart Circ Physiol 2005;288:H424–35.10.1152/ajpheart.00482.2003Search in Google Scholar PubMed
29. Clifford GD, Tarassenko L. Quantifying errors in spectral estimates of HRV due to beat replacement and resampling. IEEE Trans Biomed Eng 2005;52:630–8.10.1109/TBME.2005.844028Search in Google Scholar PubMed
30. Holland A, Aboy M. A novel recursive Fourier transform for nonuniform sampled signals: application to heart rate variability spectrum estimation. Med Biol Eng Comput 2009;47:697–707.10.1007/s11517-009-0461-0Search in Google Scholar PubMed
31. Krause M, Rudiger H, Bald M, Nake A, Paditz E. Autonomic blood pressure control in children and adolescents with type 1 diabetes mellitus. Pediatric Diabetes 2009;10:255–63.10.1111/j.1399-5448.2008.00447.xSearch in Google Scholar PubMed
32. Reimann M, Friedrich C, Gasch J, Reichmann H, Rudiger H, Ziemssen T. Trigonometric regressive spectral analysis reliably maps dynamic changes in baroreflex sensitivity and autonomic tone: the effect of gender and age. PloS One 2010;5:e12187.10.1371/journal.pone.0012187Search in Google Scholar PubMed PubMed Central
33. Tank J, Baevsky RM, Funtova II, Diedrich A, Slepchenkova IN, Jordan J. Orthostatic heart rate responses after prolonged space flights. Clin Auton Res 2011;21:121–4.10.1007/s10286-010-0106-2Search in Google Scholar PubMed
34. Schaffer T, Hensel B, Weigand C, Schuttler J, Jeleazcov C. Evaluation of techniques for estimating the power spectral density of RR-intervals under paced respiration conditions. J Clin Monit Comput 2014;28:481–6.10.1007/s10877-013-9447-4Search in Google Scholar PubMed
35. Smith AL, Owen H, Reynolds KJ. Heart rate variability indices for very shorterm (30 beat) analysis. Part 1: survey and toolbox. J Clin Monit Comput 2013;27:569–76.10.1007/s10877-013-9471-4Search in Google Scholar PubMed
36. Smith AL, Owen H, Reynolds KJ. Heart rate variability indices for very short-term (30 beat) analysis. Part 2: validation. J Clin Monit Comput 2013;27:577–85.10.1007/s10877-013-9473-2Search in Google Scholar PubMed
40. Eleuteri A, Fisher AC, Groves D, Dewhurst CJ. An efficient time-varying filter for detrending and bandwidth limiting the heart rate variability tachogram without resampling: MATLAB open-source code and Internet web-based implementation. Comput Math Methods Med 2012:578–785.10.1155/2012/578785Search in Google Scholar PubMed PubMed Central
41. Chang KI, Monahan KJ, Griffin MP, Lake D, Moorman JR. Comparison and clinical application of frequency domain methods in analysis of neonatal heart rate time series. Ann Biomed Eng 2001;29:764–74.10.1114/1.1397791Search in Google Scholar PubMed
42. Dantas EM, Sant’anna ML, Varejao AR, Goncalves CP, Morra EA, Baldo MP, et al. Spectral analysis of heart rate variability with the autoregressive method: What model order to choose? Comput Biol Med, 2012;42:164–70.10.1016/j.compbiomed.2011.11.004Search in Google Scholar PubMed
43. Kuusela TA, Kaila TJ, Kahonen M. Fine structure of the low-frequency spectra of heart rate and blood pressure. BMC Physiol 2003;13:1–11.Search in Google Scholar
44. Piskorski J, Guzik P, Krauze T, Zurek S. Cardiopulmonary resonance at 0.1 Hz demonstrated by averaged Lomb-Scargle periodogram. Cent Eur J Phys 2010;8:386–92.Search in Google Scholar
45. Montes-Brown J, Machado A, Estevez M, Carricarte C, Velazquez-Perez L. Autonomic dysfunction in presymptomatic spinocerebellar ataxia type-2. Acta Neurol Scand 2012;125: 24–9.10.1111/j.1600-0404.2011.01494.xSearch in Google Scholar PubMed
46. Montes-Brown J, Sanchez-Cruz G, Garcia AM, Baez ME, Velazquez-Perez L. Heart rate variability in type 2 spinocerebellar ataxia. Acta Neurol Scand 2010;122:329–35.Search in Google Scholar
47. 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:1060–74.10.1093/oxfordjournals.eurheartj.a059428Search in Google Scholar PubMed
48. Tatsuoka MM. Multivariate analysis. New York, NY: Wily, 1971.Search in Google Scholar
49. Finn JD. A general model for multivariate analysis. New York: Holt, Rinehart & Winston, 1974.Search in Google Scholar
50. Scheffé H. A method for judging all possible contrasts in the analysis of variance. Biometrica 1953;40:87–104.Search in Google Scholar
51. Abhishekh HA, Nisarga P, Kisan R, Meghana A, Chandran S, Trichur R, Sathyaprabha TN. Influence of age and gender on autonomic regulation of heart. J Clin Monit Comput 2013;27:259–64.10.1007/s10877-012-9424-3Search in Google Scholar PubMed
52. Machado-Ferrer Y, Estevez M, Machado C, Hernandez-Cruz A, Carrick FR, Leisman G, et al. Heart rate variability for assessing comatose patients with different Glasgow Coma Scale scores. Clin Neurophysiol 2013;124:589–97.10.1016/j.clinph.2012.09.008Search in Google Scholar PubMed
54. Badilini F, Maison-Blanche P, Coumel P. Heart rate variability in passive tilt test: comparative evaluation of autoregressive and FFT spectral analyses. Pacing Clin Electrophysiol 1998;21:1122–32.10.1111/j.1540-8159.1998.tb00159.xSearch in Google Scholar PubMed
56. Pichon A, Roulaud M, Antoine-Jonville S, de Bisschop C, Denjean A. Spectral analysis of heart rate variability: interchangeability between autoregressive analysis and fast Fourier transform. J Electrocardiol 2006;39:31–7.10.1016/j.jelectrocard.2005.08.001Search in Google Scholar PubMed
57. Ryan ML, Ogilvie MP, Pereira BMT, Gomez-Rodriguez JC, Manning RJ, Vargas PA, et al. Heart rate variability is an independent predictor of morbidity and mortality in hemodynamically stable trauma patients. J Trauma Inj Infect Crit Care 2011;70:1371–9.10.1097/TA.0b013e31821858e6Search in Google Scholar PubMed
58. Machado C, Estévez M, Pérez-Nellar J, Gutiérrez J, Rodríguez R, Carballo M, et al. Autonomic, EEG, and behavioral arousal signs in a PVS case after zolpidem intake. Can J Neurol Sci 2011;38:341–44.10.1017/S0317167100011562Search in Google Scholar PubMed
59. Riganello F, Dolce G, Sannita WG. Heart rate variability and the central autonomic network in the severe disorder of consciousness. J Rehab Med 2012;44:495–501.10.2340/16501977-0975Search in Google Scholar PubMed
60. Machado C, Estévez M, Rodríguez R, Pérez-Nellar J, Chinchilla M, Defina P, et al. Zolpidem arousing effect in persistent vegetative state patients: Autonomic, EEG and behavioral assessment. Curr Pharm Des, 2014;20:4185–202.Search in Google Scholar
©2016 by De Gruyter