Due to the positioning of the sensing elements of cardiac measurement modalities (for example ECG and BCG), the recorded signals do not only capture the periodic reaction of the heart, but also capture respiration in form of a frequency modulation and/or amplitude modulation. Thus, they can be described by so-called “intrinsic mode functions” (IMF). An IMF is charaterised by a cosine function with a time-dependent frequency and amplitude. The BCG signal can be assumed to consist of mainly two intrinsic mode functions [22, 23], in form of an amplitude modulation, i.e. an IMF associated with the heart rate superposed by one associated with the respiratory rate, i.e.
$$\begin{aligned} s(t) = \sum _{i\in \{\text {HR,RR}\}} A_\text {i}(t) \cdot \text {cos}(2\pi f_\text {i}t) + n(t) \text {.} \end{aligned}$$
(1)
Here, A describes the amplitude signal, f describes the frequency and n(t) describes additive noise and motion artefacts. HR and RR refer to the IMF for the heart rate (HR) and the respiratory rate (RR). Similar to [22, 23], to separate the respiratory IMF from the one associated with the heart, a Butterworth bandpass filter with cut-off frequencies at \({0.1}\,\hbox {Hz}\) and \({0.5}\,\hbox {Hz}\), according to respiration rates of 6 and 30 breaths per minute (BPM), can be used. The residual signal can thus be assumed to consist only of the IMF for the respiration as well as additive filtered noise and motion artefacts
$$\begin{aligned} \tilde{s}(t) = A_\text {RR}(t) \cdot \text {cos}(2\pi f_\text {RR}(t)t) + \tilde{n}(t)\text {.} \end{aligned}$$
(2)
Applying the Hilbert transform to this IMF then leads to the so-called analytic component
$$\begin{aligned} \tilde{s}_+(t) \approx A_\text {RR}(t) \cdot \text {exp}(i2\pi f_\text {RR}(t)t)\text {,} \end{aligned}$$
(3)
where i denotes the imaginary unit. The analytic component is a complex phasor, where \(A_\text {RR}(t)\) represents the waveform of the signal and the complex exponential function its periodic repetition. Therefore, the phase signal of \(\tilde{s}_+\) can be assumed to be periodic with the respiratory rate. The respiratory rate can then be extracted by differentiation of the phase signal \(\varphi _{\tilde{\mathrm{s}}_+}(t)\):
$$\begin{aligned} f_\text {RR}(t)&= \frac{1}{2\pi } \frac{d\varphi _{\tilde{\mathrm{s}}_+}(t)}{dt} \quad \text {with} \end{aligned}$$
(4)
$$\begin{aligned} \varphi _{\tilde{\mathrm{s}}_+}(t)&= \text {arg}(\tilde{s}_+)\text {.} \end{aligned}$$
(5)
This frequency signal is further median-filtered to remove artefacts due to phase discontinuities which appear periodically with the respiratory rate. Additionally, the signal is low-pass filtered to remove BCG shape-related artefacts [19] in which additional peaks in the phase signal of one respiration appear (see Fig. 3). For filtering, a Butterworth low-pass filter with a cut-off frequency of \({0.1}\,\hbox {Hz}\) is applied since the oscillation’s frequency was found to be around \({0.2}\,\hbox {Hz}\).
Furthermore, the frequency can also be calculated with a peak detection directly on the phase signal \(\varphi _{\tilde{s}_+}\) since it has clear peaks at its discontinuities, i.e. at phase jumps from \(\pi\) to \(-\pi\). An example can be seen in the bottom plot of Fig. 4. Each of the saws represents one breathing cycle. For the peak detection, it is advantageous that the phase signal is independent of amplitude changes in the BCG signal which often appear due to a shift in position of the subject on the sensor (see Fig. 4). The complete workflow is depicted in Fig. 5.
As a baseline for comparison, two methods from literature were chosen. First, the approach from Karlen et al. [14] originally for flow signals was adapted to BCG signals. The BCG signal was resampled to \({50}\,\hbox {Hz}\) and band-pass filtered with a third-order Butterworth bandpass filter with cut-off frequencies at \({0.1}\,\hbox {Hz}\) and \({0.5}\,\hbox {Hz}\) to remove the heart’s IMF. Subsequently, the signal was windowed by a Hamming window. Each segment has a length of 2048 samples which corresponds to \({40.96}\,\hbox {s}\). The signal was then transformed into the Fourier domain. The DC-component was removed and only the frequencies smaller than \({8}\,\hbox {Hz}\) are analysed. The peak in the spectrum is chosen to be the respiratory rate. Second, the method from Paalasmaa et al. [19] was used. The BCG signal here was resampled to \({300}\,\hbox {Hz}\). To discard segments in which motion artefacts occur, the signal was split into segments of \({10}\,\hbox {s}\). Each window’s peak-to-peak value, i.e. lowest to highest amplitude, was calculated. If the peak-to-peak value was found to be larger than twice the average, 15 s before and after the segment were discarded. The signal was then low-pass filtered with four low-pass filters with different cut-off frequencies, i.e. \({0.154}\,\hbox {Hz}\), \({0.22}\,\hbox {Hz}\), \({0.33}\,\hbox {Hz}\) and \({0.5}\,\hbox {Hz}\). A peak detection was applied on all four signals. For every 3 s, a respiratory rate was chosen from the four signals by the inter-breath intervals. The rate of that signal was chosen in which the variability with respect to the amplitude of the last five cycles was smallest. The variability V for each interval was calculated by
$$\begin{aligned}&V = \max _{2 \le i \le 5} = |\text {log}(b_i)-\text {log}(b_{i-1})|, \end{aligned}$$
(6)
with \(b_i\) being the amplitudes at the start of each interval.
