IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 7, JULY 2014
1931
Streaming Updates for Heart Rate Variability Algorithms Stergios Stergiou∗ and Rajalakshmi Balakrishnan
Abstract—Heart rate variability (HRV) quantifies the fluctuations of the lengths of consecutive heart beat intervals, and is a reliable descriptor of many physiological factors modulating the normal rhythm of the heart. As the heart rate signal is nonstationary, indicators deduced from it may be present at all times, but may also occur episodically at nonpredetermined time instances. The potential for real-time feedback long-term ambulatory recordings is thus apparent. Numerous methods for measuring HRV have been standardized and are in active use, but are typically not designed to operate at real time. In this paper, we study the most popular HRV quantification methods and propose streaming algorithms that maximally utilize previously computed information without altering the output of the methods. We demonstrate speedups of more than two orders of magnitude for typical use-case scenarios. Using our algorithms on embedded systems that compute HRV leads to dramatic decreases in power consumption and in some cases allows for computation of metrics that were not previously possible at real time. Index Terms—Heart rate variability, streaming algorithms.
I. INTRODUCTION EART rate variability (HRV) quantifies the fluctuations of the lengths of consecutive heart beat intervals. HRV reflects the regularity of heart beats, with higher regularity implying lower HRV, and is a reliable descriptor of many physiological factors modulating the normal rhythm of the heart. Crucially, HRV provides a powerful means of observing the interplay between the sympathetic and parasympathetic nervous systems and is thought to reflect the heart’s ability to adapt to changing circumstances by detecting and quickly responding to unpredictable stimuli. In [1], Dagres et al. investigated the influence of Holter duration on the detection of recurrences after ablation for atrial fibrillation. During the complete seven-day recording, 30% of patients had a recurrence. A 24-h Holter would have detected 59%, a 48-h Holter 67%, and a 72-h Holter 80% of patients with recurrences, whereas a four-day recording would have detected 91% of the recurrences that were detected with the complete seven-day recording. They concluded that a Holter duration of less than four days misses a great portion of recurrences. Recent technological advances, including miniaturization, more power efficient CPUs, vast deployment of mobile
H
Manuscript received November 15, 2013; revised January 9, 2014; accepted February 11, 2014. Date of publication February 20, 2014; date of current version June 14, 2014. Asterisk indicates corresponding author. ∗ S. Stergiou is with the Yahoo! Inc., Sunnyvale, CA 94089 USA (e-mail: [email protected]). R. Balakrishnan is with Fujitsu Laboratories of America, Sunnyvale CA 94085 USA (e-mail: [email protected]). Digital Object Identifier 10.1109/TBME.2014.2307014
platforms, and new low-power RF protocols have resulted in the development of wearable ECG monitoring devices that are now capable of continuously recording and transmitting ECG signals for seven-day, 14-day, or even longer periods [2]–[6]. HRV methods and their clinical applications have been extensively studied [7], [8]. In 1996, the European Society of Cardiology and the North American Society of Pacing and Electrophysiology formed a task force with the goal of standardizing on measurements, physiological interpretation, and clinical uses of HRV methods [9]. Their methods accept as input a sequence of normal-to-normal (NN) intervals between adjacent QRS complexes resulting from sinus node depolarizations. The standardized methods were classified into time domain and frequency domain. In the time domain, they include 1) simple metrics: mean NN interval, mean heart rate, difference between the longest and shortest NN interval, difference between night and day heart rate, 2) statistical metrics (applied on longer, traditionally 24 h, periods): standard deviation of the NN intervals (SDNN), and its variations, the square root of the mean squared differences of successive NN intervals (RMSSD), the proportion of interval differences of successive NN intervals greater than 50 ms (pNN50), and 3) geometrical interpretation metrics. In the frequency domain, the standardized methods quantify the power within specific frequency bands via power spectral density (PSD) analysis. The primary metric is LF/HF, which is the ratio of the power within bands LF ([0.04 Hz, 0.15 Hz]) and HF ([0.15 Hz, 0.4 Hz]). Although the performance of these methods has been extensively analyzed in terms of quantifying HRV, their computational aspects have not been studied. As a motivational example, let us consider pNN50, one of the simpler HRV algorithms. Given a sequence of k + 1 NN intervals (x0 , . . . , xk ), the algorithm computes the percentage of successive NN interval differences whose absolute value is larger than 50 ms. Specifically: pNN50 = 1/k kj=1 Xj , where Xj is 1, if |xj − xj −1 | > 50 and 0 otherwise. Upon the detection of a new N N interval xk +1 , we would like to compute pNN50 , the updated value of pNN50 for sequence (x1 , . . . , xk +1 ). This can be performed either by recomputing the sum on the new sequence, or by reusing previously computed information as follows: pNN50’ = pNN50 + (Xk +1 − X1 )/k, thus avoiding the superfluous recomputation of the sum. In this paper, we propose and analyze a collection of algorithms that are designed to efficiently update at real time the most popular time domain and nonparametric frequency domain HRV metrics. Specifically, we will present algorithms that optimally update the SDNN, RMSSD, Median NN, NN Interval Range, Triangular Index, pNN50 and LF/HF HRV metrics. In Section II,
0018-9294 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications standards/publications/rights/index.html for more information.
1932
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 7, JULY 2014
buffer, which is a fixed-size buffer data structure that supports insertions and deletions at its two endpoints in O(1) time, along with O(1) random access to its elements. A. Computing Central Moments Over Sliding Windows The kth raw moment for a window of elements (x1 , . . . , xN ) is defined as μk =
N 1 k x N n =1 n
where mean μ = μ1 . The kth central moment is defined as we present our algorithms for the time domain HRV methods. In Section III, we present our algorithms for the frequency-domain HRV methods. In Section IV, we demonstrate the performance gains of adopting the proposed algorithms. Section V concludes this paper.
μk =
N 1 (xn − μ)k . N n =1
We will show how to update the kth central moment in O(k) time. We first observe that raw moments of any degree can be updated in constant time. For j = 1. .k, we maintain the unscaled sums
II. TIME-DOMAIN METHODS In this section, we present algorithms that enable streaming updates for the most popular time-domain HRV methods. We first present an algorithm for updating raw and central moments over sliding windows. This allows for constant-time updates for 1) the mean over a sliding window of NN intervals, 2) SDNN, which is the root of the second central moment of a window of NN intervals, 3) SDANN, which is an HRV algorithm that is computationally similar to SDNN, defined as the mean of the standard deviations of consecutive 300-s NN intervals, and 4) RMSSD, which is the root of the second raw moment of a window of NN interval differences. Second, we present an algorithm for updating the median over a sliding window in expected constant time. This allows for the efficient computation of the median over a sliding window of δNN intervals. We then analyze an amortized constant-time algorithm for updating the minimum and maximum over a sliding window. This allows for the calculation of the difference between the longest and shortest NN interval over a sliding window of NN intervals. Afterward, we present an algorithm for updating the HRV Triangular (HRVtr) index over a sliding window of NN intervals in O(log n) time. Finally, we present an algorithm for updating pNN50 over a sliding window of δNN intervals. We maintain the data points of the sliding window on a ring (or circular)
sumj =
N −1
xjn .
n =0
A ring buffer is utilized for maintaining the N most recent data points. Upon receiving a new element xN , we extract the oldest element x0 from the buffer and add xN into it. We add xjN − xj0 to sumj for j = 1. .k. The updated μj is simply sumj /N . We note that by keeping the original elements in the ring buffer, we can reuse it for updating the raw moments of any degree. Via the binomial transform, the kth central moment can be expressed as k k k −n (−1) μn μk −n . μk = n n =0
As each raw moment can be updated in O(1), all raw moments up to degree k can be updated in O(k). Therefore, μk can also be updated in O(k) time. Only a single common ring buffer containing the elements of the relevant window is required to be maintained. The algorithm is depicted in Algorithm 1. B. Computing Medians Over Sliding Windows We present an algorithm that updates the median over a rolling window in O(1) expected time. Our comparative evaluations
STERGIOU AND BALAKRISHNAN: STREAMING UPDATES FOR HEART RATE VARIABILITY ALGORITHMS
have shown that the median of a sliding window of δNN intervals is one of the most valuable HRV algorithms and yields very intuitive connections with other HRV algorithms. For this reason, we include it in this paper. Nevertheless, an evaluation of its performance as an HRV algorithm is beyond the scope of this paper. 1) Initialization Phase: The initialization phase of the algorithm is depicted in Algorithm 2. We insert the first point into the ring buffer that contains the rolling window and clear the table that holds the histogram (where M axV alue is typically 100) and also mark the existence of the first data point into it. The first point is initially also the median. We additionally initialize two auxiliary variables whose purpose will become evident shortly. 2) Update Phase: The update algorithm maintains three invariants. The first is that balance ∈ [0, 1]. The second is that 0 ≤ splitter < histogram[median]. The third is that variable
1933
median is indeed the median of all the elements in the ring buffer. Variable balance should be interpreted as follows: If positive, it denotes the number of sample points smaller than the current median that would have to be inserted into the histogram such that the current median remains valid. Similarly, if negative, it denotes the number of sample points larger than the current median that would have to be inserted such that the current median remains valid. Variable splitter is a bookkeeping variable. It denotes the exact position within all histogram[median] elements where the median should lie. More concretely, splitter elements should lie to the left, while histogram[median] − splitter − 1 elements should lie to the right. The actual median value is positioned at exact position splitter within the elements of the medianth entry of the histogram. Through the balance invariant, the splitter invariant is maintained and through that, the median invariant is maintained. We can now proceed by describing the algorithm which is depicted in Algorithm 3. Lines 1–9 handle the steady-state case where the ring buffer is full. Before entering the new point, we pop the earliest point that was entered. Depending on its location in relation to the current median, the median value may have been invalidated. We update balance to denote the number and relative position of new elements that would have to be entered into the histogram in order to keep the median valid. Lines 10–16 add the new point to the ring buffer and the histogram, also properly updating balance. At this point, we are in the position to know if median needs to be updated. At the entry point, balance can either be 0 or 1. After segments 1–16, it is either increased or decreased by at most 2 so that its new possible range is [−2, 3]. We will modify balance to bring it back to the range dictated by the invariant. If balance ∈ [2, 3], the only way to reduce it is to increase the median. Increasing the median position by 1 allows us to reduce balance by 2. This case is handled by segments 17–24. The median is moved to the right by increasing splitter. If splitter becomes large enough such that it lies outside the histogram[median] limits, then it is reset and median is updated to the next nonempty histogram entry to its right. In a similar fashion, if balance ∈ [−2, −1], the only way to increase it is by decreasing the median. This symmetric case is handled by segments 25–32.
1934
The algorithm needs to handle one final bookkeeping case, where the element that was originally popped from the ring buffer was exactly at median, thus possibly invalidating the splitter invariant. This case is handled by segments 33–38. C. Computing Minimums/Maximums Over Sliding Windows For every new element xt , we would like to update the max(xt−w +1 , . . . , xt ) over a sliding window of w elements. We maintain a ring buffer C of length w. For each incoming element xt , we remove all elements that are smaller than or equal to xt from the end of C and then insert xt into C. Elements are time stamped when entering C and are removed after w updates. The max of the sliding window is the first element of C. Each element is pushed into and out of C exactly once; therefore, updates are performed in amortized O(1) time. The algorithm is depicted in Algorithm 4. The algorithm for computing the minimum is a slight modification of Algorithm 4 (where the inequality on Line 1 is reversed.) We note that the two variations require separate ring buffers. Our empirical observations have shown that the number of elements in these buffers does not grow significantly in size compared to w. In fact, it can be seen that they maintain a monotonically decreasing (respectively increasing) subsequence of the complete rolling window whose length is bounded by the difference between the maximum and the minimum element in the window (a practically small difference compared to w.) D. Computing Triangular Indices Over Sliding Windows For computing the HRVtr, we maintain the density distribution of all NN intervals in the sliding window. The density bin width is typically 7.8125 ms per Task Force’s recommendations.1 HRVtr is defined as the ratio of the total NN intervals over the frequency of the modal bin. The total number of NN intervals within a sliding window is readily available. We will now show how to update the most dense bin of the NN distribution. The algorithm is depicted in Algorithm 5. We maintain a ring buffer C which contains the bin numbers for each element in the sliding window. We maintain tuples {f req, bin} in a maxheap H, ordered on the bin frequencies. There exists a tuple in H for each bin, whose frequency is initially 0. We addition1 7.8125 ms corresponds to a sampling frequency of 128 Hz which was standard ECG sampling rate in 1996.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 7, JULY 2014
ally maintain a direct map between the bin numbers and the corresponding tuple positions within H. Upon the arrival of a new data point pt, we remove the oldest element from C, obtain its location in H via direct map, reduce the frequency of the corresponding tuple in H by 1, and adjust its location in H via a heapify-down operation. We then push the bin number of pt into C and detect its position via direct map in H. We proceed by increasing the corresponding tuple’s frequency by 1 and readjust its position in H via a heapify-up operation. Traditional heap operations heapify-down and heapify-up are extended such that whenever two tuples in H swap positions, their corresponding locations in direct map are also appropriately updated. The algorithm’s complexity is dominated by the complexity of the heapify-up and heapify-down operations, which is O(log BinCount), where BinCount is the total number of bins in the distribution and is typically 256. The space requirements are O(w + BinCount), where w is the length of the sliding window. E. Computing pNN50 Over Sliding Windows Metric pNN50 can be updated in O(1) time in a straightforward manner. As noted previously, we maintain all δNN values in a ring buffer of length w. We maintain a counter cnt50 which counts the number of |δNN| values that are larger than 50. Upon the arrival of a new element xN , the oldest element x0 is extracted from the buffer. If |x0 | > 50, then cnt50 is decreased by 1. If |xN | > 50, then cnt50 is increased by 1. Metric pNN50 is then obtained from cnt50 /w. We note that, as the crucial information per new data point is binary, we maintain a ring buffer of Booleans instead of the actual data points, thereby decreasing memory requirements for the particular algorithm. III. FREQUENCY-DOMAIN METHODS PSD analysis is useful for understanding how power is distributed as a function of frequency. PSD can either be estimated via regression analysis or, more commonly, via fast Fourier transform (FFT) analysis, while both approaches provide comparable results. The HRV measures employed are 1) the total power P , 2) the power within the VLF/LF/HF frequency bands ([0, 0.04 Hz]/[0.04 Hz, 0.15 Hz]/[0.15 Hz, 0.4 Hz]), 3) the normalized power within the LF/HF bands [LFnorm = LF/(P − VLF), HFnorm = HF/(P − VLF)], and 4) the ratio LF/HF. Ratio LF/HF seems to be the more widely used metric. Before applying FFT, the unevenly spaced peak detection signal
STERGIOU AND BALAKRISHNAN: STREAMING UPDATES FOR HEART RATE VARIABILITY ALGORITHMS
1935
be the Fourier coefficients for window (xt−L +1 , . . . , xt ). Then, it holds Xkt =
L −1
xn +t−L · e−2π ik (n −1)/L + xt · e−2π i(L −1)/L
n =1
=
L −1
xn +t−L · e−2π ik n /L
· e2π ik /L
n =1
+ xt · e−2π i(L −1)/L = (Xkt−1 − xt−L + xt ) · e2π ik /L .
typically needs to be resampled. Many studies in the field of HRV report resampling rates between 2 and 4 Hz. In [10], a study is performed on the effects of beat replacement and resampling, and it is concluded that a resampling rate of 7 Hz is appropriate as it allows for the resolution of heart rates of up to 210 bpm. We will now present a method for updating Fourier coefficients in constant time. As will be shown, each coefficient can be updated independent of the rest. This allows for the updating of just the appropriate subset of coefficients corresponding to the desired frequency bands, as opposed to updating all of them at each iteration. We note that the problem of real-time FFT computation has been extensively analyzed in [11] and [12]. We present an algorithm for updating the power in the LF and HF bands, and therefore, the ratio LF/HF in Algorithm 7. We also show how the total power can be updated in a straightforward manner via Parseval’s theorem in Algorithm 8. We note that the proposed algorithms are applicable only when FFT is directly used on the input signal, specifically when windowing is not adopted. A. Updating FFT Coefficients Over Sliding Windows Let {x}t be a streaming sequence of numbers and L be the length of a sliding window. We would like to update the Fourier coefficients of window (xt−L +1 , . . . , xt ) for every new element xt arriving at time t. Let, Xkt =
L −1 n =0
xn +t−L +1 · e−2π ik n /L
We observe that each Fourier coefficient Xkt can be updated from Xkt−1 in constant time. More concretely, if all powers of the Lth primitive root of unity e2π i/L are precomputed then one real addition and one complex multiplication will suffice for updating Xkt . Alternatively, one additional complex multiplication per coefficient will be sufficient for generating the appropriate root power. B. Computing Power Spectral Densities Over Sliding Windows The periodogram of {x}t is defined as ⎧ 1 ⎪ ⎨ 2 |Xk |2 , k = 0, N/2 N P (k) = ⎪ ⎩ 2 |X |2 , k = 1, . . . , N/2 − 1. k N2 P (k) is an estimate of the power spectrum of the time series at frequency fk = ΔkN , where Δ is the sampling period. Since xt is real, then Xk = XN −k . By Parseval’s theorem
N /2
k =0
P (k) =
N −1 1 |xj |2 . N j =0
IV. EVALUATION A. Data Analytics Platform and Experimental Setup All algorithms were implemented in C++ and were executed on the Sprout platform. The core of the platform is the Sprout, a special-purpose built mobile data aggregation and analytics device, depicted in Fig. 1(a). It is based on a 454-MHz ARM926EJ-S processor and communicates with sensors and the user via 802.11n, Bluetooth 3.0, and USB 2.0. The Linux-based platform facilitates data collection and analysis,
1936
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 7, JULY 2014
Fig. 1. Data aggregation and analytics platform. (a) Sprout mobile device for real-time biosignals aggregation and analytics. (b) Zephyr BioHarness 3.0 ECG strap. Fig. 3. Speedup results for the HRVtr index algorithm for (a) 7.7125 ms and (b) 1 ms bin widths.
Fig. 2. Speedup results for SDNN and RMSSD algorithms based on the central and raw moments update algorithm.
and provides a real-time streaming user interface accessible via an on-board web server, while additionally transmitting information to a central server. ECG data were collected with a Zephyr Bioharness 3.0 ECG sensor, shown in Fig. 1(b). Sampling and peak detection was performed at 1 kHz. For the purposes of evaluating the proposed algorithms, ECG data were continuously collected over a 48-h period. The time-domain algorithms are individually benchmarked across 12 different sliding window lengths, specifically of 300, 600, 1200, 2400, 3600, 7200, 14400, 28800, 43200, 57600, 72000, and 86400 s. Twenty-four hours of generated data were collected from each algorithm. We present results on SDNN, RMSSD, HRVtr index with binning of 7.8125 ms (128 Hz) and 1 ms (1 kHz), pNN50, NN Median, and NN interval range. As a representative frequency domain algorithm, we present results on updating the LF/HF ratio. FFT window lengths used are 29 , . . . , 216 . The signal is linearly interpolated and resampled at rates from 1 to 10 Hz before applying FFT. Each algorithm is compared against its “classic” implementation. Specifically, for SDNN, RMSSD, NN Median, NN Interval Range, and pNN50, the metrics are computed as efficiently as possible from the data points on the sliding window without reusing previously computed information. For the HRVtr index, the classic algorithm also maintains the histogram of the NN intervals. B. Results In Fig. 2, we present the performance ratios between the proposed and classical implementations for SDNN and RMSSD. We observe speedups that are linearly increasing with the
Fig. 4.
Speedup results for the NN median update algorithm.
Fig. 5.
Speedup results for the NN interval range update algorithm.
window size, as theoretically expected. In Fig. 3, we present speedups for the HRVtr index algorithm. We observe that the algorithm is more sensitive to the input signal. We note that both the classic and proposed algorithms operate on the NN histogram, whose size depends only on the granularity of the binning and not on the size of the window. The speedups obtained are 5.8× to 13.9× for 7.8125 ms bins, and 36.5× to 79.9× for 1 ms bins. In Fig. 4, we present the speedups for the NN median algorithm. We note that the speedups obtained are significant, and comparable to the ones obtained from the moments algorithms, even though the proposed algorithm performs updates only in expected constant time. In Fig. 5, we present speedups for the NN interval range algorithm, which is based on the maximum/minimum update
STERGIOU AND BALAKRISHNAN: STREAMING UPDATES FOR HEART RATE VARIABILITY ALGORITHMS
1937
In Fig. 8, we present speedups for the LF/HF algorithm for varying window sizes and resampling rates. We note that for the classic algorithm, the dominating part is computing the FFT. For a window of length 216 , the classic algorithm requires 110 020” to complete the 24 h dataset, and thus cannot be run at real time. As the window size increases, the proposed algorithm is relatively faster. Also, we observe that for a fixed window size, when resampling rate increases, the width of the LF and HF bands decrease. As the proposed algorithm only updates the FT coefficients within those bands, it becomes faster as resampling rate increases. Fig. 6.
Speedup results for the pNN50 update algorithm.
V. CONCLUSION In this paper, we studied the most popular HRV quantification methods and proposed streaming algorithms that maximally utilize previously computed information without altering the output of the methods. We demonstrated speedups of more than two orders of magnitude for typical use-case scenarios that translate to significantly longer battery lives for mobile systems that adopt them. For the most computationally heavy frequencydomain methods, our algorithms allow for real-time updates in cases where it was not previously possible. REFERENCES
Fig. 7.
Aggregate speedup results for all time-domain algorithms.
Fig. 8. Speedup results for Fourier transform LF/HF update algorithm. The original signal was linearly interpolated and resampled at rates ranging from 1 to 10 Hz.
algorithm. In Fig. 6, we present speedups for the pNN50 algorithm. This algorithm is the most straightforward to optimize and indeed yields the largest speedups, due in large to the compressed way the sliding window is maintained. In Fig. 7, we present the overall speedup obtained when executing all the time-domain algorithms. We observe that, even though the speedups of the HRVtr algorithm are lower relative to the rest of the time-domain algorithms, the algorithm is fast in absolute execution times as it only operates on an input whose size does not depend on the sliding window length.
[1] N. Dagres, H. Kottkamp, C. Piorkowski, S. Weis, A. Arya, P. Sommer, K. Bode, J.-H. Gerds-Li, D. T. Kremastinos, and G. Hindricks, “Influence of the duration of Holter monitoring on the detection of arrhythmia recurrences after catheter ablation of atrial fibrillation: Implications for patient follow-up,” Int. J. Cardiol., vol. 139, no. 3, pp. 305–306, 2010. [2] S. S. Lobodzinski and M. M. Laks, “New devices for very long-term ECG monitoring,” Cardiol. J., vol. 19, no. 2, pp. 210–214, 2012. [3] B. Cinaz, R. La Marca, B. Arnrich, and G. Tr¨oster, “Towards continuous monitoring of mental workload,” presented at the 5th Int. Workshop Ubiquit. Health Wellness, Copenhagen, Denmark, 2010. [4] “Vital Connect.” (2013). [Online]. Available:http://www.vitalconnect. com/ [5] “CardioLeaf.” (2013). [Online]. Available: http://www. clearbridgevitalsigns.com/ [6] B. Gyselinckx, C. Van Hoof, J. Ryckaert, R. F. Yazicioglu, P. Fiorini, and V. Leonov, “Human++: Autonomous wireless sensors for body area networks,” in Proc. IEEE Custom Integr. Circuits Conf., 2005, pp. 13–19. [7] U. R. Acharya, K. P. Joseph, N. Kannathal, C. M. Lim, and J. S. Suri, “Heart rate variability: A review,” Med. Biol. Eng. Comput., vol. 44, no. 12, pp. 1031–1051, 2006. [8] M. P. Tarvainen, J.-P. Niskanen, J. Lipponen, P. Ranta-Aho, and P. Karjalainen, “Kubios HRV—A software for advanced heart rate variability analysis,” in Proc. 4th Eur. Conf. Int. Federat. Med. Biol. Eng., 2009, pp. 1022–1025. [9] A. Camm, M. Malik, J. T. Bigger, Jr, G. Breithardt, S. Cerutti, R. J. Cohen, P. Coumel, E. L. Fallen, H. L. Kennedy, R. E. Kleiger, F. Lombardi, A. Malliani, A. J. Moss, J. N. Rottman, G. Schmidt, P. J. Schwartz, and D. H. Singer, “Heart rate variability: Standards of measurement, physiological interpretation and clinical use. Task force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology,” Circulation, vol. 93, no. 5, pp. 1043–1065, 1996. [10] G. Clifford and L. Tarassenko, “Quantifying errors in spectral estimates of HRV due to beat replacement and resampling,” IEEE Trans. Biomed. Eng., vol. 52, no. 4, pp. 630–638, Apr. 2005. [11] B. Farhang-Boroujeny and Y. Lim, “A comment on the computational complexity of sliding FFT,” IEEE Trans. Circuits Syst. II: Analog Digital Signal Process., vol. 39, no. 12, pp. 875–876, Dec. 1992. [12] P.-C. Lo and Y.-Y. Lee, “Real-time implementation of the moving FFT algorithm,” Signal Process., vol. 79, no. 3, pp. 251–259, 1999.
Authors’ photographs and biographies not available at the time of publication.
1931
Streaming Updates for Heart Rate Variability Algorithms Stergios Stergiou∗ and Rajalakshmi Balakrishnan
Abstract—Heart rate variability (HRV) quantifies the fluctuations of the lengths of consecutive heart beat intervals, and is a reliable descriptor of many physiological factors modulating the normal rhythm of the heart. As the heart rate signal is nonstationary, indicators deduced from it may be present at all times, but may also occur episodically at nonpredetermined time instances. The potential for real-time feedback long-term ambulatory recordings is thus apparent. Numerous methods for measuring HRV have been standardized and are in active use, but are typically not designed to operate at real time. In this paper, we study the most popular HRV quantification methods and propose streaming algorithms that maximally utilize previously computed information without altering the output of the methods. We demonstrate speedups of more than two orders of magnitude for typical use-case scenarios. Using our algorithms on embedded systems that compute HRV leads to dramatic decreases in power consumption and in some cases allows for computation of metrics that were not previously possible at real time. Index Terms—Heart rate variability, streaming algorithms.
I. INTRODUCTION EART rate variability (HRV) quantifies the fluctuations of the lengths of consecutive heart beat intervals. HRV reflects the regularity of heart beats, with higher regularity implying lower HRV, and is a reliable descriptor of many physiological factors modulating the normal rhythm of the heart. Crucially, HRV provides a powerful means of observing the interplay between the sympathetic and parasympathetic nervous systems and is thought to reflect the heart’s ability to adapt to changing circumstances by detecting and quickly responding to unpredictable stimuli. In [1], Dagres et al. investigated the influence of Holter duration on the detection of recurrences after ablation for atrial fibrillation. During the complete seven-day recording, 30% of patients had a recurrence. A 24-h Holter would have detected 59%, a 48-h Holter 67%, and a 72-h Holter 80% of patients with recurrences, whereas a four-day recording would have detected 91% of the recurrences that were detected with the complete seven-day recording. They concluded that a Holter duration of less than four days misses a great portion of recurrences. Recent technological advances, including miniaturization, more power efficient CPUs, vast deployment of mobile
H
Manuscript received November 15, 2013; revised January 9, 2014; accepted February 11, 2014. Date of publication February 20, 2014; date of current version June 14, 2014. Asterisk indicates corresponding author. ∗ S. Stergiou is with the Yahoo! Inc., Sunnyvale, CA 94089 USA (e-mail: [email protected]). R. Balakrishnan is with Fujitsu Laboratories of America, Sunnyvale CA 94085 USA (e-mail: [email protected]). Digital Object Identifier 10.1109/TBME.2014.2307014
platforms, and new low-power RF protocols have resulted in the development of wearable ECG monitoring devices that are now capable of continuously recording and transmitting ECG signals for seven-day, 14-day, or even longer periods [2]–[6]. HRV methods and their clinical applications have been extensively studied [7], [8]. In 1996, the European Society of Cardiology and the North American Society of Pacing and Electrophysiology formed a task force with the goal of standardizing on measurements, physiological interpretation, and clinical uses of HRV methods [9]. Their methods accept as input a sequence of normal-to-normal (NN) intervals between adjacent QRS complexes resulting from sinus node depolarizations. The standardized methods were classified into time domain and frequency domain. In the time domain, they include 1) simple metrics: mean NN interval, mean heart rate, difference between the longest and shortest NN interval, difference between night and day heart rate, 2) statistical metrics (applied on longer, traditionally 24 h, periods): standard deviation of the NN intervals (SDNN), and its variations, the square root of the mean squared differences of successive NN intervals (RMSSD), the proportion of interval differences of successive NN intervals greater than 50 ms (pNN50), and 3) geometrical interpretation metrics. In the frequency domain, the standardized methods quantify the power within specific frequency bands via power spectral density (PSD) analysis. The primary metric is LF/HF, which is the ratio of the power within bands LF ([0.04 Hz, 0.15 Hz]) and HF ([0.15 Hz, 0.4 Hz]). Although the performance of these methods has been extensively analyzed in terms of quantifying HRV, their computational aspects have not been studied. As a motivational example, let us consider pNN50, one of the simpler HRV algorithms. Given a sequence of k + 1 NN intervals (x0 , . . . , xk ), the algorithm computes the percentage of successive NN interval differences whose absolute value is larger than 50 ms. Specifically: pNN50 = 1/k kj=1 Xj , where Xj is 1, if |xj − xj −1 | > 50 and 0 otherwise. Upon the detection of a new N N interval xk +1 , we would like to compute pNN50 , the updated value of pNN50 for sequence (x1 , . . . , xk +1 ). This can be performed either by recomputing the sum on the new sequence, or by reusing previously computed information as follows: pNN50’ = pNN50 + (Xk +1 − X1 )/k, thus avoiding the superfluous recomputation of the sum. In this paper, we propose and analyze a collection of algorithms that are designed to efficiently update at real time the most popular time domain and nonparametric frequency domain HRV metrics. Specifically, we will present algorithms that optimally update the SDNN, RMSSD, Median NN, NN Interval Range, Triangular Index, pNN50 and LF/HF HRV metrics. In Section II,
0018-9294 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications standards/publications/rights/index.html for more information.
1932
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 7, JULY 2014
buffer, which is a fixed-size buffer data structure that supports insertions and deletions at its two endpoints in O(1) time, along with O(1) random access to its elements. A. Computing Central Moments Over Sliding Windows The kth raw moment for a window of elements (x1 , . . . , xN ) is defined as μk =
N 1 k x N n =1 n
where mean μ = μ1 . The kth central moment is defined as we present our algorithms for the time domain HRV methods. In Section III, we present our algorithms for the frequency-domain HRV methods. In Section IV, we demonstrate the performance gains of adopting the proposed algorithms. Section V concludes this paper.
μk =
N 1 (xn − μ)k . N n =1
We will show how to update the kth central moment in O(k) time. We first observe that raw moments of any degree can be updated in constant time. For j = 1. .k, we maintain the unscaled sums
II. TIME-DOMAIN METHODS In this section, we present algorithms that enable streaming updates for the most popular time-domain HRV methods. We first present an algorithm for updating raw and central moments over sliding windows. This allows for constant-time updates for 1) the mean over a sliding window of NN intervals, 2) SDNN, which is the root of the second central moment of a window of NN intervals, 3) SDANN, which is an HRV algorithm that is computationally similar to SDNN, defined as the mean of the standard deviations of consecutive 300-s NN intervals, and 4) RMSSD, which is the root of the second raw moment of a window of NN interval differences. Second, we present an algorithm for updating the median over a sliding window in expected constant time. This allows for the efficient computation of the median over a sliding window of δNN intervals. We then analyze an amortized constant-time algorithm for updating the minimum and maximum over a sliding window. This allows for the calculation of the difference between the longest and shortest NN interval over a sliding window of NN intervals. Afterward, we present an algorithm for updating the HRV Triangular (HRVtr) index over a sliding window of NN intervals in O(log n) time. Finally, we present an algorithm for updating pNN50 over a sliding window of δNN intervals. We maintain the data points of the sliding window on a ring (or circular)
sumj =
N −1
xjn .
n =0
A ring buffer is utilized for maintaining the N most recent data points. Upon receiving a new element xN , we extract the oldest element x0 from the buffer and add xN into it. We add xjN − xj0 to sumj for j = 1. .k. The updated μj is simply sumj /N . We note that by keeping the original elements in the ring buffer, we can reuse it for updating the raw moments of any degree. Via the binomial transform, the kth central moment can be expressed as k k k −n (−1) μn μk −n . μk = n n =0
As each raw moment can be updated in O(1), all raw moments up to degree k can be updated in O(k). Therefore, μk can also be updated in O(k) time. Only a single common ring buffer containing the elements of the relevant window is required to be maintained. The algorithm is depicted in Algorithm 1. B. Computing Medians Over Sliding Windows We present an algorithm that updates the median over a rolling window in O(1) expected time. Our comparative evaluations
STERGIOU AND BALAKRISHNAN: STREAMING UPDATES FOR HEART RATE VARIABILITY ALGORITHMS
have shown that the median of a sliding window of δNN intervals is one of the most valuable HRV algorithms and yields very intuitive connections with other HRV algorithms. For this reason, we include it in this paper. Nevertheless, an evaluation of its performance as an HRV algorithm is beyond the scope of this paper. 1) Initialization Phase: The initialization phase of the algorithm is depicted in Algorithm 2. We insert the first point into the ring buffer that contains the rolling window and clear the table that holds the histogram (where M axV alue is typically 100) and also mark the existence of the first data point into it. The first point is initially also the median. We additionally initialize two auxiliary variables whose purpose will become evident shortly. 2) Update Phase: The update algorithm maintains three invariants. The first is that balance ∈ [0, 1]. The second is that 0 ≤ splitter < histogram[median]. The third is that variable
1933
median is indeed the median of all the elements in the ring buffer. Variable balance should be interpreted as follows: If positive, it denotes the number of sample points smaller than the current median that would have to be inserted into the histogram such that the current median remains valid. Similarly, if negative, it denotes the number of sample points larger than the current median that would have to be inserted such that the current median remains valid. Variable splitter is a bookkeeping variable. It denotes the exact position within all histogram[median] elements where the median should lie. More concretely, splitter elements should lie to the left, while histogram[median] − splitter − 1 elements should lie to the right. The actual median value is positioned at exact position splitter within the elements of the medianth entry of the histogram. Through the balance invariant, the splitter invariant is maintained and through that, the median invariant is maintained. We can now proceed by describing the algorithm which is depicted in Algorithm 3. Lines 1–9 handle the steady-state case where the ring buffer is full. Before entering the new point, we pop the earliest point that was entered. Depending on its location in relation to the current median, the median value may have been invalidated. We update balance to denote the number and relative position of new elements that would have to be entered into the histogram in order to keep the median valid. Lines 10–16 add the new point to the ring buffer and the histogram, also properly updating balance. At this point, we are in the position to know if median needs to be updated. At the entry point, balance can either be 0 or 1. After segments 1–16, it is either increased or decreased by at most 2 so that its new possible range is [−2, 3]. We will modify balance to bring it back to the range dictated by the invariant. If balance ∈ [2, 3], the only way to reduce it is to increase the median. Increasing the median position by 1 allows us to reduce balance by 2. This case is handled by segments 17–24. The median is moved to the right by increasing splitter. If splitter becomes large enough such that it lies outside the histogram[median] limits, then it is reset and median is updated to the next nonempty histogram entry to its right. In a similar fashion, if balance ∈ [−2, −1], the only way to increase it is by decreasing the median. This symmetric case is handled by segments 25–32.
1934
The algorithm needs to handle one final bookkeeping case, where the element that was originally popped from the ring buffer was exactly at median, thus possibly invalidating the splitter invariant. This case is handled by segments 33–38. C. Computing Minimums/Maximums Over Sliding Windows For every new element xt , we would like to update the max(xt−w +1 , . . . , xt ) over a sliding window of w elements. We maintain a ring buffer C of length w. For each incoming element xt , we remove all elements that are smaller than or equal to xt from the end of C and then insert xt into C. Elements are time stamped when entering C and are removed after w updates. The max of the sliding window is the first element of C. Each element is pushed into and out of C exactly once; therefore, updates are performed in amortized O(1) time. The algorithm is depicted in Algorithm 4. The algorithm for computing the minimum is a slight modification of Algorithm 4 (where the inequality on Line 1 is reversed.) We note that the two variations require separate ring buffers. Our empirical observations have shown that the number of elements in these buffers does not grow significantly in size compared to w. In fact, it can be seen that they maintain a monotonically decreasing (respectively increasing) subsequence of the complete rolling window whose length is bounded by the difference between the maximum and the minimum element in the window (a practically small difference compared to w.) D. Computing Triangular Indices Over Sliding Windows For computing the HRVtr, we maintain the density distribution of all NN intervals in the sliding window. The density bin width is typically 7.8125 ms per Task Force’s recommendations.1 HRVtr is defined as the ratio of the total NN intervals over the frequency of the modal bin. The total number of NN intervals within a sliding window is readily available. We will now show how to update the most dense bin of the NN distribution. The algorithm is depicted in Algorithm 5. We maintain a ring buffer C which contains the bin numbers for each element in the sliding window. We maintain tuples {f req, bin} in a maxheap H, ordered on the bin frequencies. There exists a tuple in H for each bin, whose frequency is initially 0. We addition1 7.8125 ms corresponds to a sampling frequency of 128 Hz which was standard ECG sampling rate in 1996.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 7, JULY 2014
ally maintain a direct map between the bin numbers and the corresponding tuple positions within H. Upon the arrival of a new data point pt, we remove the oldest element from C, obtain its location in H via direct map, reduce the frequency of the corresponding tuple in H by 1, and adjust its location in H via a heapify-down operation. We then push the bin number of pt into C and detect its position via direct map in H. We proceed by increasing the corresponding tuple’s frequency by 1 and readjust its position in H via a heapify-up operation. Traditional heap operations heapify-down and heapify-up are extended such that whenever two tuples in H swap positions, their corresponding locations in direct map are also appropriately updated. The algorithm’s complexity is dominated by the complexity of the heapify-up and heapify-down operations, which is O(log BinCount), where BinCount is the total number of bins in the distribution and is typically 256. The space requirements are O(w + BinCount), where w is the length of the sliding window. E. Computing pNN50 Over Sliding Windows Metric pNN50 can be updated in O(1) time in a straightforward manner. As noted previously, we maintain all δNN values in a ring buffer of length w. We maintain a counter cnt50 which counts the number of |δNN| values that are larger than 50. Upon the arrival of a new element xN , the oldest element x0 is extracted from the buffer. If |x0 | > 50, then cnt50 is decreased by 1. If |xN | > 50, then cnt50 is increased by 1. Metric pNN50 is then obtained from cnt50 /w. We note that, as the crucial information per new data point is binary, we maintain a ring buffer of Booleans instead of the actual data points, thereby decreasing memory requirements for the particular algorithm. III. FREQUENCY-DOMAIN METHODS PSD analysis is useful for understanding how power is distributed as a function of frequency. PSD can either be estimated via regression analysis or, more commonly, via fast Fourier transform (FFT) analysis, while both approaches provide comparable results. The HRV measures employed are 1) the total power P , 2) the power within the VLF/LF/HF frequency bands ([0, 0.04 Hz]/[0.04 Hz, 0.15 Hz]/[0.15 Hz, 0.4 Hz]), 3) the normalized power within the LF/HF bands [LFnorm = LF/(P − VLF), HFnorm = HF/(P − VLF)], and 4) the ratio LF/HF. Ratio LF/HF seems to be the more widely used metric. Before applying FFT, the unevenly spaced peak detection signal
STERGIOU AND BALAKRISHNAN: STREAMING UPDATES FOR HEART RATE VARIABILITY ALGORITHMS
1935
be the Fourier coefficients for window (xt−L +1 , . . . , xt ). Then, it holds Xkt =
L −1
xn +t−L · e−2π ik (n −1)/L + xt · e−2π i(L −1)/L
n =1
=
L −1
xn +t−L · e−2π ik n /L
· e2π ik /L
n =1
+ xt · e−2π i(L −1)/L = (Xkt−1 − xt−L + xt ) · e2π ik /L .
typically needs to be resampled. Many studies in the field of HRV report resampling rates between 2 and 4 Hz. In [10], a study is performed on the effects of beat replacement and resampling, and it is concluded that a resampling rate of 7 Hz is appropriate as it allows for the resolution of heart rates of up to 210 bpm. We will now present a method for updating Fourier coefficients in constant time. As will be shown, each coefficient can be updated independent of the rest. This allows for the updating of just the appropriate subset of coefficients corresponding to the desired frequency bands, as opposed to updating all of them at each iteration. We note that the problem of real-time FFT computation has been extensively analyzed in [11] and [12]. We present an algorithm for updating the power in the LF and HF bands, and therefore, the ratio LF/HF in Algorithm 7. We also show how the total power can be updated in a straightforward manner via Parseval’s theorem in Algorithm 8. We note that the proposed algorithms are applicable only when FFT is directly used on the input signal, specifically when windowing is not adopted. A. Updating FFT Coefficients Over Sliding Windows Let {x}t be a streaming sequence of numbers and L be the length of a sliding window. We would like to update the Fourier coefficients of window (xt−L +1 , . . . , xt ) for every new element xt arriving at time t. Let, Xkt =
L −1 n =0
xn +t−L +1 · e−2π ik n /L
We observe that each Fourier coefficient Xkt can be updated from Xkt−1 in constant time. More concretely, if all powers of the Lth primitive root of unity e2π i/L are precomputed then one real addition and one complex multiplication will suffice for updating Xkt . Alternatively, one additional complex multiplication per coefficient will be sufficient for generating the appropriate root power. B. Computing Power Spectral Densities Over Sliding Windows The periodogram of {x}t is defined as ⎧ 1 ⎪ ⎨ 2 |Xk |2 , k = 0, N/2 N P (k) = ⎪ ⎩ 2 |X |2 , k = 1, . . . , N/2 − 1. k N2 P (k) is an estimate of the power spectrum of the time series at frequency fk = ΔkN , where Δ is the sampling period. Since xt is real, then Xk = XN −k . By Parseval’s theorem
N /2
k =0
P (k) =
N −1 1 |xj |2 . N j =0
IV. EVALUATION A. Data Analytics Platform and Experimental Setup All algorithms were implemented in C++ and were executed on the Sprout platform. The core of the platform is the Sprout, a special-purpose built mobile data aggregation and analytics device, depicted in Fig. 1(a). It is based on a 454-MHz ARM926EJ-S processor and communicates with sensors and the user via 802.11n, Bluetooth 3.0, and USB 2.0. The Linux-based platform facilitates data collection and analysis,
1936
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 7, JULY 2014
Fig. 1. Data aggregation and analytics platform. (a) Sprout mobile device for real-time biosignals aggregation and analytics. (b) Zephyr BioHarness 3.0 ECG strap. Fig. 3. Speedup results for the HRVtr index algorithm for (a) 7.7125 ms and (b) 1 ms bin widths.
Fig. 2. Speedup results for SDNN and RMSSD algorithms based on the central and raw moments update algorithm.
and provides a real-time streaming user interface accessible via an on-board web server, while additionally transmitting information to a central server. ECG data were collected with a Zephyr Bioharness 3.0 ECG sensor, shown in Fig. 1(b). Sampling and peak detection was performed at 1 kHz. For the purposes of evaluating the proposed algorithms, ECG data were continuously collected over a 48-h period. The time-domain algorithms are individually benchmarked across 12 different sliding window lengths, specifically of 300, 600, 1200, 2400, 3600, 7200, 14400, 28800, 43200, 57600, 72000, and 86400 s. Twenty-four hours of generated data were collected from each algorithm. We present results on SDNN, RMSSD, HRVtr index with binning of 7.8125 ms (128 Hz) and 1 ms (1 kHz), pNN50, NN Median, and NN interval range. As a representative frequency domain algorithm, we present results on updating the LF/HF ratio. FFT window lengths used are 29 , . . . , 216 . The signal is linearly interpolated and resampled at rates from 1 to 10 Hz before applying FFT. Each algorithm is compared against its “classic” implementation. Specifically, for SDNN, RMSSD, NN Median, NN Interval Range, and pNN50, the metrics are computed as efficiently as possible from the data points on the sliding window without reusing previously computed information. For the HRVtr index, the classic algorithm also maintains the histogram of the NN intervals. B. Results In Fig. 2, we present the performance ratios between the proposed and classical implementations for SDNN and RMSSD. We observe speedups that are linearly increasing with the
Fig. 4.
Speedup results for the NN median update algorithm.
Fig. 5.
Speedup results for the NN interval range update algorithm.
window size, as theoretically expected. In Fig. 3, we present speedups for the HRVtr index algorithm. We observe that the algorithm is more sensitive to the input signal. We note that both the classic and proposed algorithms operate on the NN histogram, whose size depends only on the granularity of the binning and not on the size of the window. The speedups obtained are 5.8× to 13.9× for 7.8125 ms bins, and 36.5× to 79.9× for 1 ms bins. In Fig. 4, we present the speedups for the NN median algorithm. We note that the speedups obtained are significant, and comparable to the ones obtained from the moments algorithms, even though the proposed algorithm performs updates only in expected constant time. In Fig. 5, we present speedups for the NN interval range algorithm, which is based on the maximum/minimum update
STERGIOU AND BALAKRISHNAN: STREAMING UPDATES FOR HEART RATE VARIABILITY ALGORITHMS
1937
In Fig. 8, we present speedups for the LF/HF algorithm for varying window sizes and resampling rates. We note that for the classic algorithm, the dominating part is computing the FFT. For a window of length 216 , the classic algorithm requires 110 020” to complete the 24 h dataset, and thus cannot be run at real time. As the window size increases, the proposed algorithm is relatively faster. Also, we observe that for a fixed window size, when resampling rate increases, the width of the LF and HF bands decrease. As the proposed algorithm only updates the FT coefficients within those bands, it becomes faster as resampling rate increases. Fig. 6.
Speedup results for the pNN50 update algorithm.
V. CONCLUSION In this paper, we studied the most popular HRV quantification methods and proposed streaming algorithms that maximally utilize previously computed information without altering the output of the methods. We demonstrated speedups of more than two orders of magnitude for typical use-case scenarios that translate to significantly longer battery lives for mobile systems that adopt them. For the most computationally heavy frequencydomain methods, our algorithms allow for real-time updates in cases where it was not previously possible. REFERENCES
Fig. 7.
Aggregate speedup results for all time-domain algorithms.
Fig. 8. Speedup results for Fourier transform LF/HF update algorithm. The original signal was linearly interpolated and resampled at rates ranging from 1 to 10 Hz.
algorithm. In Fig. 6, we present speedups for the pNN50 algorithm. This algorithm is the most straightforward to optimize and indeed yields the largest speedups, due in large to the compressed way the sliding window is maintained. In Fig. 7, we present the overall speedup obtained when executing all the time-domain algorithms. We observe that, even though the speedups of the HRVtr algorithm are lower relative to the rest of the time-domain algorithms, the algorithm is fast in absolute execution times as it only operates on an input whose size does not depend on the sliding window length.
[1] N. Dagres, H. Kottkamp, C. Piorkowski, S. Weis, A. Arya, P. Sommer, K. Bode, J.-H. Gerds-Li, D. T. Kremastinos, and G. Hindricks, “Influence of the duration of Holter monitoring on the detection of arrhythmia recurrences after catheter ablation of atrial fibrillation: Implications for patient follow-up,” Int. J. Cardiol., vol. 139, no. 3, pp. 305–306, 2010. [2] S. S. Lobodzinski and M. M. Laks, “New devices for very long-term ECG monitoring,” Cardiol. J., vol. 19, no. 2, pp. 210–214, 2012. [3] B. Cinaz, R. La Marca, B. Arnrich, and G. Tr¨oster, “Towards continuous monitoring of mental workload,” presented at the 5th Int. Workshop Ubiquit. Health Wellness, Copenhagen, Denmark, 2010. [4] “Vital Connect.” (2013). [Online]. Available:http://www.vitalconnect. com/ [5] “CardioLeaf.” (2013). [Online]. Available: http://www. clearbridgevitalsigns.com/ [6] B. Gyselinckx, C. Van Hoof, J. Ryckaert, R. F. Yazicioglu, P. Fiorini, and V. Leonov, “Human++: Autonomous wireless sensors for body area networks,” in Proc. IEEE Custom Integr. Circuits Conf., 2005, pp. 13–19. [7] U. R. Acharya, K. P. Joseph, N. Kannathal, C. M. Lim, and J. S. Suri, “Heart rate variability: A review,” Med. Biol. Eng. Comput., vol. 44, no. 12, pp. 1031–1051, 2006. [8] M. P. Tarvainen, J.-P. Niskanen, J. Lipponen, P. Ranta-Aho, and P. Karjalainen, “Kubios HRV—A software for advanced heart rate variability analysis,” in Proc. 4th Eur. Conf. Int. Federat. Med. Biol. Eng., 2009, pp. 1022–1025. [9] A. Camm, M. Malik, J. T. Bigger, Jr, G. Breithardt, S. Cerutti, R. J. Cohen, P. Coumel, E. L. Fallen, H. L. Kennedy, R. E. Kleiger, F. Lombardi, A. Malliani, A. J. Moss, J. N. Rottman, G. Schmidt, P. J. Schwartz, and D. H. Singer, “Heart rate variability: Standards of measurement, physiological interpretation and clinical use. Task force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology,” Circulation, vol. 93, no. 5, pp. 1043–1065, 1996. [10] G. Clifford and L. Tarassenko, “Quantifying errors in spectral estimates of HRV due to beat replacement and resampling,” IEEE Trans. Biomed. Eng., vol. 52, no. 4, pp. 630–638, Apr. 2005. [11] B. Farhang-Boroujeny and Y. Lim, “A comment on the computational complexity of sliding FFT,” IEEE Trans. Circuits Syst. II: Analog Digital Signal Process., vol. 39, no. 12, pp. 875–876, Dec. 1992. [12] P.-C. Lo and Y.-Y. Lee, “Real-time implementation of the moving FFT algorithm,” Signal Process., vol. 79, no. 3, pp. 251–259, 1999.
Authors’ photographs and biographies not available at the time of publication.
