Data Fusion for Improved Respiration Rate Estimation

We present an application of a modified Kalman-Filter (KF) framework for data fusion to the estimation of respiratory rate from multiple physiological sources which is robust to background noise. A novel index of the underlying signal quality of respiratory signals is presented and then used to modify the noise covariance matrix of the KF which discounts the effect of noisy data. The signal quality index, together with the KF innovation sequence, is also used to weight multiple independent estimates of the respiratory rate from independent KFs. The approach is evaluated on both a realistic artificial ECG model (with real additive noise), and on real data taken from 30 subjects with overnight polysomnograms, containing ECG, respiration and peripheral tonometry waveforms from which respiration rates were estimated. Results indicate that our automated voting system can out-perform any individual respiration rate estimation technique at all levels of noise and respiration rates exhibited in our data. We also demonstrate that even the addition of a noisier extra signal leads to an improved estimate using our framework. Moreover, our simulations demonstrate that different ECG respiration extraction techniques have different error profiles with respect to the respiration rate, and therefore a respiration rate-related modification of any fusion algorithm may be appropriate.


I. INTRODUCTION
ESTIMATION of respiratory rate from waveforms recorded from passively breathing subjects is notoriously difficult, due in part to the nonstationary nature of the signals and in part the frequent nonstationary noise [1]. Methods for recoding a time series of respiratory effort include impedance pneumography (differential changes in capacitance recorded at high frequencies), impedance plethysmograpy (stretch sensors on the chest wall) and flow thermography which measures the changes in temperature of air flow as it moves in and out of the mouth and/or nose over a thermistor. Of these methods, the impedance pneumogram (IP) is the most common method employed in hospitals [2,3], where a current is passed between two ECG electrodes and the differential change in capacitance due to air volume changes is measured. It is also possible to automatically record respiratory activity from accelerometers, laser-reflectivity, ultrasound, or by audio or video processing. However, such signals are not commonly recorded in most patient monitoring scenarios.
Another class of respiratory signal sources comes from measurement of indirect effects on cardiovascular physiology. Respiratory information is present in other commonly monitored physiological signals, such as the electrocardiogram (ECG) [4], photoplethysmogram (PPG) [5], arterial blood pressure (ABP) [6] and the peripheral arterial tonometry (PAT) waveforms.
Amplitude-based ECG derived respiration (EDR) algorithms have been reported to perform satisfactorily when only single-lead ECGs are available, as is usually the case in sleep apnea monitoring. When multilead ECGs are available, EDR algorithms based on either multilead QRS area or QRS-VCG loop alignment are preferable [7]. The reason is that due to thorax anisotropy and its inter-subject variability together with the inter-subject electrical axis variability, respiration influences ECG leads in different ways; the direction of the electrical axis, containing multilead information, is likely to better reflect the effect of respiration than wave amplitudes of a single lead.
Although the IP waveform is usually more representative of the air flow (but not volume), the IP is often unusable for accurate respiration rate evaluation around 37% to 61% of the time [8,9].
Literature on deriving a respiratory signal from other related signals is dense in the case of the electrocardiogram, but relatively sparse for other signals (such as the photoplethysmogram and blood pressure waveform [10]). The field of respiration rate estimation from respiratory signals (whether derived or not) is scantily covered in the public literature, particularly with respect to large data sets. Key works by O'Brien et al. [11], de Chazal et al. [12], Ishijima et al. [13], Park et al. [14] and Tarassenko et al. [1] highlight the importance of modeling noise and combining information from multiple ECG leads and sensor modalities to compensate for noisy measurements. To our knowledge, very little work has been published in the domain of respiration signal quality estimation.
The approach detailed in this paper is based on our earlier work related in the context of robust heart rate and blood pressure estimation [20][21]. Figure 1 summarizes the proposed robust respiration rate estimation technique using a signal fusion framework and signal quality indices. Fusion is the processes of combining signals from multiple instruments and sources in order to reduce measurement noise and improve overall signal quality. The field of data fusion in the context of physiological signals is described elsewhere [20]. When applied to respiration in particular, the key works of Mason [15][16] form the basis for our approach. However, our key innovation is the inclusion of a signal quality metric (as well as the past behavior of each signal) to control the KF noise covariance estimate and decide automatically how to weight each source of information. This is in contrast to standard industry approaches such as the work of Park et al. [14], who developed a system for deciding which single channel of ECG-derived respiration was the most informative. In this article we present a method for combining all the available information from every channel, even if it is noisy, to produce a superior estimate.
Given that artifacts in the ECG, IP and PPG or PAT do not always manifest simultaneously, it is likely that fusing respiratory rate data from each signal will result in an overall improved estimate of this parameter [15][16]. The difficulty lies in determining which signal to trust, for if we add parameter estimates from noisy signals to those from clean signals, we may in fact degrade any one single parameter estimate. We present a solution to this problem for accurate robust estimation of the respiration rate using signal quality indicies (SQI's) and a modified Kalman Filter (KF) fusion framework which uses the SQI's to adaptively update the KF noise covariance estimate. The SQI's are derived in real time and therefore no assumptions concerning the signal-to-noise ratio (SNR) are required.
The paper is organized as follows: section II-A describes the datasets used in this study. Section II-B introduces the methods for deriving respiratory waveforms from the ECG and PAT signals, section II-C provides an overview of a previously introduced ECG signal quality index and expounds on the applicability of a newly developed signal quality measure to assess the quality of the derived respiratory waveforms, section II-D describes the utilized respiratory rate extraction algorithm, followed by a discussion of the proposed KF-based fusion framework in section II-E. In sections III-A and III-B we present the results of applying the proposed respiratory rate estimation algorithm to simulated data and real recordings from 30 subjects during an overnight polysomnographic study.

A. Simulations and Data Collection
Two datasets were used for the analysis; one set of computer simulated ECG recordings with known respiratory signal modulation, and a set of real recordings from 30 subjects during an overnight polysomnographic study.

Figure 1
The proposed robust respiration rate estimation technique using a signal fusion framework and signal quality indices.

1) Simulated Data
Our simulated data were based on the synthetic ECG generation framework described in Clifford et al. [17], where we presented generalizations of our previously published artificial models for generating multi-channel ECG to provide simulations of cardiac rhythms. Using a three dimensional vector-cardiogram (VCG) formulation, we generated the normal cardiac dipole for a patient using a sum of Gaussian kernels, fitted to real VCG recordings. The RR interval time series were generated using our previously described model whereby time-and frequency-domain heart rate (HR) and heart rate variability (HRV) characteristics could be specified. Furthermore, following Astrom et al. [18] we incorporated a model of respiratory sinus arrhythmia (RSA) and RS-amplitude modulation to reflect influence of respiration on ECG. All the ECG signals were generated with a sampling frequency of 500 Hz and 16 bit amplitude resolution. Finally, realistic noise consisting of a combination of white noise (with SNRs of 10, 20, and 40 dB), baseline wander, muscle artifacts, and electrode motion (obtained from the MIT-BIH Noise Stress Test Database [19]), were separately generated and added to the simulated ECG. This simulation study has been used to compare breathing rate estimates from individual respiratory waveforms, as well, the combined breathing rate estimates using the proposed KF-based fusion framework.

2) Real Data; Overnight polysomnography
To evaluate our methods on real data, we chose a database of 30 subjects undergoing overnight polysomnogram analysis. Recorded signals included: one lead of ECG (V6), 4 channels of respiratory waveform data (chest and abdomen plethysmograph, nasal and oral thermistor, and nasal pressure), and one channel of PAT (Itamar Medical, Israel). The PAT signal is a pulsatile waveform recorded at the peripheral artery (fingertip), which is much like the ABP or PPG in morphological appearance, and reflects the rapid changes in blood pressure at the periphery from beat to beat. In particular, the PAT waveform, rather like the ABP or PPG waveform, exhibits amplitude fluctuations due to respiration.
All channels of data were recorded at 500Hz, 16 bits, with the exception of the PAT which was recorded at 100Hz. The length of each recording varied between 6 hours and 8 hours. Subjects were enrolled in the study for screening for sleep apnea, with an apnea-hypopnea index (AHI) ranging from 0 to 69.3 events/hour with a mean AHI of 14.0 events/hour. Since the respiration rate was not scored by humans, no gold standard measure of respiration rate was available. To provide a gold standard, we used the highest quality signal (the nasal and oral thermistor) and used the same AR spectral estimation technique for rate estimation, and the signal purity for noisy section rejection.

B. Deriving Respiratory Waveforms
ECG beat detection was performed using a combination of two open source QRS detectors; Hamilton & Tompkins' eplimited QRS detector [22] and Zong's wqrs algorithm [23]. Beat detection for the PAT waveform was achieved using Zong's wabp blood pressure onset detector [31]. We considered three categories of derived respiratory waveforms:

2) RSA-derived respiratory waveform
Respiratory Sinus Arrhythmia (EDR_RSA) is known to contain a strong respiratory component, amongst other components [4].

3) PAT-derived respiratory waveform
The pulse amplitude is known to be modulated by respiration. Therefore, if we pick the onset and peak of each pulse (using an open source algorithm wabp [31]) we may derive an oscillatory signal indicative of the respiration effort.
Since normal range of breathing rates in adult humans is between zero (apnea) and 60 breath/min (extreme hyperventilation) each respiratory waveform was resampled to 4Hz using a cubic spline method. Each respiratory waveform was segmented into 20s windows with 15s overlap from which the breathing rates were estimatedtherefore, after an initial delay of 20s the breathing rate estimates were updated every 5s. These estimated values of the breathing rate were the measurement inputs to the Kalman filter.

C. Signal Quality Metrics
We have described our approach to determining the quality of the ECG previously [20,21]. Briefly, we combine measures of abnormal statistics and power spectral density distribution with measures of QRS-detection miss-matches to provide an overall quality estimate (between 0 being poor and 1 being excellent) for any given segment of ECG. The estimation of PPG signal quality is described in Gil et al. [24,25]. We have extended these ideas to PAT and respiratory signal quality. Fidelity of respiration rate extraction is directly related to the periodicity of respiration waveform. A regular breathing pattern produces a highly sinusoidal reparation waveform with the dominant frequency at the frequency of breathing rate. A power spectral based respiration rate extraction can be accomplished by identifying the most prominent peaks in any of the respiration waveforms (PAT derived respiration, or ECG derived respiration waveform). Each peak is characterized by its frequency, amplitude, or its relative coherence with respect to similar spectral peaks in other respiratory waveforms. One measure of characterizing spectral characteristics of a signal is through Hjorth descriptors [25]. The n th -order spectral moment n ω is defined as: is the power spectrum of the signal as a function of angular frequency: , with f being in units of cycles/second. A particularly useful descriptor in the context of estimating the dominant frequency and assessing the quality of a signal with periodic components (such as respiratory waveform) is the so called spectral purity waveform and is defined as [26]: Here the term "purity" refers to the presence of a single signal frequency, as we would expect in an ideal respiratory waveform. In the case of a periodic signal with a single dominant frequency, s Γ takes the value of one and approaches zero for non-sinusoidal noisy signals. One of the attractive features of Hjorth descriptors is the feasibility of their calculation in time domain with low computational cost.

D. Deriving respiratory rates
The breathing rate extraction method was based on the work of Mason et al. [16], who utilized autoregressive (AR) modeling, a parametric spectral analysis technique. One advantage of AR modeling of spectral analysis over the traditional Fourier transform based methods is its superior performance when the number of available data points is small (<100 points). The steps involved in extracting breathing rates from the respiratory waveforms are: • Fit an all-pole model to each 20s segment (use Akiake Criteria [30] to decide on the right model order). The window was advanced by 5s and the above process repeated so that a respiration rate was available every 5s. Figure 2 is an example of the respiration rate estimation applied to a 300s long record.

Figure 2
Respiratory rate estimation; from top to bottom: nasal thermistor (taken as the gold-standard), EDR_RS (amplitude modulation-based method), and EDR_RSA (respiratory sinus arrhythmia-based method). The right panels are the estimated respiratory rates. The left panels are the respiratory waveforms corresponding to the shaded area on the right (note: the original recording was 300s long. Only the respiratory waveforms corresponding to the shaded area on the right is presented on the left to enhance clarity). Each estimate (closed-circle marks on the right) correspond to a 20s long window with 15s overlaps between consecutive windows (thus estimates are updated every 5s). For this example, within the shaded area the RSA method provided a better respiratory estimate than the RS method.

E. Kalman Filter Framework and Data Fusion
The KF is an optimal state estimation method for a stochastic signal [27] that estimates the state of a discretetime controlled process, x , with measurement data z , where x and are governed by the linear stochastic difference equations: Previously we employed the KF to estimate the systolic, mean and diastolic blood pressure derived from ABP and HR from the ECG [20][21]. In order to more heavily weight estimate derived from cleaner data, the SQI is used to adjust the measurement noise covariance, R. When the SQI is low, k z should be trusted less, and hence we force R to be s large. Thi is achieved by modifying at the m as follows: where k SQI is the signal quality of the k th segment of data and may y various measures of the underlying signal quality, such as the purity index ned in Eq. (5), that is, This nonlinear weighting function therefore tends to unity as the value of k SQI tends to unity, (at which point the asurement noise covariance is no longer affected by the SQI) forcing the KF to trust the current measurement k z given the baseline measurement noise covariance matrix k R . At low values of k SQI , k R tends to infinity (but in practice is limited to a large value to deal with issues of convergence) and therefore forces the KF to trust current measurements less. This is the key factor in the modified KF framework; we allow the KF to make a varying estimate of the noise covariance u me sing independent sig nd PAT features and SQI were used by the KF to btain the optimal breathing rate estimation on a 5s-by-5s ly equal can then ew measurement is a new estimate of the respiratory nal quality estimates based upon domain knowledge of the underlying ECG, and not the KF itself.
The SQI of each heart beat was calculated ±5s around each beat. Second-by-second ECG and PAT SQI were acquired by calculating the median values of these beats within a moving 20s window with 50% overlap. Then, the ECG a o basis.

1) Kalman Filter Initialization and Operation
Following Tarassenko and Townsend [1,15], we pick the simplest form of the KF, and set the state to be a scalar. We assumed that the breathing rate at each moment is approximate to the breathing rate at the next moment. After neglecting the control input u , Eq. (8) then reduces to . In order to initialize the KF, one must estimate Q, the state noise covariance matrix, and R, the measurement noise covariance. R was similarly initialized to unity, noting that it is immediately modified by the SQI to reflect our trust in the data. Q was empirically adjusted to have an initial value of Q=5 (±5 breath/minute). Values of Q<5 lead to the KF trusting the initial state estimate too little and not adapting to the new initial observations. Values of Q>5 lead to the KF trusting the new observations too much, and simply following the new values too closely. The filter run online with only a few iterations for convergence. The Kalman residual is then given as here a n be rate from one of the sensors as described in section II-B and Townsend et al. [1] such that the final estimate at the k th time-step is given by

2) Merging of Multiple Kalman Filter Estimates
In order to calculate a single estimate of the respiratory rate, estimates from individual Kalman filters must be fused in a manner that takes into account the uncertainty associated with each estimate. In general, it is possible to fuse any number of independent Kalman filtered estimates, s k x , , using the technique of Mason [15]    PAT-derived data.
iii. Fusion2 algorithm used the signal purity index (for all derived respiration signals).
In the case of the real data, we report the performance of all seven respiration algorithms, while for the simulation studies only the ECG based algorithms are considered. Note that, fusion0 a et al. [15] which weigh imates proportional to the inverse of the square innovations.

A. Simulation Results
In order to evaluate the respiration rate estimation algorithm, we used the simulated ECG data with varying respiration rates and under different heart rates (60-100beats/min), and signal-to-noise ratio scenarios (10, 20, and 40dB). The heart rate dependence was negligible; therefore the results that follow are the average performances over all heart rates. Note that, all values reported in this work are in root mean square (RMS) breaths per minute (BPM).
Results of these simulations (see Figure 3) indicate that the RS method is the best estimator for breathing rates in the range of 16-24 BPM, the RSA method is best for breathing rat s in the range of 8-12 e best for rates in the range of 16-24 BPM. At the lowest SNR, the KF fusion algorithm provides a good estimate for the rates of 8-24 BPM. Table 1 summarizes the overall performance of each algorithm individually, as well as, the fusion results using the purity SQI. A comparison across a wide range of SNRs indicated that the fusion re   (indicated by a black arrow). Figure 4 is an example of the KF-based fusion algorithm using four different respiration rate signals. The bottom plot is the fusion result using the purity signal quality index, while the proceeding plot (one to the bottom panel) shows the results of using no signal quality (thus, using a fixed measurement noise covariance matrix). In the 140-280s region, where the quality of the estimation is poor, the signal quality/KF-based fusion method shows clear improvements in the estimated respiratory rate. Table 2 presents results of the average (RMS) error for all 30 patients in the overnight polysomnographic database. Four single channel algorithms and three fusion algorithms were evaluated. The four single channel algorithms were RS amplitude EDR_RS, QRS area (EDR_G) and RSA-based methods (EDR_RS) of deriving respiratory waveforms from the ECG, and the pulse amplitude modulation method of deriving respiratory waveform from the PAT signal (PDR). Respiratory rates and purity-based signal quality indices were extracted from these waveforms using the AR modeling method of section II-D and signal quality metrics of section II-C.

B. Overnight Polysomnography Results
Note that the three fusion algorithms consistently perform better than any one single algorithm. In particular, when the SQI is used, a greater improvement (lower RMS error) is seen. Interestingly, results from our real data are better than our simulations, even at high SNR, indicating that even though the PAT derived respiration is particularly poor, the inclusion of this signal still improves the overall performance. This is because the PAT signal was occasionally good quality at times when the ECG was noisy.

IV. DISCUSSION AND CONCLUSIONS
In this work several respiratory rate estimation algorithms have been presented which estimate a respiratory rate from the respiratory waveforms derived from ECG or PAT. They have been divided into three categories: 1) EDR algorithms based on beat morphology, namely, those based on ECG wave amplitude or QRS area (EDR_G).

2) EDR algorithms based on HR information; i.e., respiratory sinus arrhythmia (EDR_RSA) 3) Algorithms based on pulse or PAT amplitude variations (PDR).
The choice of a particular EDR algorithm depends on the application. In general, EDR algorithms based on beat morphology are more accurate than those based on HR information, particularly at high respiration rates (>16 BPM) since the modulation of ECG by respiration is sometimes too small or embedded in other parasympathetic interactions. At all other rates the KF-based fusion approach is superior.
Electrocardiogram-derived respiration algorithms based on both beat morphology and HR may be appropriate when only a single-lead ECG is available and the respiration effect on that lead is not pronounced. Although not detailed here, the power spectra of the EDR signals based on morphology and HR can also be cross-correlated to reduce spurious peaks and enhance the respiratory frequency [28]. However, the likelihood of having an EDR signal with pronounced respiration modulation is better when the signal is derived from multilead ECGs; cross-correlation with the HR power spectrum may in those situations worsen the results due to poor respiratory HR modulation. This is particularly true during active wakefulness [28].
The median value of the fusion based estimation only diminishes by 1.7 BPM from an SNR of 40 dB down to an SNR of 10 dB. The reason for such robust performance is not only due to the effectiveness of the fusion technique but also due to the good performance of the EDR_RS and EDR_RSA in the presence of noise. This is because they only rely on accurate detection of fiducial points on the ECG and therefore are not affected by the presence of noise in the other segments of the ECG (note that the R-peak location has the largest SNR on the ECG trace and the S-point detection can be fairly reliable given the location of the Rpeak). It should also be noted that the QRS area method's effectiveness in allowing respiration rate estimation is sensitive to the size of the integration window, which was not optimized in this study. Further improvements are therefore possible. Specifically, the QRS integration window could be shortened at higher heart rates, (which is often correlated with higher respiration rates) to compensate for shortening of the QRS interval at these higher heart rates and elevated sympathetic tone.
Our results on real subjects' data recorded during sleep indicate that although the PAT is often noiser that the ECG and provides a poorer estimation of respiratory rate, if included in our KF-based fusion framework, it still adds value to the estimation process, particularly when an SQI is used. Furthermore, our respiration signal quality index, Γ, is superior to the use of the ECG-and PAT-specific SQI metric in providing a trust metric to automatically discount noisy data sources. In other words, testing the derived signal for information content provides a better trust index of the derived respiration rate than does an estimate of the quality of the underlying data from which the respiration signal is derived. This may be because a good quality waveform does not always carry respiration information, and may therefore degrade the KF-based fusion algorithm if a post-processing quality metric is not used. It may of course be productive to look at using both signal quality metrics.
Although our gold standard respiration rate is derived from nasal thermistor, it is unlikely that the nasal thermistor always provides a clean and representative signal, being susceptible to movement artifact at least. If the nasal thermistor provided a perfect evaluation of the respiration rate, we would expect our results to be even better, with lower error rates for all algorithms. However, we do not expect any of our individual algorithms to improve substantially in such a scenario, because the nasal thermistor is independent of the other measurement methods. The improved results of our KF-based fusion approach would not be invalidated in this case, and perhaps improved.
We can use the results of our investigation to inform a more intelligent combination of signals. By studying tables 1 and 2 we see that the RSA and QRS area based algorithms perform relatively weakly at higher respiration rates (≥16 BPM), which suggests that we should weight the RS amplitude more strongly when the latter algorithm indicates high respiration rates. This can be incorporated by making the trust factor, λ k , respiration rate-dependent. The exact dependence could be calibrated for a particular population or recording environment if sufficient data were available.
There are still certain topics in the field of derivedrespiration which deserve further study. One is the robustness of the derived-respiration algorithms in different physiological conditions; robustness to long-term beat morphology variations due to, for example, ischemia. The study of multimodal respiratory patterns should be considered when estimating the respiratory frequency from physiological signals by techniques like, for example, spectral coherence, particularly with wavelets to deal with nonstationarities. Weighting of the individual algorithms based upon respiration rate (through adjustment of λ k in Eq. (12)) also require further exploration given our simulation results. Also, λ k may be adjusted to weight certain sensor measurements over others which are less trustworthy by nature. For instance, our results on the real data indicate that in general the PAT-derived respiration may not be as reliable as the ECG-derived respiration for our population. This may not be true (or may be reversed) for other populations, such as neonates for example. Therefore, demographics may inform the values of λ k .
Although our approach automatically rejects noisy sections of data, it may be further improved by using statistical tests which can reject spurious frequencies that could have appeared just by mere chance. This may be accomplished in a nonparametric and nonstationary manner using a surrogate analysis approach [29], for example.
The added computational burden of our fusion step (equation 11) is negligible, being essentially a weighted average it involves only a very few divides and adds. Generally one would expect the derived respiratory waveforms and respiratory rates to be already calculated. In fact, signal qualities are generally computed in most monitors, although are often not available to the standard user. Even in the case that a signal quality is not available, the purity-based signal quality method utilized in this work is calculated in time domain using a finite-differencing approach (see [26]) and is therefore computationally very efficient. The majority of the computation is in the standard ECG analysis, which must be performed regardless.
The most important point to emphasize about our approach is that it does not require any a priori knowledge of the (changing SNR) of any input signal. Therefore our approach is robust in a wide sense, weighting the best estimators at any given epoch to provide a consistently superior estimate to any single given technique.
It should also be noted that the method presented in this paper is quite general and can be extended to any set of sources that provide a respiratory-related oscillatory data, such as the photoplethysmogram (processed in an almost identical manner as PAT), an airway thermistor, a flow meter, an impedance pneumogram, an accelerometer attached to a patient's chest, or an infra-red camera pointing at the patient's mouth and nose. Moreover, this robust KFbased fusion framework is extensible to any set of independent (or nearly independent) observations, providing a suitable signal quality parameter can be defined. Table 1. Overall performance summary of each algorithm at SNRs of 10, 20, and 40 individually and the fusion results using the purity signal quality indices. All figures are RMS error in BPM. Note that, both EDR_RS and the QRS-Area (EDR_G) methods perform poorly at low breathing rates, while EDR_RSA performance degrades at higher respiratory rates.