An assessment of algorithms to estimate respiratory rate from the electrocardiogram and photoplethysmogram

The extraction stage derives a respiratory signal (a time series dominated by respiratory modulation) from the original signal (ECG or PPG). The initial step, common to all techniques, was elimination of very low frequencies. This was achieved using a high-pass filter with  −3 dB cutoff frequency of 4 breaths per minute (bpm). Intermediate steps were determined by whether a filter- or feature-based extraction technique was being used, as described below. The final step, again common to all techniques, was elimination of non-respiratory frequencies from the extracted signal using a band-pass filter (with  −3 dB cutoff frequencies of 4 and 60 bpm).

Filter-based extraction consisted of filtering the signal using one of the four techniques in table 1, ${{X}_{\text{A}1,\ldots,4}}$. Feature-based extraction consisted of extracting a time series of beat-by-beat feature measurements through several steps as follows. Very high frequencies were eliminated using low-pass filters with  −3 dB cutoffs of 100 and 35 Hz for the ECG and PPG respectively. An additional 50 Hz notch filter was used to eliminate mains interference in the ECG. Beat detection was performed on the ECG using a QRS detector based upon the algorithm of Pan and Tompkins (1985) and Hamilton and Tompkins (1986), and on the PPG using the incremental-merge segmentation (IMS) algorithm (Karlen et al 2012). R-waves and pulse peaks were detected as the maxima at or between detected beats. QRS troughs were detected as the minima within the 0.10 s prior to R-waves (Ruangsuwana et al 2010), and pulse troughs as the minima between pulse peaks (Johansson 2003). One of the beat-by-beat features, ${{X}_{\text{B}1,\ldots,10}}$, was obtained as described in table 1. Features derived from ectopic beats were eliminated using the algorithm described in Mateo and Laguna (2003). The irregularly sampled signal of features was resampled at 5 Hz using linear interpolation (Karlen et al 2013). Finally, very low frequencies were eliminated using a high-pass filter with  −3 dB cutoff of 4 bpm.

Estimation of RR from a respiratory signal was performed using each of the techniques ${{E}_{\text{F}1,\ldots,7}}$ and ${{E}_{\text{T}1,\ldots,5}}$ in table 2. When using spectral analysis (${{E}_{\text{F}1,2,3,7}}$) the RR was identified as corresponding to the frequency of the spectral peak with greatest magnitude between 4 and 60 bpm. When using breath detection methods (${{E}_{\text{T}1,\ldots,5}}$) the estimated RR was calculated as the mean breath duration.

The fusion stage has received particular attention recently due to the improvements in algorithm performance observed with its use (Karlen et al 2013). Four modulation fusion methods, ${{F}_{\text{M}1,\ldots,4}}$ in table 3, were optionally used to fuse simultaneous RR estimates corresponding to each modulation, extracted using ${{X}_{\text{B}1,2,3}}$. Temporal fusion, F_T1, was optionally used to smooth successive RR estimates obtained from the same individual. It may be used with or without preceding modulation fusion.

We verified that techniques had been implemented successfully by testing on simulated data. Techniques were combined into algorithms as described above. All algorithms were assessed on both simulated ECG and PPG, except where the respiratory extraction technique was signal-specific: X_B7 and X_B8 only operate on ECG, and X_B10 only operates on PPG.

Simulated data, such as those shown in figure 1, were generated as follows. Single ECG and PPG beats were replicated to generate simulated signals consisting of trains of ECG and PPG beats, each of 210 s duration sampled at 500 Hz. Each of these two simulated signals was modulated separately to mimic each respiratory modulation (BW, AM, and FM), leading to six modulated signals. The modulations were repeated for a range of physiologically plausible heart rates (HRs) and RRs. A first set of 35 simulated signals was generated for each modulation by holding RR constant at 18 bpm whilst the HR was varied from 30 to 200 bpm at intervals of 5 bpm. A second set of 29 simulated signals was simulated by holding the HR constant at 80 bpm whilst RR was varied from 4 to 60 bpm at intervals of 2 bpm. Hence, a total of 192 simulated ECG and PPG signals were generated with which to verify algorithm implementations.

A technique was considered to be successfully implemented if over half of the algorithms containing that technique were acceptably accurate. 'Acceptably accurate' was defined as an absolute error of  ⩽1 bpm for at least 50% of the simulated signals for any of the three modulations, on either the ECG or PPG. The threshold of 50% was chosen to account for algorithms which perform well for subsets of the trialled HR and RR combinations, but fail under extreme physiological conditions. The correlation between performance on the simulated data and real data was unknown at the time of verification. However, it seemed reasonable to assume that algorithms which performed extremely poorly on simulated data would not perform well on real data.

A total of 370 algorithms were implemented. 314 algorithms (85%) met our acceptability criteria. Algorithms containing either E_F5 or X_B10 were excluded. E_F5 performed poorly at RRs outside of 12–20 bpm, which may be due to its bias towards identifying lower frequencies as the RR. X_B10 variably detected the end of the PPG pulse as either the time of the minimum immediately before the diastolic peak, or the time of the diastolic peak, causing inaccuracies.

Healthy adults aged between 18 and 40 years of age participated as part of the VORTAL study (National Clinical Trial 01472133). Ethical approval was obtained from the London Westminster Research Ethics Committee (11/LO/1667). Exclusion criteria were the presence of co-morbidities or medications that might significantly affect the functioning of the cardiac, respiratory and autonomic nervous systems.

Lead II ECG, PPG recorded from a finger probe, and oral-nasal pressure signals were acquired using a 1902 amplifier, a Power 1401 analogue-to-digital converter and Spike2 v.7.09 acquisition software (all Cambridge Electronic Design, Cambridge, UK). PPG was transduced using an MLT1020FC infrared reflection plethysmograph (AD Instruments, CO Springs, New Zealand). Oral-nasal pressure was transduced using an Ultima Dual Airflow differential pressure transducer (Braebon Medical Corporation, Kantata, ON, Canada) connected to a P1300 Pro-Flow oronasal cannula (Philips Respironics, Murrysville, PA, USA). Signals were sampled at 500 Hz.

Thoracic IP RRs were acquired at 1 Hz using an IntelliVue MP30 monitor (Philips Medical Systems, Boeblingen, Germany) and ixTrend acquisition software (v.2.0.0 Express, Ixellence GmbH, Wildau, Germany).

We prepared the subject's skin prior to application of ECG electrodes in the Mason-Likar configuration. We shaved excess chest hair and cleansed the skin with an alcohol wipe before removing dead skin using 10–15 light strokes with abrasive skin preparation tape. The PPG probe was placed on the right hand.

Two 10 min recordings were acquired from each subject laid supine in a room maintained between 22 and 25 °C. After the first recording subjects were asked to run on a treadmill until their HR reached 85% of their age-predicted maximum, as described in Gibbons et al (1997). The second recording was acquired immediately following exercise. The post-exercise recordings contained a wider range of HRs and RRs than the pre-exercise recordings.

Respiratory signals derived from the ECG or PPG were segmented into adjacent windows of 32 s. Signal quality was assessed using the algorithm described in Orphanidou et al (2015). Firstly, the algorithm detects beats in the signal. Secondly, beat-to-beat intervals and HRs are measured. If either are physiologically implausible, then the window is deemed to be of low quality. For each high quality window a template beat is constructed. The correlation between individual beats and the template is assessed. If the correlation is below an empirically determined threshold, then the window is deemed to be of low quality, as shown in figure 3. Windows labelled as low quality were excluded from the analysis.

Zoom In Zoom Out Reset image size

For each window RR was estimated using all possible combinations of algorithms and input signal. IP RRs were also calculated, based on the readings provided by the monitor. They were estimated as the mean of all IP RR values during that window.

Reference RRs were calculated from the oral-nasal pressure signal using a custom breath-detection algorithm. The signal was band-pass filtered with  −3 dB cutoff frequencies of 4 and 60 bpm. It was segmented into the same 32s windows as used for the ECG and PPG. The signal was normalised within each window to have a mean of 0 and standard deviation of 1. Inhalations were identified as positive-gradient crossings of a specified threshold, k. A threshold of k  =  0.42 was determined by minimising the difference between the value of 2SD (explained below) derived from the comparison of the reference RRs to expert breath annotations on the first 10 subjects. The performance of the reference RR calculator was then assessed against a second expert's breath annotations on the subsequent 10 subjects. It had a bias of 0.0 bpm and 2SD of 1.3 bpm. Windows in which the pressure signal had a low signal-to-noise ratio were excluded from the analysis. The threshold for exclusion was chosen to eliminate windows in which breaths could not be identified visually.

Agreement between each algorithm and the reference was assessed using the limits of agreement (LOA) technique. The bias, 95% LOAs, within subject variability (WSV) and between subject variability (BSV) were estimated using a random effects model to account for repeated measures in the data. Following Carstensen et al (2008) the random effects model used was

$$\begin{eqnarray}\begin{array}{*{20}{l}}&y_{mir} = \alpha_m + \mu_i + a_{ir} + c_{mi} + e_{mir} \quad \\ & \quad a_{ir} \sim \mathcal{N}(0,\omega^2) \quad c_{mi} \sim \mathcal{N}(0,\tau^2) \quad e_{mir} \sim \mathcal{N}(0,\sigma_m^2) \quad ,\end{array}\end{eqnarray} \tag{ 1 }$$

where y is the measurement, m is the method, i is the subject and r is the replicate. ${{\alpha}_{m}}$ is a fixed effect representing a method's constant underlying bias, ${{\mu}_{i}}$ is the mean bias for a particular subject, a_ir is a random effect accounting for the similarity of measurements within a subject, c_mi is a random effect related to the BSV (${{\tau}^{2}}$), and e_mir is the residual error related to the WSV of the method ($\sigma _{m}^{2}$). The bias was calculated for each algorithm as ${{\alpha}_{\text{reference}}}-{{\alpha}_{\text{algorithm}}}$, and the 95% LOAs as mean bias  ±2 standard deviations (2SD) where

$$\begin{eqnarray}&&2\text{SD}=2\sqrt{2{{{\hat{\tau}}}^{2}}+\hat{\sigma}_{\text{ref}}^{2}+\hat{\sigma}_{\text{alg}}^{2}}.\end{eqnarray} \tag{ 2 }$$

A second measure of agreement, the coverage probability (CP$_{\delta}$) is also reported. The CP$_{\delta}$ is the probability of measurement error falling within pre-defined bounds, δ. The non-parametric form of CP, expressed as a percentage, was calculated using the empirical cumulative distribution of the absolute error (Barnhart et al 2007) with δ set at 2 bpm.

The effect of the input signal on 2SD and bias was assessed using the Wilcoxon signed-rank test for those algorithms which could operate on both ECG and PPG.

Methods were ranked by sorting first by 2SD and then by absolute bias. Values of 2SD and absolute bias were rounded to one decimal place for sorting. 2SD, a statistic related to precision, was prioritised over bias when ranking since one is commonly more concerned with tracking trends in RR than absolute values when using wearable sensors for continuous monitoring. Furthermore, bias can be corrected by calibration if an instrument is precise.

Data were acquired from 45 subjects. Five subjects were excluded since their recordings were incomplete, and one due to the development of an arrhythmia following exercise. 39 subjects' data were analysed. The median (lower, upper quartiles) age of analysed subjects was 29 (26, 32) years. Their median BMI was 23 (21, 26) kgm⁻². 21 (54%) were female. Data from each subject contained a median of 40 (37, 41) windows. The oral-nasal pressure signal was of high quality in 40 (37, 40) windows. Of these, ECG and PPG were of high quality in 37 (34, 40) and 36 (34, 39) windows respectively. The ranges of RR and HR in the dataset were 5–32 bpm and 41–111 beats per minute respectively. The distributions of RR and HR are shown in figure 4.

Zoom In Zoom Out Reset image size

Of the 314 algorithms assessed, 270 could operate on both ECG and PPG and 44 were specific to the ECG. Therefore 584 algorithm-signal combinations were tested. The random effects model did not converge for 34 of these, all of which used E_F4. Inspection revealed minimal correlation between the reference and estimated RRs in these instances. These combinations were excluded from further analysis.

The values of 2SD of the remaining algorithm-signal combinations ranged from 4.7 bpm at best to 50.1 bpm at worst, as shown in figure 5(a). The top ranked method was ${{X}_{\text{B}1,2,3}}{{E}_{\text{T}4}}{{F}_{\text{M}1}}\left(\text{ECG}\right)$, whose measurements had a bias of 0.0 bpm, 95% LOA of  −4.7 to 4.7 bpm and CP₂ of 80.5%. IP was ranked 5th with a bias of  −0.2 bpm, 95% LOA of  −5.6 to 5.2 bpm and CP₂ of 76.0%. Performance metrics for the top 10 algorithms using ECG and PPG are shown in table 4. Metrics for all algorithm-signal combinations are reported in the supplementary material (stacks.iop.org/PM/37/610/mmedia).

Zoom In Zoom Out Reset image size

Table 4. Performances of the ten highest ranked algorithms for the ECG and PPG, and of IP-Derived RR: Ranked by 2SD followed by absolute bias.

Signal	Algorithm	Overall rank	2SD (bpm)	Bias (bpm)	95% LOA (bpm)	Proportion of windows with RR estimate (%)	CP2 (%)
IP	Clinical monitor	5	5.4	−0.2	−5.6 to 5.2	100.0	76.0
	${{X}_{\text{B}1,2,3}}$ET4FM1	1	4.7	0.0	−4.7 to 4.7	73.8	80.5
	${{X}_{\text{B}1,2,3}}$ ${{E}_{\text{T}2}}$ FM1	2	5.2	1.4	−3.8 to 6.4	72.3	72.6
	${{X}_{\text{B}1,2,3}}$${{E}_{\text{T}5}}$ ${{F}_{\text{M}1}}$	3	5.2	2.0	−3.3 to 7.2	75.4	69.1
ECG	${{X}_{\text{B}1,2,3}}$${{E}_{\text{T}3}}$ ${{F}_{\text{M}1}}$	4	5.3	1.4	−3.8 to 6.7	72.5	73.0
	XB2${{E}_{\text{T}2}}$	6	5.6	−0.2	−5.8 to 5.4	100.0	75.2
	XB2${{E}_{\text{T}3}}$	7	5.7	−0.2	−5.9 to 5.4	100.0	74.3
	XB2${{E}_{\text{T}2}}$ ${{F}_{\text{T}1}}$	8	5.7	−0.2	−6 to 5.5	100.0	69.3
	XB2${{E}_{\text{T}5}}$	9	5.7	0.5	−5.2 to 6.3	100.0	74.9
	XB2${{E}_{\text{T}3}}$ ${{F}_{\text{T}1}}$	10	5.8	−0.2	−6.0 to 5.6	100.0	69.8
	${{X}_{\text{B}1,2,3}}$${{E}_{\text{T}4}}$${{F}_{\text{M}1}}$${{F}_{\text{T}1}}$	11	5.9	0.0	−5.9 to 6.0	100.0	66.6
	${{X}_{\text{B}1,2,3}}$${{E}_{\text{T}4}}$ ${{F}_{\text{M}1}}$	15	6.2	1.0	−5.1 to 7.2	54.2	71.5
	${{X}_{\text{B}1,2,3}}$${{E}_{\text{T}1}}$ ${{F}_{\text{M}1}}$	17	6.5	−1.0	−7.5 to 5.5	62.1	62.1
	${{X}_{\text{B}1,2,3}}$${{E}_{\text{T}1}}$ FM1${{F}_{\text{T}1}}$	35	7.0	−1.3	−8.3 to 5.7	100.0	54.2
PPG	XB2${{E}_{\text{T}5}}$ ${{F}_{\text{T}1}}$	46	7.5	3.0	−4.5 to 10.5	100.0	44.3
	XB5${{E}_{\text{T}1}}$ ${{F}_{\text{T}1}}$	48	7.6	0.7	−6.9 to 8.3	100.0	57.0
	XB2${{E}_{\text{T}2}}$ ${{F}_{\text{T}1}}$	53	7.6	2.7	−4.9 to 10.3	100.0	47.2
	${{X}_{\text{B}1,2,3}}$ ${{E}_{\text{T}4}}$ ${{F}_{\text{M}1}}$ ${{F}_{\text{T}1}}$	54	7.8	1.1	−6.8 to 8.9	97.3	57.2
	${{X}_{\text{B}1,2,3}}$ ${{E}_{\text{T}5}}$ ${{F}_{\text{M}1}}$	55	7.8	3.8	−4.0 to 11.5	70.9	49.8
	XB2${{E}_{\text{T}4}}$ ${{F}_{\text{T}1}}$	56	7.9	0.3	−7.7 to 8.2	100.0	60.5
	${{X}_{\text{B}1,2,3}}$ ${{E}_{\text{T}2}}$ ${{F}_{\text{M}1}}$ ${{F}_{\text{T}1}}$	58	7.9	3.7	−4.2 to 11.5	100.0	60.5

Note: Definitions of the techniques are provided in tables 1–3.

Four algorithm-signal combinations were more precise than IP. These combinations all used a time-domain RR estimation technique and a modulation fusion technique, and operated on the ECG. They provided RR estimates for between 72.3 and 75.4% of windows. In contrast, IP and the next best-performing combinations provided estimates for all windows.

There was a sharp increase in values of 2SD for algorithm-signal combinations ranked in the lowest quintile, as shown in figure 5(a). Algorithms within this quintile contained the X_A1, X_A4 or E_F4 components suggesting that these components were responsible for poorer performance (figure 5(b)). Figure 5(c) shows that more of the better-performing combinations used a time-domain technique, and more of the worse-performing combinations used a frequency-domain technique.

Of the 283 algorithms which could operate on both ECG and PPG, 30 were excluded because the random effects model did not converge when operating on at least one of the signals. Of the remaining 253 algorithms, 161 (63.6%) were more precise when using ECG and 92 (36.4%) were more precise when using PPG. There was a statistically significant median decrease in the 2SD of 0.8 bpm when algorithms operated on ECG (11.6 bpm) compared to when they operated on PPG (12.4 bpm), z  =  −4.024, p  =  0.001.

${{X}_{\text{B}1,2,3}}$ ${{E}_{\text{T}4}}$ F_M1 gave the most precise estimates for both ECG and PPG-based methods. When operating on PPG it ranked 15th overall. The algorithms which were ranked in the top 10 when using ECG and those in the top 10 when using PPG shared common features. All used feature-based extraction techniques (${{X}_{\text{B}1,\ldots,5}}$) and time-domain breath-detection estimation techniques (${{E}_{\text{T}1,\ldots,5}}$). All except three used either smart fusion (F_M1, a type of modulation fusion) or temporal fusion (F_T1).

We cannot claim that all our implementations match those described previously. In some cases this is because algorithms have not been described in sufficient detail to be replicated exactly (Addison et al 2015). In other cases we have omitted algorithm-specific tuning steps which did not readily fit into our framework. However, the implementations of each technique were verified using simulated data, ensuring that poor performance was not due to incorrect implementation. By providing the source code of our algorithms we have provided the opportunity for others, including the authors of the original implementations, to verify our work and suggest improvements.

Although we have attempted to be comprehensive, we cannot claim to have implemented all techniques reported in the literature. This was not feasible due to the number of techniques previously proposed. We welcome contributions of new techniques for inclusion in future versions of the toolkit.

Performance rankings are approximate. In common with most statistical measures of agreement, calculation of the LOA depends on the distribution of measurement differences being normal and constant across the range of RRs. The distribution of RR measurement differences varies between algorithms, with some algorithms violating these assumptions to a greater degree than others. Therefore the LOA may be more accurately estimated for some than others. However, the estimation error is unlikely to be large enough to change the conclusions of our study.

While the 2SD provides a measure of the trustworthiness of a single measurement in isolation, the CP₂ provides a descriptor of the algorithms' typical performance. Algorithms with a high CP₂ but relatively large value of 2SD may merit further investigation. Such algorithms performed very well for a large proportion of the data and extremely poorly for the remainder. Modified versions of these algorithms might perform very well if it is possible to mitigate against the factors causing occasional poor performance.

The results cannot be directly extrapolated to other monitoring scenarios, such as during exercise, nor other populations, such as elderly or unwell subjects. We designed the study to demonstrate the best possible results which can be obtained from RR algorithms using ECG and PPG signals. To do so we recorded data from young healthy volunteers at rest using high-fidelity recording equipment and pre-processed signals to exclude artifact. In a clinical setting, for instance, the patient population and signal acquisition equipment are different. These differences may significantly affect algorithm performance.

Future studies should examine the effects of the subject population, artifact, and recording equipment on algorithm performance. This will determine the generalisability of the conclusions of this study. It will also be valuable to assess the utility of RR algorithms in a wider range of monitoring scenarios, including in the clinical setting.

This study indicates promising directions for future algorithm development. The smart fusion technique (F_M1) improved algorithm performance in this study. This technique weights RR estimates derived from each respiratory modulation equally. However, age and comorbidities may influence the strength of modulations. For instance, respiratory sinus arrhythmia (which causes FM) and chest wall excursion (which is linked to BW and AM) both diminish with age (Moll and Wright 1972, O'Brien et al 1986). Weighting the fusion according to the strength of each modulation's manifestation in a particular subject's data may improve performance.

In this study we have not investigated the computational requirements of algorithms. This may be important when deciding whether algorithms are suitable for use in real-time, power constrained settings, such as in wearable sensors in the wellbeing sector. In these settings the choice of algorithm may be strongly determined by computational capacity as well as precision.