Exploration of analytical tools for heart rate variability analysis: Insights into differences

INTRODUCTION

Heart rate variability (HRV) provides insight into the physiological underpinnings of the modulation of the heart. The autonomic nervous system is amongst the major influencers of perturbations in resting heart rate (HR). In addition, respiration, thermoregulation and the neuroendocrine system also have a bearing on it.^[1,2] The variability of the seemingly chaotic HR signal can be computed using multiple mathematical indices that come under the common umbrella of HRV. These indices are categorised as time domain, frequency domain and nonlinear indices. Recommendations for standardisation and computation of these indices were described in the seminal Task Force guidelines, described by the European Society of Cardiology in 1996.^[3]

Analysis of 5-min-long epochs of inter-beat interval data, known as ‘short-term’ HRV, is commonly used in the clinical physiology laboratory for evaluation of cardiac autonomic tone. Furthermore, short-term HRV values have been widely reported in published literature. There have been multiple studies globally in different populations and age groups with an objective to define normative cut-off values for different populations.^[4-8] In addition, multiple studies have been undertaken in our country to describe normative data for different geographical regions.^[9-13] Such data are particularly relevant for the purpose of defining age- and gender-specific cut-offs for different populations to strengthen the use of HRV as a tool for autonomic neuropathy screening.

Common metrics derived from time domain and frequency domain analysis of HRV have been defined to describe the variability of inter-beat interval data. Time domain indices include standard deviation of normal-to-normal intervals (SDNN), root mean square of successive differences (RMSSD) of normal-to-normal intervals, standard deviation (SD) of successive differences, number of normal-to-normal intervals varying by more than 50 ms from preceding intervals (NN50) and their percentage (pNN50%). Frequency domain metrics are derived from Fast Fourier Transform (FFT)/Autoregressive analysis of the inter-beat interval time series and include total power (TP), power in low-frequency (LF) band and high-frequency (HF) band and LF/HF ratio. In addition, frequency domain indices are expressed as normalised power in LF band (LF-nu) and HF band (HF-nu) to negate the effect of very low-frequency (VLF) component, which is primarily respiration driven.

A careful scrutiny of HRV literature reveals a plethora of software solutions for analytical purposes. Popular software tools for HRV measurement include Nevrokard^™,^[14,15] Kubios^™[16] and LabChart^™.^[17,18] We propose that time domain indices (SDNN, RMSSD, nn50 and pNN50) employ simple mathematical calculations and therefore should be similar irrespective of the software used for analysis. However, frequency domain parameters involve resampling of data followed by FFT using different windowing and data padding techniques. Therefore, while time domain indices may be similar across software tools, frequency domain indices should be different and dependent on the analytical tool. However, we observed a paucity of literature comparing popular HRV analysis tools. Therefore, we undertook the present work to compare HRV indices derived using popular analytical tools and explore the differences, if any.

MATERIALS AND METHODS

The protocol was approved by the Institutional Ethics Committee of our institution. Apparently healthy subjects of either sex were recruited. The sample size was calculated using HRV values in the Indian population in a previous study by Rajalakshmi et al.^[19] wherein three HRV software were compared. We assumed a power of 80% and Type I error to be 5% and used a comparison of means approach to calculate the sample size. The sample size obtained was 30.

Written informed consent was obtained after a detailed explanation of the study protocol. Participants were requested to report to the autonomic function laboratory of the department of physiology of our institute, 2 h after a light meal. Abstinence from tea, coffee and nicotine was ensured on the morning of the test. All recordings were done in the forenoon between 9 am and 1 pm. Female subjects were requested to report to the laboratory between days 2 and 6 of the cycle. All study participants were requested to void their bladders and were provided a supine rest of 10 min before commencement of data acquisition.

Disposable adhesive ECG electrodes were applied to the skin after thorough cleaning of the application sites with alcohol-based swabs. The electrodes were connected to a wireless BioNomandix^™ wireless module of Biopac MP150^™ (Biopac Inc., USA) system. Lead II ECG data were acquired in the supine position for 6 min. The subjects were requested to breathe normally and avoid movement to prevent motion artefacts. The investigator ensured that the data acquired was free from artefacts. Data were acquired using Acqknowledge^™ software version 4.4.0 (Biopac Inc., USA) and stored for offline analysis.

We exported a clean noise and artefact-free data segment of 5 min length as Acqknowledge 3^™ (*.acq3) format for the purpose of analysis for every study subject. We utilised Nevrokard aHRV^™ version 13.2.2 (Nevrokard Kiauta, Slovenia), LabChart Pro^™ version 8.1.30 (AD Instruments, Australia) and Kubios Scientific Lite^™ version 4.1.2.1 (Kubios Oy, Finland) for HRV analysis. These software tools are referred to by acronyms ‘N’, ‘L’ and ‘K’ henceforth. ‘N’ and ‘L’ can directly import ‘acq3^™’ files for HRV analysis. However, the free version of ‘K’ used in the study can only import American Standard Code for Information Interchange (ASCII)/text files containing RR intervals. Therefore, RR intervals were extracted from the entire length of ECG data using ‘L’ software and imported into ‘K’ for HRV analysis. All the data segments analysed were exactly 5 min long and were performed by a single observer to prevent inter-observer variability. Visual inspection of the data was done to ensure that all peaks were selected for analysis.

Time domain indices – SDNN, RMSSD and pNN50 – were compared for all three software tools. Similarly, frequency domain indices – TP, LF power, HF power, Normalised LF and HF power (LF-nu and HF-nu respectively) – were compared for the three software. We would like to emphasise that the individual tools provide a wide array of HRV indices, in addition to the aforementioned time domain and frequency domain indices. However, we chose only those indices common to the three software. Furthermore, these metrics are widely reported in HRV research.

Postscripts ‘_K’, ‘_L’ and ‘_N’ were added to the metrics to denote analysis done using Kubios Scientific Lite^™ version 4.1.2.1, LabChart Pro^™ version 8.1.30 and Nevrokard aHRV^™ version 13.2.2 software, respectively.

Values were described as Mean ± SD or Median (interquartile range) based on the distribution, which was assessed using Shapiro–Wilk test. Repeated measures analysis of variance or Friedman’s test with appropriate post hoc tests was used for comparison of the parameters, based on the distribution of parameters. Bland–Altman Analysis was used to evaluate the agreement between the methods, employing a ‘Difference versus Average’ approach. We arbitrarily chose Nevrokard^™ software as the reference, since we observed that it was older than the other tools. In addition, the program has been used in multiple peer-reviewed publications. A P < 0.05 was considered statistically significant. MedCalc^™ Statistical Software version 23.1.2 (MedCalc Software Ltd, Ostend, Belgium; https://www.medcalc.org; 2025)^[20] was used for statistical analysis.

RESULTS

We recruited 53 healthy adults (40 males and 13 females, mean age = 29.89 ± 6.94 years) in the present study. Mean height, weight and BMI of the subjects were 167.74 ± 9.37 cm, 68.33 ± 15.08 kg and 24.23 ± 4.68 kg/m^2, respectively.

While there was statistically significant difference between time domain indices (SDNN, RMSSD and pNN50, P = 0.00052, 0.0066 and 0.021 respectively), median values were very close to each other [Figure 1]. In frequency domain indices, there was no significant difference between TP and LF power (P = 0.121 and 0.090, respectively). However, HF power was significantly different between the software tools (P < 0.00001), with median values exhibiting an increasing trend from HF_K, HF_L and HF_N indices. The values appeared to be within normal range. We also observed statistically significant difference between LF-nu and HFnu indices derived using the software (P < 0.0001 and 0.004 respectively, [Figure 2]).

Close

Figure 1:

Comparison of time domain parameters calculated using different software tools. There was a significant difference between (a) standard deviation of normal-to-normal intervals (SDNN), (b) root mean square of successive differences (RMSSD) and (c) percentage of RR intervals varying by more than 50 ms than preceding intervals (pNN50)estimated using the three tools (P = 0.0005, 0.0066 and 0.021 respectively, Friedman’s test with Dunn’s post hoc test). Parameters represented as box and whisker plots. Individual data points also depicted. ‘_K’, ‘_L’ and ‘_N’ represent parameters derived using Kubios^™, LabChart^™ and Nevrokard^™ software respectively. SDNN: Standard deviation of normal to normal intervals, RMSSD: Root mean square of successive differences, pNN50: percentage of RR intervals varying by more than 50 ms than preceeding intervals. P values from intergroup comparisons are mentioned in-line.

Export to PPT

Close

Figure 2:

Comparison of frequency domain parameters calculated using different software tools. (a) Total Power, (b) and Power in low frequency band (LF Power) were comparable between the software (P = 0.121 and 0.0902, respectively). There was a significant difference between (c) power in High frequency band (HF Power), (d), normalized power in Low frequency band (LFnu), (e) and normalized power in High frequency band (HF-nu) (P < 0.00001, <0.0001 and 0.0040, respectively, Friedman’s test with Dunn’s post hoc test). Parameters represented as box and whisker plots. Individual data points are also depicted. ‘_K’, ‘_L’ and ‘_N’ represent parameters derived using Kubios^™, LabChart^™ and Nevrokard^™ software, respectively. p values from intergroup comparisons are mentioned in-line.

Export to PPT

Bias between the software tools and limits of agreement were analysed using Bland–Altman analysis. Bland–Altman analysis compares the differences between an established and a test method and provides information regarding bias between the parameters. The values are expressed as difference versus arithmetic mean of methods,^[21] difference versus geometric mean or difference versus sample rank-based approaches. As discussed previously, we chose Nevrokard^™ as the reference tool since we found it to be the oldest amongst the three being evaluated. We performed the analysis of difference versus average of the parameter shown in the two tools under comparison. We observed that SDNN was underestimated by ‘K’ and ‘L’ software as compared to NK (Mean difference = −0.13 and −0.20, respectively). Similar trend was observed in RMSSD and pNN50 values (Mean difference = −0.4 and −0.6 and −0.4 and −0.3 respectively, [Figure 3]).

Close

Figure 3:

Comparison of bias between Kubios^™ and LabChart^™ software versus Nevrokard^™ software using Bland–Altman analysis: Time domain indices. Comparison of bias between time domain indices standard deviation of normal-to-normal intervals (SDNN), root mean square of successive differences (RMSSD) and percentage of RR intervals varying by more than 50 ms than preceding intervals (pNN50) calculated using Kubios^™ (_K) and LabChart^™ (_L) software versus Nevrokard^™ (_N) software for heart rate variability analysis. The central solid line represents mean bias, whereas the two flanking dotted lines represent limits of agreement. Whiskers on the limits of agreement depict 95% confidence interval. The dotted line at 0 mark represents the line of equality. SDNN: Standard deviation of normal to normal intervals, RMSSD: Root mean square of successive differences, pNN50: percentage of RR intervals varying by more than 50 ms than preceeding intervals

Export to PPT

Amongst frequency domain indices, ‘K’ and ‘L’ overestimated TP (Mean difference = 77.2 and 205.7 respectively). Both ‘K’ and ‘L’ underestimated LF-power (Mean difference = −32.8 and −2.9 respectively) and HF power (Mean difference = −134.8 and −79.7 respectively). Normalised powers of LF and HF (LF-nu and HF-nu) were overestimated by ‘K’ and ‘L’ software when compared with ‘N’ (mean difference = 5.5 and 4.4 and 2.2 and 0.1, respectively). The mean differences and limits of agreement for the aforementioned parameters are shown in Figures 4 and 5, respectively.

Close

Figure 4:

Comparison of bias between Kubios^™ and LabChart^™ software versus Nevrokard^™ software using Bland–Altman analysis: Frequency domain indices. Comparison of bias between frequency domain indices (Total power, Power in low frequency band (LF) and high frequency band (HF)) using Kubios^™ (_K) and LabChart^™ (_L) software versus Nevrokard^™ (_N) software for heart rate variability analysis. The central solid line represents mean bias, whereas the two flanking dotted lines represent limits of agreement. Whiskers on the limits of agreement depict 95% confidence interval. Dotted line at 0 mark represents the line of equality.

Export to PPT

Close

Figure 5:

Comparison of bias between Kubios^™ and LabChart^™ software versus Nevrokard^™ software using Bland–Altman analysis: Normalised frequency domain indices. Comparison of bias between normalised frequency domain indices (low frequency-nu power and high frequency-nu power) using Kubios^™ (_K) and LabChart^™ (_L) software versus Nevrokard^™ (_N) software for heart rate variability analysis. The central solid line represents mean bias, whereas the two flanking dotted lines represent limits of agreement. Whiskers on the limits of agreement depict 95% confidence interval. Dotted line at 0 mark represents the line of equality. LFNU and HFNU represent normalized powers in low frequency and high frequency bands respectively. SD represents Standard Deviation.

Export to PPT

DISCUSSION

HRV has emerged as an integral component for the evaluation of autonomic neuropathy in various disorders such as Diabetes mellitus, Heart failure, Myocardial Infarction and Epilepsy.^[22-26] It has been shown to serve as a marker of morbidity, mortality and long-term survival in aforementioned disorders. In addition, HRV has found applications in diverse physiological domains such as sports physiology, cognitive physiology and wellness.^[27-29] The rise of wearable technology^[30,31] has invigorated the interest in HRV measurement from the wrist. Furthermore, an upsurge has been observed in ultra-short HRV metrics involving computation of HRV measures from small segments of inter-beat interval data. The segments typically range from 30 s to 2 min.^[32,33]

The ECG signal provides a robust and reliable method for quantification of inter-beat intervals, corresponding to RR intervals in this signal. Pan and Tompkins algorithm is a commonly used QRS detector algorithm used to identify the R waves and compute RR interval data.^[34] The algorithms for processing the ECG signal of the different tools evaluated in the present work are built on the same principle. After the extraction of inter-beat interval data, the tools proceed to calculate time and frequency domain indices. While calculation of time domain indices is relatively straightforward, that of frequency domain indices varies according to different tools. FFT involves resampling of data to account for unevenly placed RR intervals, followed by using a windowing function such as Hanning or Hamming.^[35] This is followed by calculation of spectral powers in well-established frequency bands described previously, i.e. VLF, LF and HF. The frequency bands are classified based on previously well-established works.^[36-39] The commonly defined ranges of power spectral bands are as follows: VLF: 0–0.04 Hz, LF: 0.04–0.15 Hz and HF: 0.15–0.40 Hz.

HRV involves the quantification of the scatter in inter-beat interval data using various metrics. These metrics are clubbed together under time domain and frequency domain indices. Time domain indices include SDNN, RMSSD, NN50 and pNN50. SDNN denotes the SDNNs, while RMSSD denotes the root mean square of the difference between the successive intervals. NN50 is the number of RR intervals varying by more than 50 ms with respect to previous intervals, and pNN50 denotes the percentage of such intervals in the total signal. Frequency domain indices, such as TP, LF-Power and HF-Power, denote the spectral powers in the frequency bands from 0 to 0.04 Hz, 0.04 to 0.15 Hz and 0.15 to 0.4 Hz, respectively. These spectral powers are derived from FFT processing of the time series data. The overall variability in the frequency domain is indicated by TP. A detailed summary of these indices has been discussed in various sources previously.^[40,41]

In the present work, we observed a significant difference between time and frequency domain indices derived using the three tools [Figures 1 and 2]. While time domain indices (SDNN, RMSSD and pNN50) exhibited significant differences between the median values, Bland–Altman analysis demonstrated very small differences between them. This supports our hypothesis that since time domain indices are relatively straightforward to compute; therefore, they are consistent across various software.

The differences emerge when frequency domain parameters are computed using different tools. In the present work, while TP and Low frequency power were statistically comparable, HF power, HF-nu and LF-nu demonstrated statistically significant differences. There was a bias between the aforementioned indices when comparing ‘K’ and ‘L’ with ‘N’ [Figures 4 and 5].

We attempted to investigate the cause of the differences in the values reported by the tools used for assessment. One of the causes of the difference in values could have been the differences in beat detection in the three tools. Therefore, we compared the beat detection methods of the three tools. ‘L’ uses two parameters for defining a normal beat – RR interval and ‘Complexity’ factor, based on the morphology of the ECG wave. QRS detection is based on an algorithm proposed by Hamilton and Tompkins.^[42] The software allows the user to define a range of normal RR intervals and thus ensure selection of maximum peaks. In addition, ‘beat classifier’ view of the software allows changes in the thresholds for RR intervals as well as ‘complexity’ by the user to allow detection of all normal peaks. ‘N’ software allows modifications to baseline at the time of signal import. In addition, the software also allows the user to scroll through the signal and manually insert/removal of missed/spuriously detected R peaks, respectively. Thereafter, the user can opt for ‘deletion’ or ‘digital interpolation’ of the missed peaks. ‘K’ platform allows introduction of ‘beat correction’ algorithm after the signal import, whose threshold may be changed from ‘very low’ to very strong’, enabling optimal detection of all peaks.

In the present work, beat detection was optimal for all the software tools. As discussed previously, the analysis was performed by a single observer to negate inter-observer bias. The observer manually scrolled through the selected ECG data in both ‘L’ and ‘N’ software to ensure that all R peaks were detected. The ECG data segments used in the study were all in Lead II configuration and free from ectopics/artefacts. In addition, the free tier of ‘K’ software used in the present work imported R-R interval data extracted from ‘L’. Therefore, there was little, if any, chance of misdetection of R peaks in ‘K’ software.

‘K’ software has a pre-processing ‘detrending’ algorithm applied to the data at the time of import, switched on by default. [16,43] The ‘detrending’ option automatically removes low-frequency components from the RR interval series. Since it was likely to influence the HRV metrics, we turned it off for all data imports. Furthermore, previous literature has shown that ‘very strong’ threshold correction in ‘K’ software may lead to an unacceptable level of interpolation, and lower levels of threshold correction are recommended.^[44-47] Since our ECG data were free from ectopics and artefacts, we did not require application of threshold correction for any study subject.

The probable reason for differences in frequency domain indices was the underlying analytical methodology adopted by the software. The Scientific Lite^™ version of ‘K’ software uses a Welch’s periodogram for FFT, wherein the window width is 300 s with 50% window overlap.^[16] ‘L’ software uses a Lomb’s periodogram method to perform spectral analysis. In addition, the default value of the upper cut-off of HF band in ‘L’ was 0.45, which was modified to 0.40 by us in the present work to ensure comparability across tools.

A systematic understanding of the different methodologies adopted by commercially available analytical software that may affect HRV parameters is still lacking. In this context, there are very few studies that have compared popular HRV analysis solutions used globally. Rajalakshmi et al.^[19] compared Kubios^™, Nevrokard^™ and HRV Soft^™ software tools in healthy adults using correlation and intra-class correlation coefficient (ICC) as a measure of agreement. However, we believe that ICC may not be a good statistical measure to compare bias between the tools. Therefore, we have used Bland–Altman analysis for comparison and found differences between the values reported by different software.

There are multiple strengths in the present work. Our study compares popular HRV tools that have widespread use in the scientific community. Hence, the findings have widespread applicability. We performed acquisition and analysis of artefact- and ectopic-free data by a single observer; therefore, there is little, if any, likelihood of values being affected by aberrant EGC data or interobserver bias. Based on our data, we may infer that the same data analysed by a single observer using different software tools may lead to different results, especially with reference to frequency domain parameters. We recommend that researchers looking to interpret their HRV data should look for comparable data acquired using the same software tool they are using. Comparison with HRV data acquired and processed using different software may lead to erroneous interpretation of results. This is particularly important in the context of normative data for HRV across different age and population subgroups.^[48-51] Researchers looking to access and interpret these data should carefully sift through the software tool and analytical parameters defined to ensure comparability.

There are some limitations to our present work. We did not include data with artefacts and/or ectopics. Processing of such data may require deletion, interpolation or other necessary measures (‘threshold’ application in ‘K’) and therefore is likely to have a bearing on the agreement between software tools.^[52-54] Another potential limitation of the present work is the use of the free-tier of ‘K’ (Kubios Scientific Lite^™) in the present work, which used exported RR intervals rather than raw data files for analysis. We believe this step bypasses the R wave detection process of the said software and is likely to improve concordance between the tools. Therefore, this is not a potential limitation. For Kubios Scientific Lite^™ software, the default methodology for spectral analysis was the Welch’s periodogram. A subscription-based version of the software would have allowed us to use Lomb’s periodogram for this analytical purpose, which may have improved the agreement between the software. However, we could not obtain access to the premium version due to budget constraints.

We would like to emphasise that the present work did not intend to, nor is sufficient to, establish superiority of one software over another. We have merely compared these tools and attempted to document the differences between them. Therefore, the findings may not be construed as one software tool being better than the other two.

CONCLUSION

Notwithstanding the limitations of the present work, we may infer that frequency domain HRV indices are often different when reported by popular analytical tools. Therefore, interpretation of HRV values should be done in light of the analytical parameters used by the software solution. In addition, we recommend that researchers looking for reference values of HRV should refer to the data analysed and reported by the same platform that they intend to use for their work. We believe such an approach will ensure comparability and robustness of HRV data.