# Human Breathing Rate Estimation from Radar Returns Using Harmonically Related Filters

Authors: Mohamed Mabrouk; Sreeraman Rajan; Miodrag Bolic; Mohamad Forouzanfar; Hilmi R. Dajani; Izmail Batkin
Journal: Journal of Sensors
 (2016)

---

## Abstract

Radar-based noncontact sensing of life sign signals is often used in safety and rescue missions during disasters such as earthquakes and avalanches and for home care applications. The radar returns obtained from a human target contain the breathing frequency along with its strong higher harmonics depending on the target’s posture. As a consequence, well understood, computationally efficient, and the most popular traditional FFT-based estimators that rely only on the strongest peak for estimates of breathing rates may be inaccurate. The paper proposes a solution for correcting the estimation errors of such single peak-based algorithms. The proposed method is based on using harmonically related comb filters over a set of all possible breathing frequencies. The method is tested on three subjects for different postures, for different distances between the radar and the subject, and for two different radar platforms: PN-UWB and phase modulated-CW (PM-CW) radars. Simplified algorithms more suitable for real-time implementation have also been proposed and compared using accuracy and computational complexity. The proposed breathing rate estimation algorithms provide a reduction of about 81% and 80% in the mean absolute error of breathing rates in comparison to the traditional FFT-based methods using strongest peak detection, for PN-UWB and PM-CW radars, respectively.


---

## Body

## 1. Introduction and Motivation

Measurement of vital signals using ultra-wideband (UWB) and both phase modulated and frequency-modulated continuous wave (CW) radars is gaining popularity as a noncontact, noninvasive methodology [1–3]. In many applications, it is necessary to acquire and monitor a subject’s vital signs remotely and in a noncontact, noninvasive manner [4–6] as in the case of sleep apnea patients [7] and people trapped under collapsed buildings [8] and for fall detection [9, 10]. Also such noncontact, noninvasive vital sign measurement is useful for long-term monitoring of elderly people [9]. A noncontact method of measuring the human breathing signal [11] may be an economical way of monitoring the health of people [12, 13] and of tracking of subjects of interest, in the case of military operations [2].Breathing signals are used not only as a vital sign indicator but also in various applications ranging from tumor targeting during radiation treatment to speech monitoring of people [12, 13]. In almost all these applications, a reliable breathing frequency estimate needs to be obtained from the breathing signal. However, the return signal received at a radar is a superposition of signals due to the movements of the torso and limbs and expansion and contraction of the chest cavity in addition to noise and interference. Even for stationary people, detecting the contributions of the small displacement due to breathing in the radar returns is challenging due to strong clutter from objects in the vicinity of the stationary people. Physiological radar signal processing, therefore, requires the following steps to be performed: noise and clutter suppression, range estimation, classification of activities, and breathing and heartbeat signal extraction for stationary subjects.Traditional methods such as fast Fourier transform (FFT) and short time Fourier transform (STFT) have been used to detect the human breathing from radar returns in a specific range bin [2, 7, 8] and to estimate the human breathing rate [13–15], once the location of the target is determined. Recently, several FFT-based estimators were proposed for breathing rate estimation. These methods estimate the breathing rate using techniques such as detecting the peak of the FFT magnitude and the peak of the FFT magnitude squared [14, 16]. As the breathing signal is not invariant to motion [12, 13], these estimators perform poorly in the presence of movements. Several techniques to reduce the effects of movements on breathing rate estimation such as moving average method and ellipse fitting method are presented in the literature [13]. Clutter suppression in the presence of body movements during the detection of breathing activity through UWB radars was also attempted and then a Lomb periodogram was used to estimate the breathing rate [17]. Recently, a Multiple Signal Classification (MUSIC) algorithm was proposed for breathing frequency estimation [18]. This estimation method is computationally burdensome and its performance for different postures is unknown. Harmonics of the breathing frequency appear depending on the posture. A strong first harmonic of the breathing frequency impacts the breathing rate estimation as seen in [17]. However, the impact of posture on the breathing frequency has not been studied extensively. In this work, we present a methodology to estimate breathing frequency under different postures. Further, we will assume that the detection of activities has already been performed and the paper will address only the issue of breathing rate estimation from the radar returns.Most research in sensing of breathing signal using radar considers the breathing signal to be sinusoidal as the motion of only thorax is assumed. Recently, the displacement of the abdomen was observed to be on average four times that of the thorax during respiratory motion analysis [19]. The contribution of abdomen was incorporated through a two-cylinder model in [20] which led to a nonsinusoidal model for breathing signal in [21] and was further exploited in [22] to study the effects of posture on breathing signal. Furthermore, in [22], it was shown that the abdomen also contributes to the radar returns and to accommodate this contribution, a multiharmonic breathing model was considered. In the frequency domain of the radar returns in the model, the effects of phase modulations due to both chest and abdomen can be observed. As movements of abdomen may dominate over that of the chest in certain postures, especially while sitting [22], the higher harmonics of breathing may be stronger than the fundamental. Based on this fact, it is difficult to decide which one of the peaks in the spectrum of the breathing signal is the actual breathing rate because many harmonics may fall in the human breathing frequency range.As there are multiple sources, namely, the thorax and the abdomen, that contribute to the micro-Doppler in the radar return signal, we propose a harmonic filter-based approach to estimate the breathing rate accurately regardless of the posture. As we have recognized the similarity between the problem at hand and the pitch estimation problem in speech processing, we propose a maximum likelihood approach that almost mimics the harmonically related filtering approach proposed for pitch estimation [23]. Our proposed method achieves more accurate breathing rate estimation than the FFT-based techniques that rely only on a single strongest peak. Henceforth in this paper, we will refer to the FFT-based techniques that rely on a single strongest peak for breathing rate estimation as single strong peak-based FFT algorithms (SSP-FFT). A comprehensive set of experiments were done to analyze the sensitivity of the SSP-FFT to different postures, the improvement in breathing rate estimation using the proposed algorithm for different postures, and breathing rate estimation at different distances. Also, two more variants of the proposed algorithms that present a trade-off between computational complexity and accuracy are also presented and evaluated.Experiments were performed on two different types of radars: (1) PN-UWB with center frequency of 4.3 GHz, bandwidth of 2.2 GHz, and range bin of 0.91 cm and (2) phase modulated-CW (PM-CW) surveillance radar operating at 24.125 GHz with 500 MHz bandwidth and 3 m zone (range bin). These two radars have very different characteristics and require different preprocessing steps; however, breathing rate estimation is performed in the same way. Therefore, it was important to observe if multiple breathing harmonics appear in signals acquired from both radars and affect the breathing rate estimation of SSP-FFT and our proposed algorithms in the same way. Our results show that, for both radars, SSP-FFT algorithm often incorrectly estimates the breathing rate resulting in doubling the estimated breathing rate (the amplitude of the first harmonic is higher than the amplitude at the fundamental frequency) in sitting and standing positions or picking some other frequency as a breathing frequency. The proposed approach more accurately detects the fundamental breathing frequency for both radars with any posture.


## 2. Methodology

### 2.1. Experimental Setup

In this work, the reflected signals from three human subjects in various postures in an office environment using PN-UWB and PM-CW radars were obtained. Both these radars are monostatic as the transmitter and the receiver antennas are colocated from the perspective of the subjects. The advantages of PN-UWB over PM-CW radar are that it has high spatial resolution, good immunity against multipath interference, ability to reject clutter, and high spatial resolution and discriminates between targets close to one another. On the other hand, CW radars are less complex and cheap and require less preprocessing [24]. PN-UWB radar is a pulsed radar while the PM-CW radar is a nonpulsed radar. The intent of using two different radar technologies is to demonstrate that the proposed method is applicable to any type of radar. Table 1 presents the key parameters of the PN-UWB and PM-CW radars. The subjects, who participated in the experiments, were facing the transmitter and the receiver antennas. Postures that were used in this work were standing, sitting, and lying down at different ranges, 1 m, 2 m, and 3 m, from the transmitter and the receiver planes as shown in Figure 1. A total of 27 measurements were obtained with each radar. Subjects’ reference breathing rate was measured using a piezoelectric breathing belt produced by CleveMed™, Cleveland, USA. The subjects were asked to breathe when a tone generated from a software program was heard. The tone was repeated every 3.33 seconds.Table 1
Key parameters of PN-UWB and PM-CW radars used in this work.
Parameter
PN-UWB
CW
Centre frequency
4.3 GHz
24.125 GHz
Frequency band
2.2 GHz
500 MHz
Transmitted power
1 mW
17 mW
Pulse repetition rate
10.1 MHz
Not applicable
Dynamic range
75 dB
120 dBFigure 1
Experimental setup for one subject in three different postures at three different ranges using the PN-UWB radar.


#### 2.1.1. PN-UWB Radar Experimental Measurements

We used a PN-UWB (PulsON 410) radar of TIME DOMAIN™, Huntsville, USA, to study the breathing signal from 3 subjects in various postures and at varying distance from the radar. This radar provided the cross-correlation between the transmitted radar signal and the radar returns as the output. The sampled correlator output data, y(t,τ), was arranged as a two-dimensional matrix K×L [25], where each row represents correlator output at instants in slow time t=kTs, k=1,…,K, while the correlator output was sampled at discrete fast sample time τ=lTf, l=1,…,L, where K is the number of scans and L is the number of range bins. The algorithm in [25] was applied to estimate the target’s range bin which related to the fast time τ. The correlator output at the target range bin along the slow time direction is preprocessed by applying first-order high-pass filter in each range bin, followed by a bandpass filter in the range from 0.1 Hz to 1 Hz in each range bin. The output of the bandpass filter is x(t), whose Fourier transform, X(f), is the input to the proposed algorithms in this paper.


#### 2.1.2. Phase Modulated-CW (PM-CW) Radar Experimental Measurements

We used a millimeter PM-CW (SR4505) radar of K&G Spectrum, Gatineau, Quebec, Canada, to conduct measurements similar to those of the PN-UWB radar. The transmitted signal is biphase shift keying (BPSK) modulated signal with a 24.125 GHz center frequency and 500 MHz bandwidth. The receiver output is a sampled signalx(t). This signal is downsampled from 1750 Hz to 50 Hz sampling rate. A high-pass filter is applied to the downsampled signal to remove clutter. After that, a 6th-order IIR Butterworth bandpass filter in the range from 0.1 Hz to 1 Hz is applied to remove the frequencies out of the band of the human breathing range. Finally, the Fourier transform Xf is computed for this bandpass filtered signal x(t). It should be noted that the signals at the output of the preprocessing stages of both radars are labeled the same because they are used as the input to the same breathing rate estimation algorithm presented in Sections 2.2 and 2.3.


### 2.2. Proposed Harmonically Related Filter for Breathing Rate Estimation

One of the problems in detecting the breathing rate using SSP-FFT is that the first harmonic may have a magnitude higher than the magnitude at the fundamental frequency in some human postures. Figures2 and 3 present the time domain and the (normalized) FFT of a breathing signal, obtained from a human subject who was breathing at 0.3 Hz and sitting 1.5 m from the PN-UWB and PM-CW radars in an anechoic chamber at the University of Ottawa, Ontario, Canada, respectively. It is observed that the first harmonic has the maximum magnitude in both the figures. The SSP-FFT method fails to provide an accurate breathing rate estimation in such situations as SSP-FFT relies on the strongest peak for frequency estimation. In this section, a harmonically related filter based on the maximum likelihood approach is designed to address the problem of breathing rate estimation in the presence of stronger harmonics.Figure 2
The time domain and the discrete time Fourier transform of the breathing signal reflected of a human breathing with 0.3 Hz breathing rate and sitting 1.5 m from PN-UWB radar.Figure 3
The time domain and the discrete time Fourier transform of the breathing signal reflected of a human breathing with 0.3 Hz breathing rate and sitting 1.5 m from PM-CW radar.Consider an IIR comb filter with notches located at the normalized frequencies,θ,2θ,…,Nθ, with parameter r=0.99, where θ=2·fb/fs, where fb is the breathing frequency in Hz which lies in the range of θbl=0.2Hz/fs to θbh=2Hz/fs and fs is the sampling frequency in Hz, is applied to x(t). The parameter r controls the bandwidth of the filter. Values of r closer to 1 provide very narrow notches.The comb filter transfer function forN second-order sections is as shown in [26]:(1)Hθ,z=Yθ,zXz=∏m=1NHmθ,z,where(2)Hmθ,z=1-2z-1cos⁡mθ+z-21-2rz-1cos⁡mθ+r2z-2.The output of the comb filter in the frequency domain at a chosen breathing frequency, fb, is given by(3)Yθ,f=Hθ,fXf,where Hθ,f is the frequency response of the comb filter for a chosen set of N harmonic frequencies and X(f) is the Fourier transform of xt.It is clear thatYθ,f is minimum when the comb filter is adjusted to have N harmonic frequencies equal to the breathing frequency and its harmonics. The maximum likelihood breathing frequency estimation is obtained by(4)f^ML=fs2×arg⁡min︸θ∈θbl:θbh⁡∑fYθ,f2.In this work, the number of second-order filters in (4) was selected to be N=3 as third and higher harmonics were found to be small.Figures4 and 5 illustrate the normalized summed power of the breathing signal after applying the comb filter with r=0.99, at each θ, and the minimum power is detected at the fundamental harmonic using (4).Figure 4
The power of the breathing signal illustrated in Figure2 after applying the comb filter and computing the power at every frequency fs×θ/2.Figure 5
The power of the breathing signal illustrated in Figure3 after applying the comb filter and computing the power at every frequency fs×θ/2.


### 2.3. Simplification of the Algorithm for Breathing Rate Estimation

Execution time of the algorithm presented using (3) and (4) is long as it requires computing Fourier transform of the signal x(t) and then multiplying the spectrum of the signal K times with Hθ,f, where K is the number of selected frequency points of interest in the range [fs×θbl/2,fs×θbh/2]. In addition, sum of squares of each of K spectra Yθ,f needs to be calculated as well as the minimum of all the points in (4).Modified algorithms are introduced below to address the issue of computational complexity.


#### 2.3.1. Modified Algorithm  1 (FFT + Moving Windows)

Three narrow windows centered on the fundamental frequency and the first two harmonics of the fundamental breathing frequency were chosen. Then, the magnitude value at the fundamental frequency was replaced by the summation of the magnitudes after applying the three windows. Finally, the three windows were slid on every possible human breathing frequency (0.1~1 Hz) and the magnitude values were updated at these frequencies. Using this method, it can be expected that the magnitude at the fundamental breathing frequency will be higher than the magnitude of the first harmonic.


#### 2.3.2. Modified Algorithm  2 (FFT + Sum Squares)

A more simplified solution based on the following formula was used for extracting the breathing frequency:(5)f^ML=fs2×arg⁡max︸f∈fbl:fbh⁡∑k=13Xkf2.This algorithm requires only computation of Fourier transform of the signal xt and then the computation of the sum of squares of the absolute values at the fundamental and the first two harmonics. Therefore, this algorithm does not significantly increase the complexity in comparison with the original FFT-based algorithm. Here, fbl and fbh present minimum and maximum frequency of interest.


## 2.1. Experimental Setup

In this work, the reflected signals from three human subjects in various postures in an office environment using PN-UWB and PM-CW radars were obtained. Both these radars are monostatic as the transmitter and the receiver antennas are colocated from the perspective of the subjects. The advantages of PN-UWB over PM-CW radar are that it has high spatial resolution, good immunity against multipath interference, ability to reject clutter, and high spatial resolution and discriminates between targets close to one another. On the other hand, CW radars are less complex and cheap and require less preprocessing [24]. PN-UWB radar is a pulsed radar while the PM-CW radar is a nonpulsed radar. The intent of using two different radar technologies is to demonstrate that the proposed method is applicable to any type of radar. Table 1 presents the key parameters of the PN-UWB and PM-CW radars. The subjects, who participated in the experiments, were facing the transmitter and the receiver antennas. Postures that were used in this work were standing, sitting, and lying down at different ranges, 1 m, 2 m, and 3 m, from the transmitter and the receiver planes as shown in Figure 1. A total of 27 measurements were obtained with each radar. Subjects’ reference breathing rate was measured using a piezoelectric breathing belt produced by CleveMed™, Cleveland, USA. The subjects were asked to breathe when a tone generated from a software program was heard. The tone was repeated every 3.33 seconds.Table 1
Key parameters of PN-UWB and PM-CW radars used in this work.
Parameter
PN-UWB
CW
Centre frequency
4.3 GHz
24.125 GHz
Frequency band
2.2 GHz
500 MHz
Transmitted power
1 mW
17 mW
Pulse repetition rate
10.1 MHz
Not applicable
Dynamic range
75 dB
120 dBFigure 1
Experimental setup for one subject in three different postures at three different ranges using the PN-UWB radar.


### 2.1.1. PN-UWB Radar Experimental Measurements

We used a PN-UWB (PulsON 410) radar of TIME DOMAIN™, Huntsville, USA, to study the breathing signal from 3 subjects in various postures and at varying distance from the radar. This radar provided the cross-correlation between the transmitted radar signal and the radar returns as the output. The sampled correlator output data, y(t,τ), was arranged as a two-dimensional matrix K×L [25], where each row represents correlator output at instants in slow time t=kTs, k=1,…,K, while the correlator output was sampled at discrete fast sample time τ=lTf, l=1,…,L, where K is the number of scans and L is the number of range bins. The algorithm in [25] was applied to estimate the target’s range bin which related to the fast time τ. The correlator output at the target range bin along the slow time direction is preprocessed by applying first-order high-pass filter in each range bin, followed by a bandpass filter in the range from 0.1 Hz to 1 Hz in each range bin. The output of the bandpass filter is x(t), whose Fourier transform, X(f), is the input to the proposed algorithms in this paper.


### 2.1.2. Phase Modulated-CW (PM-CW) Radar Experimental Measurements

We used a millimeter PM-CW (SR4505) radar of K&G Spectrum, Gatineau, Quebec, Canada, to conduct measurements similar to those of the PN-UWB radar. The transmitted signal is biphase shift keying (BPSK) modulated signal with a 24.125 GHz center frequency and 500 MHz bandwidth. The receiver output is a sampled signalx(t). This signal is downsampled from 1750 Hz to 50 Hz sampling rate. A high-pass filter is applied to the downsampled signal to remove clutter. After that, a 6th-order IIR Butterworth bandpass filter in the range from 0.1 Hz to 1 Hz is applied to remove the frequencies out of the band of the human breathing range. Finally, the Fourier transform Xf is computed for this bandpass filtered signal x(t). It should be noted that the signals at the output of the preprocessing stages of both radars are labeled the same because they are used as the input to the same breathing rate estimation algorithm presented in Sections 2.2 and 2.3.


## 2.1.1. PN-UWB Radar Experimental Measurements

We used a PN-UWB (PulsON 410) radar of TIME DOMAIN™, Huntsville, USA, to study the breathing signal from 3 subjects in various postures and at varying distance from the radar. This radar provided the cross-correlation between the transmitted radar signal and the radar returns as the output. The sampled correlator output data, y(t,τ), was arranged as a two-dimensional matrix K×L [25], where each row represents correlator output at instants in slow time t=kTs, k=1,…,K, while the correlator output was sampled at discrete fast sample time τ=lTf, l=1,…,L, where K is the number of scans and L is the number of range bins. The algorithm in [25] was applied to estimate the target’s range bin which related to the fast time τ. The correlator output at the target range bin along the slow time direction is preprocessed by applying first-order high-pass filter in each range bin, followed by a bandpass filter in the range from 0.1 Hz to 1 Hz in each range bin. The output of the bandpass filter is x(t), whose Fourier transform, X(f), is the input to the proposed algorithms in this paper.


## 2.1.2. Phase Modulated-CW (PM-CW) Radar Experimental Measurements

We used a millimeter PM-CW (SR4505) radar of K&G Spectrum, Gatineau, Quebec, Canada, to conduct measurements similar to those of the PN-UWB radar. The transmitted signal is biphase shift keying (BPSK) modulated signal with a 24.125 GHz center frequency and 500 MHz bandwidth. The receiver output is a sampled signalx(t). This signal is downsampled from 1750 Hz to 50 Hz sampling rate. A high-pass filter is applied to the downsampled signal to remove clutter. After that, a 6th-order IIR Butterworth bandpass filter in the range from 0.1 Hz to 1 Hz is applied to remove the frequencies out of the band of the human breathing range. Finally, the Fourier transform Xf is computed for this bandpass filtered signal x(t). It should be noted that the signals at the output of the preprocessing stages of both radars are labeled the same because they are used as the input to the same breathing rate estimation algorithm presented in Sections 2.2 and 2.3.


## 2.2. Proposed Harmonically Related Filter for Breathing Rate Estimation

One of the problems in detecting the breathing rate using SSP-FFT is that the first harmonic may have a magnitude higher than the magnitude at the fundamental frequency in some human postures. Figures2 and 3 present the time domain and the (normalized) FFT of a breathing signal, obtained from a human subject who was breathing at 0.3 Hz and sitting 1.5 m from the PN-UWB and PM-CW radars in an anechoic chamber at the University of Ottawa, Ontario, Canada, respectively. It is observed that the first harmonic has the maximum magnitude in both the figures. The SSP-FFT method fails to provide an accurate breathing rate estimation in such situations as SSP-FFT relies on the strongest peak for frequency estimation. In this section, a harmonically related filter based on the maximum likelihood approach is designed to address the problem of breathing rate estimation in the presence of stronger harmonics.Figure 2
The time domain and the discrete time Fourier transform of the breathing signal reflected of a human breathing with 0.3 Hz breathing rate and sitting 1.5 m from PN-UWB radar.Figure 3
The time domain and the discrete time Fourier transform of the breathing signal reflected of a human breathing with 0.3 Hz breathing rate and sitting 1.5 m from PM-CW radar.Consider an IIR comb filter with notches located at the normalized frequencies,θ,2θ,…,Nθ, with parameter r=0.99, where θ=2·fb/fs, where fb is the breathing frequency in Hz which lies in the range of θbl=0.2Hz/fs to θbh=2Hz/fs and fs is the sampling frequency in Hz, is applied to x(t). The parameter r controls the bandwidth of the filter. Values of r closer to 1 provide very narrow notches.The comb filter transfer function forN second-order sections is as shown in [26]:(1)Hθ,z=Yθ,zXz=∏m=1NHmθ,z,where(2)Hmθ,z=1-2z-1cos⁡mθ+z-21-2rz-1cos⁡mθ+r2z-2.The output of the comb filter in the frequency domain at a chosen breathing frequency, fb, is given by(3)Yθ,f=Hθ,fXf,where Hθ,f is the frequency response of the comb filter for a chosen set of N harmonic frequencies and X(f) is the Fourier transform of xt.It is clear thatYθ,f is minimum when the comb filter is adjusted to have N harmonic frequencies equal to the breathing frequency and its harmonics. The maximum likelihood breathing frequency estimation is obtained by(4)f^ML=fs2×arg⁡min︸θ∈θbl:θbh⁡∑fYθ,f2.In this work, the number of second-order filters in (4) was selected to be N=3 as third and higher harmonics were found to be small.Figures4 and 5 illustrate the normalized summed power of the breathing signal after applying the comb filter with r=0.99, at each θ, and the minimum power is detected at the fundamental harmonic using (4).Figure 4
The power of the breathing signal illustrated in Figure2 after applying the comb filter and computing the power at every frequency fs×θ/2.Figure 5
The power of the breathing signal illustrated in Figure3 after applying the comb filter and computing the power at every frequency fs×θ/2.


## 2.3. Simplification of the Algorithm for Breathing Rate Estimation

Execution time of the algorithm presented using (3) and (4) is long as it requires computing Fourier transform of the signal x(t) and then multiplying the spectrum of the signal K times with Hθ,f, where K is the number of selected frequency points of interest in the range [fs×θbl/2,fs×θbh/2]. In addition, sum of squares of each of K spectra Yθ,f needs to be calculated as well as the minimum of all the points in (4).Modified algorithms are introduced below to address the issue of computational complexity.


### 2.3.1. Modified Algorithm  1 (FFT + Moving Windows)

Three narrow windows centered on the fundamental frequency and the first two harmonics of the fundamental breathing frequency were chosen. Then, the magnitude value at the fundamental frequency was replaced by the summation of the magnitudes after applying the three windows. Finally, the three windows were slid on every possible human breathing frequency (0.1~1 Hz) and the magnitude values were updated at these frequencies. Using this method, it can be expected that the magnitude at the fundamental breathing frequency will be higher than the magnitude of the first harmonic.


### 2.3.2. Modified Algorithm  2 (FFT + Sum Squares)

A more simplified solution based on the following formula was used for extracting the breathing frequency:(5)f^ML=fs2×arg⁡max︸f∈fbl:fbh⁡∑k=13Xkf2.This algorithm requires only computation of Fourier transform of the signal xt and then the computation of the sum of squares of the absolute values at the fundamental and the first two harmonics. Therefore, this algorithm does not significantly increase the complexity in comparison with the original FFT-based algorithm. Here, fbl and fbh present minimum and maximum frequency of interest.


## 2.3.1. Modified Algorithm  1 (FFT + Moving Windows)

Three narrow windows centered on the fundamental frequency and the first two harmonics of the fundamental breathing frequency were chosen. Then, the magnitude value at the fundamental frequency was replaced by the summation of the magnitudes after applying the three windows. Finally, the three windows were slid on every possible human breathing frequency (0.1~1 Hz) and the magnitude values were updated at these frequencies. Using this method, it can be expected that the magnitude at the fundamental breathing frequency will be higher than the magnitude of the first harmonic.


## 2.3.2. Modified Algorithm  2 (FFT + Sum Squares)

A more simplified solution based on the following formula was used for extracting the breathing frequency:(5)f^ML=fs2×arg⁡max︸f∈fbl:fbh⁡∑k=13Xkf2.This algorithm requires only computation of Fourier transform of the signal xt and then the computation of the sum of squares of the absolute values at the fundamental and the first two harmonics. Therefore, this algorithm does not significantly increase the complexity in comparison with the original FFT-based algorithm. Here, fbl and fbh present minimum and maximum frequency of interest.


## 3. Experimental Results

This work is compared against the FFT-based breathing rate estimator (SSP-FFT) used in [13, 14] to estimate the breathing rate of a target breathing at a known distance from PN-UWB and PM-CW radars. Breathing rate estimation performance of several algorithms has been compared including the algorithm given by (3) and (4) that we refer to as FFT + comb, the algorithm in Section 2.3.1 that is based on three moving windows that we refer to as FFT + moving windows, and the algorithm presented in (5) that we refer to as FFT + sum squares. The mean absolute error was used as the evaluation metric. In addition, the execution time of each algorithm without including preprocessing steps is discussed below.An FFT with number of frequency points equal to 216 was used for the comparison between the proposed estimator and the SSP-FFT estimator. In FFT + comb algorithm, the comb filter −3 dB bandwidth was fixed at 0.003 Hz. The first breathing harmonic had a higher magnitude than the fundamental harmonic in ten and thirteen measurements with the PN-UWB and PM-CW radars, respectively, where the SSP-FFT estimator failed to obtain the correct breathing frequency estimate.For FFT + moving windows estimator, it was found that a moving average window length equal to 1% of the breathing frequency (5 points in the case of breathing rate ~0.3 Hz) is enough to eliminate the effect of the first harmonic in the human breathing spectrum in most cases. The error in the estimated breathing rate increased with longer windows and lower SNR. Different types of windows (rectangular, Gaussian, Hann, and Hamming windows) were applied and it was found that the rectangular window outperformed the other widows.Figure6 shows the mean absolute error and standard deviation of all proposed algorithms studied in this paper for different postures. As expected, larger errors are observed in standing and sitting posture than in the lying posture for the traditional SSP-FFT algorithm because the magnitude at the first breathing harmonic is often higher than the magnitude at the fundamental frequency at sitting and standing postures. Standard deviation is very large because the SSP-FFT algorithm is making large errors, sometimes even 100% when the first harmonic is estimated as the breathing frequency. The proposed algorithms have lower errors. As seen in Figure 3, the SNR of PM-CW radar is much lower and its range bin (often referred to as zone since CW radars, in general, have poor range resolution) is very large (3 m) so that the radar picks up some other signals (nonbreathing signals) and the algorithms make significant errors when the magnitude of the breathing signal is much smaller than the magnitude of these nonbreathing signals. This accounts for worse results for PM-CW radar than for PN-UWB radar.Figure 6
The mean absolute estimation error in Hz of the breathing frequency of the SSP-FFT (blue) and the proposed methods over three different ranges for (a) PN-UWB radar and (b) PM-CW radar for 3 subjects in lying, sitting, and standing postures.
(a)
(b)As can be seen in Figure6, posture contributes to the error and our proposed algorithms are a step forward towards reducing the error. In addition, it can be concluded that the SSP-FFT estimator is more accurate when using PN-UWB radar than PM-CW radar (0.16 Hz versus 0.33 Hz mean absolute error for all postures and distances). The error in breathing rate is higher in standing than in sitting position when SSP-FFT estimator is applied on the data collected using PM-CW radar. As it is not possible to stand perfectly still, standing posture produces the largest error in the case of PM-CW radar. In general, PM-CW radars are more sensitive to small movements than PN-UWB radars. The error in the breathing estimates is the largest for the sitting position when the proposed estimators are applied to the data collected using the PN-UWB radar.While using a PM-CW radar, as shown in Table2, an overall error reduction versus conventional SSP-FFT algorithm of about 80% can be expected and with PN-UWB radar, an error overall reduction of 81% can be expected by using the proposed methodology.Table 2
Error reduction of our proposed methods versus the conventional SSP-FFT algorithm for all distances and postures using PN-UWB and PM-CW radars.
Algorithm/radar
PN-UWB
PM-CW
FFT + comb
83.5%
78.8%
FFT + moving windows
80.3%
81.7%
FFT  + sum squares
80.3%
79.2%Execution time of the algorithms is measured in MATLAB. Please note that MATLAB functions were not modified to minimize the execution time. In order to compare execution times, we did not take into account the preprocessing steps described in Sections2.1.1 and 2.1.2. The FFT + comb algorithm is very slow as it computes FFT many times. The total execution times for the methods are in relative terms: SSP-FFT: 1, FFT + comb: 332, FFT + moving windows: 65, and FFT + sum squares: 8. Even though the execution time of the FFT + sum squares is 8 times slower than the execution time of FFT itself, this algorithm (as well as FFT + moving windows) will allow for very efficient parallel implementation because there are no dependencies in the processing step (after FFT). Therefore, we believe that the conventional SSP-FFT-based algorithms need to be augmented with the proposed solution in order to improve the estimation while still supporting parallel processing and the efficient implementation.


## 4. Conclusion

The movements of the chest and the abdomen while breathing result in multiple frequency harmonics in the human breathing signal received by PN-UWB and PM-CW radars. The main sources of error in the human breathing rate estimators are the first breathing harmonic and any human large-scale displacements. The proposed estimator using harmonically related comb filters and its simplified versions are able to eliminate the effect of the higher harmonics of the breathing signal, thus making the estimation of breathing frequency from the radar return measurements more robust. The proposed method can be also applied to breathing estimation of human targets hidden behind obstacles.Future work involves performing new experiments that will include more subjects as well as subjects at different orientations relative to the radars. The fact that the magnitude at the first harmonic is larger than the magnitude of the fundamental frequency with some postures can be utilized for classification of postures. Our future work will be directed towards researching the robustness of classification algorithms for posture estimation based on the estimated breathing signal.



---

*Source: 9891852-2016-02-03.xml*