Online Algorithm for Deriving Heart Rate Variability Components and Their Time–Frequency Analysis

Abstract

Featured Application

This work makes possible to extract heart rate variability components online in order to monitor the underlying human body systems, in particular to determine the activity of the sympathetic and parasympathetic branches of the autonomic nervous system (ANS) as well as the balance between them, and indirectly in detecting and monitoring many common diseases related to the cardiovascular system and ANS. Such a solution can, for example, be directly embedded into Holter devices. A more precise determination of the components’ properties can give the opportunity to link them to specific physiological processes, especially those of very low and ultra-low frequencies, which has not yet been fully achieved, increasing the practical importance of this research.

Abstract

Heart rate variability (HRV) containing four components of high (HF), low (LF), very low (VLF), and ultra-low (ULF) frequencies provides insight into the cardiovascular and autonomic nervous system functions. Classical spectral analysis is most often used in research on HRV and its components. The aim of this work was to develop and validate an online HRV decomposition algorithm for monitoring the associated physiological processes. The online algorithm was developed based on variational mode decomposition (VMD), validated on synthetic HRV with known properties and compared with its offline adaptive version AVMD, standard VMD, continuous wavelet transform (CWT), and wavelet package decomposition (WPD). Finally, it was used to decompose 36 real all-night HRVs from two datasets to analyze the properties of the four extracted components using the Hilbert transform. The statistical tests confirmed that the online VMD (VMDon) algorithm returned results of comparable quality to AVMD and CWT, and outperformed standard VMD and WPD. VMDon, AVMD, and CWT extracted four components from the real HRV with frequency content slightly exceeding the previously recognized ranges, suggesting the possibility of their modes mixing. Their ranges of variability were assessed as follows: HF: 0.11–0.40 Hz; LF: 0.029–0.14 Hz; VLF: 4.7–31 mHz; and ULF: 0.002–3.0 mHz.

1. Introduction

Heart rate variability (HRV) describes the temporal changes between successive R peaks in the QRS complexes of the electrocardiogram (ECG). HRV provides insight into the function of both the cardiovascular system and autonomic nervous system (ANS) through the activity of the sympathetic (SNS) and parasympathetic (PNS) nervous subsystems [1,2,3,4,5,6,7,8,9,10]. The influence of neural drive on the recorded HRV is described by increasingly complex integral pulse frequency modulation (IPFM) models [11,12]. A well-studied relationship is the change in heart rate (HR) due to increased SNS or decreased PNS activities [13,14,15]. However, much more information is hidden in HRV spectral properties and, more specifically, in its four frequency bands. The high-frequency (HF) component (0.15–0.4 Hz) is mostly influenced by the PNS activity via the vagus nerve, being synchronized with the respiratory cycle [13,14,16,17]. The low-frequency (LF) component (0.04–0.15 Hz) is modulated by the SNS, PNS, and baroreflex system [13,14,16,17]. Both very-low-frequency (VLF) (0.0033–0.04 Hz) and ultra-low-frequency (ULF) components (<0.0033 Hz) are not well understood, and not as often used as HF or LF. VLF and ULF are generally thought to be associated with the circadian rhythm [14,16,18], temperature regulation [14,17,19], renin–angiotensin system [17,20], physical activity [21], and vasomotor activity [14]. In summary, these properties are not well known (especially for VLF and ULF), and their derivation may contribute to further research in the area of network physiology focusing on physiological interactions within the human organism.

HRV is a non-stationary, irregularly time-sampled signal and is mostly processed after interpolation to the evenly spaced samples, usually using cubic spline interpolation with a 4 Hz sampling rate [22,23,24]. The short-term analysis can minimize the impact of signal non-stationarity using fragments of 5 min or less [23,25]. This is a common approach but does not allow for the separation of lower-frequency components like VLF and ULF. To extract ULF, a long-term analysis is needed, a minimum of several hours (typically 24 h) [16], but then the problem of stationarity must be considered. Another issue is the presence of outliers in the signal. The typical approach for short-term analysis contains the removal of the whole corrupted fragments [4,26,27] and modification of erroneous samples [28,29,30], but this can result in the modification of the phase of the frequency components, especially VLF and ULF. Another method, rarely used, is online analysis, often performed similarly to short-term analysis, extracting only HF and LF with sliding windows of a few or tens of seconds [29,31].

There are many methods by which the HRV spectrum can be derived [32]. Among the most popular are autoregressive modeling, which requires finding the optimal model order [5,9,24,33] and the short-time Fourier transform (STFT) that can determine local frequency components from short quasi-stationary fragments of a signal [1,29,33,34]. Empirical mode decomposition (EMD), on the other hand, is a method designed to extract oscillating components from nonlinear and non-stationary signals [35] and is often followed by the Hilbert transform (HT) to derive instantaneous amplitudes and frequencies [4,23,29,34,36,37]. Multiband filtering (MBF) is an old and seldom used method nowadays, requiring the knowledge of frequency bands for filtering [38], but a recent study has shown that it can give comparable or better results than STFT or EMD [39]. Spectral filtering effects can also be obtained using the wavelet transform, which is a relatively new tool for analyzing non-stationary signals. This is achieved due to the local time properties of the analyzing wavelets and the well-defined frequency properties of the components obtained by the inverse transform. The wavelet transform (WT) was recently used to decompose HRV into three variants: continuous (CWT), discrete (DWT), and wavelet package decomposition (WPD) [40,41,42,43]. One of the newest methods for nonlinear analysis of non-stationary signals is variational mode decomposition (VMD), splitting a signal into several intrinsic mode functions (IMFs), similarly to EMD [44]. While the original VMD procedure requires specifying the number of searched IMFs [45,46], which is a disadvantage compared to EMD, recently an adaptive and data-driven version of this algorithm (AVMD) has been proposed that automatically determines the optimal value of the IMFs [47]. This solution was tested on overnight HRV signals and was able to successfully extract all four components in the offline analysis, outperforming MBF, STFT, and EMD.

In summary, the existing current algorithms operate offline for the vast majority of time, which does not allow for monitoring changes in the frequency components continuously. Meanwhile, the courses of individual components indicate the current processes in various systems of the human body. In addition, older decomposition methods are used (like STFT, EMD, WT, or classic VMD), and VLF and ULF are mostly ignored, still leaving many questions about their nature unanswered. Moreover, there is still a large uncertainty as to the instantaneous properties of the VLF and ULF waveforms, which prevents further research on their internal connections with the physiological processes of various body systems. The solution to these problems would be to develop an algorithm allowing for the online HRV decomposition and thus the monitoring of the underlying processes.

The aim of this work was to develop and validate an online HRV decomposition algorithm yielding simultaneously all four frequency components and to use it to analyze the real HRV signals in a continuous mode. By design, this algorithm should not impose sharp boundaries on the frequency content of the components, and thus, the obtained results can be used to more precisely determine their amplitude–frequency variability.

In the following sections, the data used and the methodology applied are explained, and the results obtained are presented and discussed. Finally, the main conclusion regarding the suitability of the proposed method is provided.

2. Materials and Methods

2.1. Data and Preprocessing

The real data come from two different databases available on the PhysioNet platform [48]. The first one, from St. Vincent’s University Hospital Sleep Disorders Clinic in Dublin, contains recordings from overnight polysomnographic study of 25 adults (age range: 28–68 years) with suspected sleep-disordered breathing. Subjects in the database were randomly selected from patients of the clinic over a half-year period. Additionally, the selected subjects had no known heart disease or autonomic dysfunction and were not taking medications that affect their heart rate. All polysomnographic signals, including ECG sampled at 128 Hz (modified lead V2), were recorded during one night of sleep using the Jaeger–Toennies system. Seven subjects were rejected (#5, 6, 8, 15, 21, 23, 26) due to frequent and long discontinuities in the ECG signal recordings. Because of the lack of ECG signal, it was also necessary to remove the first 54 min of Subject #13 and the last 2 min of Subject #19, resulting in a total of 18 ECG signals from the first database. The second Apnea-ECG Database was also based on sleep apnea-oriented polysomnographic studies. It contains 70 ECG signals lasted around 7–10 h recorded during sleep with a sampling rate of 100 Hz. This database was created from two sets of studies between 1993 and 1995 and 1998 and 1999 [49]. Of all the recordings, it was decided to use 18 ECGs of the best possible quality, especially without signal breaks or strong artifacts (#a04, a05, a09, a11, a18, a20, c09, c10, x01, x03, x05, x07, x08, x12, x16, x21, x23, x26). Even then, due to the lack of signal or its total distortion, it was necessary to remove the first 17 min of c09, 12 min of c10, 11 min of x01, 2 min of x03, 1 min of x26, and the last 10 min of x12. As a result, the recorded ECG signals cover natural, undisturbed by physical activity heartbeat and are long enough to allow the derivation of VLF and ULF components. Ultimately, data from 36 subjects, 18 from each database, were analyzed.

The suggested optimal ECG sampling rate for the HRV spectral analysis is at least 250–500 Hz [22]. Therefore, each ECG signal was resampled from the original 128 or 100 Hz to 1000 Hz to find R peaks with a resolution of 1 ms. The exact location of the R peak was specified as the signal maximum near the QRS complex found using the PhysioToolkit software package (WFDB version 10.6.2) [48]. This software has ready-to-use algorithms and functions designed specifically for medical data processing and has been widely used by researchers since its public release, especially for the analysis of data from the PhysioNet database. However, these methods require parameter tuning, sometimes for each signal individually, and are prone to errors when the signal quality is poor. Then, the HRV signals were formed from the time intervals between consecutive R peaks. Finally, local modifications were made to the outliers by taking care not to modify the original phase of the signal (preserving the lower-frequency components of HRV) [39]: RRs lower than 0.61 s were added to the larger neighbor, and then, values greater than 1.22 s were split into two equal ones (related to the heart rate of 49–98 bpm when sleeping). Then, the timeline was reconstructed as a cumulative sum of RR intervals. It should be emphasized that the purpose of this work is to present quantitative information about HRV components; therefore, ECG signals should be of high-quality, without many artifacts or discontinuities, and sample modifications should be minimal and only local.

In addition to 36 real signals, one synthetic HRV signal with features mimicking the real data was also used to assess the effectiveness of the decomposition algorithm by comparing the obtained results with the known parameters used to generate the synthetic waveform (ground truth). This non-stationary 6 h HRV was generated in the same way as in [39]. It consists four components x_c, each simulated as amplitude (AM)- and frequency (FM)-modulated continuous signals:

(1)$\begin{array}{c}{x}_c(t)=A_c(t)\sin \left({\phi }_c(t)\right),\\ A_c(t)=A_{am0}+A_{am}\sin \left(2\mathsf{\pi }{f}_{am}t\right),\\ {\phi }_c\left(t\right)=2\mathsf{\pi }{A}_{fm0}t+A_{fm}/f_{fm}\sin \left(2\mathsf{\pi }{f}_{fm}t\right),\end{array}$

where A_am0, A_am and A_fm0, A_fm represent the range of magnitude variation in amplitudes and frequencies and f_am, f_fm are their modulation frequencies (Table 1). The values of these parameters were estimated from the study of three real HRV signals [39]. Then, the continuous HRV, sampled at 1 kHz, was created as the sum of all four simulated components (HF, LF, VLF, and ULF), shifted up by 0.95 s. Then, irregularly spread HRV samples were derived from it, taking into account the relationship between the HRV values and the timeline. Next, the signal was disturbed by adding Gaussian noise with a standard deviation of 0.01 s. To obtain the final uniformly sampled HRV signals, a necessary condition for using the intended signal processing procedures, cubic spline interpolation at 2 Hz was used, and an antialiasing zero-phase filter with a cutoff frequency of 0.5 Hz was applied to minimize HRV spectrum distortions.

In the rest of this work, we assume that the HRV signal is ready for use, and new samples appear successively as a result of preprocessing.

2.2. Variational Mode Decomposition

As mentioned in the Introduction, VMD-based HRV decomposition procedures outperform other approaches such as MBF, STFT, or EMD; therefore, this research focuses on adapting VMD to online HRV decomposition.

VMD is a relatively new method enabling the optimal decomposition of a non-stationary signal into several intrinsic mode functions (IMFs) of different center frequencies [44]. It uses adaptive Wiener filters, and optimization is performed with the alternating-direction multiplier method. As a non-recursive approach, all searched modes are extracted simultaneously, and the final IMFs result from minimizing the sum of the spectral widths of all modes [44]. The original VMD procedure requires the predetermined number of searched modes (maxIMFs), and the incorrect setting of this parameter may result in their omission or mixing [50]. It is especially important since each of the non-stationary component of HRV can potentially contain several dominant frequencies in its characteristic range. Therefore, an adaptive VMD (AVMD) that can find the optimal number of IMFs has been recently proposed for the offline analysis of HRV [47]. The phenomenon underlying this method is the significant reduction in residuals with a justified increase in maxIMFs. At the same time, this parameter should be as low as possible due to the increasing complexity of AVMD with less and less improvement in the quality of decomposition for larger values. An additional condition, and in this work the only change from the original AVMD [47], is that at least one IMF is assigned to each HRV component based on its frequency range. If this is not the case, the algorithm continues to run. Assigning an IMF to a specific HRV component involves the use of the HT, commonly applied in conjunction with VMD [44]: the instantaneous frequencies are calculated, and their histogram is created; all frequencies are counted within the standard ranges corresponding to the HRV components, and finally, the IMF is assigned to the one whose range is most represented—this prevents the imposition of rigid frequency limits on individual components. The final step after finding the optimal maxIMFs is to create the HRV components as the sum of all IMFs assigned to them [47].

2.3. Online Algorithm

The proposed online algorithm based on VMD (VMDon) differs significantly from the offline one. HRV is analyzed in windows moving by one sample, containing the latest measured value and several older ones. The window length for each component is different and, by design, is equal to the period corresponding to the lower limit of the given frequency range. Such a window length allows to capture all the frequencies of interest. In practice, after preprocessing, HRV was resampled to 2 Hz, so the window lengths were rounded by 0.5 s. Another issue is the lack of a defined lower-frequency limit for ULF; in this case, the length of the window was chosen empirically as the shortest one, allowing the effective extraction of the lowest synthetic frequencies. In summary, the window lengths for the HF, LF, VLF and ULF components were, respectively, T_HF = 6.5 s, T_LF = 25 s, T_VLF = 303 s, and T_ULF = 1800 s.

In the analyzed window, samples are affected by lower frequencies from outside the range of the component. For this reason, the linear trend within the window is first removed to eliminate this influence. Since only one component is expected in the remaining part of the signal, VMD returning only one IMF is employed, directly being the estimated frequency component for that particular window. Then, the analyzing window moves by one sample with a new HRV value, and the calculations are repeated, which means that for a given moment of time, many different estimates of a given component sample are calculated. The final estimate is the average of the estimates for that sample calculated within all window positions covering it. Averaging is performed at the very beginning of the window (the oldest sample) because, only for that sample, all the estimates from the previous locations of the window are accessible (Figure 1).

Each time the higher-frequency component is obtained, it is subtracted from the original HRV to form the input signal for extracting the lower-frequency component. Thus, only the searched component is expected in the window since the lower frequencies are eliminated by trend removal, and the higher frequencies have already been extracted. This causes additional delay because the extraction of the lower frequencies requires ready-made results for the higher ones (Figure 1).

The entire procedure for one VMDon iteration is shown as a pseudocode in Algorithm 1. The extraction of lower-frequency components is carried out in the same way as in the case of HF. Additionally, ULF calculated as a last waveform is increased by the value of the trend previously subtracted before its extraction, as it represents the basal HR. The additional flowchart of the VMDon is presented in Figure 2.

Algorithm 1 Variational Mode Decomposition Online (VMDon)

1:.       while do         2:.               collect new $HRV[t]$ sample         3:.               ${HRV}_{HF}=HRV[t-T_{HF}: t]$ ← shift HF window         4:.               ${HRV}_{HF}={HRV}_{HF}-\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{d}\left({HRV}_{HF}\right)$ ← remove trend         5:.               ${IMF}_t=\mathrm{V}\mathrm{M}\mathrm{D}\left({HRV}_{HF}\right)$ ← extract one IMF from current HF window         6:.               $HF[t-T_{HF}]={\displaystyle \frac{1}{T_{HF}}}{\sum }_{i=t-T_{HF}}^t{IMF}_i[t-i]$ ← extract a single HF sample as the average of all IMFs for that specific point in time         7:.               ${HRV}_{LF}=HRV-HF$ ← remove HF from HRV before LF extraction         8:.               extract LF, VLF and ULF repeating steps 4–7 in analogous way         9:.               $ULF=ULF+\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{d}\left(ULF\right)$ ← restore basal HRV       10:.               wait until the next HRV sample is measured       11:.       end while

2.4. Wavelet Transforms

The essence of the wavelet transform (WT) is the use of a local mother wavelet with a shape similar to the sought features in the non-stationary signal together with daughter wavelets of increasingly smaller scale (closely related to the frequency band), which are moved along the signal and high- and low-pass filter it, producing the subsequent levels of decomposition. The effects are the approximation and detail coefficients obtained for each of the levels, from which, using the inverse transform, the components of the analyzed signal with a specific frequency content can be obtained. The WT is applied in a continuous (CWT) or discrete (DWT) form. A special case of the second one is the wavelet package decomposition (WPD). In this work context, WPD is the better option than the original DWT due to the ability to precisely determine the frequency ranges of the extracted components.

In this work, the CWT was used to decompose HRV in a similar manner as in [42], using the Morlet wavelet. The resulting scalogram (2D) was then transformed into the HRV components by applying the inverse transform (ICWT) and summing its samples over four standard frequency ranges (see Introduction) for each time instant.

Similarly, in the case of WPD, an analogous approach was used as in the compared works [40,41]. Since the analyzed HRV is sampled at a frequency of 2 Hz, a 10-level decomposition with the biorthogonal wavelet (bior4.4) was applied, which for the spectrum down to 1 Hz allowed obtaining band boundaries with an accuracy of 1/2¹⁰ ≈ 1 mHz. Finally, the individual HRV components were obtained by summing the inverse transform (IWPD) results for the levels: 1–3 for ULF (including the HRV mean value), 4–40 for VLF, 41–153 for LF, and 154–410 for HF.

2.5. AM and FM Extraction

The amplitude (AM) and frequency (FM) modulation ranges of HRV components are still not well explored, so special attention has been paid to this aspect. To calculate instantaneous amplitudes and phases, the Hilbert transform (HT) was used for each extracted component. This part of the analysis was carried out after all components had been previously extracted and was not part of online processing. The instantaneous amplitudes obtained in this way may contain disturbances caused by online component extraction, and like post-processing, they were smoothed using the Savitzky–Golay filter (SGF) [51]. The SGF is a least-squares smoothing filter that fits a low-degree polynomial to a fragment of the signal in a window moving over it. Another application of SGF is the computation of the derivative that was applied here to obtain instantaneous frequencies directly from instantaneous phases, reducing the problem of negative frequencies resulting from numerical differentiation after HT [52]. Calculated amplitudes and frequencies were smoothed by SGF in the same way. In each case, the order of the polynomial was 2, and the length of the filter corresponded to the length of the window related to the analyzed component.

2.6. Validation of the Method

To quantitatively evaluate the extraction efficiency of individual components, the synthetic signal with exactly known properties values was used. The following similarity measures between ground truth (synthetic components) and the algorithms’ results were calculated, i.e., Spearman’s coefficient r_s of rank correlation between the signals’ samples, the medians of absolute residuals |$\widetilde{{\epsilon }_c}$|, and the relative error of extraction:

(2)${\delta }_{ex}=\frac{{‖x_c-x_s‖}_2}{{‖x_s‖}_2}100\%=\frac{{‖{\epsilon }_c‖}_2}{{‖x_s‖}_2}100\%,$

where x_c and x_s represent the extracted and synthetic components, respectively. All metrics were calculated for VMDon and AVMD, as well as for the compared methods used in the literature: standard VMD (set to decompose the signal into four expected IMFs), CWT, and WPD. Additionally, the one-sample Kolmogorov–Smirnov test was performed to check whether the distributions of |$\widetilde{{\epsilon }_c}$| values were consistent with the normal distribution, and the one-sided Wilcoxon signed-rank test (α = 0.05) to verify if |$\widetilde{{\epsilon }_c}$| for one method was statistically significantly lower than for another. Validation was performed for each of the derived components as well as for their AM and FM.

In a similar way, the similarities between components obtained using VMDon, AVMD, and CWT (which turned out to be the best of the previously used methods) were evaluated for all real HRV signals. Although, in this case, it is not possible to assess the correctness of extraction because the true values are not known, it is still possible to check the consistency level of the results returned by both methods.

2.7. Analysis of Real HRVs

The VMDon, AVMD, and CWT algorithms were run on the overnight HRV signals from two datasets, each comprising 18 subjects, to split them into the frequency components and calculate their AM and FM. On this basis, AM and FM variability ranges were determined for each subject using HT. The results obtained when analyzing HRVs from these two different datasets were compared by performing the two-sample nonparametric Kolmogorov–Smirnov test (α = 0.05), checking the equality of the probability distributions of the AM and FM samples of components extracted separately using each of the three decomposition methods. For deeper insight, the population AM and FM values for the combined data from all subjects were also presented in histograms, and their ranges were assessed based on the empirical cumulative distribution function (99th percentile). Moreover, it was calculated what percentage of the FM samples were inside the generally accepted literature frequency ranges of the corresponding components.

3. Results

3.1. Validation of Methods

The efficiency of the newly proposed online algorithm VMDon and modified AVMD was assessed based on the known properties of the synthetic signal and additionally compared to those obtained with the previously used offline procedures: standard VMD, CWT, and WPD (Figure 3, Figure 4 and Figure 5 and Table 2).

Figure 3, Figure 4 and Figure 5 show that the results of decomposition by classic VMD and WPD are clearly inferior to the other algorithms. VMDon and CWT perform best with lower-frequency components, where the extracted signals almost match the true values, although there are temporary deviations visible in ULF AM. At higher frequencies, more noise and fluctuations are observed, but the overall reproduction of the signals is proper.

All the numerical results of the validation of the proposed methods against those previously used in the literature are summarized in Table 2. They confirm the inferiority of the standard VMD method in comparison to its online and adaptive versions. Similarly, WPD also turns out to be worse than them. Nevertheless, the CWT proves to be a method giving similar results or even slightly better in terms of the VLF component. AVMD is superior in the derivation of HF and LF components, but VMDon yields slightly better results in case of VLF and ULF. In general, these three algorithms (VMDon, AVMD, and CWT) present good alignment to the synthetic signal, with the minor exception of HF FM and LF FM in the case of VMDon. The result of the Wilcoxon test for ULF AM is not obvious, as it indicates a significantly smaller median for VMDon where |$\widetilde{{\epsilon }_c}$| is larger than for AVMD. It is also worth noting that VMDon, AVMD, CWT, and WPD always have lower |$\widetilde{{\epsilon }_c}$| values than VMD, also according to the Wilcoxon test.

In the next step, the best performing three algorithms, VMDon, AVMD and CWT, were tested on real overnight HRVs from the two datasets. The correlation coefficients (see Supplementary Material for details) show overall consistency between the methods, with most results having r_s greater than 0.8. Nevertheless, there is noticeable discrepancy in the frequency modulation of HF, LF, and VLF, where r_s is mostly below 0.5, and the values are more diverse among subjects.

Results of VMDon, AVMD, and CWT were also compared by visualizing the extracted components, as well as their AM and FM (representative results for Subject #2 are shown in Figure 6). The results correspond well, although VMDon tends to slightly underestimate the magnitude range of FM. It is also visible that, while most of the signals propagate in some range of magnitude, there are well-rapid spikes out of the usual boundary, especially in case of HF, and these occurrences are detected regardless of the derivation method.

3.2. Derived HRV Components and Their Properties

Figure 7 shows the empirical probability density of instantaneous amplitude values based on the combined data from all subjects, however, separately for the first and the second dataset. This points out the overall distribution of the instantaneous amplitudes of HRV components. In general, the histograms overlap quite well for the three methods, which confirms the correctness of the found distribution in the population. There is also a great similarity between the corresponding histograms representing both datasets. AM distributions for higher-frequency components are more concentrated at lower amplitudes.

The extracted real instantaneous frequencies are illustrated in Figure 8 for all components and the three decomposition methods, for both datasets separately. Although the separation between them is evident, adjacent bands partially overlap in all cases.

4. Discussion

The aim of this study was to derive and validate an online HRV decomposition method and determine the characteristics of the HF, LF, VLF, and ULF components based on real data. For the latter purpose, it was compared with the modified AVMD and other offline decomposition methods used so far. Moreover, although the preprocessing of ECG signals had to be adapted to their database-dependent form, all algorithms compared in this work operated on the same HRV, so the preprocessing itself had no influence on the obtained results.

In the study presented here, only a moderate-level random noise in the synthetic HRV, simulating R peak-positioning errors in ECG, was considered. On the other hand, it would be interesting to comprehensively evaluate the influence of typical ECG disturbances on HRV-decomposing results by the proposed algorithms. In this case, it would be necessary to use a standard pure ECG signal and numerical models simulating individual disturbances, which, however, is beyond the scope and purpose of this work. Therefore, solving this problem requires new research, which the authors plan for the future.

It has been previously shown that AVMD is better than the use of STFT and EMD and comparable to MBF (while still being better for extracting HF and LF) [47]. With this in mind, it can be now stated that the evaluated own-elaborated algorithms, VMDon and AVMD, as well as standard CWT, proved to be better than others previously used for HRV decomposition. In addition, the modification of the aforementioned AVMD algorithm used in this work enabled the extraction of ULF from all analyzed HRVs (the original version did not always succeed).

The analyzed algorithms require uniform sampling. A resampling to 2 Hz was chosen after testing the impact of lower and higher sampling rates. The lower ones led to low HRV resolution and the need for an upper limit on the HRV spectrum, while higher sampling frequencies did not change the results but required more computing power to, for example, be able to update the data every 0.25 s for 4 Hz. It was also noted that removing the window-averaged residuals after each extracted component does not significantly affect the results. It turns out, however, that the residuals after the VMDon algorithm are an order of magnitude larger than after AVMD. Using the synthetic signal, it is possible to observe a part of the HF component in them, while at the same time, VMDon performs noticeably worse with HF than AVMD (Table 2). This leads to the conclusion that, in the future, one can try to improve the extraction of the HF component by a modified version of VMDon, considering the residuals.

Figure 3, Figure 4 and Figure 5 and more accurately Table 2 show, using synthetic data representing ground truth, that both proposed methods, VMDon and modified AVMD, perform well in decomposing HRV into the four known frequency components. The distortions seen in the extracted components, especially HF and LF (Figure 3), are due to the operation of the VMDon procedure using a sliding window and determining only one component sample for a given window position. Although the averaging of such obtained values for a given moment covered by many consecutive window positions and SGF filtering reduce the level of interference, nevertheless, the essence of the implemented procedure is the discrete nature of the extracted components. The situation is similar when summing the corresponding IMFs forming a given component in the AVMD procedure. This effect is enhanced when the FM of the components is determined (Figure 5), since this procedure requires the numerical calculation of the derivative of the instantaneous phase returned by the Hilbert transform, and the differentiation, having the nature of high-pass filtering, amplifies the energy of high-frequency interference. Nonetheless, using SGF to directly calculate the derivative from the instantaneous phase after HHT gave much better results compared to the previously used method based on the discrete first-order derivative and subsequent median filtering [47].

In the case of VMDon, the efficiency of component extraction is lower at the beginning and end of a finite signal due to a reduced number of averaged samples. On the other hand, during online analysis, older samples already exist and the ECG recording continues. Another attribute of the online algorithm is that it operates with a certain delay—component extraction is delayed by the length of the associated window and the operating rule used. This delay is 6.5, 25, and 303 s for HF, LF and VLF, respectively, and 30 min for ULF. Such a window length for ULF was adopted despite the lack of verified information about its lower frequency limit; however, no difference in results was noticed for trials with windows as long as 90 min, while the 30 min window has shorter delay and much lower computational complexity. In addition, it should be noted that the 30 min delay applies only to ULF, which is inherently slow to change, while the delays of the most used HF or LF components are much shorter, which should enable, for example, the real-time detection of sought-after episodes related to the ANS.

The Kolmogorov–Smirnov tests showed that the distributions of values in all analyzed data groups were strongly non-normal, justifying the use of Spearman’s correlation coefficient and the Wilcoxon test. In the case of AM of ULF, the indication of the test in Table 2 is inconsistent with the value of |$\widetilde{{\epsilon }_c}$|, where the median value of the residuals is smaller for AVMD, but the Wilcoxon test shows that VMDon has a statistically lower median in the population. With the Wilcoxon test, this can indeed happen when discrete and extremely skew data are analyzed [53,54]. This is the case in this particular instance because the Wilcoxon test calculates the differences between the corresponding elements from both samples from AVMD and VMDon, and they are mostly positive for AM ULF, having also a positive median.

The different nature of the procedures has resulted in a higher accuracy of AVMD in extracting HF and LF and an advantage of VMDon in extracting VLF and ULF; although for AM and FM of ULF, the results of both are comparable. Additionally, CWT proves to be the best for extracting the VLF component (Table 2). The importance of this observation is highlighted by the results of the real data decomposition analysis, in particular the lower correlation between the FM samples of components extracted for HF, LF, and VLF (see Supplementary Material). This effect is also seen in Figure 6.

Temporary peaks in the values observed in the components and thus in the AM (Figure 6, best seen in the HF) occur for all three presented algorithms at the same time. This is therefore a feature of the components themselves, not the method. Their explanation, however, requires the further analysis of the temporary cardiac rhythm controlled by the ANS.

The differences between the methods seen in Figure 6 (a slight shift in results from VMDon toward smaller values, especially for FM) when decomposing real HRVs, and not visible for the synthetic data, are evidence of imperfections in the HRV synthesizing procedure. They may come from the narrower range of variability of the component parameters used in the HRV simulation, previously estimated using data from only three subjects [39]. This work presents a more accurate estimation of AM and FM variability based on data from 18 individuals. This is a good starting point for the generation of a more reliable synthetic HRV [12,32,39], which will be helpful for understanding the reasons for the observed difference in the results of the two algorithms.

Important information on the characteristics of the HRV components comes from the analysis of the distributions of AM and FM values (Figure 7 and Figure 8). AM most often takes on small values and mostly does not exceed about 230 ms in all four cases; however, the higher the component frequencies (shorter instantaneous periods), the smaller the range of their variability (Figure 7). Similar conclusions emerge regarding the ranges of FM variability (Figure 8). This observation is in agreement with the power laws (scale invariance) observed in nature [55].

Figure 7 and Figure 8, as well as Table 3, allow the comparison of results obtained with the same methods from two datasets of different origin, which allows for cross-dataset validation. The obtained outcomes show a high similarity between them and allow one to conclude about the reliability of components extracted with significantly different algorithms from two separate datasets. At the same time, the two-sample Kolmogorov–Smirnov tests performed for the corresponding results from these two datasets rejected the hypothesis that they come from the same distributions. This indicates, however, a statistical difference in the results obtained from HRVs originating from these different databases. Nevertheless, additional checks showed that the main reason for this outcome is the sensitivity of the test itself to small differences in the data; it was enough to remove data of 1 patient from a given set of 18 patients, and the test returned a statistical difference between these two cases.

The above results, which validate the HRV component waveforms obtained using the proposed algorithms, are of particular importance for further analyses of the VLF and ULF components, which are currently the least understood. Thanks to the reliable determination of the frequency ranges of these components and the possibility of isolating their waveforms obtained in these studies, it will be possible in the future to test hypotheses regarding their origin. However, the verification of such hypotheses will require both appropriate HRV decomposition algorithms and simultaneous measurement of other physiological quantities. It is already suggested that they can support diagnosis, monitoring, and/or therapy in cases of aortic regurgitation, coronary artery bypass grafting, vasovagal syncope, postural tachycardia syndrome, or heart failure [56,57,58,59,60].

Using the FM distribution in the combined two populations (shown separately in Figure 8), it was checked what percentage of the actual instantaneous frequencies are within the accepted literature frequency range of a given component. For VLF (based on CWT) and ULF (based on VMDon), it is more than 99.9%, and for HF and LF (based on AVMD), it is 94.7% and 95.0%, respectively; thus, the obtained results are in general agreement with the ranges reported earlier. Nevertheless, the results of the analysis of real data show an overlap of these ranges (Figure 8, Table 3), which is known as the mode-mixing problem (particularly for EMD [61]) and poses a significant challenge for all signal decomposition algorithms [32,62,63]. The feasibility of observing band overlapping in this work is due to the component extraction procedures used, which do not assume rigid frequency limits, unlike using other approaches [39].

Data from the population-based AM and FM distributions for two datasets (results from AVMD for HF and LF, CWT for VLF, and VMDon for ULF methods gathered together) were used to quantify their ranges for all HRV components. Table 3 contains these ranges defined by the 99th percentile of the empirical cumulative distribution function.

There were also attempts to identify a possible additional very-high-frequency component (higher than 0.4 Hz) using a 2.5 s window in VMDon to extract this range first, but the results were significantly worse, and even such a short window still captured primarily the usual HF component.

5. Conclusions

This study presents and validates two new HRV signal decomposition methods, VMDon and AVMD, compared to the algorithms previously used for this purpose. To the best of the authors’ knowledge, the online algorithm decomposing HRV into HF, LF, VLF, and VLF waveforms has been proposed for the first time. The results of the work also indicate which decomposition approaches are the best, in particular which method (including new and previously used ones) is the most precise for a given component: AVMD is the more accurate for HF and LF components extraction and VMDon for VLF and ULF ones, the latter of which can also be used offline in this form. Furthermore, the suitability of CWT for VFL component derivation is also demonstrated. In particular, they can be applied to the same data, allowing for the accurate extraction of all four HRV components.

HRV decomposition results from 36 subjects suggest slightly different frequency ranges of the components than those usually proposed in the literature and the possibility of mixing their modes, which should be considered in the future work on this topic. Consistent results in this respect were obtained using the three best decomposition methods, which confirm the reliability of the component characteristics assessed in this work. Now, for the convenience of use in the future work, they are gathered together in Table 4.

There are still a few other aspects to work on. First, one can try to develop a comprehensive online algorithm to extract HRV components directly from ECG signal, possibly to be embedded in, e.g., Holter devices. Having isolated HRV components, it will be easier to identify the relationship of low-frequency waveforms (VLF, ULF) with specific physiological processes through ECG/HRV analysis for cases when other biosignals were recorded at the same time. In addition, the correct extraction of HRV components also creates a new possibility to infer sinus rhythm stimulation by attempting to invert the integral pulse frequency modulation (IPFM) models.

Author Contributions

Conceptualization, K.A. and A.G.P.; methodology, K.A. and A.G.P.; software, K.A.; validation, K.A. and A.G.P.; formal analysis, K.A. and A.G.P.; investigation, K.A.; resources, K.A.; data curation, K.A.; writing—original draft preparation, K.A.; writing—review and editing, A.G.P.; visualization, K.A.; supervision, A.G.P. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Datasets used in this study are available on the Physionet website: https://physionet.org/content/ucddb/1.0.0/ (accessed on 11 January 2025); https://physionet.org/content/apnea-ecg/1.0.0/ (accessed on 11 January 2025).

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The principle of operation of the VMDon algorithm (window lengths of virtual components A (3 samples) and B (5 samples) are set only for demonstration of the algorithm).

Figure 2. The flowchart of the VMDon algorithm including the main steps for one iteration (new HRV sample).

Figure 3. Results of extraction of HRV components using VMDon (blue), AVMD (orange), classic VMD (purple), CWT (green), and WPD (red) algorithms compared to the synthetic signal (ground truth) components (black).

Figure 4. Results of amplitude modulation (AM) of HRV components extracted using VMDon (blue), AVMD (orange), classic VMD (purple), CWT (green), and WPD (red) algorithms compared to the synthetic signal (ground truth) components (black).

Figure 5. Results of frequency modulation (FM) of HRV components extracted using VMDon (blue), AVMD (orange), classic VMD (purple), CWT (green), and WPD (red) algorithms compared to the synthetic signal (ground truth) components (black).

Figure 6. Subject #2’s real HRV components extracted using VMDon (blue), AVMD (orange), and CWT (green) with their amplitude and frequency modulations.

Figure 7. Distributions of the real amplitude values of HRV components in the population of subjects from the first (St. Vincent University Hospital) and second (Apnea-ECG) databases obtained using VMDon (blue), AVMD (orange), and CWT (light green). Additional colors indicate shared result areas of the algorithms: olive: all algorithms; turquoise: VMDon and CWT; red: VMDon and AVMD; light brown: AVMD and CWT.

Figure 8. Distributions of the real frequency values of HRV components in the population of subjects from the first (St. Vincent University Hospital) and second (Apnea-ECG) databases obtained using VMDon, AVMD, and CWT.

Table 1

Core parameters for synthetic HRV signal generation.

Component	Aam0(ms)	Aam(ms)	Afm0(mHz)	Afm(mHz)	fam(mHz)	ffm(mHz)
HF	25	15	240	60	1.05	0.96
LF	35	25	80	20	0.67	0.81
VLF	50	30	17.5	7.5	0.37	0.45
ULF	75	35	0.8	0.6	0.22	0.27

Table 2

The accuracy of the methods: relative error δ, Spearman’s correlation coefficient r_s, and median of absolute residuals |$\widetilde{{\epsilon }_c}$| (the medians of absolute residuals are in ms for components and AMs and in mHz for FMs).

Signal	VMDon			AVMD			VMD			CWT			WPD
	δ (%)	rs	|$\widetilde{{\epsilon }_c}$|	δ (%)	rs	|$\widetilde{{\epsilon }_c}$|	δ (%)	rs	|$\widetilde{{\epsilon }_c}$|	δ (%)	rs	|$\widetilde{{\epsilon }_c}$|	δ (%)	rs	|$\widetilde{{\epsilon }_c}$|
HF comp.	42.7	0.915	5.6 &	35.5	0.941	4.6 #,*,&	65.4	0.782	7.6	38.8	0.929	5.0 #,&	52.6	0.878	6.3
HF AM	31.8	0.841	6.5	22.1	0.852	4.1 #,*,&	53.9	0.390	9.6	22.9	0.842	4.2 #,&	30.4	0.741	5.3 #
HF FM	23.5	0.676	28.0	16.0	0.898	26.1 #,*,&	36.6	0.243	46.7	18.1	0.874	25.8 #,&	20.5	0.696	26.8 #
LF comp.	25.6	0.975	4.6 &	12.7	0.991	2.4 #,*,&	69.3	0.694	7.4	16.5	0.985	3.1 #,&	34.3	0.945	5.7
LF AM	19.7	0.960	5.9	8.4	0.974	2.2 #,*,&	59.5	0.346	7.2	8.8	0.973	2.3 #,&	15.3	0.934	3.7 #
LF FM	16.9	0.851	0.9 &	3.7	0.976	0.7 #,*,&	65.9	−0.321	2.9	5.0	0.961	0.8 #,&	11.5	0.867	1.4
VLF comp.	11.7	0.997	2.5 &	19.6	0.994	2.2 #,&	100.3	0.332	26.4	2.6	1.000	0.7 #,$,&	16.8	0.988	3.7
VLF AM	9.9	0.981	2.9 $	19.3	0.932	3.5	71.0	0.216	30.5	0.8	0.999	0.2 #,$,&	3.9	0.993	1.1 #,$
VLF FM	0.5	1.000	0.1 $,&	0.6	1.000	0.1	300.7	−0.468	46.3	0.5	1.000	0.1 #,$,&	0.6	1.000	0.1 $
ULF comp.	0.2	1.000	0.9 $,*,&	0.7	0.990	1.4 *	3.4	0.856	21.3	0.4	0.998	3.6	0.3	0.998	1.6 $,*
ULF AM	11.6	0.939	3.9 $,*,&	11.3	0.936	3.5 *	18.2	0.845	10.0	12.1	0.927	3.4	11.4	0.936	3.3 *
ULF FM	12.1	0.999	0.1	12.0	0.999	0.1 #	402.5	0.190	0.8	11.8	0.998	0.1 #,$,&	11.9	0.999	0.1 #,$

Based on the Wilcoxon signed-rank test (α = 0.05): ^#—value significantly lower than for VMDon; ^$—value significantly lower than for AVMD; *—value significantly lower than for CWT; and ^&—value significantly lower than for WPD; |$\widetilde{{\epsilon }_c}$| for VMDon, AVMD, CWT, and WPD are always significantly lower than that for VMD; bold values indicate the best result of a given measure among all algorithms.

Table 3

The AM and FM ranges in the population of all subjects, separately for both datasets.

	First Database		Second Database
Component	AM (ms)	FM (mHz)	AM (ms)	FM (mHz)
HF	1.1–233	112–398	1.1–118	110–398
LF	2.1–188	31–142	3.2–139	28–132
VLF	7.2–195	4.7–30	8.2–147	4.8–31
ULF	8.2–221	0.002–3.4	7.2–134	0.002–2.9

Table 4

The AM and FM ranges (99th percentile) in the population of 36 subjects.

Component	AM (ms)	FM (mHz)
HF	1.1–202	111–398
LF	2.5–169	29–137
VLF	7.7–162	4.7–31
ULF	7.7–175	0.002–3.0

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/app15031210/s1, Table S1: Spearman’s correlation r_s between results of VMDon and AVMD for real HRVs from St. Vincent University Hospital database; Table S2: Spearman’s correlation r_s between results of VMDon and AVMD for real HRVs from Apnea-ECG database; Table S3: Spearman’s correlation r_s between results of VMDon and CWT for real HRVs from St. Vincent University Hospital database; Table S4: Spearman’s correlation r_s between results of VMDon and CWT for real HRVs from Apnea-ECG database; Table S5: The AM and FM ranges in the population of 18 subjects from St. Vincent University Hospital database; Table S6: The AM and FM ranges in the population of 18 subjects from Apnea-ECG database; Table S7: The AM and FM ranges in the population of 36 subjects from combined databases.

References

1. Acharya, U.R.; Joseph, K.P.; Kannathal, N.; Lim, C.M.; Suri, J.S. Heart rate variability: A review. Med. Biol. Eng. Comput.; 2006; 44, pp. 1031-1051. [DOI: https://dx.doi.org/10.1007/s11517-006-0119-0]

2. Cerutti, S.; Goldberger, A.L.; Yamamoto, Y. Recent advances in heart rate variability signal processing and interpretation. IEEE Trans. Biomed. Eng.; 2006; 53, pp. 1-3. [DOI: https://dx.doi.org/10.1109/TBME.2005.859777]

3. Fujiwara, K.; Abe, E.; Kamata, K.; Nakayama, C.; Suzuki, Y.; Yamakawa, T.; Hiraoka, T.; Kano, M.; Sumi, Y.; Masuda, F. et al. Heart rate variability-based driver drowsiness detection and its validation with EEG. IEEE Trans. Biomed. Eng.; 2019; 66, pp. 1769-1778. [DOI: https://dx.doi.org/10.1109/TBME.2018.2879346] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30403616]

4. Pan, W.; He, A.; Feng, K.; Li, Y.; Wu, D.; Liu, G. Multi-frequency components entropy as novel heart rate variability indices in congestive heart failure assessment. IEEE Access; 2019; 7, pp. 37708-37717. [DOI: https://dx.doi.org/10.1109/ACCESS.2019.2896342]

5. Callara, A.L.; Sebastiani, L.; Vanello, N.; Scilingo, E.P.; Greco, A. Parasympathetic-sympathetic causal interactions assessed by time-varying multivariate autoregressive modeling of electrodermal activity and heart-rate-variability. IEEE Trans. Biomed. Eng.; 2021; 68, pp. 3019-3028. [DOI: https://dx.doi.org/10.1109/TBME.2021.3060867]

6. Patural, H.; Franco, P.; Pichot, V.; Giraud, A. Heart rate variability analysis to evaluate autonomic nervous system maturation in Neonates: An expert opinion. Front. Pediatr.; 2022; 10, 860145. [DOI: https://dx.doi.org/10.3389/fped.2022.860145]

7. Grégoire, J.M.; Gilon, C.; Carlier, S.; Bersini, H. Autonomic nervous system assessment using heart rate variability. Acta Cardiol.; 2023; 78, pp. 648-662. [DOI: https://dx.doi.org/10.1080/00015385.2023.2177371]

8. Nayak, S.K.; Pradhan, B.; Mohanty, B.; Sivaraman, J.; Ray, S.S.; Wawrzyniak, J.; Jarzębski, M.; Pal, K. A Review of Methods and Applications for a Heart Rate Variability Analysis. Algorithms; 2023; 16, 433. [DOI: https://dx.doi.org/10.3390/a16090433]

9. Pinto, H.; Antonacci, Y.; Pernice, R.; Barà, C.; Javorka, M.; Faes, L.; Rocha, A.P. Decomposing the mutual information rate of heart period and respiration variability series to assess cardiorespiratory interactions. Proceedings of the 2023 45th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC); Sydney, Australia, 24–27 July 2023; [DOI: https://dx.doi.org/10.1109/EMBC40787.2023.10341174]

10. Lu, L.; Zhu, T.; Morelli, D.; Creagh, A.; Liu, Z.; Yang, J.; Liu, F.; Zhang, Y.-T.; Clifton, D.A. Uncertainties in the Analysis of Heart Rate Variability: A Systematic Review. IEEE Rev. Biomed. Eng.; 2024; 17, pp. 180-196. [DOI: https://dx.doi.org/10.1109/RBME.2023.3271595]

11. Bayly, E.J. Spectral analysis of pulse frequency modulation in the nervous systems. IEEE Trans. Biomed. Eng.; 1968; 15, pp. 257-265. [DOI: https://dx.doi.org/10.1109/TBME.1968.4502576]

12. von Rosenberg, W.; Hoting, M.O.; Mandic, D.P. A physiology based model of heart rate variability. Biomed. Eng. Let.; 2019; 9, pp. 425-434. [DOI: https://dx.doi.org/10.1007/s13534-019-00124-w] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31799012]

13. Shaffer, F.; Ginsberg, J.P. An overview of heart rate variability metrics and norms. Front. Public Health; 2017; 5, 258. [DOI: https://dx.doi.org/10.3389/fpubh.2017.00258] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29034226]

14. Pham, T.; Lau, Z.J.; Chen, S.H.A.; Makowski, D. Heart rate variability in psychology: A review of HRV indices and an analysis tutorial. Sensors; 2021; 21, 3998. [DOI: https://dx.doi.org/10.3390/s21123998]

15. Bari, V.; De Maria, B.; Cairo, B.; Gelpi, F.; Lambert, E.; Esler, M.; Baumert, M.; Porta, A. Assessing correlation between heart rate variability markers based on laguerre expansion and direct measures of sympathetic activity during incremental head-up tilt. Proceedings of the 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC); Online, Mexico, 1–5 November 2021; [DOI: https://dx.doi.org/10.1109/EMBC46164.2021.9630873]

16. Karim, N.; Hasan, J.A.; Ali, S. Heart rate variability—A review. J. Basic Appl. Sci.; 2011; 7, pp. 71-77.

17. Noda, A.; Hayano, J.; Ito, N.; Miyata, S.; Yasuma, F.; Yasuda, Y. Very low frequency component of heart rate variability as a marker for therapeutic efficacy in patients with obstructive sleep apnea: Preliminary study. J. Res. Med. Sci.; 2019; 24, 84. [DOI: https://dx.doi.org/10.4103/jrms.JRMS_62_18]

18. Stampfer, H.G.; Dimmit, S.B. Variations in circadian heart rate in psychiatric disorders: Theoretical and practical implications. ChronoPhysiol. Ther.; 2023; 3, pp. 41-50. [DOI: https://dx.doi.org/10.2147/CPT.S43623]

19. Aoki, K.; Stephens, D.P.; Johnson, J.M. Diurnal variation in cutaneous vasodilator and vasoconstrictor systems during heat stress. Am. J. Physiol. Regul. Integr. Comp. Physiol.; 2001; 281, pp. 591-595. [DOI: https://dx.doi.org/10.1152/ajpregu.2001.281.2.R591]

20. Akselrod, S.; Gordon, D.; Ubel, F.A.; Shannon, D.C.; Berger, A.C.; Cohen, R.J. Power spectrum analysis of heart rate fluctuation: A quantitative probe of beat-to-beat cardiovascular control. Science; 1981; 213, pp. 220-222. [DOI: https://dx.doi.org/10.1126/science.6166045]

21. Bernardi, L.; Valle, F.; Coco, M.; Calciati, A.; Sleight, P. Physical activity influences heart rate variability and very-low-frequency components in Holter electrocardiograms. Cardiovasc. Res.; 1996; 32, pp. 234-237. [DOI: https://dx.doi.org/10.1016/0008-6363(96)00081-8]

22. 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. Circulation; 1996; 93, pp. 1043-1065. [DOI: https://dx.doi.org/10.1161/01.CIR.93.5.1043]

23. Li, H.; Kwong, S.; Yang, L.; Huang, D.; Xiao, D. Hilbert-Huang transform for analysis of heart rate variability in cardiac health. IEEE/ACM Trans. Comput. Biol. Bioinform.; 2011; 8, pp. 1557-1567. [DOI: https://dx.doi.org/10.1109/TCBB.2011.43] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21383423]

24. Montano, N.; Porta, A.; Cogliati, C.; Costantino, G.; Tobaldini, E.; Casali, K.R.; Iellamo, F. Heart rate variability explored in the frequency domain: A tool to investigate the link between heart and behavior. Neurosci. Biobehav. Rev.; 2009; 33, pp. 71-80. [DOI: https://dx.doi.org/10.1016/j.neubiorev.2008.07.006] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18706440]

25. Jelinek, H.F.; Huang, Z.Q.; Khandoker, A.H.; Chang, D.; Kiat, H. Cardiac rehabilitation outcomes following a 6-week program of PCI and CABG Patients. Front. Physiol.; 2013; 4, 302. [DOI: https://dx.doi.org/10.3389/fphys.2013.00302]

26. Penzel, T.; Kantelhardt, J.W.; Grote, L.; Peter, J.H.; Bunde, A. Comparison of detrended fluctuation analysis and spectral analysis for heart rate variability in sleep and sleep apnea. IEEE Trans. Biomed. Eng.; 2003; 50, pp. 1143-1151. [DOI: https://dx.doi.org/10.1109/TBME.2003.817636]

27. Yokobori, Y.; Nakane, H.; Uehara, C.; Nagasawa, T.; Mitsuyama, S.; Ohkawa, K.; Kario, K.; Ozawa, S. Temporal relationships among changes in the RR-interval and the powers of the low- and high-frequency components of heart rate variability in normal subjects. Physiol. Rep.; 2023; 11, e15557. [DOI: https://dx.doi.org/10.14814/phy2.15557]

28. Mendez, M.O.; Bianchi, A.M.; Matteucci, M.; Cerutti, S.; Penzel, T. Sleep apnea screening by autoregressive models from a single ECG lead. IEEE Trans. Biomed. Eng.; 2009; 56, pp. 2838-2850. [DOI: https://dx.doi.org/10.1109/TBME.2009.2029563]

29. Kudrynski, K.; Strumillo, P. Real-time estimation of the spectral parameters of heart rate variability. Biocybern. Biomed. Eng.; 2015; 35, pp. 304-316. [DOI: https://dx.doi.org/10.1016/j.bbe.2015.05.002]

30. Polak, A.G.; Klich, B.; Saganowski, S.; Prucnal, M.A.; Kazienko, P. Processing photoplethysmograms recorded by smartwatches to improve the quality of derived pulse rate variability. Sensors; 2022; 22, 7047. [DOI: https://dx.doi.org/10.3390/s22187047]

31. Chan, H.L.; Chou, W.S.; Chen, S.W.; Fang, S.C.; Liou, C.S.; Hwang, Y.S. Continuous and online analysis of heart rate variability. J. Med. Eng. Technol.; 2005; 29, pp. 227-234. [DOI: https://dx.doi.org/10.1080/03091900512331332587]

32. Sornmo, L.; Bailon, R.; Laguna, P. Spectral analysis of heart rate variability in time-varying conditions and in the presence of confounding factors. IEEE Rev. Biomed. Eng.; 2024; 17, pp. 322-341. [DOI: https://dx.doi.org/10.1109/RBME.2022.3220636]

33. Romano, M.; Faiella, G.; Clemente, F.; Iuppariello, L.; Bifulco, P.; Cesarelli, M. Analysis of foetal heart rate variability components by means of empirical mode decomposition. Proceedings of the XIV Mediterranean Conference on Medical and Biological Engineering and Computing 2016; Paphos, Cyprus, 31 March–2 April 2016; [DOI: https://dx.doi.org/10.1007/978-3-319-32703-7_15]

34. Chung, I.Q.; Yu, J.T.; Hu, W.C. Estimating heart rate and respiratory rate from a single lead electrocardiogram using ensemble empirical mode decomposition and spectral data fusion. Sensors; 2021; 21, 1184. [DOI: https://dx.doi.org/10.3390/s21041184] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33567575]

35. Huang, N.E.; Shen, Z.; Long, S.R.; Wu, M.C.; Shih, H.H.; Zheng, Q.; Yen, N.-C.; Tung, C.C.; Liu, H.H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. A.; 1998; 454, pp. 903-995. [DOI: https://dx.doi.org/10.1098/rspa.1998.0193]

36. Shi, M.; He, H.; Geng, W.; Wu, R.; Zhan, C.; Jin, Y.; Zhu, F.; Ren, S.; Shen, B. Early detection of sudden cardiac death by using ensemble empirical mode decomposition-based entropy and classical linear features from heart rate variability signals. Front. Physiol.; 2020; 11, 118. [DOI: https://dx.doi.org/10.3389/fphys.2020.00118] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32158399]

37. Cossul, S.; Andreis, F.R.; Favretto, M.A.; Marques, J.L.B. The Use of Empirical Mode Decomposition on Heart Rate Variability Signals to Assess Autonomic Neuropathy Progression in Type 2 Diabetes. Appl. Sci.; 2023; 13, 7824. [DOI: https://dx.doi.org/10.3390/app13137824]

38. Billman, G.E. Heart rate variability—A historical perspective. Front. Physiol.; 2011; 2, 86. [DOI: https://dx.doi.org/10.3389/fphys.2011.00086]

39. Adamczyk, K.; Polak, A.G. Comparison of multiband filtering, empirical mode decomposition and short-time Fourier transform used to extract physiological components from long-term heart rate variability. Metrol. Meas. Syst.; 2021; 28, pp. 643-660. [DOI: https://dx.doi.org/10.24425/mms.2021.137700]

40. Enamamu, T.; Otebolaku, A.; Marchang, J.; Dany, J. Continuous m-Health Data Authentication Using Wavelet Decomposition for Feature Extraction. Sensors; 2020; 20, 5690. [DOI: https://dx.doi.org/10.3390/s20195690]

41. Deka, D.; Deka, B. Stratification of high-risk hypertensive patients using hybrid heart rate variability features and boosting algorithms. IEEE Access; 2021; 9, pp. 62665-62675. [DOI: https://dx.doi.org/10.1109/ACCESS.2021.3074967]

42. Nesaragi, N.; Sharma, A.; Patidar, S.; Acharya, U.R. Automated diagnosis of coronary artery disease using scalogram-based tensor decomposition with heart rate signals. Med. Eng. Phys.; 2022; 110, 103811. [DOI: https://dx.doi.org/10.1016/j.medengphy.2022.103811]

43. Wang, L.; Hao, J.; Zhou, T.H. ECG Multi-Emotion Recognition Based on Heart Rate Variability Signal Features Mining. Sensors; 2023; 23, 8636. [DOI: https://dx.doi.org/10.3390/s23208636]

44. Dragomiretskiy, K.; Zosso, D. Variational mode decomposition. IEEE Trans. Signal Proces.; 2014; 62, pp. 531-544. [DOI: https://dx.doi.org/10.1109/TSP.2013.2288675]

45. Bae, D.; Matthews, J.J.; Chen, J.J.; Mah, L. Increased exhalation to inhalation ratio during breathing enhances high-frequency heart rate variability in healthy adults. Psychophysiology; 2021; 58, e13905. [DOI: https://dx.doi.org/10.1111/psyp.13905] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34289128]

46. Tsubota, K.; Ohhira, T.; Hashimoto, H. Heart Rate Variability and Body-Movement Extractions using Wi-Fi Channel State Information. Proceedings of the 2023 IEEE/SICE International Symposium on System Integration (SII); Atlanta, GA, USA, 17–20 January 2023; [DOI: https://dx.doi.org/10.1109/SII55687.2023.10039288]

47. Adamczyk, K.; Polak, A.G. Derivation of frequency components from overnight heart rate variability using an adaptive variational mode decomposition. Proceedings of the 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC); Online, Mexico, 1–5 November 2021; [DOI: https://dx.doi.org/10.1109/EMBC46164.2021.9630414]

48. Goldberger, A.L.; Amaral, L.A.; Glass, L.; Hausdorff, J.M.; Ivanov, P.C.; Mark, R.G.; Mietus, J.E.; Moody, G.B.; Peng, C.K.; Stanley, H.E. PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals. Circulation; 2000; 101, pp. 215-220. [DOI: https://dx.doi.org/10.1161/01.CIR.101.23.e215]

49. Pentzel, T.; Moody, G.B.; Mark, R.G.; Goldberger, A.L.; Peter, J.H. The apnea-ECG database. Proceedings of the Computers in Cardiology 2000; Cambridge, MA, USA, 27–27 September 2000; [DOI: https://dx.doi.org/10.1109/CIC.2000.898505]

50. Lian, J.; Liu, Z.; Wang, H.; Dong, X. Adaptive variational mode decomposition method for signal processing based on mode characteristic. Mech. Syst. Signal Process.; 2018; 107, pp. 53-77. [DOI: https://dx.doi.org/10.1016/j.ymssp.2018.01.019]

51. Savitzky, A.; Golay, M.J.E. Smoothing and differentiation of data by simplified least squares procedures. Anal. Chem.; 1964; 36, pp. 1627-1639. [DOI: https://dx.doi.org/10.1021/ac60214a047]

52. Hsu, M.-K.; Sheu, J.-C.; Hsue, C. Overcoming the negative frequencies-instantaneous frequency and amplitude estimation using osculating circle method. J. Mar. Sci. Technol.; 2011; 19, pp. 514-521. [DOI: https://dx.doi.org/10.51400/2709-6998.2165]

53. McElduff, F.; Cortina-Borja, M.; Chan, S.K.; Wade, A. When t-tests or Wilcoxon-Mann-Whitney tests won’t do. Adv. Physiol. Educ.; 2010; 34, pp. 128-133. [DOI: https://dx.doi.org/10.1152/advan.00017.2010]

54. Divine, G.W.; Norton, H.J.; Barón, A.E.; Juarez-Colunga, E. The Wilcoxon–Mann–Whitney procedure fails as a test of medians. Am. Stat.; 2018; 72, pp. 278-286. [DOI: https://dx.doi.org/10.1080/00031305.2017.1305291]

55. Sethna, J.P. Power laws in physics. Nat. Rev. Phys.; 2022; 4, pp. 501-503. [DOI: https://dx.doi.org/10.1038/s42254-022-00491-x]

56. Freed, L.A.; Stein, K.M.; Borer, J.S.; Hochreiter, C.; Supino, P.; Devereux, R.B.; Roman, M.J.; Kligfield, P. Relation of ultra-low frequency heart rate variability to the clinical course of chronic aortic regurgitation. Am. J. Cardiol.; 1997; 79, pp. 1482-1487. [DOI: https://dx.doi.org/10.1016/S0002-9149(97)00175-6]

57. Suda, Y.; Otsuka, K.; Niinami, H.; Ichikawa, S.; Ban, T.; Higashita, R.; Takeuchi, Y. Changes in ultra-low and very low frequency heart rate variability after coronary artery bypass grafting. Biomed. Pharmacother.; 2000; 55, pp. s110-s114. [DOI: https://dx.doi.org/10.1016/S0753-3322(01)90015-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11774857]

58. Wang, Y.; Zhang, C.; Chen, S.; Li, X.; Jin, H.; Du, J. Frequency domain indices of heart rate variability are useful for differentiating vasovagal syncope and postural tachycardia syndrome in children. J. Pediatr.; 2019; 207, pp. 59-63. [DOI: https://dx.doi.org/10.1016/j.jpeds.2018.11.054] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30626483]

59. Hadase, M.; Azuma, A.; Zen, K.; Asada, S.; Kawasaki, T.; Kamitani, T.; Kawasaki, S.; Sugihara, H.; Matsubara, H. Very low frequency power of heart rate variability is a powerful predictor of clinical prognosis in patients with congestive heart failure. Circ. J.; 2004; 68, pp. 343-347. [DOI: https://dx.doi.org/10.1253/circj.68.343] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15056832]

60. Zeid, S.; Buch, G.; Velmeden, D.; Söhne, J.; Schulz, A.; Schuch, A.; Tröbs, S.-O.; Heidorn, M.W.; Müller, F.; Strauch, K. et al. Heart rate variability: Reference values and role for clinical profile and mortality in individuals with heart failure. Clin. Res. Cardiol.; 2024; 113, pp. 1317-1330. [DOI: https://dx.doi.org/10.1007/s00392-023-02248-7]

61. Huang, H.-P.; Wei, S.-Y.; Chao, H.-H.; Hsu, C.F.; Hsu, L.; Chi, S. An investigation study on mode mixing separation in empirical mode decomposition. IEEE Access; 2019; 7, pp. 100684-100691. [DOI: https://dx.doi.org/10.1109/ACCESS.2019.2930543]

62. Zhou, W.; Feng, Z.; Xu, Y.F.; Wang, X.; Lv, H. Empirical Fourier decomposition: An accurate signal decomposition method for nonlinear and non-stationary time series analysis. Mech. Syst. Signal Process.; 2022; 163, 108511. [DOI: https://dx.doi.org/10.1016/j.ymssp.2021.108155]

63. Bazghandi, R.; Hoseintabar Marzebali, M.; Abolghasemi, V.; Hedayati Kia, S. A Novel Mode Un-Mixing Approach in Variational Mode Decomposition for Fault Detection in Wound Rotor Induction Machines. Energies; 2023; 16, 5551. [DOI: https://dx.doi.org/10.3390/en16145551]