Extracting Instantaneous Respiratory Rate From Multiple Photoplethysmogram Respiratory-Induced Variations

In this study, we proposed a novel method for extracting the instantaneous respiratory rate (IRR) from the pulse oximeter photoplethysmogram (PPG). The method was performed in three main steps: (1) a time-frequency transform called synchrosqueezing transform (SST) was used to extract the respiratory-induced intensity, amplitude and frequency variation signals from PPG, (2) the second SST was applied to each extracted respiratory-induced variation signal to estimate the corresponding IRR, and (3) the proposed peak-conditioned fusion method then combined the IRR estimates to calculate the final IRR. We validated the implemented method with capnography and nasal/oral airflow as the reference RR using the limits of agreement (LOA) approach. Compared to simple fusion and single respiratory-induced variation estimations, peak-conditioned fusion shows better performance. It provided a bias of 0.28 bpm with the 95% LOAs ranging from −3.62 to 4.17, validated against capnography and a bias of 0.04 bpm with the 95% LOAs ranging from −5.74 to 5.82, validated against nasal/oral airflow. This algorithm would expand the functionality of a conventional pulse oximetry beyond the measurement of heart rate and oxygen saturation to measure the respiratory rate continuously and instantly.


INTRODUCTION
Respiratory rate (RR), along with other vital signs like heart rate (HR) and blood pressure, is monitored for primary or continuous assessment of patient wellness. There is significant evidence that an abnormal respiratory rate is an important predictor of serious illness. For example, in children aged 1-5 years old, an elevated RR (>40 breaths/min) is an important criterion for the diagnosis of pneumonia (WHO, 2013). Furthermore, Fieselmann et al. analyzed the measurements of vital signs during the 72 h prior to cardiac arrest and showed that a high respiratory rate (>27 breaths/min) was a significant predictor of cardiac arrest in hospitals (Fieselmann et al., 1993). In addition, Subbe et al. showed that relative changes in respiratory rate are much more significant than changes in HR or systolic blood pressure in unstable patients and therefore the respiratory rate is more likely to be a better predictor for identifying the patient at risk (Subbe et al., 2003). RR can be measured by a nurse counting the number of times the chest rises in 1 min (Lovett et al., 2005). Continuous monitoring of RR, though, needs a monitoring device and can be performed using capnography, transthoracic impedance pneumography, nasal/oral pressure transducers and abdominal/thoracic respiratory inductance plethysmography belts, among others. However, recent studies have found that neither the nurses nor the monitoring devices provide accurate and reliable measurements of RR (Lovett et al., 2005). Therefore, there is a clear need for a robust, automatic, reliable and noninvasive measure of RR for performing a spot-check and for continuous monitoring.
Analysis of the photoplethysmogram (PPG) recorded using a pulse oximeter could offer an alternative method for monitoring RR. The PPG waveform contains information about a wide range of physiological parameters such as HR, heart rate variability (HRV), oxygen saturation (SpO 2 ), vascular tone, blood pressure, cardiac output and respiration (Shelley, 2007). However, most conventional pulse oximeters only provide information about HR and SpO 2 . In this study, we have presented an algorithm for robust estimation of instantaneous respiratory rate (IRR) from PPG with the aim of developing a portable solution based on pulse oximetry, suitable for both continuous monitoring and spot-check applications.

Background
A pulse oximeter measures the blood volume changes or PPG, based on the light absorption characteristics of the blood at the measuring site on the body (e.g., finger, forehead, and earlobe). A conventional transmission pulse oximeter sensor consists of two LEDs, red and infra-red, and a photo-detector. The LEDs emit light at both red and infrared wavelengths (650 and 950 nm, respectively) and the photo-detector captures the transmitted light. Based on the Beer-Lambert's law, the transmitted light density decreases during systole when the peripheral arterial blood volume is at its maximum value and increases during diastole. The PPG signal generated in the photo-detector, then, has a pulsatile waveform (AC) whose peaks and troughs reflect light transmitted through the tissue when blood volume is minimal and maximal, respectively. PPG also has a large baseline component (DC) mainly rises because of constant absorption and scattering of light traveling through skin, bones and tissues. The small variation observed in DC is mostly due to venous blood variation which changes the intensity of the light captured by the photo-detector (Webster et al., 1997). Respiration may induce variation in PPG in three different ways (Meredith et al., 2012) (Figure 1): (1) Respiratory-Induced Intensity Variation (RIIV): Changes in venous return due to changes in intra-thoracic pressure throughout the respiratory cycle cause a baseline (DC) modulation of the PPG signal. During inspiration, decreases in intra-thoracic pressure result in a small decrease in central venous pressure increasing venous return. The opposite occurs during expiration. As the venous bed at probing site cyclically fills and drains, the baseline is modulated accordingly. (2) Respiratory-Induced Amplitude Variation (RIAV): During inspiration, left ventricular stroke volume decreases due to changes in intra-thoracic pressure leading to the decreased pulse amplitude. The opposite happens during expiration. (3) Respiratory-Induced Frequency Variation (RIFV): HR varies throughout the respiratory cycle; HR increases during inspiration and decreases during expiration. This phenomenon well-known as respiratory sinus arrhythmia (RSA) is mainly due to the autonomic regulation of HR during respiration.
Respiration may induce variation in PPG differently among different individuals in health and disease. For instance, RIFV, as an indicator of autonomic activity, maybe affected by diseases and disorders (e.g., myocardial infarction, diabetic neuropathy or sleep breathing disorders Dehkordi et al., 2016). RIAV and RIIV are also very sensitive to dehydration and hypovolemia. In addition, respiratory-induce variations are different for women and men (Li et al., 2010). Estimation of IRR by combining the information from three respiratory-induce waveform variations, then, improves the algorithm performance and increases the robustness of results . Many algorithms have been proposed to estimate RR from PPG. Auto-regression (Thayer et al., 2002), Fourier transform analysis , correntropy spectral density (Garde et al., 2014), digital filters (Nakajima et al., 1996) and empirical mode decomposition (Garde et al., 2013) were successfully used, among others. These algorithms have mostly focused on estimating average RR from a window of PPG. For example, Karlen et al. (2013) and Garde et al. (2014) estimated RR every second using 16, 32, 64-s windows of PPG data.
In this study, we have proposed a novel method for extracting the instantaneous respiratory rate (IRR) from PPG. The method was performed in three main steps: extraction of RIIV, RIAV, and RIFV signals from PPG, estimation of IRR from each extracted respiratory-induced variation signals and fusion of IRR estimates. A time-frequency transform called synchrosqueezing transform (SST) (Daubechies et al., 2011) was used to extract RIIV, RIAV, and RIFV from PPG. Later, a second SST was applied to estimate IIR from respiratory-induced variation signals (Addison and Watson, 2004). To combine IRR estimates corresponding to each respiratory-induced variation signal, a novel method, called peak-conditioned fusion algorithm was proposed.

Instantaneous Frequency (IF)
The instantaneous frequency is the frequency at a given time.
Consider a multicomponent signal f that can be modeled as where A k (t) and φ k (t) are the time-varying amplitude and phase of kth frequency component, respectively.
In theory, the instantaneous frequency (IF) is defined as the derivative of the phase function with respect to time as

Synchrosqueezing Transform
Synchrosqueezing Transform (SST) was first introduced by Daubechies et al. (2011) in 1996 and then implemented by Thakur et al. (2013). SST is a combination of wavelet analysis and a reallocation method which sharpens a time-frequency representation by allocating its points to another locations in the time-frequency plane. SST can provide an accurate estimation of IF.
As defined in Daubechies et al. (2011), SST involves three steps: Step 1: Estimation of the continuous wavelet transform (CWT) The CWT of f is calculated as where ψ is a wavelet withψ(ξ ) = 0 for ξ ≤ 0 and a and b are scale and location variables, respectively. ψ(ξ ) is the complex conjugate of ψ(ξ ) andψ(ξ ) is the Fourier transform of ψ(ξ ) estimated asψ Step 2: Estimation of the instantaneous frequency Ifψ(ξ ) is concentrated around ξ = ω 0 , then W f (a, b) will be spread out around the horizontal line a = ω 0 / ω on the time-scale presentation for a given frequency of ω. However, Daubechies et al. (2011) showed that the oscillation of W f (a, b) around b tends to the original frequency ω, irrespective of the value of a. Therefore, for any (a, b) where W f (a, b) = 0, the instantaneous frequency ω f (a, b) for signal f can be defined as Step 3: Transfer to the time-frequency plane In this step, each point on the time-scale plane is allocated to a point on the time-frequency plane using the map (a, b) → (ω f (a, b), b). The frequency variable ω and the scale variable a are both binned: W f (a, b) is computed only at discrete values a k , with a k − a k−1 = ( a) k and its SST, T f (ω, b), is estimated only at the centers ω l of the successive bins [ω l − 1 2 , ω l + 1 2 ], with ω l − ω l−1 = ω, by summing different points:

Capnobase Data Set
The Capnobase contains test and calibration data sets (Karlen et al., 2010). Test data set contains forty-two 8-min segments of recordings obtained from 29 pediatric and 13 adults receiving general anesthesia at the British Columbia Children's Hospital and St. Paul's Hospital, Vancouver, BC, respectively. Calibration data set contains one hundred twenty-four 2-min segments of recordings used for tuning the parameters of the proposed algorithm.
In both data sets, the recordings included ECG, capnometry, and PPG (sampled at 300, 300, and 100 Hz, respectively) obtained with S/5 collect software (Datex-Ohmeda, Finland). The capnography waveform was used as the reference gold standard recording for RR. A research assistant manually labeled each breath in the capnogram. The beginning and end of all artifacts in the PPG waveforms were also manually labeled. Both datasets can be downloaded from the on-line database, CapnoBase.org.

Sleep Data Set
Sleep database contains forty-three 20-min segments of recording from 43 children (age = 9.1 ± 4.1, Apnea/Hypopnea Index = 8.9 ± 17.9) referred to the British Columbia Children's Hospital for overnight standard polysomnography (PSG). The children had been recruited following approval by the University of British Columbia Clinic Research Ethics Board (H11-01769) and informed parental consent. Parental/guardian written informed consent was obtained for all children. Children with a cardiac arrhythmia or abnormal hemoglobin were excluded from the study.
In addition to PSG, the PPG was recorded simultaneously using the Phone Oximeter TM (Karlen et al., 2011) sampled at 62.5 Hz with 32-bit resolution.
For each subject, the 20 min of recordings were selected after 2 h of sleep. Apnea/hypopnea events were annotated as the artifacts and no specific sleep stage was selected.
The nasal/oral airflow waveform was used as the reference gold standard recording for RR. Two experts manually labeled each breath in nasal/oral airflow waveform. The beginning and end of all artifacts in the oral/nasal waveforms were also manually labeled.

Estimation of IRR From PPG
To perform IRR estimation, after a preprocessing stage, a first SST was applied to PPG to extract RIIV, RIAV and RIFV. Later, a second SST was performed to estimate IIR from a respiratory-induced variation signals. The peak-conditioned fusion algorithm was then used to fuse simultaneous IRR estimates. This procedure, inspired by the method known as secondary wavelet feature decoupling (SWFD) (Addison and Watson, 2004), involves the following steps (Figure 2): (1) The first SST is applied to the PPG signal.
In this study, reference ranges of cardiac and respiratory bands were extracted from a review of observational studies that used HR from 143,346 children and RR data from 3,881 children (from 6 months to 18 years old) (Fleming et al., 2011). Based on 99th and 1st centiles for children and young adults, the HR could range from 30 to 180 beats/min (0.50 to 3 Hz, respectively) and RR from 8 to 54 breaths/min (0.14 to 0.9 Hz, respectively). The range in adults is much more restricted, thus it would be included in this range.
(3) The respiratory component in SST surface shows RIIV and its ridge in the frequency-time plane represents RIIV-derived IRR (IRR riiv ) (Figure 3). (4) The ridge of cardiac component is followed either in the amplitude-time plane to get RIAV or in the frequency-time plane to get RIFV. This is done by projecting the cardiac ridge points onto the amplitude-time or frequency-time planes, respectively. (5) The second SST applied to RIAV results in a dominant single component in the respiratory band (0.14-0.9 Hz, 8-54 breaths/min) whose ridge represents RIAV-derived IRR (IRR riav ). (6) A second SST is applied to the RIFV signal as well to get a dominant single component in the respiratory band whose ridge represents RIFV-derived IRR (IRR rifv ). (7) Estimation of final IRR (IRR ppg ) is performed using a proposed peak frequency tracking method (so-called peakconditioning fusion) which combines the instantaneous frequency information from (IRR riiv ), (IRR riav ) and (IRR rifv ).

Preprocessing
The PPG signals were lowpass filtered by a lowpass Chebyshev Type I IIR filter of order 8 and down sampled to 10 Hz.
FIGURE 2 | To extract IRR from PPG, the first SST was applied to PPG to extract RIIV, RIAV, and RIFV. Later, the second SST was performed to estimate IIR from a respiratory-induced variation signals. The peak-conditioned fusion algorithm was then used to fuse simultaneous IRR estimates.

Estimation of IRR riiv
Consider a PPG signal as a vector ppg ∈ R n , n = 2 L+1 where L is a nonnegative integer. The CWT of ppg, W ppg , was calculated using the Morlet wavelet, ψ, where its Fourier transform was concentrated around 1.25 Hz. The W ppg was sampled at the location (a j , b), where a j = 2 j/n v , j = 1, ..., Ln v , n v = 32 and b = 1, ..., n. The result is a Ln v × n matrix denotedW ppg . WhenW ppg > 0,ω ppg was implemented as follow where D bWppg was the finite differences ofW ppg with respect to b.
The SST of PPG was calculated as T ppg over time shows both cardiac and respiratory bands (Figure 3). A ridge fitting the dominant area of T ppg in the respiratory band (0.14-1 Hz) represented IRR riiv and was extracted by tracking the local maximum values in this region.

Estimation of IRR riav
Consider RIAV as a vector riav ∈ R n , where n is the length of ppg. In the amplitude-time plane of T ppg , riav estimated as a ridge fitting the dominant area of T ppg in the cardiac band (0.5-3 Hz, 30-180 beats/min). The ridge extracted by finding the local maximum values which minimize the following cost function (Abid et al., 2007): (9) The SST of riav, T riav was calculated using the same implementation described in the previous section.
A ridge fitting the dominant area of T riav in the respiratory band (0.14-0.9 Hz) represented the RIAV-derived IRR (IRR riav ) and can be extracted by tracking the local maximum values in this region.

Estimation of IRR rifv
Consider RIFV as a vector rifv ∈ R n , where n is the length of ppg. In the frequency-time plane of T ppg , rifv estimated as a ridge fitting the dominant area of T ppg in the cardiac band (0.5-3 Hz, 30-180 beats/min). The ridge extracted by finding the local maximum values which minimize the following cost function (Abid et al., 2007): (10) The SST of riav, T riav was calculated using the same implementation described in the section 3.2.2.
A ridge fitting the dominant area of T riav in the respiratory band (0.14-0.9 Hz) represented the RIFV-derived IRR (IRR rifv ) and can be extracted by tracking the local maximum values in this region.

Peak-Conditioned Fusion
The peak-conditioned fusion method, inspired by Lázaro Plaza (2015), was proposed to combine the IRR estimates from three respiratory-induced variations to provide the final IRR ppg .
The calculated T ppg , T riav and T rifv are two-dimensional matrices ∈ R Ln v n , n = 2 L+1 where L is a nonnegative integer and n v = 32. Each column of T ppg , T riav and T rifv matrices shows the frequency distribution of PPG, RIAV and RIFV signals at each time instance, respectively. To reduce the variance, each matrix is averaged in time dimension using a moving window of length T m = 16 s every t s = 5 s. The averaged matrix is denoted asT k , where k refers to ppg, riav, or rifv (Figure 4).
At instant b, the location of the largest peak in respiratory band of eachT k (:, b) column (for k = ppg, riav or rifv) is detected and denoted as ir I k (b). Then, a reference frequency interval, k (b), was defined as where (b − 1) was a respiratory rate reference estimated from the b − 1 previous step. All peaks larger than 85% of ir I k (b) inside (b) were detected and ir II k (b) was chosen as the nearest to f (b−1) . By reaching to this point, ir II riiv (b), ir II riav (b) and ir II rifv (b) were available simultaneously. The final respiratory peak at instant b, IIR ppg ((b)), was then chosen among ir II riiv (b), ir II riav (b) and ir II rifv (b) estimates with the largest P k . P k is a measure of the peakness and was defined as the ratio of power contained in an interval centered around the largest peak to the power of k (b). P mathematically calculated as Estimation of respiratory rate as the largest peak in the respiratory band would increase the risk of choosing the location of false peaks. To decrease this risk, the search for the largest peak was limited to the reference frequency interval, k (b) (Lázaro Plaza, 2015). This is an asymmetric interval of 3δ centered around a reference frequency. At each step the respiratory rate reference was updated using where f (b) = arg max(T k (:, 1)) in the frequency band of [0.2Hz, 0.7Hz].
Value of δ was set as 0.1 and the value of a was tuned as 0.6 over the calibration data set.

Algorithm Evaluation
To evaluate the performance of SST-based algorithms, agreement between reference IRR and estimated IRR (using peakconditioned fusion, simple fusion, single respiratory-induce variation) was assessed using the limits of agreement (LOA) technique. The bias and 95% LOA were estimated using the Bland-Altman plot. Since for each subject multiple measurement were observed, the Bland-Altman method for multiple observations per individual (Zou, 2013) was used instead of the standard Bland-Altman method. The bias was calculated as mean of IRR est -IRR ref and the 95% LOAs as mean bias ± 1.95 standard deviations. Two standard deviations (2SD) were also estimated in the purpose of ranking the proposed algorithm in this study based on the statistical analysis reported by Charlton et al. (2016).
The coverage probability (CP 2 ) was also reported as the probability of measurement error falling within pre-defined bounds, set as 2 breaths per minute (bpm) in this study (Barnhart et al., 2007).
In addition, the performance of algorithms was assessed using the root mean square error (RMSE) (breaths/min) defined as where n is the number of observations and IRR ref and IRR est are the reference and the algorithm estimates, respectively.

Capnobase Data Base
IRR extracted from the capnography waveform (IRR CO2 ) was used as the reference gold standard. The distribution of the respiratory rates contained 3,542 data points estimated every 5 s from IRR CO2 for the 16 s moving windows over the whole dataset. The respiratory rates ranged from the lowest value of 3.6521 bpm to the highest value of 44.22 bpm. The mean rate was 15.02 bpm with standard deviation of 7.66 bpm. About 7.7% of the data points were excluded from the further analysis due to to poor signal quality of the capnography signals. For each algorithm, the measures of agreement between the estimated IRR from PPG (IRR est ) and IRR CO2 and also RMSE were estimated (Table 1). For peak selection algorithm, bias was estimated as 0.28 bpm with the 95% LOAs from −3.62 to 4.17 (Figure 5). The value of 2SD was estimated as 3.97 bpm. The value of CP 2 showed that for 89% of the IRR estimates, the measured error was less than 2 breaths/min. RMSE was estimated as 1.8 bpm.
The values of 2SD of the other algorithms ranged from 8.32 bpm to 16.00 bpm.

Sleep Database
IRR extracted from the nasal/oral airflow waveform (IRR nas ) was used as the reference gold standard in the sleep dataset. The distribution of the respiratory rates contained 10,553 data points estimated every 5 s from IRR nas over the 16 s moving window for all subjects. The respiratory rates ranged from the lowest value of 9.561 bpm to the highest value of 50.85 bpm. The mean rate was 18.64 bpm with standard deviation of 5.66 bpm. About 0.66% of the data points were excluded from the further analysis due to to poor signal quality of the nasal/oral airflow signals.
The measures of agreement between the estimated IRR from PPG (IRR est ) and IRR nas also RMSE were estimated for each algorithm ( Table 1). For peak selection algorithm, bias was estimated as 0.04 bpm with the 95% LOAs from −5.74 to 5.82 (Figure 6). The value of 2SD was estimated as 5.90 bpm. The value of CP 2 showed that for 85% of the IRR estimates, the measured error was less than 2 breaths/min. RMSE was estimated as 2.3 bpm.
The values of 2SD of the other algorithms ranged from 8.4 bpm to 21.3 bpm.

DISCUSSION AND CONCLUSION
In this study, we presented an algorithm to extract IRR from PPG. We extracted RIIV, RIAV, and RIFV from PPG using Synchrosqueezing Transform (SST), a sharpening timefrequency method which provides instantaneous frequency rate. The peak-conditioned fusion was proposed to combine the extracted information from three respiratory induced variations waveforms to estimate respiratory rate at each instance. We validated the implemented method with capnography and nasal/oral airflow as the reference RR. Compared to simple fusion and single respiratory-induced variation estimations, peak-conditioned fusion shows better performance ( Table 1). It  provided a bias of 0.28 bpm with the 95% LOAs ranging from −3.62 to 4.17, validated against capnography (in the Capnobase dataset) ( Figure 5) and a bias of 0.04 bpm with the 95% LOAs ranging from −5.74 to 5.82, validated against nasal/oral airflow (in the Sleep dataset) (Figure 6). In this study, the proposed method estimated IRR from three sources of respiratory-induced variation and fused the estimated rates to measure the final IRR. Our findings showed that fusion of estimation rates would increase the accuracy and robustness of RR estimation. Even the simple fusion compared to single respiratory-induced variation estimations showed higher rank (narrower 2SD and greater CP2). It is consistent with the findings of Li et al. (2010) that respiratory activity may induce variation in PPG differently in different individuals. As discussed by Karlen et al. (2013), ventilatory conditions (spontaneous or mechanical ventilation) can change the behavior of respiratory induced variations.
In this study, we applied the proposed algorithm to two different data sets to include a broad range of subjects into the study. The Capnobase data set includes children adults, under controlled ventilation or spontaneously breathing over a wide RR range. The subjects were under general anesthesia and were continuously monitored. The sleep dataset includes children from 1-month to 17 years old spontaneously breathing during the overnight sleep in a sleep lab. During the overnight recordings, respiratory rates might change significantly during different sleep stages (light sleep, deep sleep or rapid eye movement (REM) sleep). In addition, some of the children may have experienced periods of breathing cessation, or obstructive sleep apnea. Since the apnea/hypopnea events induce changes on the morphology of nasal/oral airflow, i.e., flat line during apnea or very fast and lowamplitude oscillations during hypopnea, to extract the accurate reference RR from the nasal/oral airflow waveform, we labeled the apnea/hypopnea events as the artifact on nasal/oral airflow waveform.
A recent study (Charlton et al., 2016) presented a very complete assessment of RR estimation using PPG. A wide range of available techniques for estimation of respiratory-induced variations from PPG, estimation of RR from respiratory-induced variations, and fusion of RR estimates were identified and then more than 300 algorithms were implemented by assembling all possible combinations of available techniques. The algorithms were applied to the Vortal data set of 39 resting adults, with RR ranging between 5 and 32 bpm and ranked based on 2SD. The first ten top-ranked algorithms had the 2SD values ranging from 6.2 to 7.9. To be able to compare the performance of our algorithm with the ten top-ranked algorithms, we applied our proposed method on the same data set. For peak conditioned algorithm, the value of 2SD was estimated as 6.83 bpm. It ranked our algorithm among the top ten algorithms; however, the performance of our algorithm slightly decreased when RR was lower than 8 breaths/min. In addition, for the top ranked algorithms, the value of CP 2 was reported as 71.5 while we obtained a CP 2 of 75 applying our proposed algorithm.
It is important to note that all top ranked algorithm reported in Charlton et al. (2016) estimated RR using 32-s windows while our method can estimate RR instantaneously. It suggests that our algorithm shows better performance compared to methods that extract IRR based on time-frequency analysis (Li et al., 2010;Addison et al., 2012).
In Charlton et al. (2016), the methods for extracting RR from ECG were assessed as well. The findings of that study showed that algorithms performed better when using ECG than PPG. The best algorithm had 95% LOAs of −4.7 to 4.7 bpm and a bias of 0.0 bpm when using the ECG.
In Charlton et al. (2016), the performance of thoracic Impedance Pneumography (IP) were assessed as well providing a bias of −0.2 bpm with 95% LOAs of −5.6 to 5.2 bpm. Thoracic IP is a commonly-used technique for continuous monitoring of RR that measures changes in the electrical impedance of the person's chest during respiration. Our results showed that the performance of our algorithm is comparable with the performance of thoracic IP.
Several studies based on the continuous wavelet transform (CWT) (Clifton et al., 2007;Addison et al., 2012), the short-time Fourier transform (STFT) (Shelley et al., 2006), and empirical mode decomposition (EMD) (Garde et al., 2014) have been proposed to detect RR from PPG. The results of a study conducted by Thakur et al. (2013) to compare SST to CWT, STFT and EMD showed the superior precision of SST at identifying components of complicated oscillatory signals. Moreover, the study showed that time-varying instantaneous frequencies could be clearly distinguished in the SST while there is much more smearing and distortion in the CWT and STFT.
In this study, the SST algorithm has been implemented with O(nlog(n)) computations per scale and peak-conditioned fusion algorithm has been implemented with O(n) computations where n is the number of samples. However, there are some fast algorithms to reduce the number of computations of SST implementation to O(n) per scale Rioul and Duhamel (1992).
This study introduces a new method to estimate IRR from pulse oximetry. This would expand the functionality of a conventional pulse oximetry beyond the measurement of HR and SpO 2 to measure the respiratory rate continuously and instantly in the clinical setting and at home. Importantly, these are all achievable with a simple, cheap, single-sensor solution.

AUTHOR CONTRIBUTIONS
All authors made substantial contributions to the conception and design of the paper. PD contributed to data acquisition and processed the obtained data, designed, developed and evaluated the algorithm, and drafted the review. AG contributed to data acquisition and the design of algorithm and revised the paper critically for content. BM contributed to the design of algorithm and revised the paper critically for content. JA and GD revised the paper critically for content. AG and BM have contributed equally to this work. All authors approved the final submission of the document and agreed to be accountable for all aspects of the work.