 Research
 Open Access
Spectral fusionbased breathing frequency estimation; experiment on activities of daily living
 Iman Alikhani^{1}Email authorView ORCID ID profile,
 Kai Noponen^{1},
 Arto Hautala^{1},
 Rahel Ammann^{2} and
 Tapio Seppänen^{1}
 Received: 21 February 2018
 Accepted: 21 July 2018
 Published: 27 July 2018
Abstract
Background
We study the estimation of breathing frequency (BF) derived from wearable singlechannel ECG signal in the context of mobile daily life activities. Although respiration effects on heart rate variability and ECG morphology have been well established, studies on ECGderived respiration in daily living settings are scarce; possibly due to considerable amount of disturbances in such data. Yet, unobtrusive BF estimation during everyday activities can provide vital information for both disease management and athletic performance optimization.
Method and data
For robust ECGderived BF estimation, we combine the respiratory information derived from R–R interval (RRI) variability and morphological scale variation of QRS complexes (MSV), acquired from ECG signals. Two different fusion techniques are applied on MSV and RRI signals: crosspower spectral density (CPSD) estimation and power spectrum multiplication (PSM). The algorithms were tested on large sets of data collected from 67 participants during office, household and sport activities, simulating daily living activities. We use spirometer reference BF to evaluate and compare our estimations made by different models.
Results and conclusion
PSM acquires the least average error of BF estimation, \(\%D^{2\sigma }=9.86\) and \(\%E = 9.45\), compared to the reference spirometer values. PSM offers approximately 25 and 75% less error in comparison with the CPSD fusion estimation and the estimation by those two exclusive sources, respectively. Our results demonstrate the superiority of both of the fusion approaches, compared to the estimation derived from either of RRI or MSV signals exclusively.
Keywords
 Heart rate variability
 Time–frequency analysis
 Singlechannel ECG
 Breathing rate estimation
 ECG morphology
Background
Breathing frequency (BF) is a vital biomarker utilized for diagnostics, and sport physiology applications. However, measuring BF using respiratory sensors over longterm monitoring sessions can be uncomfortable. Indirect monitoring of respiratory frequency can be conducted in different modalities, including videobased [1], electrical impedance pneumographybased [2] or wearable accelerometerbased respiration reconstruction [3]. Single or multiple channel electrocardiography (ECG) is one of the wellestablished modalities explored for BF and respiratory pattern reconstruction, known as ECGderived respiration (EDR).
EDR was introduced during the 80s in [4, 5]. The idea was based on the fact that respiratory sinus arrhythmia (RSA)—which is obtainable from heart rate variability (HRV)—correlates with the respiratory pattern. Additionally, the beattobeat morphological variation of ECG signal, e.g. modulation of Rpeak amplitudes or areas under Rwave, is largely a result of inhale and exhale. Thus, it is feasible to derive respiratory information such as the BF indirectly by analyzing a singlechannel ECG signal.
The development of wearable devices has made it practical and inexpensive to monitor the biosignals of subjects who can benefit from continuous monitoring, including patients suffering from sleep apnea [6–8], and professional athletes managing their exercise regimes according to the biosignal feedback [9]. Consequently, there is an increasing demand on the biosignal processing algorithm development to enhance the capabilities of wearables and/or to introduce new features.
In order to rely on BF reconstruction from singlechannel wearable ECG signal, the algorithm employed and its performance should be validated across a variety of daily activities. We found a few studies about BF estimation during physical activities using singlechannel ECG signals. However, most of the research in this discipline has contributed to the resting state EDR [10–12]. What is more, severe challenges are introduced to ECG processing in physical activity contexts, including variable mean heart rate (HR), high level of movement artifacts and introduction of cardio–locomotion coupling (CLC) components [13, 14].
In this paper, we hypothesize that the combined information taken from different existing sources of respiratory signal in ECG indices yields a more accurate and reliable BF estimation. To this end, we tested two frequencydomain data fusion techniques on two respiratorily modulated indices. Namely, we apply crosspower spectral density (CPSD) estimation and power spectrum multiplication (PSM) on R–R interval variability (RRI) and morphological scale variation of QRS complexes (MSV). We tested our hypothesis on a data set collected from 67 subjects in reallife activities such as office, households and sports.
Methods
Preprocessing
This preprocessing step is followed by signal conditioning wherein anomalies from RRI are detected and replaced with linear interpolation to keep the number of beats unaltered. This ectopic beat detection and editing is explained in more details in the following subsection.
Ectopic beat detection and editing
In practice, the RRI signal derivation must be followed by ectopic and anomaly beat detection and editing [17, 18]. Especially when the data is recorded during physical activities via wearables, the signal quality of ECG is generally lower, which can cause problems in the Rpeak detection. During such activities, motion artifacts and ectopic beats are also more abundant. Hence, ectopic beat detection and editing is an important preprocessing step in HRV analysis.
 (a)
Detect the evident outliers and edit them using linear interpolation. Although the HR ranges from 60 beats/min (bpm) to 200(bpm) in our dataset, we define a wider range for healthy intervals between 250 and 1500 ms. This is to ensure that healthy controlled beats are not discarded and to make the outlier detector more general for any dataset.
 (b)Construct an nbythree dimensional feature matrix describing local HRV around each beat. More precisely, the following three descriptors are calculated for each beat:

\(RR_i,\)

\(\mid RR_iRR_{i1}\mid +\mid RR_iRR_{i2}\mid\)

\(\mid RR_iRR_{i+1}\mid +\mid RR_iRR_{i+2}\mid\)
At the borders, use symmetric padding by mirroring the samples.
(\(\{RR_1, RR_2, RR_3,\ldots , RR_{n2}, RR_{n1}, RR_n\} \Longrightarrow\)
\(\{RR_3, RR_2, RR_1, RR_2, RR_3,\ldots , RR_{n2}, RR_{n1}, RR_n, RR_{n1}, RR_{n2}\}\)).
Duration unit of millisecond is used.

 (c)
Normalize the feature matrix columnwise for zeromean and unitstandard deviation.
 (d)A multivariate Gaussian distribution is fitted on the feature matrix using \(\mu \in \mathrm{I\!R}^3\) and the covariance matrix of features \(\sigma \in \mathrm{I\!R}^{3\times 3}\). Given the feature set \(\{ x^{(1)}, x^{(2)},\ldots , x^{(n)} \}\), where every \(x^{(i)} \in \mathrm{I\!R}^3\),The model can be fitted by computing$$\begin{aligned} p(x;\mu ,\sigma ) = \frac{1}{(2\pi )^{\frac{n}{2}}\mid \sigma \mid ^{\frac{n}{2}}}e^{\frac{1}{2}(x\mu )^T\sigma ^{1}(x\mu )}. \end{aligned}$$(1)and$$\begin{aligned} \mu = \frac{1}{n}\sum ^n_{i=1}x^{(i)} \end{aligned}$$(2)$$\begin{aligned} \sigma =\frac{1}{n}\sum ^n_{i=1}(x^{(i)}\mu )(x^{(i)}\mu )^T \end{aligned}$$(3)
 (e)
Detect as anomalous beats those with \(p(x^{(i)})\le \epsilon\) where \(\epsilon\) is a threshold value. The value of \(\epsilon\) tunes the strictness of the anomaly detection and is adjusted according to the HR level. In this work \(\epsilon\) is set in the range of \([10^{11},10^{5}]\), where lower range was used for high exercise intensity data, and upper range for nonphysical context.
 (f)
Edit the detected anomalies using linear interpolation by nearby controlled beats to keep the same number of beats. Linear interpolation is selected because it is a widely used method.
 (g)
Steps (b) to (f) are iterated as long as there is a newly detected anomaly, otherwise the procedure is terminated.
Spectral analysis
Naturally, the heart is beating irregularly. Hence, the constructed signals (RRI and MSV) are not evenly sampled. Therefore, we resample the signals to make them equidistant which also prepares them for spectral analysis, as suggested in [19]. We interpolate the signals by 8 Hz sampling frequency and compute their baseline using a 5s window moving average smoothing function. We apply the standardization procedure in [20] for correction of mean HR, induced by physical activities as follows: The baseline is subtracted from the signal and divided up by the baseline, i.e. (signalbaseline)/baseline.
The HRV spectral components during physical activities differs from those of resting condition [14]. During steady state resting condition, high frequency (HF) components of HRV spans the range of [0.15 Hz 0.4 Hz] [21], but during exercise this range is redefined to be [0.15 Hz \(\frac{HR_m}{2}],\) in which \(\frac{HR_m}{2}\) is the half of local mean HR^{1} in Hertz.
Considering the heart as a natural signal sampler with the sampling frequency of \(HR_m,\) \(\frac{HR_m}{2}\) is the intrinsic Nyquist frequency of HRV signal, meaning that interpretable physiological components should fall within the range of HF. Aliasing phenomenon of the components over this intrinsic Nyquist frequency of HRV signal might lead to misinterpretation of autonomic nervous system (ANS) activity. For instance, let’s assume an abnormal BF just a bit over the intrinsic Nyquist frequency, e.g. \(\frac{HR_m}{2}+\Delta F\). Due to the sampling theorem, it will fold back into the HF range of interest at the frequency of \(\frac{HR_m}{2}\Delta F\). In the spectral analysis this could lead into misguided detection of BF at this lower frequency, instead of the higher actual one. Speaking of which, a prevalent factor that might influence the spectral interpretation of HRV signal is CLC components that arise from cadence during walking or running; or pedaling frequency during cycling [13, 14]. Because of the mentioned aliasing phenomenon, these components will fold back to the HF range when they exceed the Nyquist frequency. The folded components might possess significant energy level compared to the energy at BF at some time instants. Our proposed fusion models aim to highlight the joint BF components in both signals and attenuate the influence of existing disturbances exclusively in either of those signals.
Spectral fusion
The existence of energy around the BF range in HRV and morphologicallyderived signals is justified in the literature as noted in the Introduction. Moreover, some of the possible challenges related to these signals were also described in the previous subsection.
To enhance the estimation of BF particularly during a nonstationary recording situation (e.g. physical activities), we propose to combine the spectral information of BF components of RSA and morphological variation (situated in RRI and MSV signals, respectively). We address the spectral fusion in two different ways: by PSM and CPSD estimation of MSV and RRI signals. The key point is that potentially there are mutual joint energies in the spectrum of these two signals, supposedly greater than the background energy and corresponding to the BF at each time instant. Thus, attempting to find the significant joint energy content present in both signals at each time instant makes sense, since the influence of unwanted distortions (e.g. CLC components) between those two sources might vary.
Elementwise PSM basically enhances the shared joint energy bands and diminishes the unshared ones [23]. This technique is advantageous in strengthening shared energy bands, potentially also the BF component. Assume the time–frequency representation of RRI and MSV respectively as \(S_{xx}(t_i,f_j)\) and \(S_{yy}(t_i,f_j)\), where i and j represent the time and frequency indices. The spectral multiplication at each time and frequency index can be computed by elementwise multiplication of the amplitude spectrum (\(S_{xx}(t_i,f_j) \times S_{yy}(t_i,f_j)\)) and is expressed as \(\mid S_{xx} \cdot S_{yy}\mid\).
We quantify the coupling between RRI and MSV signals using squared magnitude of \(S_{xy}(t,f)\), assuming that BF is a joint frequency component among those signals. The second last subfigure of Fig. 2 exemplifies the normalized (columnwise Zscore normalization) squared CSPD magnitude for sample sets of RRI and MSV depicted in the first and second subfigures. Similarly, the normalized spectral multiplication is depicted in the last subfigure. The coupling of the time series in the spectrums (high energy narrow band) is closely aligned with the reference BF (depicted as dashed lines) in the last two subfigures. We derive the frequency which owns the maximum power at each time instant as the BF estimation.
Results
Data
Since our goal is to evaluate the proposed method in reallife context, data was collected from measurements in different activity protocols. In total, there are 67 subjects (30 female and 37 male) aged from 18 to 60 years who participated in eight different protocol phases, including office, household and sport activities. General physical characteristics of the participants are summarized in Table 1.
General characteristics of participants
Characteristic  Mean  Min  Max 

Height (cm)  175  160  195 
Weight (kg)  75.4  45.6  122.8 
Age (years)  37.9  18  60 
BMI (kg/m^{2})  24.51  14.72  35.5 
Table of exercise intensity
Activity protocol  Mean  Min  Max 

Office work  45  28  66 
Emotional stress  43  26  73 
Floor sweeping  52  35  73 
Tidying up  54  38  75 
Table cleaning  50  34  77 
Walking  52  36  74 
Cycling  66  48  83 
Running  75  48  91 
Evaluation metrics
Table of results
Activity protocol  \(\mid S_{xy}\mid ^2\)  \(\mid S_{xx} \cdot S_{yy}\mid\)  \(S_{xx}\)  \(S_{yy}\)  

\(\%D^{2\sigma }\)  \(\%E\)  \(\%D^{2\sigma }\)  \(\%E\)  \(\%D^{2\sigma }\)  \(\%E\)  \(\%D^{2\sigma }\)  \(\%E\)  
Office work  13.09  7.82  8.65  5.61  16.48  11.56  14.08  12.21 
Emotional stress  13.62  10.98  9.03  8.68  16.47  11.67  16.70  14.42 
Floor sweeping  13.43  11.08  10.74  9.03  22.17  15.77  18.75  11.57 
Tidying up  17.96  13.74  14.05  11.50  25.37  17.94  21.25  14.18 
Table cleaning  16.91  13.01  9.93  9.05  22.21  15.82  18.51  12.33 
Walking  11.31  11.31  8.51  8.74  17.35  15.99  23.00  16.95 
Cycling  8.10  9.44  5.40  6.96  21.35  16.97  10.10  11.71 
Running  15.59  19.17  12.58  16.05  28.20  24.06  25.16  25.30 
Average  13.75  12.06  9.86  9.45  21.20  16.22  18.44  14.83 
Quantitative results
In this subsection, we evaluate the performance of estimated BF by the proposed method. The fusion methods, including squared magnitude of CPSD and the PSMbased BF estimation were compared with the BF estimation from either of RRI or MSV spectrograms [10, 11], using the metrics introduced in the previous subsection. All the software implementation and processing are done in MATLAB R2017a.
Table 3 summarizes the average figures computed for different BF estimation methods. The numbers show that both of the fusion methods outperform the BF estimation from either of the RRI or MSV spectrums exclusively. Among the two fusion methods, the PSMbased BF estimation (\(\mid S_{xx} \cdot S_{yy}\mid\)) outperforms the other one. Among all the protocol phases, cycling acquires the lowest BF estimation error in the chosen metrics. In contrast, the running and tidying up phases are the most erroneous BF estimation.
Discussion and conclusion
Computationally efficient methods to boost the performance of BF estimation via ECG surrogate signal processing could enlarge the scope of BF monitoring applications, increase useracceptance, and provide the users more accurate data. Singlechannel ECGderived BF estimation was investigated in this paper. Unlike many studies in the literature of this discipline, the purpose of our practical approach was to examine the performance of methods in a reallife like context. Thus, our database was comprised accordingly of different reallife settings, including office, households and sport activities.
The existence of breathing component in ECG and feasibility of deriving BF from the RSA and morphological variation of the signal is wellrecognized in the literature. Nevertheless, the examination of BF estimation during daily activities measured by a singlechannel wearable ECG recorder has not been largely studied. This is most likely due to the challenges involved in such a context such as noteworthy movement artifacts, variable and nonstationary HR, and CLC component’s introduction to ECG.
We proposed to fuse two sources of existing respiratory components in ECG, since those sources might be differently influenced by noise, movements, physiological factors such as age and health, as well as aliasing artifacts. Spectraldomain fusion methods, including CPSD and PSM were applied on RRI and MSV signals, constructed and derived from ECG. Table 3 shows that the performance of fusion methods in all the activity protocols in a dailylife situation is superior to the BF estimation derived from a single source, whether RRI or MSVderived BF.
Among the fusion methods, the PSM offers more accurate estimation, compared to the CPSDbased estimation. Both fusion methods are computationally efficient, while relatively PSM technique is slightly more demanding than the CPSD as two separate spectrograms are computed, and then BF is estimated from the multiplied spectrum. In CPSD, the computational cost is reduced by avoiding another transform computation and then the subsequent multiplication of the spectrograms. Quantitatively in a nonoptimized implementation, for two time series with 130 s of data, the elapsed time for computation of PSM fusion is 0.61 s versus 0.19 s for CPSD fusion. These numbers were acquired with an Intel(R) Core(TM)i53570 processor @ 3.40 GHz and 8 GB of RAM on a 64bit operating Windows 7.
It should also be noted that whereas the PSM only requires sufficient energies to exist at the same time for a high reading, the CPSD requires also phase difference stability between the two signals [26]. Based on our results provided in Table 3, the superiority of spectral multiplication could hence be explained by phase dispersion between the two derived time series.
In this study, local mean HR is derived from the computed baseline of RRI using a 5s window moving average smoothing function.
Declarations
Authors' contributions
IA implemented the algorithms and analyzed the data with close supervision of KN and TS in designing the research procedure and physiological advice of AH. RH was responsible for the experiments and data collection. IA,KN and TS were interpreting the performance of the methods. IA and KN were the major contributers in writing the manuscript. All authors read and approved the final manuscript.
Acknowledgements
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Availability of data and materials
The datasets analyzed during the current study are not publicly available and it was used under license of the owner for the current study. Data are however available from the authors upon reasonable request and with permission of the owner.
Funding
There is no funding source for this study.
Statement of human rights
All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Informed consent was obtained from all individual participants included in the study.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Wiede C, Richter J, Manuel M, Hirtz G. Remote respiration rate determination in video data—vital parameter extraction based on optical flow and principal component analysis. In: Proceedings of the 12th international joint conference on computer vision, imaging and computer graphics theory and applications, vol 4. VISAPP, (VISIGRAPP 2017). SciTePress; 2017. https://doi.org/10.5220/0006095003260333. INSTICC. p. 326–33
 Jeyhani V, Vuorinen T, Mäntysalo M, Vehkaoja A. Comparison of simple algorithms for estimating respiration rate from electrical impedance pneumography signals in wearable devices. Health Technol. 2017;7(1):21–31.View ArticleGoogle Scholar
 Vertens J, Fischer F, Heyde C, Hoeflinger F, Zhang R, Reindl L, Gollhofer A. Measuring respiration and heart rate using two acceleration sensors on a fully embedded platform. In: Proceedings of the 3rd international congress on sport sciences research and technology support, Vol 1. icSPORTS. SciTePress; 2015. https://doi.org/10.5220/0005604000150023. INSTICC. p. 15–23.
 Moody GB, Mark RG, Zoccola A, Mantero S. Derivation of respiratory signals from multilead ECGs. Comput Cardiol. 1985;12(1985):113–6.Google Scholar
 Moody GB, Mark RG, Bump MA, Weinstein JS, Berman AD, Mietus JE, Goldberger AL. Clinical validation of the ECGderived respiration (edr) technique, Group 1. Comput Cardiol. 1986;3:507–10.Google Scholar
 Oliver N, FloresMangas F. Healthgear: a realtime wearable system for monitoring and analyzing physiological signals. In: BSN 2006. International workshop on wearable and implantable body sensor networks, 2006. IEEE. 2006. p. 4.Google Scholar
 Ishida R, Yonezawa Y, Maki H, Ogawa H, Ninomiya I, Sada K, Hamada S, Hahn AW, Caldwell WM. A wearable, mobile phonebased respiration monitoring system for sleep apnea syndrome detection. Biomed Sci Instrum. 2004;41:289–93.Google Scholar
 Sani SN. Using respiratory sinus arrhythmia to detect obstructive sleep apnea. Health Technol. 2017;7:1–6.View ArticleGoogle Scholar
 Glaros C, Fotiadis D, Likas A, Stafylopatis A. A wearable intelligent system for monitoring health condition and rehabilitation of running athletes. In: 4th International IEEE EMBS special topic conference on information technology applications in biomedicine, 2003. IEEE; 2003. p. 276–79Google Scholar
 Cysarz D, Zerm R, Bettermann H, Frühwirth M, Moser M, Kröz M. Comparison of respiratory rates derived from heart rate variability, ECG amplitude, and nasal/oral airflow. Ann Biomed Eng. 2008;36(12):2085–94.View ArticleGoogle Scholar
 Schäfer A, Kratky KW. Estimation of breathing rate from respiratory sinus arrhythmia: comparison of various methods. Ann Biom Eng. 2008;36(3):476–85.View ArticleGoogle Scholar
 Tiinanen S, Noponen K, Tulppo M, Kiviniemi A, Seppänen T. ECGderived respiration methods: adapted ica and pca. Med Eng Phys. 2015;37(5):512–7.View ArticleGoogle Scholar
 Blain G, Meste O, Blain A, Bermon S. Time–frequency analysis of heart rate variability reveals cardiolocomotor coupling during dynamic cycling exercise in humans. Am J Physiol Heart Circ Physiol. 2009;296(5):1651–9.View ArticleGoogle Scholar
 Bailón R, Garatachea N, de la Iglesia I, Casajús JA, Laguna P. Influence of running stride frequency in heart rate variability analysis during treadmill exercise testing. IEEE Trans Biomed Eng. 2013;60(7):1796–805.View ArticleGoogle Scholar
 Pan J, Tompkins WJ. A realtime qrs detection algorithm. IEEE Trans Biomed Eng. 1985;3:230–6.View ArticleGoogle Scholar
 Noponen K, Tiinanen S, Seppänen T. Deriving respiration from the electrocardiogram by serial comparison with statistical mean shape. In: 2012 computing in cardiology. IEEE; 2012. p. 809–12.Google Scholar
 Salo MA, Huikuri HV, Seppanen T. Ectopic beats in heart rate variability analysis: effects of editing on time and frequency domain measures. Ann Noninvasive Electrocardiol. 2001;6(1):5–17.View ArticleGoogle Scholar
 Peltola M. Role of editing of RR intervals in the analysis of heart rate variability. Front Physiol. 2012;3:148.View ArticleGoogle Scholar
 T.F. of the European Society of Cardiology. Heart rate variability standards of measurement, physiological interpretation, and clinical use. Eur heart J. 1996;17:354–81.View ArticleGoogle Scholar
 Bailón R, Laouini G, Grao C, Orini M, Laguna P, Meste O. The integral pulse frequency modulation model with timevarying threshold: application to heart rate variability analysis during exercise stress testing. IEEE Trans Biomed Eng. 2011;58(3):642–52.View ArticleGoogle Scholar
 Camm AJ, Malik M, Bigger J, Breithardt G, Cerutti S, Cohen RJ, Coumel P, Fallen EL, Kennedy HL, Kleiger RE. Heart rate variability: standards of measurement, physiological interpretation and clinical use. task force of the european society of cardiology and the north american society of pacing and electrophysiology. Circulation. 1996;93(5):1043–65.View ArticleGoogle Scholar
 Martin W, Flandrin P. Wigner–Ville spectral analysis of nonstationary processes. IEEE Trans Acoustics Speech Signal Process. 1985;33(6):1461–70.View ArticleGoogle Scholar
 Alikhani I, Noponen K, Hautala A, Ammann R, Seppänen T. Spectral data fusion for robust ECGderived respiration with experiments in different physical activity levels. In: Proceedings of the 10th international joint conference on biomedical engineering systems and technologies, Vol 5. HEALTHINF, (BIOSTEC 2017). SciTePress; 2017. https://doi.org/10.5220/0006144100880095. INSTICC. p. 88–95.
 Flandrin P, et al. Timefrequency/timescale analysis, vol. 10. Cambridge: Academic press; 1998.MATHGoogle Scholar
 Orini M, Bailón R, Mainardi LT, Laguna P, Flandrin P. Characterization of dynamic interactions between cardiovascular signals by time–frequency coherence. IEEE Trans Biomed Eng. 2012;59(3):663–73.View ArticleGoogle Scholar
 Friston KJ, Bastos A, Litvak V, Stephan KE, Fries P, Moran RJ. Dcm for complexvalued data: crossspectra, coherence and phasedelays. Neuroimage. 2012;59(1):439–55.View ArticleGoogle Scholar