Beat detection and HRV time series

For heart rate variability (HRV) analysis, inter-beat interval (IBI) data is required. Such data can be extracted from electrocardiogram (ECG) recording as time intervals between successive ECG R-waves, i.e. RR time intervals. Similarly IBI data can be extracted from photoplethysmography (PPG) recording as time intervals between successive pulsations. The detection of beats, beat detection accuracy and formation of HRV time series is discussed below.

The aim in heart rate variability analysis is to examine the sinus rhythm modulated by the autonomic nervous system. Therefore, one should technically detect the occurrence times of the SA-node action potentials, which initiate every heart beat. In practical applications, this is not however possible. Instead, an the electrocardiogram (ECG) is recorded by placing two or more electrodes on skin contact and heart beats are detected from the ECG. The nearest observable activity in the ECG compared to SA-node firing is the P-wave resulting from atrial depolarization (see Fig. 1). The signal-to-noise ratio of the P-wave is, however, clearly lower than that of the strong QRS complex which results primarily from ventricular depolarization. Therefore, the heart beat period is commonly evaluated as the time difference between the easily detectable QRS complexes. A typical QRS detector consists of a preprocessing part followed by a decision rule. Several different QRS detectors have been proposed within last decades^{Thakor et al. 1983, Pahlm et al. 1985, Pan & Tompkins 1985, Hamilton et al. 1984, Friesen et al. 1990}.

Figure 1: Electrophysiology of the heart (redrawn from ^{Malmivuo & Plonsey 1995}). The different waveforms for each of the specialized cells found in the heart are shown. The latency shown approximates that normally found in the healthy heart.

The accuracy of the R-wave occurrence time estimates is often required to be 1–2 ms, and thus, the sampling frequency of ECG should be at least 500–1000 Hz ^{Task Force 1996}. If the sampling frequency of ECG is less than 500 Hz, the errors in R-wave occurrence times can cause critical distortion to HRV analysis results, especially to spectrum estimates ^{Merri et al. 1990}. The distortion of the spectrum is even bigger if the overall variability in heart rate is small ^{Pinna et al. 1994}. The estimation accuracy can however be improved by interpolating the QRS complex e.g. by using a cubic spline interpolation ^{Daskalov et al. 1997} or some model based approach ^{Bragge et al. 2005}.

Figure 2: Normal PPG end ECG signal and definitions of pulse transmit time (PTT) and pulse to pulse interval (PP).

Depending on the pulse wave velocity and the vascular path from the heart, there is a delay between each QRS complex and the onset of its corresponding pulse wave. The delay is called pulse transit time (PTT) and is negatively correlated with blood pressure, arterial stiffness, and age. Physiological variability in PTT causes deviation between the PP-intervals and the RR-intervals. Since PP and RR intervals are not equal, it is always better to use the term pulse rate variability (PRV) rather than heart rate variability (HRV) when PPG measurement is used. Usability and accuracy of PRV as an estimate of HRV has been widely studied. A good review on the topic can be found from ^{Schafer et al 2013}, with the following concluding remarks:

• PRV as an estimate of HRV has been proved to be sufficiently accurate only for healthy (and mostly younger) subjects at rest.
• Moderate physical or mental stress tends to diminish agreement between PRV and HRV to an extent that is or is not acceptable.
• Physically more active states, such as walking or physical exercising, the agreement between PRV and HRV often becomes insufficient, mostly due to motion artifacts.

Pulse wave detector of Kubios HRV software is based on the matched filtering approach. Firstly, maximum of 1st derivative representing the steepest part of the pulse wave is used for initial pulse location estimation. Secondly, a template for the pulse wave (and matched filter) is constructed using the initial pulses. Decision of the final pulse wave locations are defined by comparing the filtered signal against varying threshold and comparing normalized error between the template and the PPG signal. Allowed normalized error between template and pulse wave under inspection can be adjusted in software preferences. That is, the smaller the acceptance threshold percent is the more similar the pulse wave have to be with the template in order to be accepted. The accuracy of the pulse wave detection algorithm is shown in Fig. 3. The left panel showing the Bland-Altman plot illustrating the agreement between detected PP intervals and corresponding RR intervals during a resting measurement. The right panel shows error percentages of commonly used heart rate variability parameters estimated from PP interval compared to RR interval time series. Used dataset contains 20 healthy volunteers with wide age scale (20 to 50 years). The error between the RR and PP intervals was -0.01±5.16 ms (mean ± SD). This ±5 ms error in heart beat detection produces approximately ±10 % maximum errors to the HRV parameters.

Figure 3: Bland-Altman plot for the difference between PP intervals and RR intervals extracted from resting measurements (left axes). The errors in standard HRV parameters when computed from PRV compared to HRV time series (right axes); where red lines show the median, the box shows the 25-75 percentiles, and the whiskers show the most extreme error values.

After the beat occurrence times (i.e. QRS complex or pulse wave fiducial points) have been detected, an inter-beat interval (IBI) or HRV time series can be derived. Instead of IBI, we will here use the term RR interval, which refers to time interval between successive ECG R-wave occurrence times. The n‘th RR interval is obtained as the difference between the R-wave occurrence times RR_n=t_n − t_n-1. In some context, normal-to-normal (NN) may also be used when referring to these intervals indicating strictly intervals between successive QRS complexes resulting from SA-node depolarization ^{Task Force 1996}. In practice when analysing normal sinus rhythm, the NN and RR intervals are the same, and thus, the term RR is preferred here.

The time series constructed from all available RR intervals is not equidistantly sampled, but has to be presented as a function of time, i.e. as values (t_n, RR_n). This fact has to be taken into account before frequency-domain analysis. In general, three different approaches have been used to get around this issue ^{Task Force 1996}. The simplest approach that have been adopted e.g. in ^{Baselli et al. 1987} is to assume equidistant sampling and calculate the spectrum directly from the RR interval tachogram (RR intervals as a function of beat number), see Fig. 4 middle panel. This assumption can however cause distortion into the spectrum ^{Mateo et al. 2000}. This distortion becomes substantial when the variability is large in comparison to the mean RR interval length. Furthermore, the spectrum can not be considered to be a function of frequency, but rather a function of cycles per beat ^{DeBoer et al. 1984}, see Fig. 5 bottom panel. One choice for the interpolation method is the cubic spline interpolation ^{Mateo et al. 2000}. After interpolation, regular spectrum estimation methods can be applied. The third general approach is to apply methodology, which are designed for analysing non-equidistantly sampled data. Such a method is for example the Lomb-Scargle periodogram, which computes the periodogram spectrum estimate for non-equidistantly sampled data ^{van Dongen et al. 1999}.

Figure 4: Derivation of two HRV time series from ECG: the interval tachogram (middle panel) and interpolated RR interval series (bottom panel).