RSS-based respiratory rate monitoring using periodic Gaussian processes and Kalman filtering

In this paper, we propose a method for respiratory rate estimation based on the received signal strength of narrowband radio frequency transceivers. We employ a state-space formulation of periodic Gaussian processes to model the observed variations in the signal strength. This is then used in a Rao-Blackwellized unscented Kalman filter which exploits the linear substructure of the proposed model and thereby greatly improves computational efficiency. The proposed method is evaluated on measurement data from commercially available off the shelf transceivers. It is found that the proposed method accurately estimates the respiratory rate and provides a systematic way of fusing the measurements of asynchronous frequency channels.


I. INTRODUCTION
In the past years, interest toward ubiquitous health monitoring has increased and it has been envisioned that the advances in health monitoring technologies enable future smart homes that would continuously monitor the vital signs of the residents.This information could be used for: i) increasing our healthawareness; ii) diagnostic purposes by medical professionals; and iii) detecting aberrations that need immediate clinical attention.The respiration rate is one essential indicator of a person's health and hence, monitoring it is important [1], [2].
Traditional approaches for measuring the respiratory rate require body contact in which an instrument is attached to the subject's body to measure for example respiratory airflow, chest/abdominal movement, CO 2 concentration, or blood oxygen saturation [2].Such sensors are unsuitable in home healthcare applications since the person's mobility is restricted, an elderly suffering from dementia may forget to wear a sensor, and infants may remove them.An alternative is to use non-invasive methods that do not rely on body contact and popular approaches include vision-based techniques [3] and radio frequency (RF)-based methods [4]- [6], a research topic we explore in this work.
The propagation delay of the RF signal that reflects from the breathing person's moving chest cavity is quasi-periodic in time [4].Measuring the propagation delay is possible for example with ultra-wideband radios [4] and frequency modulated carrier wave radars [6].These works show that respiration rate can be estimated with high accuracy and the technologies can also be used to estimate the heart rate.As a drawback, these systems rely on customized technology and require a large bandwidth.
Narrowband receivers on the other hand can not resolve individual multipath components.However, such devices are cheap and received signal strength (RSS) measurements are ubiquitously available in nearly all radios.In [5], it is shown that sinusoidal changes in the propagation delay also cause periodic variations in the RSS and breathing rate estimation is possible if many transceivers [5] or frequency channels [7] are used.These systems measure the RSS between different wireless links, then the power spectral density (PSD) is calculated over a 15 s-60 s time window for each link, the PSDs are averaged over all links, and the maximum of the mean PSD is used as the breathing rate estimate.This approach is brute force and does not take into account that the RSS measurements for each link are sampled asynchronously and that breathing is not perfectly sinusoidal but also contains higher order harmonics.
In this paper, we develop a recursive, non-parametric approach to respiratory rate tracking.Specifically, we model the oscillatory nature of breathing using a state-space formulation of periodic Gaussian processes and subsequently apply a Rao-Blackwellized unscented Kalman filter for estimation.The method can easily handle asynchronous and missing (e.g.due to packet loss) measurements from multiple sources such as multiple radio channels and does not rely on spectral analysis.The proposed method is similar to the ones proposed in [8], [9] in the context of electrocardiography and oximetry-based respiratory rate monitoring, which use batch-based Gaussian processes.However, the method proposed herein distinguishes itself by using a state-space formulation and Kalman filtering which has the beneficial properties mentioned earlier.

II. MODEL
The chest's motion due to breathing causes changes in the properties of the radio path, which reflects as a change in the RSS.Due to the quasi-periodic nature of this motion, the changes in RSS are quasi-periodic too.These changes can be modeled using, for example, Fourier series or (quasi-)periodic Gaussian processes (GP).In this work, the latter approach is chosen due to its flexibility and robustness [10], [11].
The measured RSS for channel i (i = 1, . . ., C) at time t n can thus be written as where r i,n denotes the measurement noise that is modeled as independent, identically distributed white Gaussian noise as r i,n ∼ N (0, R i ) [12].Furthermore, g i (t) is a GP according to where GP(µ(x), k(x, x )) denotes a Gaussian process with mean function µ(x) and covariance kernel k(x, x ), and τ = t − t .For periodic GPs, a suitable covariance function is the canonical periodic covariance function [11], [13] given by where σ 2 is the variance magnitude, ω the frequency of the periodicity, and l the GPs length scale.
It has been shown that such periodic GPs can be written in the form of an infinite-dimensional state-space system as [11] where γ j (t) are 2 × 1 vectors, γ 0 (t) is a scalar, H j = 1 0 , and η γj (t) is a white noise process with spectral density S γj .
In practice, the infinite sum over j is not realizable and the series has to be truncated at some upper bound J.
The matrices A j , B j , and S γj are determined by the covariance kernel k(τ ).For the kernel in (3) and j > 0, these are where I j (x) is the jth order modified Bessel function of the first kind.Discretization of ( 4)-( 5) yields the discrete-time model with q γj ,n ∼ N (0, Q γj ), and ∆t = t n − t n−1 for j > 0.
For the DC component γ 0 (t), we have that and the discretized version is with q γ0,n ∼ N (0, Q γ0 ), and Q γ0 as in (7).
The frequency ω = 2πf of the GP is normally constant and known (or optimized as one of the GP's hyperparameters).However, here it corresponds to the breathing rate which is the time-varying quantity of interest.Thus, we introduce an additional state for the breathing frequency f (t).Since f (t) ≥ 0, we model this as geometric Brownian motion with zero drift.Let ν(t) = log(f (t)), then, the stochastic differential equation for ν(t) is given by [14] where η f (t) is a white noise process with spectral density S f .Discretization of (10) yields and Note that by introducing a time-varying frequency, the covariance function becomes time-dependent and thus, the GP is non-stationary.Collecting all the coefficients for the ith channel as T , with where ⊗ denotes the Kronecker product and the dependency of F i on ω n−1 = 2π exp(ν n−1 ) is implicit.Furthermore, we can define the covariance matrix of the complete process noise vector as Note that while the logarithm of the breathing rate ν n is a shared state for all channels there is one individual GP for each channel, accounting for the different magnitudes and phases.Also, the measurements y i,n may be completely asynchronous, depending on the communication protocol.

A. State Estimation
The state-space model (12) can readily be used together with any standard state estimation algorithm such as Kalman or particle filters (and smoothers).However, the model is linear in all states but ν n for the dynamics and completely linear in the measurements.Thus, in order to not waste computational resources, a Rao-Blackwellized variant should be used to efficiently exploit the linear substructure.We use a Rao-Blackwellized unscented Kalman filter here (see [15], [16]).
The idea is to partition the state vector x n into states that enter the model non-linearly (x n n ) and linearly (x l n ) such that T .This structure can now be exploited to obtain a reduced-complexity filtering algorithm summarized in Algorithm 1 (see, e.g., [15], [16] for details).Here, the M = 2N + 1 sigma-points (N = 1 is the dimension of the non-linear subspace) and their weights as obtained from the unscented transform are given by [17] x where [•] m denotes the mth column of the matrix argument, δx = (N + λ)P n n|n , w m M and w m C are the weights of the mean and covariance, respectively, λ = α 2 (N + κ) − N , and α, β, and κ are the parameters of the unscented transform.
It is important to point out that the Kalman filter together with the model ( 12) easily accommodate asynchronous measurement updates.This is especially useful for the problem considered here as the RSS measurements arrive asynchronously due to the channel hopping nature of the communication protocol.
Also note that by using Rao-Blackwellization, the number of sigma-points is reduced to 3 as compared to using an unscented Kalman filter without Rao-Blackwellization which would require 2[C(2J + 1) + 1] + 1 sigma-points.

B. Preprocessing
Because of the limited dynamic range of low-cost, off-theshelf radio transceivers and the small amplitude of the RSS variations, significant quantization is sometimes observed.In particular, the RSS is often only reported in terms of integer values and the signal has to be preprocessed before the Kalman filter introduced above can be employed.A pragmatic solution to this is to recover the underlying signal (12b) by applying a low-pass filter to the quantized signal.For this purpose, a 5th order elliptical FIR low-pass filter with cutoff and stop-band frequencies f p = 1 Hz and f s = 1.2 Hz, respectively, passband ripple R p = 0.05 dB and stop-band attenuation R s = 40 dB has been used in this work.Strictly speaking, this does not exactly recover (12b) (e.g. the measurement noise will be correlated), but we found this to generally work well.

C. Initialization
A well initialized filter is key to achieve good performance.The crucial states here are the frequency f n (or rather its ) and the DC components γ 0,n .If the frequency is initialized poorly, there is a big risk that the filter will not converge to the correct frequency -which is a common challenge in recursive frequency estimation.Fortunately, the range of the breathing rate for healthy humans is relatively narrow, approximately between 12 bpm and 18 bpm (0.2 Hz-0.3 Hz) [18].Thus, the initial breathing rate is set to the mean of these two, that is, ν0|0 = log(15/60) with covariance such that one-sigma covers the whole range, that is, Cov{ν 0 } = 1 4 (log(18/60)−log(12/60)) 2 .However, a suitable initial breathing rate could also be determined from a small batch of initialization data.
Proper initialization of the DC component γ 0,n helps to speed up convergence.Since the amplitudes of the oscillations are generally relatively small compared to the mean RSS, a suitable initialization is simply the first RSS measurement y i,1 .
The remaining states (the coefficients of the higher orders) are not critical.Good initialization will help to speed up convergence, but initializing them to zero will not affect the performance significantly.Thus, these are set to zero.

A. Experiment Setup
The breathing rate of a person is monitored using two IEEE 802.15.4 compliant narrowband transceivers that are deployed ) and the spectrum estimation-based approach ( ) (middle) and both methods' absolute error (bottom) for 14 bpm.
on opposite sides of a bed, two meters apart from each other.In the experiments, the person is lying still in between the transceivers and breathing at a predefined rate.The transceivers communicate on C = 16 frequency channels covering a band between 2.405 GHz and 2.480 GHz and the sampling rate for each channel is on average 32 ms.The transceivers are equipped with low-resolution A/D converters and the RSS is measured with ±1 dB resolution.The experimental procedure follows that of experiment no. 1 in [7] and the reader is referred to that work for more details about the experimental setup.
In Section IV-B, the model parameters are set to S f = σ 2 = (1 × 10 −3 ) 2 , l = 0.1, J = 4, and R i = 0.25 2 .The filter parameters are α = 1, β = 0, and κ = 1.Furthermore, the prefilter discussed in the previous section is applied to each channel.The proposed method is benchmarked against a Kalman filter-based spectrum estimation method [19] which is applied for each radio channel individually and yields similar performance as the commonly used batch-based PSD estimation methods.If multiple channels are used, the estimated spectra are averaged before the person's breathing rate is estimated as done in [7].The breathing rate is estimated by taking the frequency of the highest peak in the range of 0.1 Hz to 1 Hz, excluding the strong DC term.The frequency grid spacing was chosen to be ∆f = 0.01 Hz.

B. Results and Discussion
Fig. 1 and Fig. 2 show two examples of the RSS signal together with the estimated breathing rate as well as the absolute value of the estimation error e i,n = fn − f n of the proposed method and the spectrum estimation-based method when using a single channel (channels 1 and 3, respectively).For channel 1, the RSS shows clear oscillations and both methods converge in approximately the same time to the true rate of 14 bpm. ) and the spectrum estimation-based approach ( ) (middle) and both methods' absolute error (bottom) for 14 bpm.
However, the spectrum-based method shows a jumpy behavior during convergence which is due to how the rate is estimated in this approach (selecting the peak frequency).For the proposed method, the estimate converges more smoothly.
For channel 3 (Fig. 2), the measured RSS exhibits much less clear oscillations of the base frequency.This causes problems in the spectrum-based method since the higher harmonics are stronger in this case.Using this approach, it is not straightforward to resolve this and the method does not converge.The GP-based method on the other hand copes with this since there is no requirement for the base frequency to exhibit the highest magnitude.Thus, the estimate converges as quickly as for channel 1.
Fig. 3 shows the estimated breathing rate and the absolute estimation error for the case when all channels are used simultaneously where the GP-based method simply performs asynchronous updates while the spectra for the different channels are averaged as described above.Again, both methods converge almost equally quickly.An important difference between the methods is that the spectrum-based approach actually exhibits a small and constant steady-state error.This is because of this method being bound to the discrete frequency bins which in this case are separated by 0.01 Hz.However, 14 bpm correspond to 0.233 Hz and thus, the closest bin is at 0.23 Hz.The GP method does not suffer from this limitation.
Finally, Fig. 4 illustrates the time-averaged root mean squared error (RMSE) e RMSE = 1 T n e 2 i,n for different breathing rates using all channels (similar figures could be shown for each of the 16 individual channels but that would not provide additional insight).During the transient period (0 s-30 s), the proposed method attains a significantly lower RMSE compared to the reference method.The main reason for this is the large error due to the jumpy behavior as observed in Fig. 3, which inflates the RMSE.It is also worth pointing out that the GP approach has higher RMSEs for 12 bpm, 18 bpm, and 20 bpm compared to 14 bpm and 16 bpm.This is due to the initialization of the breathing rate (15 bpm) which affects the convergence time and results in the higher RMSE for these cases.After convergence, the RMSE is much lower and of the same magnitude for both methods.The main factor contributing to the higher RMSE is how close the closest frequency bin is for the spectrum-based method, which determines the steady-state error.Note that for 12 bpm (0.2 Hz), the closest frequency bin coincides with the true breathing rate exactly and zero RMSE is achieved for the spectrum method.

V. CONCLUSIONS
In this paper, we presented a novel method for respiratory rate monitoring based on periodic Gaussian processes and Kalman filtering motivated by the variations in RSS measurements of narrowband radio transceivers.It has been shown that the method is capable of estimating the breathing rate recursively and optimally fusing the different channels, without requiring windowing or batch-processing.Furthermore, missing measurements do not pose a problem for the proposed method since there is no requirement of equispaced sampling.Despite the fact that the method has been introduced in the light of RSS measurements, it can readily be generalized to any type of measurements where the respiratory rate manifests in quasi-periodic signals (e.g.ultra-wideband transceivers or chest-mounted accelerometers).

Fig. 1 .
Fig. 1.Example of the measured RSS signal for channel i = 1 (top), the estimated breathing rates for the proposed method () and the spectrum estimation-based approach ( ) (middle) and both methods' absolute error (bottom) for 14 bpm.

Fig. 2 .
Fig. 2. Example of the measured RSS signal for channel i = 3 (top), the estimated breathing rates for the proposed method () and the spectrum estimation-based approach ( ) (middle) and both methods' absolute error (bottom) for 14 bpm.

Fig. 3 .Fig. 4 .
Fig. 3. Example of the estimated breathing rate (top) and the absolute errors (bottom) using all C = 16 channels for the proposed method () and the spectrum estimation-based approach ( ) for 14 bpm.