Preprocessing of HRV data

Careful preprocessing of HRV data is very important before applying any HRV analysis methods. Any artefact in the RR interval time series may significantly interfere the analysis of these signals. The artefacts within HRV signals can be divided into technical and physiological artefacts. The technical artifacts include missing, extra or misaligned beat detections. These artefacts may be due to measurement noise or due to inaccuracy of the detection algorithm. Physiological artefacts include ectopic beats and arrhythmic events. In order to avoid the interference of such artefacts, the ECG recording and the corresponding event series should always be manually checked for artefacts and only artefact-free sections should be included in the analysis ^{Task Force 1996}. Alternatively, if the amount of artefact-free data is insufficient, proper interpolation methods can be used to reduce these artefacts, see e.g. ^{Lippman et al. 1993, Lippman et al. 1994, Mateo et al. 2003}.

In addition to artefacts, slow changes in mean HR during the recording can have undesirable effects on certain heart rate variability analysis parameters. Such slow nonstationarities are characteristic for HRV signals and should be considered before the analysis. The origin of nonstationarities in HRV are discussed e.g. in ^{Berntson et al. 1997}. The slow nonstationary trends can be easily removed from the HRV time series before analysis. Detrending can be performed using first order ^{Litvack et al. 1995, Mitov et al. 1998} or higher order polynomial ^{Porges et al. 1990, Mitov et al. 1998} models. In addition, Kubios HRV software includes an advanced detrending procedure originally presented in ^{Tarvainen et al. 2002}. This approach is based on smoothness priors regularization and is described in detail below.

Noise detection

The first HRV preprocessing step in Kubios HRV software is noise detection, which automatically identifies noise segments based on the raw ECG data (if available) and from the interbeat interval data (RR or pulse-to-pulse intervals). Automatic noise detection is available for example in the Kubios HRV Scientific and Sport and Exercise Analytics products.

An example of automatic noise detection for ECG data is shown in the figure below (a screenshot from Kubios HRV Scientific software). Automatic noise detection is extremely helpful for those analyzing long-term ambulatory ECG, PPG or RR recordings, where signal quality is not good throughout the recording. Please note that noise detection aims to identify noise segments which distort several consecutive beat detections, and thus, would affect on HRV analysis accuracy. Individual intermittent abnormal beat intervals (e.g. ectopic beats) are not typically marked as noise, since these can be reliably corrected by the automatic beat correction which is described next.

Figure 1: Illustration of automatic noise detection (screenshot taken from Kubios HRV Scientific software).

Beat correction

Kubios HRV software provides excellent tools for HRV preprocessing. First, beat detections can be corrected manually when raw ECG (or PPG) data is available. Such beat correction provides true beat intervals but is reasonable only when when data includes a limited number of misdetections. Secondly, the software includes two beat correction methods: 1) Threshold based beat correction, and 2) Automatic beat correction. Details of these beat correction methods are given below.

Threshold based beat correction algorithm

The threshold based beat correction algorithm compares every RR interval value against a local average interval. The local average is obtained by median filtering the RR interval time series, and thus, the local average is not affected by single outliers in RR interval time series. If an RR interval differs from the local average more than a specified threshold value, the interval is identified as an artefact and is marked for correction. The threshold value can be selected from:

Very low: 0.45 sec (threshold in seconds)
Low: 0.35 sec
Medium: 0.25 sec
Strong: 0.15 sec
Very strong: 0.05 sec
Custom, for setting a custom threshold in seconds

For example, the “Medium” correction level will identify all RR intervals that are larger/smaller than 0.25 seconds compared to the local average. The correction is made by replacing the identified artefacts with interpolated values using a cubic spline interpolation. It should be noted, that these threshold are adjusted with mean heart rate. That is, thresholds shown above are for HR of 60 beats/min, but for higher HR the thresholds are smaller (because the variability is expected to decrease when HR increases) and vice versa for lower heart rates.

Because the artefacts are identified by simple thresholding and normal heart rate variability level is highly individual, the correction level should be adjusted individually as follows. First, identify if there are any beat intervals that should be corrected. If there are such intervals, then select the lowest possible correction level, which corrects the abnormal beats but does not over correct the data.

Automatic beat correction algorithm

The automatic beat correction algorithm has been described in detail in Lipponen & Tarvainen 2019. In automatic beat correction algorithm, artefacts are detected from dRR series, which is a time series consisting of differences between successive RR intervals. The dRR series provides a robust way to separate ectopic and misplaced beats from the normal sinus rhythm. To separate ectopic and normal beats, time varying threshold (Th) is used. To ensure adaptation to different HRV levels, the threshold is estimated from the time-varying distribution of dRR series. For each beat, quartile deviation of the 90 surrounding beats is calculated and multiplied by factor 5.2. Beats within this range cover 99.95% of all beats if RR series is normally distributed. However, RR interval series is not often normally distributed, and thus, also some of the normal beats exceed the threshold. Therefore, decision algorithm is needed to detect artefact beats.

Table 1: Performance of the automatic artefact correction algorithm in detecting normal sinus rhythm beats, and simulated missed, extra and misaligned (with three different levels of displacement) beats.

Kubios HRV beat correction algorithm accuracy

Ectopic beats form negative-positive-negative (NPN) or positive-negative-positive (PNP) patterns to the dRR series. Similarly long beats form positive-negative (PN) and short beats negative-positive (NP) patterns to the dRR series. Only dRR segments containing these patterns are classified as artefact beats. Missed or extra beats are detected by comparing current RR value with median of the surrounding 10 RR interval values (medRR). A missed beat is detected if current RR interval (RR(i)) satisfies condition

$\begin{equation*} \left|\frac{RR(i)}{2}-medRR(i)\right|<2 Th \end{equation*}$

and an extra beat is detected if two successive RR intervals (RR(i) and RR(i+1)) satisfies condition

$\begin{equation*} \left|RR(i)+RR(i+1)-medRR(i)\right|<2 Th. \end{equation*}$

Correction of detected artefacts: Detected ectopic beats are corrected by replacing corrupted RR times by interpolated RR values. Similarly too long and short beats are corrected by interpolating new values to the RR time series. Missed beats are corrected by adding new R-wave occurrence time and extra beats are simply corrected by removing extra R-wave detection and recalculating RR interval series.

Trend removal

To avoid the effect of slow nonstationary trends (slow changes in mean HR basically) over HRV analysis results, such trends by default are removed in Kubios HRV software using a smoothness priors method. The theory behind the smoothness priors detrending method is described shortly here. For more details, see ^{Tarvainen et al. 2002}.

Let $z\in\R^N$ denote the RR interval time series which can be considered to consist of two components

$\begin{equation*} z = z_{\rm{stat}}+z_{\rm{trend}} \end{equation*}$

where $z_{\rm{stat}}$ is the nearly stationary RR interval series of interest, $z_{\rm{trend}}$ is the low frequency aperiodic trend component, and $N$ is the number of RR intervals. Suppose that the trend component can be modeled with a linear observation model as

$\begin{equation*} z_{\rm{trend}} = H\theta + e \end{equation*}$

where $H\in\R^{N\times p}$ is the observation matrix, $\theta\in\R^p$ are the regression parameters, and $e$ is the observation error. The task is then to estimate the parameters by some fitting procedure so that $\hat{z}_{\rm{trend}} = H\hat{\theta}$ can be used as the estimate of the trend. The properties of the estimate depend strongly on the properties of the basis vectors (columns of the matrix $H$ ) in the fitting. A widely used method for the solution of the estimate $\hat{\theta}$ is the least squares method. However, a more general approach for the estimation of $\hat{\theta}$ is used here. That is, the so-called regularized least squares solution

$\begin{equation*} \hat{\theta}_\lambda = \arg \min_\theta \left\{ \Vert z-H\theta \Vert^2 + \lambda^2\Vert D_2 (H\theta)\Vert^2 \right\} \end{equation*}$

where $\lambda$ is the regularization parameter and $D_2$ indicates the discrete approximation of 2nd order derivative. This is a modification of the ordinary least squares solution, dragging the solution to the direction in which the norm $\Vert D_2(H\theta)\Vert$ gets smaller. In this way, prior information about the predicted trend $H\theta$ can be implemented to the estimation. The solution of previous equation can be written in the form

$\begin{equation*} \hat{\theta}_\lambda = (H^TH + \lambda^2H^TD_2^TD_2H)^{-1} H^T z \end{equation*}$

and the estimate for the trend which is to be removed can be written

$\begin{equation*} \hat{z}_{\rm{trend}} = H\hat{\theta}_\lambda. \end{equation*}$

The selection of the observation matrix $H$ can be implemented according to some known properties of the data $z$ . However, here we use a trivial choice of identity matrix $H = I\in\R^{N\times N}$ . In this case, the regularization part of the minimized function simply draws the solution towards a curve for which the 2nd order difference is zero, i.e. towards first order linear curve. With these specific choices, the detrended nearly stationary RR series can be written as

$\begin{equation*} \hat{z}_{\rm{stat}} = z- H\hat{\theta}_\lambda = (I-(I+\lambda^2D_2^TD_2)^{-1}) z. \end{equation*}$

This detrending method works like a time-varying highpass filter, where the cutoff frequency can be adjusted by changing the smoothing parameter $\lambda$ . In Kubios HRV Premium software, the cutoff frequency related to the given smoothing parameter is presented. The smoothing parameter is by default selected to remove only VLF frequency components from the RR time series.

References

G.G. Berntson, J.T. Bigger Jr., D.L. Eckberg, P. Grossman, P.G. Kaufmann, M. Malik, H.N. Nagaraja, S.W. Porges, J.P. Saul, P.H. Stone, and M.W. Van Der Molen. Heart rate variability: Origins, methods, and interpretive caveats. Psychophysiol, 34:623–648, 1997.
N. Lippman, K.M. Stein, and B.B. Lerman. Nonlinear predictive interpolation: a new method for the correction of ectopic beats for heart rate variability analysis. J Electrocardiol, 26:S14–S19, 1993.
N. Lippman, K.M. Stein, and B.B. Lerman. Comparison of methods for removal of ectopy in measurement of heart rate variability. Am J Physiol, 267(1):H411–H418, July 1994.
J.A. Lipponen and M.P. Tarvainen. A robust algorithm for heart rate variability time series artefact correction using novel beat classification. J Med Eng Technol, 43(3):173-181, 2019.
D.A. Litvack, T.F. Oberlander, L.H. Carney, and J.P. Saul. Time and frequency domain methods for heart rate variability analysis: a methodological comparison. Psychophysiol, 32:492–504, 1995.
J. Mateo and P. Laguna. Analysis of heart rate variability in the presence of ectopic beats using the heart timing signal. IEEE Trans Biomed Eng, 50(3):334–343, March 2003.
I.P. Mitov. A method for assessment and processing of biomedical signals containing trend and periodic components. Med Eng Phys, 20(9):660–668, November-December 1998.
S.W. Porges and R.E. Bohrer. The analysis of periodic processes in psychophysiological research. In J.T. Cacioppo and L.G. Tassinary, editors, Principles of Psychophysiology: Physical Social and Inferential Elements, pages 708–753. Cambridge University Press, 1990.
M.P. Tarvainen, P.O. Ranta-aho, and P.A. Karjalainen. An advanced detrending method with application to HRV analysis. IEEE Trans Biomed Eng, 49(2):172–175, February 2002.
Task force of the European society of cardiology and the North American society of pacing and electrophysiology. Heart rate variability – standards of measurement, physiological interpretation, and clinical use. Circulation, 93(5):1043–1065, March 1996.