- Email: [email protected]
Wavelet analysis of heart rate variability: Impact of wavelet selection
Biomedical Signal Processing and Control 40 (2018) 220–225
Contents lists available at ScienceDirect
Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc
Research Paper
Wavelet analysis of heart rate variability: Impact of wavelet selection Alexander Tzabazis a,b,∗,1 , Andreas Eisenried a,c,1 , David C. Yeomans a , Hyatt Moore IV d a
Dept. of Anesthesiology, Perioperative and Pain Medicine, Stanford University, School of Medicine, Stanford, CA, USA Dept. of Anesthesiology and Critical Care, Universitätsklinikum Schleswig-Holstein, Campus Lübeck, University of Lübeck, Lübeck, Germany c Dept. of Anesthesiology, Universitätsklinikum Erlangen, Friedrich-Alexander Universität Erlangen-Nürnberg, Erlangen, Germany d Dept. of Mechatronics, Monterey Peninsula College, Monterey, CA, USA b
a r t i c l e
i n f o
Article history: Received 6 March 2017 Received in revised form 12 September 2017 Accepted 30 September 2017 Keywords: Electrocardiography Heart rate variability Autonomic nervous system Wavelet analysis Haar Daubechies Low frequency High frequency
a b s t r a c t Background: Wavelet transform based analysis of heart rate variability is increasingly being used for a wide variety of clinical applications. There is no gold standard as to which wavelet to use and the correlation between results obtained by using different wavelets is unknown. Methods: Heart rate variability in electrocardiograms from healthy volunteers was analyzed using the following wavelets (maximum overlap discrete wavelet packet transform): Haar, Daubechies 2, 4, and 8, least asymmetric Daubechies 4 and 8, and best localized Daubechies 7 using the RHRV package in R. Correlation of power in the different frequency bands (ultra low frequency (ULF), very low frequency (VLF), low frequency (LF), high frequency (HF)) as well as total power and LF:HF ratio were calculated. Bland-Altman comparisons were also made for selected wavelets to test for agreement. Findings: Correlations between results obtained by different wavelets were all statistically significant. Most correlation coefficients were moderate (0.3 ≤ r ≤ 0.7). They were, however, generally lower for the LF:HF ratio, which is commonly used to assess balance of the autonomic nervous system. Conclusion: It is necessary to report which wavelet is used when performing wavelet transform based heart rate variability analysis and depending on whether one is interested in detecting onset or intensity of changes performance of wavelets will vary. © 2017 Published by Elsevier Ltd.
1. Introduction Analysis of heart rate variability (HRV) is increasingly being used [1] as a measure for assessing autonomic nervous system balance. Clinical applications include responses to drug administration [2], objective measurement of pain [3], measuring depth of anesthesia[4], assessing hemorrhagic shock compensation [5], depression [6], sleep staging [7], driver drowsiness [8], and many more. One of the mathematical models used to analyze HRV – the Fourier transform (FT) – has certain limitations that make it a nonideal candidate: it is limited to stationary (or periodic) signals, which assumes that the mean heart rate is constant over time, and lacks time resolution. The short time Fourier transform (STFT) has been used to overcome some of these problems, but it basically represents a compromise between either high time resolution or high frequency resolution [9]. For many of the above-mentioned
∗ Corresponding author at: Department of Anesthesiology and Critical Care, Universitätsklinikum Schleswig-Holstein, Campus Lübeck, University of Lübeck, Ratzeburger Allee 160, 23538 Lübeck, Germany. E-mail address: [email protected] (A. Tzabazis). 1 These authors contributed equally. https://doi.org/10.1016/j.bspc.2017.09.027 1746-8094/© 2017 Published by Elsevier Ltd.
clinical applications a high time resolution is critical in order to detect sudden changes without major delays. At the same time, it is imperative to maintain frequency resolution in order to isolate spectral ranges of interest. The wavelet transform (WT) has been suggested, and increasingly used, to overcome this dilemma in HRV analysis where a temporal localized sliding analysis of the electrocardiogram (ECG) signal is performed. In addition, the shape of the WT-analyzing equation can be chosen to better fit the biomedical signal than the sinusoid shapes used by the FT or STFT. Researchers are frequently offered a wide variety of these wavelets or kernels to choose from: Haar [3], Daubechies [2,6,10–12], biorthogonal [13], Morlet [5], etc. For example, Tejman-Yarden et al. used the Haar wavelet to detect responses to cold pain stimuli in volunteers. They found significant changes in the WT coefficients’ density during increasing, stable and decreasing pain [3]. Deschamps et al. [14] used the Daubechies 4 kernel to investigate changes in the HRV frequency bands and total power in laboring women receiving epidural analgesia. High frequency power increased and the low frequency to high frequency ratio decreased significantly in their study population. However, there is no gold standard as to which kernel should be used for HRV analysis and it is not clear how wavelet selection impacts results of WT based HRV analysis [15], or how comparable results obtained by different wavelets are.
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
This study therefore investigated the correlations between power in the different frequency bands obtained by maximum overlap discrete wavelet packet transform (MODWPT) based HRV analysis using several commonly used kernels in ECGs obtained from healthy volunteers. We hypothesized that wavelet selection does not have a major impact on outcome of MODWPT based HRV analysis secondary to high correlations between results obtained by different wavelets.
2. Methods 2.1. Participants and ECG acquisition After obtaining IRB (Stanford University IRB3, #28503) approval and written, informed consent, five healthy volunteers (two male, three female, age range: 28–33 years) were connected to a Propaq monitor (Welch Alynn, Skaneateles Falls, NY, USA) and a 3-lead ECG was acquired. Volunteers were seated on a comfortable, reclining
221
chair and instructed to hold still in order to minimize movement artifacts. Data was recorded in ASCII format for offline-analysis using the real-time ECG output connector of the Propaq monitor and a custom-built, Atmel ATmega2560-based device with a sampling rate of 500 Hz. In addition, we analyzed ECG obtained in identical fashion as described above from an anesthetized patient that was exposed to a painful stimulus (clamping the head for a craniotomy). After obtaining written informed consent the patient was connected to the standard anesthesia monitor (Intellivue MP60, Philips Healthcare, Andover, MA, USA). Data was recorded in ASCII format for offline-analysis using the real-time ECG output connector of the monitor and a custom-built, Atmel ATmega2560based device with a sampling rate of 500 Hz.
2.2. HRV analysis Raw data was imported to Kubios [16,17](“Kubios HRV premium software 3.0.0”, Finland), artifacts were automatically corrected
Table 1 Spearman rank correlation coefficients for spectral power obtained by wavelet transform in the ultra low frequency (ULF), very low frequency (VLF), low frequency (LF), high frequency (HF) power bands, as well as the total power (TP) and the LF:HF ratio (LFHF). Data presented as average (minimum to maximum range). All correlations were significant (p < 0.01). ULF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.90 (0.88–0.92) 0.72 (0.70–0.76) 0.41 (0.38–0.44) 0.93 (0.91–0.94) 0.94 (0.93–0.95) 0.95 (0.94–0.96)
– 0.88 (0.87–0.90) 0.55 (0.52–0.59) 0.97 (0.97–0.97) 0.86 (0.84–0.87) 0.95 (0.95–0.96)
– – 0.76 (0.74–0.79) 0.82 (0.81–0.85) 0.65 (0.63–0.68) 0.78 (0.77–0.81)
– – – 0.48 (0.46–0.52) 0.34 (0.30–0.37) 0.44 (0.41–0.48)
– – – – 0.92 (0.92–0.93) 0.98 (0.98–0.98)
– – – – – 0.95 (0.95–0.96)
VLF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.17 (0.11–0.23) 0.26 (0.13–0.33) 0.23 (0.05–0.31) 0.39 (0.19–0.52) 0.40 (0.21–0.53) 0.29 (0.19–0.37)
– 0.66 (0.57–0.88) 0.47 (0.41–0.52) 0.42 (0.25–0.97) 0.27 (0.10–0.85) 0.29 (0.07–0.95)
– – 0.75 (0.71–0.80) 0.53 (0.41–0.82) 0.33 (0.21–0.64) 0.19 (0.02–0.41)
– – – 0.52 (0.46–0.58) 0.35 (0.28–0.41) 0.12 (0.02–0.41)
– – – – 0.81 (0.77–0.92) 0.35 (0.16–0.98)
– – – – – 0.53 (0.40–0.95)
LF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.56 (0.48–0.59) 0.63 (0.56–0.68) 0.43 (0.35–0.48) 0.45 (0.37–0.50) 0.40 (0.31–0.48) 0.48 (0.40–0.54)
– 0.51 (0.43–0.57) 0.49 (0.42–0.55) 0.87 (0.85–0.90) 0.63 (0.32–0.75) 0.60 (0.53–0.65)
– – 0.78 (0.75–0.81) 0.53 (0.47–0.58) 0.58 (0.32–0.68) 0.77 (0.70–0.82)
– – – 0.57 (0.52–0.62) 0.60 (0.30–0.75) 0.76 (0.69–0.82)
– – – – 0.76 (0.32–0.88) 0.74 (0.70–0.79)
– – – – – 0.79 (0.32–0.91)
HF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.59 (0.49–0.67) 0.58 (0.51–0.68) 0.48 (0.37–0.59) 0.52 (0.44–0.62) 0.47 (0.38–0.59) 0.45 (0.34–0.58)
– 0.67 (0.56–0.81) 0.60 (0.51–0.72) 0.63 (0.59–0.69) 0.56 (0.50–0.62) 0.48 (0.38–0.61)
– – 0.59 (0.51–0.69) 0.56 (0.50–0.60) 0.52 (0.45–0.59) 0.36 (0.26–0.53)
– – – 0.62 (0.56–0.70) 0.63 (0.55–0.68) 0.52 (0.45–0.64)
– – – – 0.87 (0.85–0.90) 0.50 (0.43–0.58)
– – – – – 0.54 (0.48–0.60)
TP
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.87 (0.86–0.89) 0.80 (0.78–0.82) 0.59 (0.53–0.63) 0.84 (0.82–0.86) 0.84 (0.81–0.86) 0.86 (0.84–0.88)
– 0.85 (0.83–0.87) 0.66 (0.61–0.69) 0.94 (0.92–0.96) 0.85 (0.82–0.87) 0.87 (0.81–0.90)
– – 0.84 (0.82–0.87) 0.81 (0.77–0.83) 0.74 (0.70–0.79) 0.82 (0.76–0.87)
– – – 0.63 (0.56–0.68) 0.55 (0.48–0.61) 0.62 (0.57–0.69)
– – – – 0.94 (0.94–0.95) 0.91 (0.86–0.93)
– – – – – 0.94 (0.89–0.96)
LFHF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.26 (0.17–0.41) 0.42 (0.32–0.55) 0.28 (0.22–0.39) 0.25 (0.17–0.38) 0.21 (0.12–0.35) 0.22 (0.11–0.37)
– 0.44 (0.32–0.57) 0.43 (0.34–0.52) 0.70 (0.67–0.74) 0.55 (0.51–0.61) 0.43 (0.38–0.50)
– – 0.62 (0.58–0.68) 0.41 (0.34–0.49) 0.49 (0.43–0.52) 0.48 (0.42–0.52)
– – – 0.50 (0.45–0.57) 0.58 (0.51–0.63) 0.57 (0.50–0.62)
– – – – 0.83 (0.79–0.87) 0.56 (0.49–0.62)
– – – – – 0.69 (0.67–0.72)
222
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
[18,19], and then visually inspected for potential residual artifacts and accurate R-spike detection. Cumulative R–R peak times obtained from Kubios were imported to R [20] and analyzed using its RHRV package [21]. The package uses an adaptive threshold filter, which was used with default settings, and interpolates the input data at 4 Hz prior to calculating spectral power using MODWPT. The wavelet based spectral method implemented in the RHRV package, calculates the power in specified spectral bands (see below) using a given tolerance to account for the boundaries between bands. The RHRV method is based on MODWPT, but works more efficiently by observing how decomposition nodes at higher levels cover the bands of interest and pruning nodes with redundant or overlapping coverage. This approach reduces computation time and optimizes placement of the wavelet coefficients. See Garcia et al. [15] for a complete description. Power in the ultra-low frequency (ULF, 0–<0.03 Hz), very-low frequency (VLF, 0.03–<0.05 Hz), low frequency (LF, 0.05–<0.15 Hz), and high frequency (HF, 0.15–0.4 Hz) band was calculated, as well as total power (TP) and the LF:HF ratio (LFHF). This analysis was repeated for each volunteer with the following, commonly used, wavelets (RHRV name): Haar (haar), Daubechies 2 (d4), 4 (d8), and 8 (d16), least asymmetric Daubechies 4 (la8), least asymmetric Daubechies 8 (la16) and best localized Daubechies 7 (bl14). Correlation coefficients between wavelets were calculated for each power band, TP and the LFHF for each volunteer. 2.3. Statistics Data is presented as average ± standard deviation (SD) unless indicated otherwise. Data was tested for normal distribution using the D’Agostino and Pearson omnibus normality test [22]. For normally distributed data Pearson correlation coefficients (r) were calculated, while Spearman rank test coefficients (r) were used to determine correlations of non-normally distributed data. Correlation coefficients were categorized, in both cases, as follows: low (<0.3), moderate (≥0.3, ≤0.7), and high (>0.7). Bland-Altman plots were used to test agreement between selected results obtained
from using different wavelets for MODWPT HRV analysis. Data was normalized prior to Bland-Altman analysis by log transformation. All statistical analyses were done using GraphPad Prism 6 for Mac (GraphPad Software, La Jolla, California USA). The significance level was set to 0.01. 3. Results Total ECG length captured was 233 min with an average of 46.6 ± 4.3 min per volunteer. Data quality was excellent with a RHRV beat acceptance rate of 99.7% (range: 98.9-100%). None of the data was normally distributed and therefore Spearman rank correlation coefficients were calculated. Table 1 displays the mean Spearman rank correlation coefficients for spectral power obtained by wavelet transform using different kernels in the ULF, VLF, LF, HF power bands, as well as TP and the LFHF. All correlations were significant at a p level of less than the 0.01 significance level initially set. Average correlations for ULF were very high with only a few exceptions (e.g., d16 vs. la16: r = 0.34). VLF band average correlations were mixed, being generally lower, especially for bl14 vs. d16 (r = 0.12) and haar vs. d4 (r = 0.17), but there were also a few high average correlations, e.g. d16 vs. d8 (r = 0.75) or la16 vs. la8 (r = 0.81). It is interesting to note that there was a wide range for Spearman rank correlations for VLF. Correlations in the LF band were also mixed with some high average correlations (e.g., bl14 vs. la16: r = 0.79, or la8 vs. d4: r = 0.87) and some moderate correlations. It is interesting to note that the correlations between bl14 vs. d8 and bl14 vs. d16 were high, 0.77 and 0.76 respectively, while the VLF correlations between these two kernels were low (0.19 and 0.12, respectively). Average correlations for both HF and TP were in general moderate to high. For the LFHF ratio, a frequently used indicator of autonomic nervous system balance, correlation coefficients for results obtained using different kernels were in general low to moderate with some exceptions. The highest average correlation was found for la16 vs. la8 (r = 0.83) and the lowest average for la16 vs. haar (r = 0.21). Fig. 1 displays the correlation plots between LF, HF and LFHF, respectively, from one volunteer as obtained by la8
Fig. 1. Top panel: Correlation plot of the volunteer with the best correlation between wavelets (least asymmetric Daubechies 4 (la8) and least asymmetric Daubechies 8 (la16) for LFHF ratio (right graph) and corresponding correlation plots for LF (left graph) and HF (middle graph). Bottom panel: Correlation plots of a different volunteer with the worst correlation between wavelets (best localized Daubechies 7 (bl14) and Haar (haar)) for LFHF and corresponding correlation plots for LF and HF.
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
223
Fig. 2. A) Bland-Altman plot for the data presented in Fig. 1, top right panel. Difference between log of least asymmetric 4 (la8) and log of least asymmetric 8 (la16) divided by the average. B) Bland-Altman plot for the data presented in Fig. 1, bottom right panel. Difference between log of best localized Daubechie 7 (bl14) and log of Haar (haar) divided by the average. The dotted horizontal lines in both panels represent the 95% limits of agreement.
and la16 (overall best LFHF correlation and corresponding LF and HF correlation plots, 3 top graphs). The bottom panel of Fig. 1 shows correlation plots from a different volunteer that had the overall worst LFHF correlation between results – obtained by using bl14 and haar – and corresponding LF and HF correlation plots. Fig. 2 displays the Bland-Altman plots for the LFHF ratios for the two volunteers with the overall best (Fig. 2A, r = 0.87 between la8 and la16, corresponding to 3 top graphs in Fig. 1) and worst (Fig. 2B, r = 0.11 between bl14 and haar, corresponding to 3 bottom graphs in Fig. 1) correlation, respectively. Fig. 3 shows the time course of the LFHF ratios as obtained by using different kernels from an anesthetized patient before, during, and after an extremely painful stimulus. 4. Discussion The present work studied for the first time correlations between power in the different frequency bands obtained by maximum overlap discrete wavelet transform (MODWPT) based HRV analysis using several commonly used kernels in ECGs obtained from healthy volunteers. While all of the correlations between the inves-
Fig. 3. Time course for LFHF ratios as obtained by the investigated RHRV kernels before, during, and after exposing an anesthetized patient to an extremely painful stimulus.
tigated wavelets reached statistical significance, it is important to notice that the abundance of data does present a potential concern. With an average recording length of approximately 46 min and 1 value per second, 2760 observation were compared between MODWPT derived results on average. A simple power analysis reveals that a relatively low correlation of 0.09 reaches – conservative – statistical significance of 0.01 with a power of 0.9. Therefore we believe that it is better to use correlation coefficients to make assumptions about clinical usefulness. For example, while the correlation for LFHF between results obtained by using the Haar and the best localized Daubechies 7 wavelet, respectively, was 0.11, yet significant, the low value of the correlation coefficient as well as visual inspection of the correlation plot (Fig. 1, bottom right panel) suggests an unsatisfactory correlation between results. Since even a high correlation coefficient does not necessarily imply good agreement between two analyzing methods, we
224
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
also performed Bland-Altman analysis for selected comparisons. Depending on which kernels were used, agreements were within reasonable limits. For example, the agreement between results obtained by using la8 and la16, respectively, were clinically acceptable with a bias of 0.004 and limits of agreement of +/−0.74 (Fig. 2A). On the contrary, the agreement of results obtained by haar and bl14, respectively, was lower with a bias of 0.086 and limits of agreement ranging from −1.33 to +1.5 (Fig. 2B). Again, one has to keep in mind though that for Bland-Altman analysis the clinical application dictates what limits of agreement are acceptable. Correlation coefficients between LFHF ratios, which are commonly used for assessment of autonomic nervous system balance, were moderate. When correlating the results obtained by the Haar wavelet with all other investigated wavelets, correlations were in general low. Depending on the clinical application, those moderate, and when using the Haar wavelet even low, correlations could lead to misinterpretation of data, if different wavelets have been used for MODWPT based HRV analysis by other groups. We have deliberately chosen relatively short recordings for our study since one of the big advantages of MODWPT based HRV analysis is its better time resolution. For example, while more traditional methods of spectral analysis such as the STFT or Burg’s methods may need a minimum of 32-s windows to determine fluctuation of power in these bands, MODWPT analysis is able to do so with 16-s windows. While our data does not allow us to give a generalized recommendation as to which kernel to use for MODWPT based HRV analysis for a specific application, it emphasizes the importance of reporting the kernel used in order to make results from different groups more comparable. Garcia et al. [15] found that the mother wavelet used in MODWPT based HRV analysis influences the timefrequency resolution because it determines the filter shape: shorter wavelets have greater time resolution than the larger wavelets, whereas the latter have a better filtering behavior than the first ones [15]. It would be desirable to standardize MODWPT based HRV analysis by selecting a “state of the art” kernel. A similar approach has been used by Jonckheere et al. [23] for analysis of surface electromyography along the paravertebral muscles. In this study the authors used a trial and error approach to determine, which specific wavelet yielded the best distinction between “noise” and the actual signal of interest. This group found the Daubechies DB3 wavelet to be the best choice. It is not clear though, whether those findings can be directly transferred to HRV analysis. In addition, depending on the clinical application, different kernels might be better suited for HRV data analysis than others. It is unclear whether other analysis methods, such as Poincare plots [24] could be useful in determining which kernel to use for further wavelet based analysis. The data presented from the anesthetized patient in Fig. 3 corroborates this hypothesis as the performance of kernels varies depending on whether one is interested in detecting onset or intensity of changes. Fig. 3 suggests that MODWPT using the Haar kernel indicates nociceptive response rapidly after stimulus onset, which is in line with the findings of Tejman-Yarden et al. [3]. A benefit of looking at data while stimulating an anesthetized patient is that potential confounders such as respiratory arrhythmia, anxiety, and others can be excluded. It would be beneficial to agree on a standard terminology for abbreviating wavelet names. For example, the RHRV package nonintuitively uses the abbreviation “d4” for the Daubechies 2 wavelet, which is commonly abbreviated as “DB2” or “db2”, e.g. in MATLAB. Additional studies are warranted to investigate potential advantages and disadvantages of the wide variety of kernels that have been used with MODWPT based HRV analysis. While correlations between results obtained by different wavelets were significant for any of the wavelet comparisons investigated in our study, low,
i.e. clinically insignificant, correlation coefficients and/or visual inspection of correlation plots might be a better indicator for unsatisfactory correlation. Especially LFHF results obtained by MODWPT based HRV analysis using different wavelets, which are commonly used to assess sympathovagal balance need to be interpreted with caution. Conflict of interest The authors have no conflicts of interest or financial disclosures to report. References [1] V. Pichot, F. Roche, S. Celle, J.-C. Barthélémy, F. Chouchou, HRVanalysis: a free software for analyzing cardiac autonomic activity, Front. Physiol. [Internet] 7 (November) (2016) 557 [cited 2017 Jul 6]. Available from: http://www.ncbi. nlm.nih.gov/pubmed/27920726. [2] V. Pichot, J.M. Gaspoz, S. Molliex, A. Antoniadis, T. Busso, F. Roche, et al., Wavelet transform to quantify heart rate variability and to assess its instantaneous changes, J. Appl. Physiol. [Internet] 86 (March (3)) (1999) 1081–1091 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/ pubmed/10066727. [3] S. Tejman-Yarden, O. Levi, A. Beizerov, Y. Parmet, T. Nguyen, M. Saunders, et al., Heart rate analysis by sparse representation for acute pain detection, Med. Biol. Eng. Comput. [Internet] 54 (April (4)) (2016) 595–606 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/26264057. [4] D.A. Kenwright, A. Bernjak, T. Draegni, S. Dzeroski, M. Entwistle, M. Horvat, et al., The discriminatory value of cardiorespiratory interactions in distinguishing awake from anaesthetised states: a randomised observational study, Anaesthesia [Internet] 70 (December (12)) (2015) 1356–1368 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/ 26350998. [5] C.G. Scully, G.C. Kramer, D.G. Strauss, Evaluation of heart rate and blood pressure variability as indicators of physiological compensation to hemorrhage before shock, Shock [Internet] 43 (May (5)) (2015) 463–469 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/ 25692248. [6] S. Akar, S. Kara, V. Bilgic¸, Investigation of heart rate variability in major depression patients using wavelet packet transform, Psychiatry Res. 30 (2016) 326–332. [7] F. Ebrahimi, Setarehdan S-KS, J. Ayala-Moyeda, H. Nazeran, Automatic sleep staging using empirical mode decomposition, discrete wavelet transform, time-domain, and nonlinear dynamics features of heart rate variability signals, Comput. Methods Programs Biomed. [Internet] 112 (October (1)) (2013) 47–57 [cited 2016 Nov 28]. Available from: http://www.ncbi.nlm.nih. gov/pubmed/23895941. [8] G. Li, W.-Y. Chung, Detection of driver drowsiness using wavelet analysis of heart rate variability and a support vector machine classifier, Sensors (Basel) [Internet] 13 (12) (2013) 16494–16511 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/24316564. [9] A.V. Oppenheim, R.W. Schafer, Discrete-time Signal Processing, Pearson Higher Education, 2010. [10] J.L. Ducla-Soares, M. Santos-Bento, S. Laranjo, A. Andrade, E. Ducla-Soares, J.P. Boto, et al., Wavelet analysis of autonomic outflow of normal subjects on head-up tilt, cold pressor test, Valsalva manoeuvre and deep breathing, Exp. Physiol. [Internet] 92 (July (4)) (2007) 677–686 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/17468200. [11] F. Ebrahimi, S. Setarehdan, J. Ayala-Moyeda, H. Nazeran, Automatic sleep staging using empirical mode decomposition, discrete wavelet transform, time-domain, and nonlinear dynamics features of heart rate variability signals, Comput. Methods Programs Biomed. 112 (1) (2013) 47–57. [12] F. Chouchou, V. Pichot, C. Perchet, V. Legrain, L. Garcia-Larrea, F. Roche, et al., Autonomic pain responses during sleep: a study of heart rate variability, Eur. J. Pain 15 (6) (2011) 554–560. [13] A. Hossen, A. Barhoum, D. Jaju, V. Gowri, K. Al-Hashmi, M.O. Hassan, et al., Identification of patients with preeclampsia from normal subjects using wavelet-based spectral analysis of heart rate variability, Technol. Health Care [Internet] 25 (August (4)) (2017) 641–649 [cited 2017 Sep 12]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/28436399. [14] A. Deschamps, I. Kaufman, S.B. Backman, G. Plourde, Autonomic nervous system response to epidural analgesia in laboring patients by wavelet transform of heart rate and blood pressure variability, Anesthesiology [Internet] 101 (July (1)) (2004) 21–27 [cited 2017 Jul 6]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/15220767. [15] C.A. García, A. Otero, X. Vila, D.G. Márquez, A new algorithm for wavelet-based heart rate variability analysis, Biomed. Signal. Process. Control [Internet] 8 (November (6)) (2013) 542–550 [cited 2017 Mar 1]. Available from: http://linkinghub.elsevier.com/retrieve/pii/S1746809413000803. [16] Kubios HRV Software 2.2 [Internet], 2015, Available from: http://kubios.uef.fi. [17] M.P. Tarvainen, J.-P. Niskanen, J.A. Lipponen, P.O. Ranta-Aho, P.A. Karjalainen, G. Berntson, et al., Kubios HRV–heart rate variability analysis software,
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
[18]
[19]
[20]
[21]
Comput. Methods Programs Biomed. [Internet]. Elsevier 113 (1) (2014) 210–220 [cited 2016 Jul 20]. Available from: http://www.ncbi.nlm.nih.gov/ pubmed/24054542. M.P. Tarvainen, J. Lipponen, J.-P. Niskanen, P.O. Ranta-aho, Kubios HRV User’s Guide [Internet], 2017 [cited 2017 Jul 13]. Available from: http://www.kubios. com/downloads/Kubios HRV Users Guide.pdf. M.P. Tarvainen, P.O. Ranta-Aho, P.A. Karjalainen, An advanced detrending method with application to HRV analysis, IEEE Trans. Biomed. Eng. [Internet] 49 (February (2)) (2002) 172–175 [cited 2017 Jul 13]. Available from: http:// ieeexplore.ieee.org/document/979357/. R Core Team, R: A Language and Environment for Statistical Computing. [Internet], R Foundation for Statistical Computing, Vienna, Austria, 2013, Available from: http://www.r-project.org/. ˜ L. Rodríguez-Linares, A. Méndez, M. Lado, D. Olivieri, X. Vila, I. Gómez-Conde,
225
An open source tool for heart rate variability spectral analysis, Comput. Methods Programs Biomed. 103 (1) (2011) 39–50. [22] R.B. D’Agostino, Tests for normal distribution, in: R.B. D’Agostino, M.A. Stephens (Eds.), Goodness-Of-Fit Techniques, Marcel Dekker, 1986. [23] E. Jonckheere, P. Lohsoonthorn, S. Musuvathy, V. Mahajan, M. Stefanovic, On a standing wave central pattern generator and the coherence problem, Biomed. Signal. Process. Control [Internet] 5 (2010) 336–347 [cited 2017 Jun 30]. Available from: http://ac.els-cdn.com/S1746809410000273/1-s2.0S1746809410000273-main.pdf? tid=a0214b5a-5dcf-11e7-a71a00000aacb361&acdnat=1498853379 9e8e2d3aaabea471b645125a0a6f8303. [24] H.D. Esperer, C. Esperer, R.J. Cohen, Cardiac arrhythmias imprint specific signatures on Lorenz plots, Ann. Noninvasive Electrocardiol. [Internet] 13 (January (1)) (2008) 44–60 [cited 2017 Jun 30]. Available from: http://doi. wiley.com/10.1111/j.1542-474X.2007.00200.x.
Contents lists available at ScienceDirect
Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc
Research Paper
Wavelet analysis of heart rate variability: Impact of wavelet selection Alexander Tzabazis a,b,∗,1 , Andreas Eisenried a,c,1 , David C. Yeomans a , Hyatt Moore IV d a
Dept. of Anesthesiology, Perioperative and Pain Medicine, Stanford University, School of Medicine, Stanford, CA, USA Dept. of Anesthesiology and Critical Care, Universitätsklinikum Schleswig-Holstein, Campus Lübeck, University of Lübeck, Lübeck, Germany c Dept. of Anesthesiology, Universitätsklinikum Erlangen, Friedrich-Alexander Universität Erlangen-Nürnberg, Erlangen, Germany d Dept. of Mechatronics, Monterey Peninsula College, Monterey, CA, USA b
a r t i c l e
i n f o
Article history: Received 6 March 2017 Received in revised form 12 September 2017 Accepted 30 September 2017 Keywords: Electrocardiography Heart rate variability Autonomic nervous system Wavelet analysis Haar Daubechies Low frequency High frequency
a b s t r a c t Background: Wavelet transform based analysis of heart rate variability is increasingly being used for a wide variety of clinical applications. There is no gold standard as to which wavelet to use and the correlation between results obtained by using different wavelets is unknown. Methods: Heart rate variability in electrocardiograms from healthy volunteers was analyzed using the following wavelets (maximum overlap discrete wavelet packet transform): Haar, Daubechies 2, 4, and 8, least asymmetric Daubechies 4 and 8, and best localized Daubechies 7 using the RHRV package in R. Correlation of power in the different frequency bands (ultra low frequency (ULF), very low frequency (VLF), low frequency (LF), high frequency (HF)) as well as total power and LF:HF ratio were calculated. Bland-Altman comparisons were also made for selected wavelets to test for agreement. Findings: Correlations between results obtained by different wavelets were all statistically significant. Most correlation coefficients were moderate (0.3 ≤ r ≤ 0.7). They were, however, generally lower for the LF:HF ratio, which is commonly used to assess balance of the autonomic nervous system. Conclusion: It is necessary to report which wavelet is used when performing wavelet transform based heart rate variability analysis and depending on whether one is interested in detecting onset or intensity of changes performance of wavelets will vary. © 2017 Published by Elsevier Ltd.
1. Introduction Analysis of heart rate variability (HRV) is increasingly being used [1] as a measure for assessing autonomic nervous system balance. Clinical applications include responses to drug administration [2], objective measurement of pain [3], measuring depth of anesthesia[4], assessing hemorrhagic shock compensation [5], depression [6], sleep staging [7], driver drowsiness [8], and many more. One of the mathematical models used to analyze HRV – the Fourier transform (FT) – has certain limitations that make it a nonideal candidate: it is limited to stationary (or periodic) signals, which assumes that the mean heart rate is constant over time, and lacks time resolution. The short time Fourier transform (STFT) has been used to overcome some of these problems, but it basically represents a compromise between either high time resolution or high frequency resolution [9]. For many of the above-mentioned
∗ Corresponding author at: Department of Anesthesiology and Critical Care, Universitätsklinikum Schleswig-Holstein, Campus Lübeck, University of Lübeck, Ratzeburger Allee 160, 23538 Lübeck, Germany. E-mail address: [email protected] (A. Tzabazis). 1 These authors contributed equally. https://doi.org/10.1016/j.bspc.2017.09.027 1746-8094/© 2017 Published by Elsevier Ltd.
clinical applications a high time resolution is critical in order to detect sudden changes without major delays. At the same time, it is imperative to maintain frequency resolution in order to isolate spectral ranges of interest. The wavelet transform (WT) has been suggested, and increasingly used, to overcome this dilemma in HRV analysis where a temporal localized sliding analysis of the electrocardiogram (ECG) signal is performed. In addition, the shape of the WT-analyzing equation can be chosen to better fit the biomedical signal than the sinusoid shapes used by the FT or STFT. Researchers are frequently offered a wide variety of these wavelets or kernels to choose from: Haar [3], Daubechies [2,6,10–12], biorthogonal [13], Morlet [5], etc. For example, Tejman-Yarden et al. used the Haar wavelet to detect responses to cold pain stimuli in volunteers. They found significant changes in the WT coefficients’ density during increasing, stable and decreasing pain [3]. Deschamps et al. [14] used the Daubechies 4 kernel to investigate changes in the HRV frequency bands and total power in laboring women receiving epidural analgesia. High frequency power increased and the low frequency to high frequency ratio decreased significantly in their study population. However, there is no gold standard as to which kernel should be used for HRV analysis and it is not clear how wavelet selection impacts results of WT based HRV analysis [15], or how comparable results obtained by different wavelets are.
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
This study therefore investigated the correlations between power in the different frequency bands obtained by maximum overlap discrete wavelet packet transform (MODWPT) based HRV analysis using several commonly used kernels in ECGs obtained from healthy volunteers. We hypothesized that wavelet selection does not have a major impact on outcome of MODWPT based HRV analysis secondary to high correlations between results obtained by different wavelets.
2. Methods 2.1. Participants and ECG acquisition After obtaining IRB (Stanford University IRB3, #28503) approval and written, informed consent, five healthy volunteers (two male, three female, age range: 28–33 years) were connected to a Propaq monitor (Welch Alynn, Skaneateles Falls, NY, USA) and a 3-lead ECG was acquired. Volunteers were seated on a comfortable, reclining
221
chair and instructed to hold still in order to minimize movement artifacts. Data was recorded in ASCII format for offline-analysis using the real-time ECG output connector of the Propaq monitor and a custom-built, Atmel ATmega2560-based device with a sampling rate of 500 Hz. In addition, we analyzed ECG obtained in identical fashion as described above from an anesthetized patient that was exposed to a painful stimulus (clamping the head for a craniotomy). After obtaining written informed consent the patient was connected to the standard anesthesia monitor (Intellivue MP60, Philips Healthcare, Andover, MA, USA). Data was recorded in ASCII format for offline-analysis using the real-time ECG output connector of the monitor and a custom-built, Atmel ATmega2560based device with a sampling rate of 500 Hz.
2.2. HRV analysis Raw data was imported to Kubios [16,17](“Kubios HRV premium software 3.0.0”, Finland), artifacts were automatically corrected
Table 1 Spearman rank correlation coefficients for spectral power obtained by wavelet transform in the ultra low frequency (ULF), very low frequency (VLF), low frequency (LF), high frequency (HF) power bands, as well as the total power (TP) and the LF:HF ratio (LFHF). Data presented as average (minimum to maximum range). All correlations were significant (p < 0.01). ULF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.90 (0.88–0.92) 0.72 (0.70–0.76) 0.41 (0.38–0.44) 0.93 (0.91–0.94) 0.94 (0.93–0.95) 0.95 (0.94–0.96)
– 0.88 (0.87–0.90) 0.55 (0.52–0.59) 0.97 (0.97–0.97) 0.86 (0.84–0.87) 0.95 (0.95–0.96)
– – 0.76 (0.74–0.79) 0.82 (0.81–0.85) 0.65 (0.63–0.68) 0.78 (0.77–0.81)
– – – 0.48 (0.46–0.52) 0.34 (0.30–0.37) 0.44 (0.41–0.48)
– – – – 0.92 (0.92–0.93) 0.98 (0.98–0.98)
– – – – – 0.95 (0.95–0.96)
VLF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.17 (0.11–0.23) 0.26 (0.13–0.33) 0.23 (0.05–0.31) 0.39 (0.19–0.52) 0.40 (0.21–0.53) 0.29 (0.19–0.37)
– 0.66 (0.57–0.88) 0.47 (0.41–0.52) 0.42 (0.25–0.97) 0.27 (0.10–0.85) 0.29 (0.07–0.95)
– – 0.75 (0.71–0.80) 0.53 (0.41–0.82) 0.33 (0.21–0.64) 0.19 (0.02–0.41)
– – – 0.52 (0.46–0.58) 0.35 (0.28–0.41) 0.12 (0.02–0.41)
– – – – 0.81 (0.77–0.92) 0.35 (0.16–0.98)
– – – – – 0.53 (0.40–0.95)
LF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.56 (0.48–0.59) 0.63 (0.56–0.68) 0.43 (0.35–0.48) 0.45 (0.37–0.50) 0.40 (0.31–0.48) 0.48 (0.40–0.54)
– 0.51 (0.43–0.57) 0.49 (0.42–0.55) 0.87 (0.85–0.90) 0.63 (0.32–0.75) 0.60 (0.53–0.65)
– – 0.78 (0.75–0.81) 0.53 (0.47–0.58) 0.58 (0.32–0.68) 0.77 (0.70–0.82)
– – – 0.57 (0.52–0.62) 0.60 (0.30–0.75) 0.76 (0.69–0.82)
– – – – 0.76 (0.32–0.88) 0.74 (0.70–0.79)
– – – – – 0.79 (0.32–0.91)
HF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.59 (0.49–0.67) 0.58 (0.51–0.68) 0.48 (0.37–0.59) 0.52 (0.44–0.62) 0.47 (0.38–0.59) 0.45 (0.34–0.58)
– 0.67 (0.56–0.81) 0.60 (0.51–0.72) 0.63 (0.59–0.69) 0.56 (0.50–0.62) 0.48 (0.38–0.61)
– – 0.59 (0.51–0.69) 0.56 (0.50–0.60) 0.52 (0.45–0.59) 0.36 (0.26–0.53)
– – – 0.62 (0.56–0.70) 0.63 (0.55–0.68) 0.52 (0.45–0.64)
– – – – 0.87 (0.85–0.90) 0.50 (0.43–0.58)
– – – – – 0.54 (0.48–0.60)
TP
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.87 (0.86–0.89) 0.80 (0.78–0.82) 0.59 (0.53–0.63) 0.84 (0.82–0.86) 0.84 (0.81–0.86) 0.86 (0.84–0.88)
– 0.85 (0.83–0.87) 0.66 (0.61–0.69) 0.94 (0.92–0.96) 0.85 (0.82–0.87) 0.87 (0.81–0.90)
– – 0.84 (0.82–0.87) 0.81 (0.77–0.83) 0.74 (0.70–0.79) 0.82 (0.76–0.87)
– – – 0.63 (0.56–0.68) 0.55 (0.48–0.61) 0.62 (0.57–0.69)
– – – – 0.94 (0.94–0.95) 0.91 (0.86–0.93)
– – – – – 0.94 (0.89–0.96)
LFHF
haar
d4
d8
d16
la8
la16
d4 d8 d16 la8 la16 bl14
0.26 (0.17–0.41) 0.42 (0.32–0.55) 0.28 (0.22–0.39) 0.25 (0.17–0.38) 0.21 (0.12–0.35) 0.22 (0.11–0.37)
– 0.44 (0.32–0.57) 0.43 (0.34–0.52) 0.70 (0.67–0.74) 0.55 (0.51–0.61) 0.43 (0.38–0.50)
– – 0.62 (0.58–0.68) 0.41 (0.34–0.49) 0.49 (0.43–0.52) 0.48 (0.42–0.52)
– – – 0.50 (0.45–0.57) 0.58 (0.51–0.63) 0.57 (0.50–0.62)
– – – – 0.83 (0.79–0.87) 0.56 (0.49–0.62)
– – – – – 0.69 (0.67–0.72)
222
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
[18,19], and then visually inspected for potential residual artifacts and accurate R-spike detection. Cumulative R–R peak times obtained from Kubios were imported to R [20] and analyzed using its RHRV package [21]. The package uses an adaptive threshold filter, which was used with default settings, and interpolates the input data at 4 Hz prior to calculating spectral power using MODWPT. The wavelet based spectral method implemented in the RHRV package, calculates the power in specified spectral bands (see below) using a given tolerance to account for the boundaries between bands. The RHRV method is based on MODWPT, but works more efficiently by observing how decomposition nodes at higher levels cover the bands of interest and pruning nodes with redundant or overlapping coverage. This approach reduces computation time and optimizes placement of the wavelet coefficients. See Garcia et al. [15] for a complete description. Power in the ultra-low frequency (ULF, 0–<0.03 Hz), very-low frequency (VLF, 0.03–<0.05 Hz), low frequency (LF, 0.05–<0.15 Hz), and high frequency (HF, 0.15–0.4 Hz) band was calculated, as well as total power (TP) and the LF:HF ratio (LFHF). This analysis was repeated for each volunteer with the following, commonly used, wavelets (RHRV name): Haar (haar), Daubechies 2 (d4), 4 (d8), and 8 (d16), least asymmetric Daubechies 4 (la8), least asymmetric Daubechies 8 (la16) and best localized Daubechies 7 (bl14). Correlation coefficients between wavelets were calculated for each power band, TP and the LFHF for each volunteer. 2.3. Statistics Data is presented as average ± standard deviation (SD) unless indicated otherwise. Data was tested for normal distribution using the D’Agostino and Pearson omnibus normality test [22]. For normally distributed data Pearson correlation coefficients (r) were calculated, while Spearman rank test coefficients (r) were used to determine correlations of non-normally distributed data. Correlation coefficients were categorized, in both cases, as follows: low (<0.3), moderate (≥0.3, ≤0.7), and high (>0.7). Bland-Altman plots were used to test agreement between selected results obtained
from using different wavelets for MODWPT HRV analysis. Data was normalized prior to Bland-Altman analysis by log transformation. All statistical analyses were done using GraphPad Prism 6 for Mac (GraphPad Software, La Jolla, California USA). The significance level was set to 0.01. 3. Results Total ECG length captured was 233 min with an average of 46.6 ± 4.3 min per volunteer. Data quality was excellent with a RHRV beat acceptance rate of 99.7% (range: 98.9-100%). None of the data was normally distributed and therefore Spearman rank correlation coefficients were calculated. Table 1 displays the mean Spearman rank correlation coefficients for spectral power obtained by wavelet transform using different kernels in the ULF, VLF, LF, HF power bands, as well as TP and the LFHF. All correlations were significant at a p level of less than the 0.01 significance level initially set. Average correlations for ULF were very high with only a few exceptions (e.g., d16 vs. la16: r = 0.34). VLF band average correlations were mixed, being generally lower, especially for bl14 vs. d16 (r = 0.12) and haar vs. d4 (r = 0.17), but there were also a few high average correlations, e.g. d16 vs. d8 (r = 0.75) or la16 vs. la8 (r = 0.81). It is interesting to note that there was a wide range for Spearman rank correlations for VLF. Correlations in the LF band were also mixed with some high average correlations (e.g., bl14 vs. la16: r = 0.79, or la8 vs. d4: r = 0.87) and some moderate correlations. It is interesting to note that the correlations between bl14 vs. d8 and bl14 vs. d16 were high, 0.77 and 0.76 respectively, while the VLF correlations between these two kernels were low (0.19 and 0.12, respectively). Average correlations for both HF and TP were in general moderate to high. For the LFHF ratio, a frequently used indicator of autonomic nervous system balance, correlation coefficients for results obtained using different kernels were in general low to moderate with some exceptions. The highest average correlation was found for la16 vs. la8 (r = 0.83) and the lowest average for la16 vs. haar (r = 0.21). Fig. 1 displays the correlation plots between LF, HF and LFHF, respectively, from one volunteer as obtained by la8
Fig. 1. Top panel: Correlation plot of the volunteer with the best correlation between wavelets (least asymmetric Daubechies 4 (la8) and least asymmetric Daubechies 8 (la16) for LFHF ratio (right graph) and corresponding correlation plots for LF (left graph) and HF (middle graph). Bottom panel: Correlation plots of a different volunteer with the worst correlation between wavelets (best localized Daubechies 7 (bl14) and Haar (haar)) for LFHF and corresponding correlation plots for LF and HF.
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
223
Fig. 2. A) Bland-Altman plot for the data presented in Fig. 1, top right panel. Difference between log of least asymmetric 4 (la8) and log of least asymmetric 8 (la16) divided by the average. B) Bland-Altman plot for the data presented in Fig. 1, bottom right panel. Difference between log of best localized Daubechie 7 (bl14) and log of Haar (haar) divided by the average. The dotted horizontal lines in both panels represent the 95% limits of agreement.
and la16 (overall best LFHF correlation and corresponding LF and HF correlation plots, 3 top graphs). The bottom panel of Fig. 1 shows correlation plots from a different volunteer that had the overall worst LFHF correlation between results – obtained by using bl14 and haar – and corresponding LF and HF correlation plots. Fig. 2 displays the Bland-Altman plots for the LFHF ratios for the two volunteers with the overall best (Fig. 2A, r = 0.87 between la8 and la16, corresponding to 3 top graphs in Fig. 1) and worst (Fig. 2B, r = 0.11 between bl14 and haar, corresponding to 3 bottom graphs in Fig. 1) correlation, respectively. Fig. 3 shows the time course of the LFHF ratios as obtained by using different kernels from an anesthetized patient before, during, and after an extremely painful stimulus. 4. Discussion The present work studied for the first time correlations between power in the different frequency bands obtained by maximum overlap discrete wavelet transform (MODWPT) based HRV analysis using several commonly used kernels in ECGs obtained from healthy volunteers. While all of the correlations between the inves-
Fig. 3. Time course for LFHF ratios as obtained by the investigated RHRV kernels before, during, and after exposing an anesthetized patient to an extremely painful stimulus.
tigated wavelets reached statistical significance, it is important to notice that the abundance of data does present a potential concern. With an average recording length of approximately 46 min and 1 value per second, 2760 observation were compared between MODWPT derived results on average. A simple power analysis reveals that a relatively low correlation of 0.09 reaches – conservative – statistical significance of 0.01 with a power of 0.9. Therefore we believe that it is better to use correlation coefficients to make assumptions about clinical usefulness. For example, while the correlation for LFHF between results obtained by using the Haar and the best localized Daubechies 7 wavelet, respectively, was 0.11, yet significant, the low value of the correlation coefficient as well as visual inspection of the correlation plot (Fig. 1, bottom right panel) suggests an unsatisfactory correlation between results. Since even a high correlation coefficient does not necessarily imply good agreement between two analyzing methods, we
224
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
also performed Bland-Altman analysis for selected comparisons. Depending on which kernels were used, agreements were within reasonable limits. For example, the agreement between results obtained by using la8 and la16, respectively, were clinically acceptable with a bias of 0.004 and limits of agreement of +/−0.74 (Fig. 2A). On the contrary, the agreement of results obtained by haar and bl14, respectively, was lower with a bias of 0.086 and limits of agreement ranging from −1.33 to +1.5 (Fig. 2B). Again, one has to keep in mind though that for Bland-Altman analysis the clinical application dictates what limits of agreement are acceptable. Correlation coefficients between LFHF ratios, which are commonly used for assessment of autonomic nervous system balance, were moderate. When correlating the results obtained by the Haar wavelet with all other investigated wavelets, correlations were in general low. Depending on the clinical application, those moderate, and when using the Haar wavelet even low, correlations could lead to misinterpretation of data, if different wavelets have been used for MODWPT based HRV analysis by other groups. We have deliberately chosen relatively short recordings for our study since one of the big advantages of MODWPT based HRV analysis is its better time resolution. For example, while more traditional methods of spectral analysis such as the STFT or Burg’s methods may need a minimum of 32-s windows to determine fluctuation of power in these bands, MODWPT analysis is able to do so with 16-s windows. While our data does not allow us to give a generalized recommendation as to which kernel to use for MODWPT based HRV analysis for a specific application, it emphasizes the importance of reporting the kernel used in order to make results from different groups more comparable. Garcia et al. [15] found that the mother wavelet used in MODWPT based HRV analysis influences the timefrequency resolution because it determines the filter shape: shorter wavelets have greater time resolution than the larger wavelets, whereas the latter have a better filtering behavior than the first ones [15]. It would be desirable to standardize MODWPT based HRV analysis by selecting a “state of the art” kernel. A similar approach has been used by Jonckheere et al. [23] for analysis of surface electromyography along the paravertebral muscles. In this study the authors used a trial and error approach to determine, which specific wavelet yielded the best distinction between “noise” and the actual signal of interest. This group found the Daubechies DB3 wavelet to be the best choice. It is not clear though, whether those findings can be directly transferred to HRV analysis. In addition, depending on the clinical application, different kernels might be better suited for HRV data analysis than others. It is unclear whether other analysis methods, such as Poincare plots [24] could be useful in determining which kernel to use for further wavelet based analysis. The data presented from the anesthetized patient in Fig. 3 corroborates this hypothesis as the performance of kernels varies depending on whether one is interested in detecting onset or intensity of changes. Fig. 3 suggests that MODWPT using the Haar kernel indicates nociceptive response rapidly after stimulus onset, which is in line with the findings of Tejman-Yarden et al. [3]. A benefit of looking at data while stimulating an anesthetized patient is that potential confounders such as respiratory arrhythmia, anxiety, and others can be excluded. It would be beneficial to agree on a standard terminology for abbreviating wavelet names. For example, the RHRV package nonintuitively uses the abbreviation “d4” for the Daubechies 2 wavelet, which is commonly abbreviated as “DB2” or “db2”, e.g. in MATLAB. Additional studies are warranted to investigate potential advantages and disadvantages of the wide variety of kernels that have been used with MODWPT based HRV analysis. While correlations between results obtained by different wavelets were significant for any of the wavelet comparisons investigated in our study, low,
i.e. clinically insignificant, correlation coefficients and/or visual inspection of correlation plots might be a better indicator for unsatisfactory correlation. Especially LFHF results obtained by MODWPT based HRV analysis using different wavelets, which are commonly used to assess sympathovagal balance need to be interpreted with caution. Conflict of interest The authors have no conflicts of interest or financial disclosures to report. References [1] V. Pichot, F. Roche, S. Celle, J.-C. Barthélémy, F. Chouchou, HRVanalysis: a free software for analyzing cardiac autonomic activity, Front. Physiol. [Internet] 7 (November) (2016) 557 [cited 2017 Jul 6]. Available from: http://www.ncbi. nlm.nih.gov/pubmed/27920726. [2] V. Pichot, J.M. Gaspoz, S. Molliex, A. Antoniadis, T. Busso, F. Roche, et al., Wavelet transform to quantify heart rate variability and to assess its instantaneous changes, J. Appl. Physiol. [Internet] 86 (March (3)) (1999) 1081–1091 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/ pubmed/10066727. [3] S. Tejman-Yarden, O. Levi, A. Beizerov, Y. Parmet, T. Nguyen, M. Saunders, et al., Heart rate analysis by sparse representation for acute pain detection, Med. Biol. Eng. Comput. [Internet] 54 (April (4)) (2016) 595–606 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/26264057. [4] D.A. Kenwright, A. Bernjak, T. Draegni, S. Dzeroski, M. Entwistle, M. Horvat, et al., The discriminatory value of cardiorespiratory interactions in distinguishing awake from anaesthetised states: a randomised observational study, Anaesthesia [Internet] 70 (December (12)) (2015) 1356–1368 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/ 26350998. [5] C.G. Scully, G.C. Kramer, D.G. Strauss, Evaluation of heart rate and blood pressure variability as indicators of physiological compensation to hemorrhage before shock, Shock [Internet] 43 (May (5)) (2015) 463–469 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/ 25692248. [6] S. Akar, S. Kara, V. Bilgic¸, Investigation of heart rate variability in major depression patients using wavelet packet transform, Psychiatry Res. 30 (2016) 326–332. [7] F. Ebrahimi, Setarehdan S-KS, J. Ayala-Moyeda, H. Nazeran, Automatic sleep staging using empirical mode decomposition, discrete wavelet transform, time-domain, and nonlinear dynamics features of heart rate variability signals, Comput. Methods Programs Biomed. [Internet] 112 (October (1)) (2013) 47–57 [cited 2016 Nov 28]. Available from: http://www.ncbi.nlm.nih. gov/pubmed/23895941. [8] G. Li, W.-Y. Chung, Detection of driver drowsiness using wavelet analysis of heart rate variability and a support vector machine classifier, Sensors (Basel) [Internet] 13 (12) (2013) 16494–16511 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/24316564. [9] A.V. Oppenheim, R.W. Schafer, Discrete-time Signal Processing, Pearson Higher Education, 2010. [10] J.L. Ducla-Soares, M. Santos-Bento, S. Laranjo, A. Andrade, E. Ducla-Soares, J.P. Boto, et al., Wavelet analysis of autonomic outflow of normal subjects on head-up tilt, cold pressor test, Valsalva manoeuvre and deep breathing, Exp. Physiol. [Internet] 92 (July (4)) (2007) 677–686 [cited 2016 Jul 14]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/17468200. [11] F. Ebrahimi, S. Setarehdan, J. Ayala-Moyeda, H. Nazeran, Automatic sleep staging using empirical mode decomposition, discrete wavelet transform, time-domain, and nonlinear dynamics features of heart rate variability signals, Comput. Methods Programs Biomed. 112 (1) (2013) 47–57. [12] F. Chouchou, V. Pichot, C. Perchet, V. Legrain, L. Garcia-Larrea, F. Roche, et al., Autonomic pain responses during sleep: a study of heart rate variability, Eur. J. Pain 15 (6) (2011) 554–560. [13] A. Hossen, A. Barhoum, D. Jaju, V. Gowri, K. Al-Hashmi, M.O. Hassan, et al., Identification of patients with preeclampsia from normal subjects using wavelet-based spectral analysis of heart rate variability, Technol. Health Care [Internet] 25 (August (4)) (2017) 641–649 [cited 2017 Sep 12]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/28436399. [14] A. Deschamps, I. Kaufman, S.B. Backman, G. Plourde, Autonomic nervous system response to epidural analgesia in laboring patients by wavelet transform of heart rate and blood pressure variability, Anesthesiology [Internet] 101 (July (1)) (2004) 21–27 [cited 2017 Jul 6]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/15220767. [15] C.A. García, A. Otero, X. Vila, D.G. Márquez, A new algorithm for wavelet-based heart rate variability analysis, Biomed. Signal. Process. Control [Internet] 8 (November (6)) (2013) 542–550 [cited 2017 Mar 1]. Available from: http://linkinghub.elsevier.com/retrieve/pii/S1746809413000803. [16] Kubios HRV Software 2.2 [Internet], 2015, Available from: http://kubios.uef.fi. [17] M.P. Tarvainen, J.-P. Niskanen, J.A. Lipponen, P.O. Ranta-Aho, P.A. Karjalainen, G. Berntson, et al., Kubios HRV–heart rate variability analysis software,
A. Tzabazis et al. / Biomedical Signal Processing and Control 40 (2018) 220–225
[18]
[19]
[20]
[21]
Comput. Methods Programs Biomed. [Internet]. Elsevier 113 (1) (2014) 210–220 [cited 2016 Jul 20]. Available from: http://www.ncbi.nlm.nih.gov/ pubmed/24054542. M.P. Tarvainen, J. Lipponen, J.-P. Niskanen, P.O. Ranta-aho, Kubios HRV User’s Guide [Internet], 2017 [cited 2017 Jul 13]. Available from: http://www.kubios. com/downloads/Kubios HRV Users Guide.pdf. M.P. Tarvainen, P.O. Ranta-Aho, P.A. Karjalainen, An advanced detrending method with application to HRV analysis, IEEE Trans. Biomed. Eng. [Internet] 49 (February (2)) (2002) 172–175 [cited 2017 Jul 13]. Available from: http:// ieeexplore.ieee.org/document/979357/. R Core Team, R: A Language and Environment for Statistical Computing. [Internet], R Foundation for Statistical Computing, Vienna, Austria, 2013, Available from: http://www.r-project.org/. ˜ L. Rodríguez-Linares, A. Méndez, M. Lado, D. Olivieri, X. Vila, I. Gómez-Conde,
225
An open source tool for heart rate variability spectral analysis, Comput. Methods Programs Biomed. 103 (1) (2011) 39–50. [22] R.B. D’Agostino, Tests for normal distribution, in: R.B. D’Agostino, M.A. Stephens (Eds.), Goodness-Of-Fit Techniques, Marcel Dekker, 1986. [23] E. Jonckheere, P. Lohsoonthorn, S. Musuvathy, V. Mahajan, M. Stefanovic, On a standing wave central pattern generator and the coherence problem, Biomed. Signal. Process. Control [Internet] 5 (2010) 336–347 [cited 2017 Jun 30]. Available from: http://ac.els-cdn.com/S1746809410000273/1-s2.0S1746809410000273-main.pdf? tid=a0214b5a-5dcf-11e7-a71a00000aacb361&acdnat=1498853379 9e8e2d3aaabea471b645125a0a6f8303. [24] H.D. Esperer, C. Esperer, R.J. Cohen, Cardiac arrhythmias imprint specific signatures on Lorenz plots, Ann. Noninvasive Electrocardiol. [Internet] 13 (January (1)) (2008) 44–60 [cited 2017 Jun 30]. Available from: http://doi. wiley.com/10.1111/j.1542-474X.2007.00200.x.