RapidHRV: an open-source toolbox for extracting heart rate and heart rate variability

Preprocessing

First, data is upsampled with cubic spline interpolation (3rd order polynomial; default = 1 kHz) to increase temporal accuracy of peak detection. To mitigate potential long-term drifts in the signal, the pipeline then applies a high pass butterworth filter (0.5 Hz) across the input data. Finally smoothing with a savitzky-golay filter (3^rd order polynomial; default = 100 ms) is applied to reduce spiking (sharp increases in the signal caused by artifacts such as motion) whilst retaining temporal precision.

Outlier detection

The last phase of the pipeline is to pass measures derived from peak extraction to outlier rejection (Fig. 3). This is applied at the level of the sliding window. If a window is declared an outlier, heart rate and heart rate variability measures are removed from the cleaned time series. By default, RapidHRV returns both the cleaned and the uncleaned time series. In addition to default parameters listed below, the package has optional arguments embedded to allow users to override these presets. Given that not all users may be entirely comfortable manually adjusting these, RapidHRV additionally contains semantically-labeled arguments as inputs for outlier constraints (‘liberal’, ‘moderate’ [default], and ‘conservative’; corresponding parameters are parenthesized under Biological Constraints and Statistical Constraints).

Biological Constraints. RapidHRV first applies restrictions to exclude data that are highly unlikely given known physiology:

1. Screening for sufficient peaks in a window (default: number of peaks > (window width/5) + 2), floored at 3; default at 10 s = 3 peaks). This is primarily for computational applicability and efficacy, screening data prior to further processing. The minimum number of peaks required to enable calculation of RMSSD is 3. As such, this is also applied to the uncleaned time series.

2. Minimal and maximal heart rate (‘moderate’ [default]: 30 > BPM > 190; ‘liberal’: 20 > BPM > 200; ‘conservative’: 40 BPM > 180). These boundaries were based on typical heart rate at rest and during exercise in the healthy population (Pierpont & Voth, 2004; Sandvik et al., 1995; Savonen et al., 2006).

3. Minimal and maximal heart rate variability (‘moderate’ [default]: 5 > RMSSD > 262; ‘liberal’: 0 > RMSSD > 300; ‘conservative’: 10 > RMSSD > 200). Default arguments correspond to the minimum/maximum 2nd/98th percentiles of resting RMSSD across ages 16–89 years (van den Berg et al., 2018).

Statistical constraints. RapidHRV next applies statistical constraints to account for noisy data that may otherwise appear to provide measures within the range of known physiology:

• 4. Median absolute deviation (MAD) of peak heights (distance from minimum value of signal in window; i.e., 0) and prominence (amplitude from baseline height; ‘moderate’ [default] = 5 MAD units; ‘liberal’ = 7; ‘conservative’ = 4). Unlike Z-scoring, this quantifies each peak’s height and prominence in a given window in terms of its deviation from the median value in the same window (for a discussion of median absolute deviation see Leys et al., 2013). Applying these constraints to height and prominence helps exclude windows with noise-driven inaccuracies in peak detection.

• 5. Median absolute deviation of interbeat intervals (‘moderate’ [default] = 5 MAD units; ‘liberal’ = 7; ‘conservative’ = 4). This was also implemented to account for inaccuracies in peak detection, either where spiking may cause detection of an irrelevant peak shortly after e.g., an R wave, or low signal to noise ratio may result in missing relevant peaks.

• 6. Time from the first peak to the last peak does not recede 50% of the fixed window width. This is to ensure the user that the actual length of time for extracting HR/HRV is not less than half of that which is specified in the window width argument. Given debates surrounding adequacy of different window lengths for HRV extraction (Baek et al., 2015; Munoz et al., 2015), this was implemented primarily as a theoretical constraint (rather than for just cleaning per se) to ensure the user is not provided data that deviated significantly from their specified window.

Analysis can be executed with one line of code, which returns a pandas DataFrame (McKinney, 2012; Reback et al., 2022) containing the analyzed data.

Visualization

To allow for selected manual inspection, we have also implemented optional interactive visualizations via matplotlib (Hunter, 2007) which allow the user to plot the time course of heart rate and heart rate variability. The user can then select and view specific data points to see the window of extraction (Fig. 4). We have provided an example pipeline below:

import rapidhrv as rhv

example_signal = rhv.get_example_data() # Load example data

preprocessed = rhv.preprocess(example_signal) # Preprocess data

analyzed = rhv.analyze(preprocessed) # Analyze data

rhv.visualize(analyzed) # Visualize data

Datasets

To validate the above pipeline we subjected it to a series of tests across both simulated and real data (Table 1). We first started by testing RapidHRV’s estimations in two sets of simulated data (PPGSynth; Tang et al., 2020). Next, we ran validation in real data across successively noisier modalities: electrocardiography (ECG), finger infrared photoplethysmography (IR PPG), and wrist PPG data. Database information and code generated in validation tests are available through the open science framework (https://osf.io/7zvn9/).

1. Simulations Across Noise and Sampling Rates. We first took the pipeline forward to validation using simulated photoplethysmography (PPG) data from PPGSynth (Tang et al., 2020) in MATLAB. This allowed us to test how well RapidHRV recovered known parameters under specified conditions, such as sampling rate and noise. We produced 5 min long synthetic datasets (1,000 Hz), each of which varied according to the following cardiac features:

• Mean heart rate (BPM range 60–120, increments of 5). These were selected to allow us to check the sensitivity of the pipeline for detecting values across typical resting heart rate, as well as elevations of these values.

• Heart rate variability (RMSSD range 0–100, increments of 5). Again, these were selected based on typical resting heart rate variability and moderate increases/decreases.

Following simulation of the data, we introduced noise via:

• White gaussian noise filtering (signal to noise ratios = 0.01, 10, 20, 30, 40, 50, and 60 dB). Here, we selected a range of values starting from near-zero to ascertain at which level of noise the pipeline could no longer recover parameters (with near-zero expected to prevent the extraction of any meaningful metrics).

• Downsampling (frequencies = 20, 50, 100, and 250 Hz). These were selected: (A) to capture the range of sampling rates currently used in photoplethysmography studies (typically >= 20 Hz); and (B) because prior work has suggested sensitivity to RMSSD deteriorates dramatically between 20–100 Hz (Choi & Shin, 2017). This offered the opportunity to test the effects of cleaning on what would otherwise be considered poor quality data.

2. Simulated ‘Anxiety’. Again, we simulated data using PPGSynth. However, simulations were based on prior anxiety literature to emulate experiment-specific effects. We wanted to validate the sensitivity of the pipeline across specified sampling rates and levels of noise but in the context of a known psychophysiological effect. For this, we used an estimated Cohen’s d of 0.384, which was derived from a previous within-subjects threat-of-shock study (RMSSD, t₂₅ = 1.96, polarity flipped for readability; Gold, Morey & McCarthy, 2015). To reflect decisions in experimental design, we ran a power calculation using the ‘pwr’ package (Champely et al., 2017) in R (R Core Team, 2013), leading us to simulate a sample of N = 171 (α = 0.001, 1 − β = 0.95). These ‘experiments’ were simulated 10 times to ascertain confidence intervals around estimated effect sizes.

For each ‘subject’, we generated a 5 min simulated time series (1 kHz) using typical resting heart rate and heart rate variability (BPM: μ = 74, σ = 13, Savonen et al., 2006; RMSSD: μ = 23, σ = 7, bounded between 5–262, Nunan, Sandercock & Brodie, 2010). These simulations were intended to emulate the ‘safe’ condition (no anxiety induction) in a threat-of-shock study. Next, we simulated another time series for each subject that deviated from their ‘safe’ RMSSD (heart rate held constant) based on our effect size estimate (d = 0.384). This emulated our ‘threat’ (anxiety induction) condition and contained general reductions in heart rate variability. Each time series was finally submitted to downsampling and noise filtering using the same parameters as in the previous simulations (white gaussian noise filtering to 0.01, 10, 20, 30, 40, 50, 60 dB; downsampling to 20, 50, 100, and 250 Hz).

3. Estimation via ECG. Following simulations, we wanted to clarify that our package was able to adequately extract measures from one of the highest standards for heart rate variability recordings, ECG. Here, we tested whether our package was able to recapitulate known age-related effects of heart rate variability. For this, we used the Fantasia database (Iyengar et al., 1996). This consisted of 40 subjects (20 ‘Young’ Age Range = 21–34 years; 20 ‘Old’ Age Range = 68–81 years) watching the movie ‘Fantasia’ (duration ≈ 2 h) whilst undergoing ECG recordings. Automated peak detection was previously run on this dataset, with every beat annotation verified by visual inspection. The full dataset and description is available on PhysioNet (https://physionet.org/content/fantasia/1.0.0/; Goldberger et al., 2000).

4. Estimation via finger PPG. We next took the pipeline forward to validation in a modality considered relatively noisier than ECG, finger photoplethysmography. Here, we used a dataset of 39 subjects (Age: Mean = 22.67; Range = 18–38 years; demographics reported prior to N = 1 exclusion) watching 2 × 28 min blocks of documentary and horror video clips undergoing finger IR PPG recording (1,000 Hz, de Groot, Kirk & Gottfried, 2020). This allowed us to contrast psychological conditions of the experiment, testing whether RapidHRV was able to detect effects of anxiety. Moreover, in the original study, data had been preprocessed and analyzed using a commercially available software (Labchart; ADInstruments, Sydney, Australia; analyzed using built-in ‘Peak Analysis’ module). This allowed us to benchmark RapidHRV against another software. The full dataset and description is available via the Open Science Framework (https://osf.io/y76p2/).

5. Estimation via wrist PPG. Finally, we analyzed the PPG-DaLiA dataset, which consists of 15 subjects (Age: Mean = 30.60 years; Range = 21–55) completing various activities whilst having wrist PPG recorded with an Empatica E4 device (Hz = 64) and simultaneous ECG measures with a RespiBAN (Hz = 700, Reiss et al., 2019). Over the course of 2.5 h, participants engaged in a range of activities designed to elicit low and high motion. These were: sitting and reading; ascending/descending stairs; 1 v 1 table soccer; cycling on pavements and gravel; driving a car; lunch break (queuing, purchasing, and eating food in a cafeteria); walking; and working at a desk (typically on a computer). This allowed us to test whether RapidHRV was able to extract heart rate and variability measures from wrist PPG and how these compared to ECG measurements. Moreover, this enabled us to highlight under what conditions estimations are optimal. For the full dataset and description, see Reiss et al. (2019).

Analyses

Unless otherwise stated, all analyses were conducted using RapidHRV’s default arguments: window width = 10 s; window movement = 10 s; outlier method = ‘moderate’ ((peak/height median absolute deviation = 5, interbeat interval median absolute deviation = 5, BPM range = 30–190, RMSSD range = 5–262); minimum window successful extraction = 5 s, minimum amplitude for peak detection = 50, minimum distance between peaks = 250 ms; for ECG data, ecg_prt_clustering = True). To assess performance across datasets, we used: visualizations; intraclass correlation coefficients (ICC; two-way mixed effects, absolute agreement, single measure); root-mean-square-error (RMSE); and sensitivity to experimental effects (Cohen’s d). For ICC values, we used the following semantic labels for interpretation: ICC < 0.5 as ‘poor’, 0.5 < ICC < 0.75 as ‘moderate’; 0.75 < ICC < 0.90 as ‘good’, and 0.90 < ICC as ‘excellent’ (Koo & Li, 2016).

In our wrist PPG dataset, we also calculated motion estimates as a proxy for the severity of noise present in PPG across conditions. We derived mean jerk magnitude (square root of the sum of squared changes in acceleration; Eager, Pendrill & Reistad, 2016), which aimed to pick up on oscillatory (‘jerk’) motions in the wrist for each condition. Across all axes (x, y, and z; m/s^2), jerk vector, j, was calculated as the change in acceleration, a, divided by the change in time, t, between samples:

$$\overrightarrow{\mathit{j}}\mathbf{(}\mathit{t}\mathbf{)}=\frac{\mathrm{\Delta }\overrightarrow{a}(t)}{\mathrm{\Delta }t}$$

Finally, we derived jerk magnitudes by calculating the square root of the sum of squared jerks, j, for axes x, y, and z:

$$|\overrightarrow{\mathit{j}}|=\sqrt{j_x^{\mathbf{2}}+j_y^{\mathbf{2}}+j_z^{\mathbf{2}}}$$

For each subject, jerk magnitudes were averaged across time for each condition.

Parameter recovery in simulated PPG

RapidHRV was able to accurately recover heart rate across most sampling frequencies and noise in our initial simulations. Accurate detection of BPM primarily started to degrade when signal to noise ratios were less than 10 dB (Fig. 5; Table 2). RapidHRV cleaning provided improvements in simulations with a signal to noise ratio of 10 dB.

Performance in recovery of heart rate variability was again primarily based on signal to noise ratio. At 20 dB RapidHRV recovery of RMSSD was good-to-excellent for higher sampling rates (>=100 Hz), whereas lower sampling rates (<100 Hz) required slightly lower levels of noise (>30 dB) for excellent recovery. RapidHRV cleaning provided clear improvements when signal to noise ratio was below 30 dB (Fig. 3).

Sensitivity to ‘Anxiety’ in simulated PPG

RapidHRV’s ability to estimate simulated effects was again primarily impacted by the level of SNR, where changes in RMSSD were not reliably detected at 0.01 and 10 dB (Fig. 6; Table 3). Effects were detected at signal to noise ratios of 20 and 30 dB, but this was estimated at around half that of the true value. Maximal effect sizes plateaued following 40 dB for uncleaned data and 30 dB for cleaned data. Moreover, reliable detection of effects in some scenarios (e.g., 50 Hz, 10 dB) appeared dependent on cleaning. This was consistent with our previous validation results, in that cleaning was beneficial at lower signal to noise ratios and sampling rates, but was not necessary (or could be relaxed) in cleaner data.

Estimation via ECG

As expected, visual inspection suggested ECG data to hold a high signal to noise ratio. In line with our simulation findings, we adjusted the outlier rejection method accordingly so as not to excessively exclude data (i.e., ‘liberal’). The dataset had previously been analyzed using an automated peak detection algorithm, with every beat annotation verified by visual inspection (Iyengar et al., 1996). Subject-specific heart rate estimates were not available from the original database. However, RapidHRV heart rate estimations were able to recapitulate sample-wide summary statistics (Fig. 7; Table 4).

The analytical method used for extraction of heart rate variability in the original database (power spectral density analysis) was inconsistent with RapidHRV’s time-domain heart rate variability measure (RMSSD). Despite this discrepancy, RapidHRV was able to capture previously reported effects of age on heart rate variability (Iyengar et al., 1996; estimated Cohen’s d of short-term heart rate variability (i.e., α_s) ≈ 1.32), such that—in cleaned data—younger participants demonstrated higher RMSSD (M = 58.94, SD = 28.82) than older participants (M = 25.54, SD = 13.02, Cohen’s d = 1.49; Fig. 6). Effects were not apparent in the uncleaned data (Cohen’s d = −0.31).

Estimation via finger PPG

Benchmarking

Overall, there was excellent agreement between RapidHRV and previous estimates (de Groot, Kirk & Gottfried, 2020; implemented using LabChart, ADInstruments, Sydney, Australia) of BPM (ICC > 0.99; Fig. 8). For heart rate variability, there was good agreement between the two when using the cleaned time series (ICC = 0.89), but poor agreement when using the uncleaned time series (ICC = 0.32).

Estimation via wrist PPG during activities

We first derived motion estimates (mean jerk magnitude; Eager, Pendrill & Reistad, 2016) as a proxy for the severity of noise present in PPG. This enabled us to split conditions into low (Reading, Working, Lunch Break, Driving) and high motion activities (Table Soccer, Stairs, Walking, Cycling; Fig. 9).

In our initial analysis using default arguments, RapidHRV was not able to produce estimates across many of the activities for both cleaned and uncleaned data (conditions estimated: M = 22.5%, SD = 5.18%). Visual inspection confirmed this was due to high levels of noise and variable peak amplitude present in the signal. Therefore, we adjusted window movement to 1 s (to increase the number of extraction windows in each condition), reduced minimum amplitude threshold for peak detection (30), and—in line with our simulation results above (Fig. 3)—tightened outlier rejection (‘conservative’). Following this, RapidHRV was able to produce estimates for all conditions in the uncleaned data and across most activities for cleaned data (conditions estimated: M = 87.5%, SD = 11.57%). Missing estimates from the latter were predominantly limited to high motion activities (1 missing estimate for a single subject in the driving condition).

Across low motion activities, RapidHRV PPG-based estimates converged with simultaneous (visually verified) ECG measures (Reiss et al., 2019). Agreement with low motion activities was excellent for BPM (ICC > 0.90) and agreement with RMSSD was moderate-to-good (0.57 < ICC < 0.82; Table 6; Fig. 10). Under high motion conditions, heart rate and heart rate variability estimates showed poor agreement (ICC < 0.32), except for heart rate within the cycling condition (ICC = 0.75).

Simulations

Within simulated data, RapidHRV was able to recover heart rate across most levels of noise (white gaussian noise filter >= 10 dB), even at relatively low sampling rates (>=20 Hz). RapidHRV’s recovery of heart rate variability was excellent at relatively low levels of signal to noise ratio (>=20 dB), though there was degradation of performance as sampling rate decreased. Additional simulations of cardiac responses to an anxiety induction demonstrated RapidHRV estimations fully captured effects at moderate levels of noise (>=30 dB) even at relatively low sampling rates (i.e., 20 Hz). RapidHRV was able to partially capture effects (~50% reduction in effect size) at very high levels of noise (>=10 dB when Hz >50). Simulations revealed RapidHRV cleaning was particularly beneficial at lower sampling rates and higher levels of noise, but was not necessary (or could be relaxed) when signal and sampling rates were high. Moreover, these simulations were able to clarify the validity of RapidHRV’s default window (10 s) for estimation of heart rate variability across a longer time period (i.e., 5 min).

Electrocardiography

In our electrocardiography analyses, we were able to recapitulate previously reported effects of age on heart rate and heart rate variability. In line with previous analyses (Iyengar et al., 1996), RapidHRV-estimates suggested older participants had lower heart rate variability than younger participants during movie-watching. Importantly, cleaning was vital to this detection.

Finger PPG

Using RapidHRV estimates, we noted effects of anxiety on heart rate and heart rate variability in a database of participants watching horror and documentary videos while undergoing finger infrared PPG recordings. Notably, the estimated effect size was analogous to that noted in threat-of-shock studies (Gold, Morey & McCarthy, 2015). Moreover, when contrasting subject-specific estimates, we found good-to-excellent agreement between RapidHRV and a previous analysis using a commercially available software. Effect sizes between conditions and convergence of estimates between softwares was significantly improved following RapidHRV cleaning.

Wrist PPG

In our wrist PPG validation, we noted RapidHRV was not able to produce estimates for the majority of conditions due to limitations of the default arguments. Despite applying scaling (0–100) at the level of a sliding window, PPG signals showed variable peaks amplitudes. Following alterations [rhv.analyze(preprocessed, outlier_detection_settings = “conservative”, amplitude_threshold = 30, window_overlap = 9)], RapidHRV-estimates during low motion activities (e.g., reading, eating lunch) demonstrated good-to-excellent agreement with a manually-verified analysis of ECG data. During high wrist-movement activities (e.g., table soccer), estimates were generally poor-to-moderate. We do note however that in one of the high motion conditions, cycling, RapidHRV-estimation of heart rate was good (though heart rate variability estimation was poor). This may reflect our proxy for motion-related artifacts (jerk magnitude), which may not always be as good an estimator of artifacts during activities which involve changes in acceleration but relatively low movements in the wrists (i.e., cycling, as hands are gripping the handle bars). Future work should seek to correlate RapidHRV quality assurance metrics with a wider range of motion estimation methods. Overall, we found adjustments to parameter arguments beneficial to PPG data, but we also note there are contextual factors and limitations (i.e., motion-related artifacts) influencing the feasibility of accurate estimation.

Overall user recommendations

For ECG data, users may find traditional, amplitude-based analyses will not work for subjects who demonstrate atypical signal morphologies (e.g., particularly prominent P and T waves). RapidHRV includes the use of k-means clustering to help discern these components of the ECG signal, though this is not enabled by default (ecg_prt_clustering = False). Additionally, the cleanliness of signal, namely the stability of peak prominences, means the data may be low in artifacts, and that minor deviations could be detected as outliers. As such, when dealing with already-clean data, users may find that outlier rejection can be omitted or relaxed (e.g., outlier method = ‘liberal’).

Results from the finger IR PPG data used in the present study did not suggest the need for alterations to default RapidHRV arguments, but did suggest that automated cleaning should be used.

Lastly, PPG data collected from naturalistic settings is typically low in signal to noise ratio, which can constrain peak detection. Consequently, lowering the minimum amplitude threshold for peak detection and decreasing window movement may help improve extraction. Furthermore, given the large amount of motion-related artifacts and the results from our simulation analyses, we recommend: (a) the use of relatively conservative cleaning (e.g., outlier method = ‘conservative’), and (b) inspection of motion across conditions as an indicator of estimation accuracy (Table 7).