High-Precision Vital Signs Monitoring Method Using a FMCW Millimeter-Wave Sensor

Abstract

The method of using millimeter-wave radar sensors to detect human vital signs, namely respiration and heart rate, has received widespread attention in non-contact monitoring. These sensors are compact, lightweight, and able to sense and detect various scenarios. However, it still faces serious problems of noisy interference in hardware, which leads to a low signal-to-noise ratio (SNR). We used a frequency-modulated continuous wave (FMCW) radar sensor operating at 77 GHz in an office environment to extract the respiration and heart rate of a person accustomed to sitting in a chair. Indeed, the proposed signal processing includes novel impulse denoising operations and the spectral estimation decision method, which are unique in terms of noise reduction and accuracy improvement. In addition, the proposed method provides high-quality, repeatable respiration and heart rates with relative errors of 1.33% and 1.96% on average compared with the reference values measured by a reliable smart bracelet.

1. Introduction

Using the millimeter wave Doppler sensor [1] to non-contact monitor the vital signs of human targets [1,2,3] can be applied to many application scenarios where contact methods are not applicable. These scenarios include special patient monitoring in medical care [4], such as burn patients [5], newborns [6], infectious disease patients, and early warning of emergencies such as home healthcare for the elderly [7] and the avoidance of driver fatigue [8], etc. Doppler radar sensors can detect the movement of human targets by measuring the low Doppler variation that affects the received backscattered signals [1]. They are focused on two different systems: ultra-wide band (UWB) radars and continuous wave (CW) radars. Compared with the UWB radar, the FMCW radar relies on phase information to obtain target displacement information with higher detection accuracy. In addition, its non-broadband characteristic achieves low power consumption.

At present, increasing the SNR level of the signal to improve the accuracy of the respiration and heart rate is a research focus of numerous researchers. Multiple channels of the MIMO radar were adopted to collect data, as covered in paper [9], and the SNR has been improved by channel diversity. This puts forward requirements for the hardware system. At the same time, the maximum ratio combination algorithm based on adaptive continuous wavelet transform was proposed to improve the precision of heart rate detection. An adaptive noise cancellation technique combined with the polynomial fitting was proposed by the authors of [10], to offset the influence of random human movement on the extraction of small vital-sign signals. Meanwhile, a new discrete cosine transform was introduced to improve detection accuracy. H.-I. Choi et al. defined the time phase coherence index as a new distance bin selection method [11]. By analyzing the minute phase shift of each range bin, they selected the range box containing the accurate heartbeat signal to extract it and then improved the estimate accuracy. A method of reconstructing heartbeat signals by their harmonics and using the Kalman filter to improve the quality of the heartbeat waveform further was proposed by the authors of [12], which avoided the influence of respiratory harmonics on heartbeat signals to a certain extent. M. Arsalan et al. introduced the Kalman filter to track the heartbeat signal. They used an adaptive band-pass filter to extract the heart rate [13], which continuously reduced the bandwidth of the band-pass filter and approximated its final center frequency as the actual heart rate to improve its accuracy.

This paper demonstrates the feasibility of measuring human vital-sign signals using a 77 GHz FMCW sensor. We extracted the respiratory and heart rate using the phase information of the IF signal, a mixture of transmitted and received waves. Then, an algorithm called Iterative VMD Wavelet-Interval-Thresholding was proposed to reduce the noise in the measurement process. The experimental results show that the proposed algorithm can improve the SNR of respiration and heartbeat signals by 1.89 dB and 1.44 dB. Section 2 of this paper describes the theoretical model of vital signs signals from the FMCW radar’s perspective. It shows the flowchart of the whole signal process, focusing on the principle and steps of the noise reduction algorithm and the rate estimation algorithm. The experimental setups are presented in Section 3, and the verification by experimental results is drawn in Section 4. Finally, Section 5 provides the discussion and the conclusion is provided in Section 6.

2. Principles

2.1. Signal Model of FMCW Radar

We analyzed the principle of the FMCW radar detecting human vital-sign signals in our previous paper [14]. The modulated mode of the transmitted wave is a sawtooth wave with the slope as K and the chirp duration as ${\mathrm{T}}_{\mathrm{c}}$. The expression of the mixing signal output by channel I and channel Q at the receiving end of the radar can be expressed in theory as follows:

$$\mathrm{I}\left(\mathrm{t}\right)={\mathrm{A}}_{\mathrm{r}}\cos \left(\frac{4\mathsf{\pi }\mathrm{x}\left(\mathrm{t}\right)}{{\mathsf{\lambda }}_{\max }}+\frac{4{\mathsf{\pi }\mathrm{R}}_0}{{\mathsf{\lambda }}_{\max }}+2\mathsf{\pi }\mathrm{K}\mathsf{\tau }\mathrm{t}+\Delta \mathsf{\varphi }\right),$$

$$\mathrm{Q}\left(\mathrm{t}\right)={\mathrm{A}}_{\mathrm{r}}\sin \left(\frac{4\mathsf{\pi }\mathrm{x}\left(\mathrm{t}\right)}{{\mathsf{\lambda }}_{\max }}+\frac{4{\mathsf{\pi }\mathrm{R}}_0}{{\mathsf{\lambda }}_{\max }}+2\mathsf{\pi }\mathrm{K}\mathsf{\tau }\mathrm{t}+\Delta {\mathsf{\varphi }}^{\prime }\right),$$

(1)

where ${\mathrm{A}}_{\mathrm{r}}$ is the amplitude of the beat signal; $\mathsf{\tau }\approx 2{\mathrm{R}}_0/\mathrm{c}$ is the time delay caused by the distance between the measured target and the radar marked as ${\mathrm{R}}_0$; ${\mathsf{\lambda }}_{\max }=\mathrm{c}/{\mathrm{f}}_{\mathrm{c}}$ is the wavelength corresponding to the chirp starting frequency ${\mathrm{f}}_{\mathrm{c}}$; $\Delta \mathsf{\varphi }$ and $\Delta \mathsf{\varphi }\prime$ are the phase noise of channel I and Q, respectively; $\mathrm{x}\left(\mathrm{t}\right)$ is the slight time-varying displacement caused by cardiopulmonary activities containing vital signs such as breathing and heartbeat or other movements. The subsequent work is to sample and process the above signals.

2.2. The Proposed Signal Processing Algorithm Chain

Figure 1 shows the proposed signal-processing flowchart, where the processing signal is the IF signal output by the I and Q channels of the FMCW millimeter wave sensor. The process comprises four steps: signal preprocessing, phase extraction in range window, impulse noise removal, and vital signs detection.

After the IF signal is sampled by ADC, the data acquisition board exports the sampled data to form a two-dimensional matrix. The sampled data in each chirp constitute the matrix’s columns, and the matrix’s rows are composed of different frames. It should be noted that we set only one chirp in each frame. The sampled data are preprocessed mainly through static signal-clutter removal to eliminate the background interference in the environment for the position information of the measured target. The fast Fourier transform (FFT) is applied to the column vectors of the matrix successively to obtain the precise range window of the measured target, and then we further extract the phase of each frame in the range window. The vital signs of the target can be obtained from the phase information. To eliminate the influence of the DC components of the I and Q channels on phase extraction, the DC offset compensation operation is carried out to obtain the original phase. The phase difference operation plays a significant role in enhancing the heartbeat signal components. Then the iterative VMD Wavelet-Interval-Thresholding algorithm is proposed to filter the impulse noise generated by the differential process. In the vital signs detection stage, respiration and heart rates are extracted successively. Among them, the respiration rate is obtained by the FFT-CZT mixed spectrum estimation algorithm, while the heart rate is estimated by FFT-CZT and peak-seeking in the time domain method separately. Then we combine the confidence degree for thresholding judgment, and one of the results obtained by the two methods above is selected as the final result. The details of each step are elaborated in the sections as follows.

2.2.1. Range FFT and Static Signal-Clutter Removal

The 1D-FFT in the distance dimension is performed on the two-dimensional matrix formed by sampling the IF signal to obtain the range slow-time matrix. Due to the complex background in the actual detection environment, the reflected static clutter signals will affect the selection of the target range bin and introduce some phase noise, which is not conducive to the extraction of the cardiopulmonary information. Therefore, the interference of static clutter should be removed. In the same test scenario, the clutter signals caused by other stationary objects basically remain unchanged with slow time, so subtracting the mean value of the slow-time dimension of each range window from the data in the slow-time dimension corresponding to each range bin can remove the background noise and retain the movement information of the measured target. Assuming that $y_0$ is the two-dimensional matrix composed of the original collected data, the matrix after the removal of static clutter signals can be calculated as follows:

$$y\left[\mathrm{m},\mathrm{n}\right]=y_0\left[\mathrm{m},\mathrm{n}\right]-\frac{1}{{\mathrm{N}}_{\mathrm{frames}}}{\displaystyle \sum }_{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{frames}}}{y}_0\left[\mathrm{m},\mathrm{n}\right],$$

(2)

where m = 1, 2, 3, …, ${\mathrm{N}}_{\mathrm{samples}}$; n = 1, 2, 3, …, ${\mathrm{N}}_{\mathrm{frames}}$. ${\mathrm{N}}_{\mathrm{samples}}$ means the number of sampling points of each chirp; and ${\mathrm{N}}_{\mathrm{frames}}$ means the number of frames. Figure 2 is the schematic diagram of the distance slow-time matrix before and after the static signal-clutter removal operation. Figure 2b shows the range bin where the measured target is located, which will be further processed later to extract the phase.

2.2.2. DC Offset Compensation

Because of the limitations of the manufacturing process of the hardware circuit in the receiver, the sensor used may have problems of insufficient mixer isolation or incomplete matching between the receiving antenna and the transmitting antenna. Therefore, a certain DC offset still exists after solving the static signal clutter caused by the fixed objects in the detection environment. DC offset will introduce more phase noise, resulting in the phase’s distortion extraction of the echo signal by the arctangent operation. The extracted phase can be expressed as:

$$\mathsf{\Psi }=\arctan \left[\frac{\mathrm{Q}\left(\mathrm{t}\right)+{\mathrm{DC}}_{\mathrm{Q}}}{\mathrm{I}\left(\mathrm{t}\right)+{\mathrm{DC}}_{\mathrm{I}}}\right]\ne \arctan \left[\frac{\mathrm{Q}\left(\mathrm{t}\right)}{\mathrm{I}\left(\mathrm{t}\right)}\right],$$

(3)

where ${\mathrm{DC}}_{\mathrm{I}}$, ${\mathrm{DC}}_{\mathrm{Q}}$ are the real and imaginary parts of the DC offset, respectively. In the complex plane, because of the existence of the DC offset component, the constellation center of the complex signal in the selected target range bin moves from the origin to (${\mathrm{DC}}_{\mathrm{I}}$, ${\mathrm{DC}}_{\mathrm{Q}}$), as shown in Figure 3a. In order to shift the constellation to the origin, we use the non-linear least square estimation (NLLS) proposed by the authors of [15] to estimate the center and radius of the constellation. The DC offset compensation can be realized by successively subtracting the estimated center point position from each sampling point corresponding to the slow-time dimension in the selected range bin. Figure 3b shows the results after compensation.

2.2.3. Phase Extraction and Difference Operation

After DC offset compensation, we need to extract the phase information of the complex data. Using the arctangent algorithm to extract the phase will have a breakpoint problem, leading to phase drift exceeding ±π. Assuming that the real and imaginary parts of the complex data above are expressed as I(n) and Q(n), respectively, we apply the extended differentiate and cross multiply (DACM) algorithm [16] to extract the phase, to avoid the phase drifts mentioned above, as follows:

$$\mathsf{\Psi }\left(\mathrm{n}\right)={\displaystyle \sum }_{\mathrm{i}=2}^{\mathrm{n}}\frac{\mathrm{I}\left(\mathrm{i}\right)·\left[\mathrm{Q}\left(\mathrm{i}\right)-\mathrm{Q}\left(\mathrm{i}-1\right)\right]-\mathrm{Q}\left(\mathrm{i}\right)·\left[\mathrm{I}\left(\mathrm{i}\right)-\mathrm{I}\left(\mathrm{i}-1\right)\right]}{\mathrm{I}{\left(\mathrm{i}\right)}^2+\mathrm{Q}{\left(\mathrm{i}\right)}^2}.$$

(4)

Figure 4 shows an example of extracting the phase of a complex signal using the extended DACM algorithm. We transform the extracted phase into the waveform of chest vibration amplitude varying with slow time and obtain its spectrum information through FFT. We found that the heartbeat frequency is covered in the spectrum in Figure 4b because the chest fluctuation caused by the heartbeat is very weak relative to breathing. The first-order phase difference operation is adopted to enhance the heartbeat signal. Figure 5 shows the time domain waveform and frequency domain spectrum of the differential signal, where the spectral amplitude in the heart rate range is increased. The generated impulse components are removed in the following steps.

2.2.4. Iterative VMD Wavelet-Interval-Thresholding

Iterative VMD Wavelet-Interval-Thresholding is mainly used to remove the impulse noise generated after the phase difference operation. This process refers to the EMD-CIIT described in paper [17]. Compared with the empirical mode decomposition (EMD) algorithm, the variational mode decomposition (VMD) algorithm also decomposes the original signal into a series of signal components, which is called the intrinsic mode function (IMF). However, the latter is supported by sufficient mathematical theory [18]. Furthermore, the defect of modal aliasing is avoided by decomposing the signal into a certain number of discrete signal components, in which each IMF has different center frequencies and limited bandwidth. Assuming that VMD decomposes the signal into K IMFs, the kth IMF is recorded as $u_k\left(t\right)$, and its center frequency is ${\mathsf{\omega }}_{\mathrm{k}}$; the goal of the algorithm is to minimize the sum of the estimated bandwidth of each mode, that is:

$$\underset{\left\{u_k\right\},\left\{{\mathsf{\omega }}_{\mathrm{k}}\right\}}{\min }\left\{{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{K}}\parallel -{\mathrm{j}\mathsf{\omega }}_{\mathrm{k}}{\left[\left(\mathsf{\delta }\left(\mathrm{t}\right)+\frac{\mathrm{j}}{\mathsf{\pi }\mathrm{t}}\right)\ast u_k\left(t\right)\right]}^{\prime }{\mathrm{e}}^{-{\mathrm{j}\mathsf{\omega }}_{\mathrm{k}}\mathrm{t}}{\parallel }_2^2\right\},$$

(5)

where $\left\{u_k\right\}$ means $\left\{u_1,\dots ,u_K\right\}$ and $\left\{{\mathsf{\omega }}_{\mathrm{k}}\right\}$ means $\left\{{\mathsf{\omega }}_1,\dots ,{\mathsf{\omega }}_{\mathrm{K}}\right\}$. In order to solve this constrained variational problem, the quadratic penalty term $\mathsf{\alpha }$ and Lagrangian multiplier $\lambda \left(t\right)$ are introduced, which can transform the problem into an unconstrained variational problem as follows:

$$\begin{gathered} L\left(\left\{u_k\right\},\left\{{\mathsf{\omega }}_{\mathrm{k}}\right\},\lambda \left(t\right)\right)=\alpha {\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{K}}\parallel -{\mathrm{j}\mathsf{\omega }}_{\mathrm{k}}{\left[\left(\mathsf{\delta }\left(\mathrm{t}\right)+\frac{\mathrm{j}}{\mathsf{\pi }\mathrm{t}}\right)\ast u_k\left(t\right)\right]}^{\prime }{\mathrm{e}}^{-{\mathrm{j}\mathsf{\omega }}_{\mathrm{k}}\mathrm{t}}{\parallel }_2^2 \\ +\parallel f\left(t\right)-{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{K}}{u}_k\left(t\right){\parallel }_2^2+\langle \lambda \left(t\right),f\left(t\right)-{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{K}}{u}_k\left(t\right)\rangle . \end{gathered}$$

(6)

Then, $\left\{{\mathsf{\omega }}_{\mathrm{k}}\right\}$ and $\left\{u_k\right\}$ are obtained by the alternate direction method of multipliers (ADMM), which can be described in the frequency domain:

$$\begin{gathered} {\mathrm{U}}_{\mathrm{k}}^{\mathrm{n}+1}\left(\mathsf{\omega }\right)=\frac{\mathrm{F}\left(\mathsf{\omega }\right)-{{\displaystyle \sum }}_{\mathrm{i}\ne \mathrm{k}}{\mathrm{U}}_{\mathrm{i}}\left(\mathsf{\omega }\right)+\frac{ \hat{\mathsf{\lambda }}\left(\mathsf{\omega }\right)}{2}}{1+2\mathsf{\alpha }{\left(\mathsf{\omega }-{\mathsf{\omega }}_{\mathrm{k}}\right)}^2}, \\ {\mathsf{\omega }}_{\mathrm{k}}^{\mathrm{n}+1}=\frac{{{\displaystyle \int }}_0^{\infty }\mathsf{\omega }{\left|{\mathrm{U}}_{\mathrm{k}}\left(\mathsf{\omega }\right)\right|}^2\mathrm{d}\mathsf{\omega }}{{{\displaystyle \int }}_0^{\infty }{\left|{\mathrm{U}}_{\mathrm{k}}\left(\mathsf{\omega }\right)\right|}^2\mathrm{d}\mathsf{\omega }}, \end{gathered}$$

$$\hat{\mathsf{\lambda }}{}^{\mathrm{n}+1}\left(\mathsf{\omega }\right)= \hat{\mathsf{\lambda }}{}^{\mathrm{n}}\left(\mathsf{\omega }\right)+\mathsf{\tau }\left(\mathrm{F}\left(\mathsf{\omega }\right)-{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{K}}{\mathrm{U}}_{\mathrm{k}}^{\mathrm{n}+1}\left(\mathsf{\omega }\right)\right).$$

(7)

The algorithm can be divided into the following steps: first, initialize $\left\{{\mathrm{U}}_{\mathrm{k}}^1\right\}$, $\left\{{\mathsf{\omega }}_{\mathrm{k}}^1\right\}$, ${\mathsf{\lambda }}^1$, and n, which is the number of loops; then update $\left\{{\mathrm{U}}_{\mathrm{k}}\right\}$, $\left\{{\mathsf{\omega }}_{\mathrm{k}}\right\}$, and $\mathsf{\lambda }$ according to (7); repeat the previous step until the following constraints are satisfied:

$${\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{K}}\frac{{\left|\left|{\mathrm{U}}_{\mathrm{k}}^{\mathrm{n}+1}-{\mathrm{U}}_{\mathrm{k}}^{\mathrm{n}}\right|\right|}_2^2}{{\left|\left|{\mathrm{U}}_{\mathrm{k}}^{\mathrm{n}}\right|\right|}_2^2}<\mathsf{\epsilon }.$$

(8)

We set $\mathsf{\alpha }$ as 2000, and τ = 0, which means that the signal is restored without introducing additional noise; at the same time, we set $\mathsf{\epsilon }$ to ${10}^{-7}$, which refers to the convergence tolerance.

According to the above theoretical inference, VMD decomposes the original signal into the sum of multiple IMFs, where each IMF is a narrow-band amplitude and frequency modulation (AM-FM) signal. After comparing the energy and spectrum of each IMF with the original signal, we find that the first IMF basically does not contain useful signal components, so it is considered to contain noise components, and the remaining IMFs are regarded as signal components. The noise reduction process is shown in Figure 6, and the specific steps are as follows:

(1) VMD expansion: Perform the VMD decomposition of the original noisy signal y to obtain ${\mathrm{IMF}}_1\left(1\right)~{\mathrm{IMF}}_1\left({\mathrm{K}}_1\right)$, in which the number of decomposition times is ${\mathrm{K}}_1$. The selection of ${\mathrm{K}}_1$ refers to the optimal number of decomposition modes in the improved VMD algorithm [19], in which the cross-correlation number is used to represent the correlation degree between the noise component and the signal component. In order to facilitate subsequent operations, the first IMF is taken as the noise component and expressed as ${\mathrm{IMF}}_1\left(1\right)$, while the remaining components are reconstructed to obtain the signal component and expressed as ${\mathrm{IMF}}_1\left(\mathrm{h}\right)$. Then, the cross-correlation coefficient of the two can be expressed as:

$${\rho }_c=\frac{R_{1,h}\left(1\right)}{\sqrt{R_{1,1}\left(0\right)R_{h,h}\left(0\right)}}=\frac{{{\displaystyle \sum }}_{n=1}^N{\mathrm{IMF}}_1\left(\mathrm{n},1\right){\mathrm{IMF}}_1\left(\mathrm{n}-1,\mathrm{h}\right)}{\sqrt{R_{1,1}\left(0\right)R_{h,h}\left(0\right)}},$$

(9)

where $R_{1,1}\left(0\right)$ and $R_{h,h}\left(0\right)$ are the autocorrelation functions of noise component ${\mathrm{IMF}}_1\left(1\right)$ and signal component ${\mathrm{IMF}}_1\left(\mathrm{h}\right)$, respectively. $R_{1,h}\left(1\right)$ is the cross-correlation function of the two components mentioned above, in which the time difference is set as 1 and then the correlation is maximum. Next, calculate the cross-correlation coefficient for each decomposition number, ranging from 2 to 8. This range is because when the number is too large, a certain modal component may exist in multiple modal components, and then the spectrum fracture phenomenon will occur. Then, select the decomposition number corresponding to the smallest cross-correlation coefficient as ${\mathrm{K}}_1$.

(2) Wavelet-thresholding: The first IMF obtained, which is called ${\mathrm{IMF}}_1\left(1\right)$, is denoised by wavelet-thresholding. The noise reduction process is summarized in the following steps:

• Select the appropriate wavelet basis function and determine the layers of the wavelet decomposition;
• Obtain a set of wavelet decomposition coefficients in the wavelet domain by wavelet decomposition of ${\mathrm{IMF}}_1\left(1\right);$

• Set the threshold, and use the soft threshold function to process the wavelet coefficients.

The wavelet coefficients larger than the threshold are retained as the components caused by the target signal. In contrast, the wavelet coefficients smaller than the threshold are discarded as the components caused by the noise signal. Assuming that the threshold is called T and the coefficient is marked as $\mathrm{x}$, it can be described as:

$$\mathrm{f}\left(\mathrm{x}\right)=\left\{\begin{array}{ll}\mathrm{sgn}\left(\mathrm{x}\right)\left(\left|\mathrm{x}\right|-\mathrm{T}\right),& \left|\mathrm{x}\right|>\mathrm{T}\\ 0,& \left|\mathrm{x}\right|\le \mathrm{T}\end{array}\right..$$

(10)

• Reconstruct the signal by the processed wavelet coefficients to obtain the denoised signal ${\widetilde{\mathrm{IMF}}}_1\left(1\right)$.

We choose the Symlets wavelet function as the basis function because it meets the compactness and regularity with good symmetry, reducing the phase distortion caused by signal analysis and reconstruction to a certain extent. This algorithm selects the sym4 wavelet, whose support length is 7 and the number of vanishing moments is 4. At the same time, we set the layers of wavelet decomposition to 4, and the threshold is determined by using the ddencmp function in MATLAB. The filtered part is taken as the redefined noise component ${\mathrm{IMF}}_1^{\mathrm{n}}\left(1\right)$, and ${\widetilde{\mathrm{IMF}}}_1\left(1\right)$ is added to the remaining IMFs to form the redefined signal component, which is marked as $y_r$;

(3) Randomly alter: Considering the randomness of the noise contained, we chose to randomly alter the noise component ${\mathrm{IMF}}_1^{\mathrm{n}}\left(1\right)$ to obtain a different version of the noise component. It means that the sequence ${\mathrm{IMF}}_1^{\mathrm{n}}\left(1\right)$ is circularly shifted randomly, that is, the number of cycle-shift-right is randomly selected;

(4) Signal reconstruction: Add the noise component ${\mathrm{IMF}}_1^{\mathrm{n}}\left(1\right)$ obtained by step (3) and the redefined signal component $y_r$ obtained by step (2) to obtain the equivalent version of the original noisy signal, which is marked as $y_n$;

(5) Interval-thresholding: Perform the VMD decomposition of $y_n$ to obtain ${\mathrm{IMF}}_2\left(1\right)~{\mathrm{IMF}}_2\left({\mathrm{K}}_2\right)$, in which the number of decomposition times is ${\mathrm{K}}_2$. Then perform Interval-thresholding on them successively to obtain ${\widetilde{\mathrm{IMF}}}_2\left(1\right)~{\widetilde{\mathrm{IMF}}}_2\left({\mathrm{K}}_2\right)$, and add all of them to obtain the noise reduction version $y_i$. The following is the noise reduction basis of Interval-thresholding.

Considering the general characteristics of the IMF obtained by VMD decomposition, we define the waveform between adjacent zero crossings as a unit of the IMF. These characteristics include that the difference between the number of extreme points and zero crossings in the whole sequence must be no more than one, and the envelope mean value in any time range must be zero. Then, the threshold detection is carried out for each unit. If the amplitude of the extreme value in one unit is greater than the set threshold, the unit is considered a valuable signal to be reserved. In comparison, if the amplitude is less than the threshold, the unit is considered a noise signal and set to zero. It can be expressed as:

$$s^{\left(i\right)}\left(n_j\right)=\left\{\begin{array}{ll}{s}^{\left(i\right)}\left(n_j\right)\times \frac{\left|s^{\left(i\right)}\left(N_j\right)\right|-{\mathrm{T}}_{\mathrm{i}}}{\left|s^{\left(i\right)}\left(N_j\right)\right|},& \left|s^{\left(i\right)}\left(N_j\right)\right|>{\mathrm{T}}_{\mathrm{i}}\\ 0,& \left|s^{\left(i\right)}\left(N_j\right)\right|\le {\mathrm{T}}_{\mathrm{i}}\end{array}\right.,$$

(11)

where $s^{\left(i\right)}\left(n_j\right)$ represents the sequence which is consist of the instant values contained in the jth unit of the ith IMF; $s^{\left(i\right)}\left(N_j\right)$ is the extreme value of the sequence above; ${\mathrm{T}}_{\mathrm{i}}$ is the threshold of the ith IMF, which can be expressed as:

$${\mathrm{T}}_{\mathrm{i}}=\mathrm{C}\sqrt{2{\mathrm{E}}_{\mathrm{i}}\mathrm{lnN}},$$

(12)

where N is the length of the IMF component and C is a constant. Assuming ${\mathrm{E}}_1$ is the energy of the first IMF, the noise-only energy ${\mathrm{E}}_{\mathrm{i}}$ of the ith IMF can be computed directly:

$${\mathrm{E}}_{\mathrm{i}}=\frac{{\mathrm{E}}_1}{{\mathsf{\beta }\mathsf{\rho }}^{\mathrm{i}}},$$

(13)

where i = 2, 3,…, and $\mathsf{\beta }$, $\mathsf{\rho }$ are parameters that mainly depend on the number of iterations. There we set C = 0.025, $\mathsf{\beta }$ = 0.719, $\mathsf{\rho }$ = 2.01 [17] for A = 30, which is the times of iterations;

(6) Iterate: Iterate A-1 times steps (3)–(5) to obtain the corresponding noise reduction versions;

(7) Results of processing: Add and average the noise reduction versions generated after each iteration. Figure 7 shows the signal obtained by the Iterative VMD Wavelet-Interval-Thresholding algorithm to process the phase difference signal shown in Figure 5a. Figure 7b shows the spectrum in the frequency domain of the denoised signal. Compared with the spectrum shown in Figure 5b, the spectral energy of the noise part decreases significantly. Because of the randomly altering operation of the noise component, the noises offset each other when averaging the iteration results, which is also the reason for the reduction in the noisy spectrum amplitude.

2.2.5. Vital Signs Detection

The respiration rate is obtained by extracting its spectral peak value through the FFT-CZT hybrid algorithm. The algorithm mainly refers to our previous work [14], and refines the spectrum based on FFT to obtain more accurate local spectrum information. Since the amplitude of chest vibration caused by breathing is much larger than that caused by the heartbeat, the extraction of the heartbeat signal is more complex. The following are the specific steps.

(1) Filter out the respiratory signal through a band-pass filter: Considering the general heart-rate range, we select an elliptical band-pass filter with a passband of 0.7–3.3 Hz. The finite transmission zeros in the stopband of the elliptic filter reduce the transition region, and then an extremely steep attenuation curve can be obtained. Thus, it has a very narrow transition bandwidth;

(2) Filter out the respiratory harmonic interference through notch filter banks: The nonlinearity of the FMCW radar will lead to the generation of respiratory harmonics, in which the second to fourth harmonic components are within the range of humans’ heart rate, which will affect the extraction of the heart rate. Therefore, the notch filter banks designed in paper [14] are used to filter them out, and the notch frequency points are set to 2–4 times the extracted respiratory frequency above;

(3) Heart rate judgment: Due to the intermodulation between the breathing and heartbeat signals and the characteristic of multiple scattering points of humans’ chests, there are still frequency components other than heartbeat frequency after removing breathing and its harmonic components. There are two extraction methods:

• The FFT-CZT hybrid algorithm mentioned above
• Peek-seeking in the time domain

If ${\mathrm{N}}_0$ is a certain sampling time at a peak in the time domain, we propose the following conditions for its amplitude $\mathrm{P}\left({\mathrm{N}}_0\right)$:

$$\begin{gathered} \mathrm{P}\left({\mathrm{N}}_0\right)>\mathrm{P}\left({\mathrm{N}}_0-1\right)>0, \mathrm{P}\left({\mathrm{N}}_0\right)>\mathrm{P}\left({\mathrm{N}}_0+1\right)>0, \\ \mathrm{P}\left({\mathrm{N}}_0-1\right)>\mathrm{P}\left({\mathrm{N}}_0-2\right), \mathrm{P}\left({\mathrm{N}}_0+1\right)>\mathrm{P}\left({\mathrm{N}}_0+2\right), \end{gathered}$$

(14)

where $\mathrm{P}\left({\mathrm{N}}_0-1\right)$ represents the amplitude of the previous sampling point, $\mathrm{P}\left({\mathrm{N}}_0+1\right)$ represents the amplitude of the next sampling point, and so on. The limitation of the above conditions can avoid the statistics of some additional spikes, which will affect the frequency estimation.

The above two methods have their advantages and limitations. Combined with the confidence of the FFT method, the threshold decision scheme of heart-rate estimation is given: if the confidence of FFT exceeds the threshold and the difference between the results of the two methods above is less than the threshold, the FFT-CZT hybrid algorithm is used for heart-rate estimation; otherwise, the peak-seeking in the time-domain is adopted. Suppose that ${\mathrm{S}}_{\mathrm{FFT}}$ is the confidence of the FFT algorithm, referring to the definition of SNR, we describe this parameter as:

$${\mathrm{S}}_{\mathrm{FFT}}=\frac{{\mathrm{E}}_{\mathrm{peak}}}{{\mathrm{E}}_{\mathrm{signal}}-{\mathrm{E}}_{\mathrm{peak}}},$$

(15)

where ${\mathrm{E}}_{\mathrm{signal}}$ refers to the total energy of the spectrum obtained by the extracted heartbeat signal after FFT; ${\mathrm{E}}_{\mathrm{peak}}$ refers to the sum of the energy contained in the peak spectral line and its adjacent two spectral lines in the spectrum. We chose the three spectral lines mentioned above as the calculation range of the energy of the useful heartbeat signal. It is because the resolution of the FFT algorithm is limited, so the actual heart rate may not correspond to the frequency at the peak.

3. Experiments

We chose the radar AWR1642 from Texas Instruments (TI) for the experiment. It is an integrated FMCW millimeter-wave sensor operating at 76–81 GHz. It integrates a high-performance DSP subsystem and ARM processor for post-processing, which can realize programming and signal processing on the chip. However, in order to further observe the experimental results and analyze the data, we hope to obtain the original data for the algorithm verification. Therefore, we used the AWR1642 evaluation board to test the target and the DCA1000 acquisition board to collect the raw radar data. Then, we imported the collected data into MATLAB for signal analysis. Figure 8a shows the sensor device we used, and the test scenario of this experiment is shown in Figure 8b. The measured person sat in front of the radar sensor while wearing a PPG+ECG smart bracelet to capture the respiratory and heartbeat rates for reference. We first verified the reliability of the smart bracelet through five groups of measurement experiments; each group lasted for 1 min. It was proved that the respiration and heartbeat rate displayed by the bracelet was no more than one bpm different from the actual value. Thus, the measured value of the smart bracelet can be taken as the reference value.

The parameter configuration of the experiment is shown in

Table 1. FMCW signal parameters for vital signs measurements.

Parameter	Value
Chirp starting frequency, ${\mathrm{f}}_{\mathrm{c}}$	77 GHz
Chirp duration, ${\mathrm{T}}_{\mathrm{c}}$	50 μs
Chirp slope, K	80 MHz/μs
Idle time, ${\mathrm{T}}_{\mathrm{i}}$	7 μs
ADC sampling rate, ${\mathrm{F}}_{\mathrm{s}}$	6.4 Msps
Frame period, ${\mathrm{T}}_{\mathrm{s}}$	50 ms
Sample points number in each chirp, ${\mathrm{N}}_{\mathrm{samples}}$	256
Transmitting antenna	${\mathrm{Tx}}_0$
Receiving antenna	${\mathrm{Rx}}_0$
Output power	12.5 dBm

, in which the parameters are described in Figure 9. Among them, idle time indicates the time from the end of the last chirp to the start of the next chirp. With the chirp duration and chirp slope setting, the sweeping bandwidth B of the chirp is calculated as 4 GHz. We collected 512 frames of echo data in each measurement, including two chirps in each frame, as shown in Figure 9. Only the sampled data of one chirp in each frame need to be imported into MATLAB and analyzed in the signal processing, while the other chirp is regarded as an alternative. We selected 256 frames for cardiopulmonary signal extraction. According to the frame period set, the observation duration is 12.8 s.

4. Results

4.1. Scenario Settings

A total of 11 subjects were selected for this experiment and divided into two groups.

Table 2. Subjects’ basic data.

Group #	Subject #	Basic Health Information			Distance (m)	Reference Value (bpm)
		Sex	Age	SD		Breathing	Heartbeat
1	1	Female	21	Normal breathing and heart rate	0.8	11	66
					1	21	55
					1.3	18	80
					1.5	26	66
	2	Female	24	Normal breathing and heart rate	0.8	20	84
					1	15	81
					1.3	25	83
					1.5	18	82
	3	Male	31	Normal breathing and heart rate	0.8	15	96
					1	17	95
					1.3	14	98
					1.5	24	100
	4	Female	24	Normal breathing and heart rate	0.8	16	75
					1	13	77
					1.3	18	80
					1.5	18	77
	5	Female	26	Normal breathing and heart rate	0.8	11	72
					1	20	81
					1.3	21	79
					1.5	18	84
	6	Female	22	Normal breathing and heart rate	0.8	13	99
					1	24	88
					1.3	24	89
					1.5	21	90
	7	Male	24	Normal breathing and heart rate	0.8	21	88
					1	19	85
					1.3	20	86
					1.5	17	95
2	1	Female	23	Rapid heartbeat with deliberately rapid breathing in some cases	0.8	22	102
					1	31	103
					1.3	18	104
					1.5	18	100
	2	Male	24	Normal heartbeat with deliberately rapid breathing in some cases	0.8	26	84
					1	23	85
					1.3	34	87
					1.5	31	82
	3	Male	24	Rapid heartbeat and deliberately rapid breathing in some cases	0.8	21	108
					1	29	94
					1.3	35	95
					1.5	29	103
	4	Male	29	Normal breathing with deliberately rapid heartbeat in some cases	0.8	15	120
					1	17	100
					1.3	25	90
					1.5	24	92

shows the basic information of each subject, where the reference values are measured by the bracelet mentioned above. The first group had seven samples measured, in which the subjects’ respiratory and heart rates were within the generally normal range, which means that the respiratory rate is within 25 times per minute, and the heart rate is within 100 times per minute. The subjects were required to rest for at least 5 min before each measurement to avoid abnormal vital signs. During the measurements, the subjects sat on the chair in front of the sensor without intentional movement. The experimental environment was an ordinary office environment, and there were no other human targets within the sensing range of the sensor, as shown in Figure 8b. In addition, the subjects’ breathing and heart rate were measured by varying the distance to the radar sensor in the range of 0.8 m, 1 m, 1.3 m, and 1.5 m, respectively. In the second group, there were four samples. The subjects deliberately accelerated the breathing rhythm or completed 15 min of aerobic exercise to increase the respiratory and heartbeat rates, respectively, exceeding the normal range above. This group verified the applicability of the proposed algorithm.

All measurement processes lasted for 25.6 s, and then the stable data were intercepted from them with a length of 12.8 s for subsequent processing.

4.2. Vital Signs Waveform Recovery

As shown in Figure 10, Figure 11 and Figure 12, we selected three different test scenarios to show the waveforms processed in each step to verify the reliability of the proposed algorithm. The scenarios are described as follows:

Scenario #1 in the first group: the subject sat on a chair at a distance of 0.8 m in front of the radar and kept the respiration and heartbeat signs stable before the test with the respiratory and heart rate at 24 and 75 times per minute, respectively, obtained by the smart bracelet;

Scenario #2 in the second group: the subject sat still on the chair at a distance of 1.3 m in front of the radar with a deliberately increased respiratory rate at 35 times per minute;

Scenario #3 in the second group: the subject sat on the chair at a distance of 1 m in front of the radar with some routine behavioral patterns and performed aerobic exercise before the measurements to obtain a deliberately increased respiratory and heart rate at 31 and 103 times per minute.

The breathing and heartbeat signals extracted from the denoised phase difference signals by VMD decomposition operation are shown in (d) and (e) in Figure 10, Figure 11 and Figure 12. We set the number of decomposition levels to 3, and α = 2000, τ = 0, and ε = ${10}^{-7}$. It is not difficult to find that the waveforms of the heartbeat signal are distorted at some peaks marked by the red ellipses, which deteriorate the spectrum estimation. This is the reason for the conditions of peak judgment in peak-seeking in the time domain we proposed. Assuming that the number of peaks in the heartbeat waveform obtained in a measurement is ${\mathrm{N}}_{\mathrm{p}}$, and considering that the length of data we process each time is 256 frames, which corresponds to 12.8 s in slow time, the heartbeat rate can be expressed as:

$${\mathrm{f}}_{\mathrm{h}}={\mathrm{N}}_{\mathrm{p}}\times \frac{1}{12.8} \left(\mathrm{Hz}\right).$$

(16)

4.3. Accuracy Analysis

Table 3. Heart rate decision result of group #1.

Subject #	Range (m)	Reference (bpm)	Measurement (Hz)		Judgment (Hz)	Relative Error
			FFT-CZT	Peek-Seeking
1	0.8	66	1.0989	1.0547	1.0989	0.1%
	1.3	80	1.8732	1.3281	1.3281	1.27%
2	0.8	84	3.5245	1.4063	1.4063	0.45%
	1.3	83	1.3657	1.4063	1.3657	0.37%
3	1.3	98	1.6165	1.6016	1.6165	1.03%
	1.5	100	1.9058	1.6797	1.6797	0.78%
4	1	77	1.2939	1.3672	1.2939	0.82%
	1.3	80	2.4911	1.3281	1.3281	0.39%
5	0.8	72	1.1765	1.25	1.1765	1.96%
	1.5	84	2.7039	1.4063	1.4063	0.45%
6	0.8	99	1.9977	1.6404	1.6404	0.58%
	1	88	1.4493	1.4453	1.4493	1.18%
7	0.8	88	1.4682	1.4844	1.4682	0.11%
	1.5	95	1.2390	1.5625	1.5625	1.32%

and

Table 4. Heart rate decision result of group #2.

Subject #	Range (m)	Reference (bpm)	Measurement (Hz)		Judgment (Hz)	Relative Error
			FFT-CZT	Peek-Seeking
1	0.8	102	2.7029	1.7188	1.7188	1.11%
	1.5	100	1.6901	1.4844	1.6901	1.41%
2	1	87	1.1084	1.4453	1.4453	0.32%
	1.5	82	1.7166	1.3672	1.3672	0.04%
3	0.8	108	0.9375	1.7969	1.7969	0.17%
	1	94	2.3318	1.5625	1.5625	0.27%
4	0.8	120	1.9373	1.4844	1.9373	3.14%
	1.3	90	2.4292	1.4844	1.4844	1.04%

show the heart rate decision process in some cases of the two groups of different measurement scenarios described in Section 4.1, respectively. The relative error is calculated as follows:

$$\mathrm{HR}=\left|1-\frac{{\mathrm{f}}_{\mathrm{h}}\times 60}{{\mathrm{N}}_{\mathrm{ref}}}\right|\times 100\%,$$

(17)

where ${\mathrm{f}}_{\mathrm{h}}$ is the judgment result of the heart rate estimation (unit: Hz); ${\mathrm{N}}_{\mathrm{ref}}$ is the reference value of the heart rate measured by smart bracelet (unit: bpm).

The effect of the Iterative VMD Wavelet-Interval-Thresholding is analyzed below. We redefined the SNR for respiratory and heartbeat signals in the frequency domain for quantitative analysis. Taking the SNR of the respiratory signal as an example, assuming that the spectral line position corresponding to the peak of its spectrum is ${\mathrm{PM}}_0$, its SNR can be expressed as:

$$\mathrm{SNR}=10\lg \left(\frac{{\mathrm{S}}^2\left({\mathrm{PM}}_0\right)+{\mathrm{S}}^2\left({\mathrm{PM}}_0-1\right)+{\mathrm{S}}^2\left({\mathrm{PM}}_0+1\right)}{{\mathrm{S}}^2\left(\mathrm{f}\right)-{\mathrm{S}}^2\left({\mathrm{PM}}_0\right)-{\mathrm{S}}^2\left({\mathrm{PM}}_0-1\right)-{\mathrm{S}}^2\left({\mathrm{PM}}_0+1\right)}\right),$$

(18)

where ${\mathrm{S}}^2\left({\mathrm{PM}}_0\right)$ refers to the energy contained at the peak of the signal spectrum, and ${\mathrm{S}}^2\left(\mathrm{f}\right)$ is the total energy of the spectrum.

Table 5. SNR comparison before and after Iterative VMD Wavelet-Interval-Thresholding.

Subject #	Breathing Signal		Heartbeat Signal
	$\mathbf{S}\mathbf{B}{\mathbf{R}}_0 \left(\mathbf{dB}\right)$	$\mathbf{S}\mathbf{B}{\mathbf{R}}_1 \left(\mathbf{dB}\right)$	$\mathbf{S}\mathbf{B}{\mathbf{R}}_0 \left(\mathbf{dB}\right)$	$\mathbf{S}\mathbf{B}{\mathbf{R}}_1 \left(\mathbf{dB}\right)$
1	−3.1784	−2.7656	−0.8765	0.3583
2	7.8689	8.8274	−2.8444	−1.0186
3	0.1020	4.4074	−3.0217	−1.7594

shows the SNR comparison of the extracted breathing and heartbeat signals of the three subjects mentioned in Section 4.2. ${\mathrm{SNR}}_0$ represents the quality of the signal extracted directly without the denoising algorithm processing, and ${\mathrm{SNR}}_1$ represents the signal quality after denoising. The graph shows that the algorithm improves the SNR of cardiopulmonary signals of the three measured subjects. The SNR of the respiratory signal increases by 1.8922 dB on average, and that of the heartbeat signal increases by 1.441 dB on average, which are better than the SNR improvement in paper [20], with an average increase of 1.11 dB and 0.99 dB, respectively.

Because the amplitude of the chest cavity vibration caused by respiration is relatively large, it is not easy to be disturbed by noise. Thus, we only analyzed the detection accuracy of the heartbeat rate here. Figure 13 shows the comparison of heartbeat rate accuracy with and without the Iterative VMD Wavelet-Interval-Thresholding algorithm. It can be found that the proposed noise reduction algorithm has good adaptability to different test objects and distances. The accuracy is obviously better than without the use of this algorithm. Considering the operation of randomly altering noise sequences in the iterative algorithm, we repeatedly processed each group of collected data 20 times and calculated the solution time each time. The average computation time of the above two algorithms is 5.78 s and 2.36 s, respectively. This time difference is acceptable, which makes it possible for high-precision real-time detection.

Figure 14 shows the Bland–Altman plots of respiratory and heart rates obtained by the proposed method against their reference values, respectively. The X axis is the reference values, and the Y axis is the relative error of the measured values relative to the reference values. The figure shows the limits of agreement at −3.6% and 2.6% for respiratory rates and at −1.7% and 2.3% for heart rates in the first group, while in the second group the limits of agreement are at −2.7% and 2.7% for respiratory rates and at −3.6% and 3.7% for heart rates. Comparing the measured value of respiratory rates with the reference value, the average relative errors of the two groups are 1.33% and 1.01%, respectively. The average relative errors of the heart rates are 0.86% and 1.96%, respectively. The average values mentioned above are different from the mean values marked in Figure 14, which are based on the absolute value and thus more referential in the degree of characterization accuracy. It can be found that the heart rate detection accuracy of the second group is decreased but can be ignored. This means that our proposed algorithm is capable of abnormal situation monitoring.

Table 6. Comparison of the proposed algorithm with other works.

Ref. No.	Output Power (mW)	Error of Breathing Rate	Error of Heartbeat Rate
[12]	17.8	Not mentioned	3.02%
[15]	15.8	6%	20%
[21]	0.5	6.89%	8.09%
This work	17.8	1.33%	1.96%

lists the comparison results with other works. It can be found that under basically the same test scenario, only one FMCW radar sensor can detect respiration and heart rate with a higher accuracy.

5. Discussion

This paper introduces a high-precision vital-signs detection algorithm based on millimeter-wave radar sensors. In practical situations, however, the monitored person does not stay completely still. So, next we discuss the applicability of the proposed algorithm in the non-stationary cases. It is mainly divided into two situations: small-scale micro-movements and large-scale and fast random body movement (RBM). The former includes coughing, sneezing, and arm swinging, while the latter refers to behaviors such as walking around. The difference between the two is whether there is a significant change in the displacement of the chest wall over time. Figure 15 shows the signal fluctuation detected when the subject’s body moves slightly. The marked parts in Figure 15a,c show the influence of cough behavior on the waveform, while the marked parts in (b) and (d) show the influence of arm-swing behavior on the waveform. It can be found that the selected distance window can still extract the phase signal containing the vital signs’ information because the body’s micro-motion behavior does not bring about significant changes in the position of the chest wall. The heartbeat signal extracted by this algorithm only has waveform distortion in the short-term window corresponding to the fretting behavior, and the other part in the entire detected time window is not affected. Therefore, the distorted part of the heartbeat waveform can be detected by the threshold and eliminated, and then spectrum estimation can be performed to obtain the heart rate. Since RBM can lead to changes in the chest wall position, the range-window selection algorithm adopted in this paper is no longer applicable. To eliminate its effect on chest wall vibrations caused by cardiorespiratory activity, two Doppler radar systems [22] or a radar system with a camera [23] were used to remove RBM. In addition, an adaptive noise cancellation technique combined with the polynomial fitting was proposed to offset the influence of RBM on the extraction of small vital sign signals [10]. We expect to implement more efficient RBM removal algorithms to reduce costs in future work.

Meanwhile, deep learning-based technology has also been studied and applied to the detection of vital signs’ information in order to adapt to complex and changeable scenarios. At present, the more mature application is mainly to distinguish whether the target’s respiratory or heart rate is abnormal with high precision. Yin. W et al. proposed using convolutional neural networks independent of each other to learn the signals collected by ECG and the UWB radar and then performed the comprehensive analysis in a cascaded way to achieve mutual complement of biological and physical characteristics of the heartbeat [24]. Li. J proposed to convert the constructed heartbeat waveform into an image and used VGG for identification [25]. For heart rate monitoring, convolutional neural network recognition technology based on the camera-based remote photoplethysmography method [26] and the fiber-optic sensor [27] has achieved high-precision detection. However, there are few related studies based on radar sensors, which may be because deep learning currently mainly relies on image training, and it is difficult to convert the electrical signals propagated by radar into images, which is also a challenge for the future.

6. Conclusions

In this paper, we discussed a systematic method for respiration and heartbeat monitoring based on a FMCW millimeter wave radar. Two scenarios with 11 human subjects at different distances were performed and studied. With the 25.6 s short-time monitoring in each experiment, a window of 12.8 s was randomly selected for signal processing and vital signs’ information extraction. The results show that the average accuracy of respiratory and heart rate estimation is not less than 98% using the FMCW sensor and the proposed processing algorithm in the test scenario where subjects are stationary. The possible noise contained in the test process and the corresponding noise removal measures were discussed. The experiments showed that the SNR of the vital signals can be increased by using the Iterative VMD Wavelet-Interval-Thresholding algorithm. In addition, the decision method based on spectrum refinement and peek-seeking in the time domain also helps to improve the accuracy of rate estimation. We envision that this approach could be applied to various practical situations, such as home health-monitoring for the elderly, medical assistance for burn victims, and driver monitoring systems.