Enhanced Vital Parameter Estimation Using Short-Range Radars with Advanced Motion Compensation and Super-Resolution Techniques

1. Introduction

Non contact vital sign estimation (VSE) techniques using radar systems have been extensively researched and employed for various applications [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. These techniques estimate vital signs, such as respiratory and cardiac rates (RR and CR), by detecting very small range variations in the skin caused by these vital signs. Because lung and heart vibrations contribute to the movement of the chest and back, the phase of the reflected radar signal contains RR and CR information [15]. Numerous types of short-range radars have been used to accurately extract these phases for VSE [1,2,3,6].

To achieve accurate and cost-effective VSE, various studies have explored the use of continuous-wave (CW), frequency-modulated continuous-wave (FMCW), and impulse-radio ultra-wideband (IR-UWB) signals. One study proposed an efficient CW radar system for RR measurement under one-dimensional body motion [16]. Another method introduced motion compensation (MOCOM) to eliminate phase errors caused by random body movements [17]. This approach assumes a CW signal and uses two echo signals from both the chest and back to cancel phase errors. Schires [18] developed an IR-UWB radar system integrated into a car seat backrest, matched to the human body. Xiang [19] presented a high-precision VSE method using FMCW radar, in which noise was considerably reduced via a variation-mode decomposition (VMD) wavelet-internal-thresholding algorithm. A recently proposed method used a fuzzy rule to accurately determine RR and remove distorted CR values due to motion [20].

Additionally, substantial research has been directed at improving both hardware and software techniques for VSE. Liu et al. [21] suggested an interferometry-based frequency-sweeping technique to achieve high ranging accuracy and a wide detection range. Gu [22] and Yavari [23] introduced methods based on digital post-distortion and compensation to enhance accuracy by mitigating undesirable interference. To further increase accuracy and suppress noise, signal decomposition methods such as wavelet transform [24] and empirical mode decomposition (EMD) [25] have been applied to VSE.

However, existing methods do not fully address the issues caused by considerable movement of an individual. Continuous motion during observation can result in substantial phase errors. The method described in [16] may fail because generally the movement is not one-dimensional. The approach in [17] assumes that the phases caused by forward and backward motions have opposite signs and equal magnitudes, allowing error removal by multiplying the spectra of the forward and backward radars. This method can fail if other motion components are present. Additionally, because the amplitude of the signal reflected from the front of the body is significantly larger than that from the back, spectrum multiplication cannot eliminate forward–backward motion artifacts.

The results in [18] may be significantly degraded if the back moves or if the matching condition between the body and the seat changes. The method in [19] did not incorporate MOCOM, and thus it could fail with body movement. The approach in [20] fails if the estimated value remains inaccurate over time, as it maintains the previous result once distortion is detected. Methods in [21,22,23] are not relevant to MOCOM and may be inaccurate in the presence of movement. While decomposition methods [24,25] can successfully separate signals in low-noise environments, their accuracy is questionable owing to spectrum blurring, even at very low velocities and accelerations. Moreover, when respiratory and cardiac activities are partially blocked or the radar line-of-sight changes rapidly, the signal-to-noise ratio (SNR) can decrease significantly, leading to poor separation results. Therefore, an efficient method to eliminate the effects of motion is required.

In this study, we propose an efficient method that achieves super-resolution vital sign estimation (VSE) using low-cost short-range radar. The proposed method involves the following steps:

• Denoising and signal separation: The radar signal is denoised using principal component analysis (PCA).
• Motion compensation (MOCOM): Based on an analysis of the phase components of the echo signal, two efficient MOCOM methods, MOCOM $\#1$ and MOCOM $\#2$, are introduced.
• Noise reduction and auto-focusing: To further reduce noise and enhance resolution, the respiratory and cardiac signals are separated again and auto-focused.
• Super-resolution spectrum estimation: The multiple signal classification (MUSIC) method [26], a super-resolution technique, is used to obtain the spectra of the separated signals with very high resolution.

The remainder of this paper is organized as follows: Section II introduces the radar signal model, analyzing the phase components and issues caused by movement. Section III details the proposed method. Section IV presents experimental results demonstrating the method’s efficiency, and Section V concludes the paper.

2. Signal Model and Problem Analysis

2.1. Radar Signal Model

To apply the proposed method to various low-cost radars, the signal models of two representative short-range radars were assumed to be IR-UWB and FMCW. Using the observation geometry shown in Figure 1, the methods employed in previous studies [20,27] were modified, where $r_{RLOS}={[1,0,0]}^T$ is the radar line-of-sight (RLOS) vector. When there is no body movement, the range vector $r_r\left(t\right)$ of the chest caused by RR is $r_r\left(t\right)={[d_{r}cos(2\pi f_r\left(t\right)t+{\phi }_{r0}),0,0]}^T$ and the range variation along the RLOS is given by

$$R_r\left(t\right)={\mathbf{r}}_{RLOS}^T·{\mathbf{r}}_r\left(t\right)=d_{r}cos(2\pi f_r\left(t\right)t+{\phi }_{r0}),$$

(1)

where · denotes the inner product operator, $d_r$ denotes the variation of chest skin, ${\phi }_{r0}$ represents the initial phase, and

$$f_r\left(t\right)=f_{r0}+▵f_{r}cos\left(2\pi f_{r0}t\right),$$

(2)

is the real-time respiration rate, where $f_{r0}$ denotes the fundamental RR frequency and $▵f_r$ represents the slight variation amplitude in the respiration rate responsible for high-order harmonics.

Similarly, the range vector ${\mathbf{r}}_c\left(t\right)$ of the chest caused by the CR is ${\mathbf{r}}_c\left(t\right)={[d_{c}cos(2\pi f_c\left(t\right)t+{\phi }_{c0}),0,0]}^T$ and its projection $R_c\left(t\right)$ along the RLOS is

$$R_c\left(t\right)={\mathbf{r}}_{RLOS}^T·{\mathbf{r}}_c\left(t\right)=d_{c}cos(2\pi f_c\left(t\right)t+{\phi }_{c0}),$$

(3)

where $d_c$ is the variation in the chest skin, ${\phi }_{c0}$ is the initial phase, and

$$f_c\left(t\right)=f_{c0}+▵f_{c}cos\left(2\pi f_{c0}t\right),$$

(4)

is the real-time CR function.

For the transmitted signal, we assumed the following functions of the IR-UWB and FMCW radars with bandwidth B:

$$s_t\left(t\right)=\sqrt{P_t}{f}_{IR}\left(t\right)\mathrm{and}{s}_t\left(t\right)=\sqrt{P_t}{f}_{FMCW}\left(t\right),$$

(5)

where

$$f_{IR}\left(t\right)=gauss\left(Bt\right)exp(-j2\pi f_{cen}t),$$

(6)

$$f_{FMCW}\left(t\right)=exp\left(j2\pi \left(f_{cen}t+{\displaystyle \frac{K_r{t}^2}{2}}\right)\right),$$

(7)

where $\sqrt{P_t}$ denotes the amplitude, $gauss\left(Bt\right)$ denotes the Gaussian pulse with B, $f_{cen}$ denotes the carrier frequency of the transmitted signal, $K_r=B/T_{chirp}$, and $T_{chirp}$ denotes the width of the transmitted chirp. The received FMCW signal is compressed into a sinc function by de-chirping and employing a Fourier transform (FT) [28].

For the chest observed with a given pulse repetition time $T_{PRI}$, the received signal of the IR-UWB and FMCW radars is expressed as

$$s_r(t,t_p)=\sqrt{P_r}psf\left(B(t-\tau \left(t_p\right))\right)exp\left(j2\pi f_{cen}(t-\tau \left(t_p\right))\right),$$

(8)

where $psf\left(Bt\right)=gauss\left(Bt\right)$ and $sinc\left(Bt\right)$ for the IR-UWB and FMCW radars, respectively. $t_p$ is the slow time sampled at $T_{PRI}$, and the time delay $\tau \left(t_p\right)$ is

$$\tau \left(t_p\right)={\displaystyle \frac{2}{c}}{R}_w\left(t_p\right)+{\displaystyle \frac{2}{c}}{R}_r\left(t_p\right)+{\displaystyle \frac{2}{c}}{R}_c\left(t_p\right).$$

(9)

c is the velocity of light, and $R_w\left(t_p\right)$ is the range variation of the body that should be removed. $\sqrt{P_r}$ is the amplitude of the received signal and can be expressed by the radar equation as follows:

$$P_r={\displaystyle \frac{P_t{G}_t{G}_r{\lambda }^2{\sigma }_{RCS}}{{\left(4\pi \right)}^3{R}^4}},$$

(10)

where $G_t$ and $G_r$ are the gains of the transmitting and receiving antennas, respectively; $\lambda$ is the wavelength, ${\sigma }_{RCS}$ is the radar cross section, and R is the distance to the chest.

$s_r(t,t_p)$ is down-converted to baseband signal $s_B(t,t_p)$ as follows:

$$\begin{array}{cc}\hfill s_B(t,t_p)& =s_r(t,t_p)exp(-j2\pi f_{cen}t)\hfill \\ \hfill & =\sqrt{P_r}psf\left(B(t-\tau \left(t_p\right))\right)exp(-j2\pi f_{cen}\tau \left(t_p\right)).\hfill \end{array}$$

(11)

Assuming $R_w\left(t_p\right)=R_0-v_0{t}_p-0.5a_0{t}_p^2$, where $R_0$ is the initial range and $v_0$ and $a_0$ are the velocity and acceleration of the body, respectively, (11) becomes

$$\begin{array}{cc}\hfill s_B(t,t_p)& =\sqrt{P_r}psf\left(B(t-\tau \left(t_p\right))\right)exp(-j2\pi f_{cen}{\displaystyle \frac{2}{c}}{R}_0)\hfill \\ \hfill & \times exp\left(j2\pi f_{cen}{\displaystyle \frac{2}{c}}\left(v_0{t}_p+{\displaystyle \frac{1}{2}}{a}_0{t}_p^2\right)\right)\hfill \\ \hfill & \times exp\left(-j2\pi f_{cen}{\displaystyle \frac{2}{c}}{R}_r\left(t_p\right)\right)exp\left(-j2\pi f_{cen}{\displaystyle \frac{2}{c}}{R}_0\right).\hfill \end{array}$$

(12)

2.2. Effect of the Motion of the Rigid Body

The clipped signal $s_{cp}\left(t_p\right)$ at range r with maximum energy is a function of $t_p$ only and contains information on $f_r$ and $f_c$. $s_{cp}\left(t_p\right)$ is given by

$$\begin{array}{cc}\hfill s_{cp}\left(t_p\right)=& \sqrt{P_r}exp\left(-j2\pi f_{cen}{\displaystyle \frac{2}{c}}{R}_0\right)\hfill \\ \hfill & \hfill \times exp\left(j2\pi f_{cen}{\displaystyle \frac{2}{c}}(v_0{t}_p+{\displaystyle \frac{1}{2}}{a}_0{t}_{p^2})\right)\\ \hfill & \times exp\left(-j2\pi f_{cen}{\displaystyle \frac{2}{c}}{R}_r\left(t_p\right)\right)\hfill \\ \hfill & \times exp\left(-j2\pi f_{cen}{\displaystyle \frac{2}{c}}{R}_c\left(t_p\right)\right).\hfill \end{array}$$

(13)

The Fourier transform of $s_{cp}\left(t_p\right)$ can be expressed as

$$\begin{array}{cc}\hfill s_{cp}\left(t_p\right)=& \sqrt{P_r}exp\left(-j2\pi f_{cen}{\displaystyle \frac{2}{c}}{R}_0\right)\hfill \\ \hfill & \otimes FT\left\{exp\left(j2\pi f_{cen}{\displaystyle \frac{2}{c}}(v_0{t}_p+{\displaystyle \frac{1}{2}}{a}_0{t}_{p^2})\right)\right\}\hfill \\ \hfill & \otimes FT\left\{exp\left(-j2\pi f_{cen}{\displaystyle \frac{2}{c}}{R}_r\left(t_p\right)\right)\right\}\hfill \\ \hfill & \otimes FT\left\{exp\left(-j2\pi f_{cen}{\displaystyle \frac{2}{c}}{R}_c\left(t_p\right)\right)\right\}.\hfill \end{array}$$

(14)

where $FT\{·\}$ denotes the FT operation, and ⊗ denotes the convolution operation.

As the first component is a constant envelope of one, its effect can be eliminated. The FT of the second exponential component in the above equation is a linear frequency modulation signal due to the acceleration shifted by the Doppler frequency $f_0=2f_{c}v/c$. Using the stationary phase method [29], it can be transformed into

$$G\left(f\right)=Rect\left({\displaystyle \frac{f-f_0-K{T}_p/2}{K{T}_p}}\right)exp\left({\displaystyle \frac{-j\pi {(f-f_0)}^2}{K}}\right),$$

(15)

where $Rect(·)$ is a rectangle function, $T_p$ is the coherent processing interval, and K is given by

$$K=-2f_c{a}_0/c.$$

(16)

The FT of the third and fourth components can be divided into in-phase and quadrature components as follows:

$$FT\{exp(j4\pi f_{cen}{R}_r\left(t_p\right)/c\}=I_r\left(t_p\right)+Q_r\left(t_p\right),$$

(17)

$$FT\{exp(j4\pi f_{cen}{R}_c\left(t_p\right)/c\}=I_c\left(t_p\right)+Q_c\left(t_p\right),$$

(18)

where $I_r\left(t_p\right)$ and $Q_r\left(t_p\right)$ are the in-phase and quadrature components of RR, respectively, and $I_c\left(t_p\right)$ and $Q_c\left(t_p\right)$ are those of CR. $I_r\left(t_p\right)$, $Q_r\left(t_p\right)$, $I_c\left(t_p\right)$, and $Q_c\left(t_p\right)$ can be approximated using the Bessel function $J_i\left(x\right)$ of the first type of order i as follows [17]:

$$I_r\left(t_p\right)=2(C_{20}cos\left(2{\omega }_r{t}_p\right)+C_{02}cos\left(2{\omega }_r{t}_p\right)+...),$$

(19)

$$Q_r\left(t_p\right)=2(C_{10}sin\left({\omega }_r{t}_p\right)+C_{01}sin\left({\omega }_r{t}_p\right)+...),$$

(20)

$$I_c\left(t_p\right)=2(C_{20}cos\left(2{\omega }_c{t}_p\right)+C_{02}cos\left(2{\omega }_c{t}_p\right)+...),$$

(21)

$$Q_r\left(t_p\right)=2(C_{10}sin\left({\omega }_c{t}_p\right)+C_{01}sin\left({\omega }_c{t}_p\right)+...),$$

(22)

where ${\omega }_r=2\pi f_r$, ${\omega }_c=2\pi f_c$, and $C_{ij}=J_i(4\pi f_{cen}{d}_r/c)J_j(4\pi f_{cen}{d}_c/c)$ for $i=0,1,...,\infty$ and $j=0,1,...,\infty$. Notably, i and j cannot be zero simultaneously, as they represent the DC component.

Equations (2), (4), and (19)-(22) represent spectra of RR and CR, composed of harmonics of $f_r$ and $f_c$; as $J_i\left(x\right)$ is higher for smaller i and RR has a larger $d_r$ (e.g., $0.001\le d_r\le 0.003m$) than $d_c$ (e.g., $0.00005\le d_c\le 0.0005m$), $f_r$ has the largest amplitude. However, the harmonics of $f_r$ may have negative effects on $f_c$ when they coincide with $f_c$.

The convolution of the discrete spectra of (17) and (18) with the FT of (15) yielded a distorted spectrum owing to the frequency modulation caused by the acceleration and shift of the modulated spectrum (Figure 1). In the case where $v_0\ne 0m/s$ and $a_0=0m/s^2$, a peak appears at a false location owing to $f_0$, and when $v_0=0m/s$ and $a_0\ne 0m/s^2$, the spectrum can be blurred because of the convolution with the Doppler bandwidth caused by $a_0$. With nonzero $v_0$ and $a_0$, the spectrum is severely distorted, as the peaks are both displaced and blurred. Therefore, MOCOM is an indispensable preprocessing step for accurate RRE and CRE.

Eliminating the effect of $v_0$ and $a_0$ using MOCOM is challenging. One probable method is to estimate $v_0$ first by identifying the largest peak in $S_{cp}\left(f\right)$ and subsequently focusing the spectrum by finding $a_0$. However, the exact value of $v_0$ can seldom be found because the shape of the spectrum caused by $a_0$ may not be flat; because of Gibbs’ phenomenon, overshot values can be larger than that at $f_0$, and the fluctuation of the Doppler bandwidth owing to $a_0$ and the interference by neighboring harmonics makes the spectrum more irregular. Therefore, $S_{cp}\left(f\right)$ should not occur in the presence of $v_0$ and $a_0$.

The problem of finding peaks in the spectrum of $S_{cp}\left(f\right)$ can be solved simply by using the phase information, that is, the arctangent demodulation $\theta \left(t_p\right)$ of $S_{cp}\left(t\right)$, given as

$$\begin{array}{cc}\hfill \theta \left(t_p\right)& =arctan\{Q_{cp}\left(t_p\right)/I_{cp}\left(t_p\right)\}\hfill \\ \hfill & \approx \{4\pi /{\lambda }_c\}\{R_0+R_w\left(t_p\right)+R_r\left(t_p\right)+R_c\left(t_p\right)\}\hfill \\ \hfill & \approx {\theta }_0+{\theta }_w\left(t_p\right)+{\theta }_r\left(t_p\right)+{\theta }_c\left(t_p\right),\hfill \end{array}$$

(23)

where $I_{cp}\left(t_p\right)$ and $Q_{cp}\left(t_p\right)$ are the in-phase and quadrature components of $s_{cp}\left(t_p\right)$, respectively and ${\lambda }_c$ is the wavelength of $f_{cen}$. $\theta \left(t_p\right)$ is significantly distorted by ${\theta }_w\left(t_p\right)$, because it is considerably larger than ${\theta }_r\left(t_p\right)$ and ${\theta }_c\left(t_p\right)$ (Figure 2). To avoid ambiguity in extracting $\theta \left(t_p\right)$, a phase unwrapping technique should be used [20].

Considering the effects of the surrounding clutter and noise, (23) is modified as follows:

$$\theta \left(t_p\right)\approx {\theta }_0+{\theta }_w\left(t_p\right)+{\theta }_r\left(t_p\right)+{\theta }_c\left(t_p\right)+{\theta }_{noise}\left(t_p\right)+{\theta }_{clutter}\left(t_p\right),$$

(24)

where ${\theta }_{noise}\left(t\right)$ and ${\theta }_{clutter}\left(t\right)$ are the phases of system noise and background clutter, respectively. As $\theta \left(t_p\right)$ includes a linear mixture of the effects of the body, noise, and clutter, RRE and CRE become more challenging.

When Fourier transformed, $FT\theta \left(t_p\right)$ is the linear component of each component, as follows:

$$\begin{array}{cc}\hfill FT\left\{\theta \left(t_p\right)\right\}=& FT\left\{{\theta }_0\right\}+FT\left\{{\theta }_w\left(t_p\right)\right\}+FT\left\{{\theta }_r\left(t_p\right)\right\}\hfill \\ & +FT\left\{{\theta }_c\left(t_p\right)\right\}+FT\left\{{\theta }_{noise}\left(t_p\right)\right\}\hfill \\ & +FT\left\{{\theta }_{clutter}\left(t_p\right)\right\},\hfill \end{array}$$

(25)

Because ${\theta }_0=4\pi /{\lambda }_c/R_0$ is constant, $FT\theta \left(t_p\right)$ has a large amplitude at $f=0Hz$. In addition, ${\theta }_r\left(t\right)$ and ${\theta }_c\left(t\right)$ are the two cosine functions of $f_r$ and $f_c$, respectively, and the peaks can appear at the corresponding frequencies. However, the effect of ${\theta }_w\left(t_p\right)$ should be removed before the FT of $\theta \left(t_p\right)$, because it has a very large amplitude (see Figure 3). In addition, the effects of ${\theta }_{noise}\left(t_p\right)$ and ${\theta }_{clutter}\left(t_p\right)$ should be minimized.

3. Proposed Method

3.1. Summary of the Proposed Method

This subsection introduces an efficient method for estimating the vital parameters $f_r$ and $f_c$ with super-resolution by compensating for motion and applying the MUSIC algorithm (Figure 4).

➀

Align the range profile (RP) history of the received radar signal to position the scatterer in an identical location during the coherent processing interval, thereby removing range migration due to movement.

➁

Clip the radar signal around the range bin with the maximum energy to estimate the vital parameters.

➂

Apply PCA to denoise the radar signal and extract the torso and vital signals.

➃

Estimate the rigid-body Doppler frequency using a fast Fourier transform (FFT) and remove the corresponding velocity.

➄

Extract the phase history using a phase unwrapping technique.

➅

Coarsely estimate the phase history via Gaussian filtering of the unwrapped phase and remove phase errors using the filtered phase history (MOCOM $♯1$).

➆

Remove residual errors using the envelope (MOCOM $♯2$).

➇

Reconstruct the complex signal using amplitude- and motion-compensated phases.

➈

Further suppress noise and obtain super-resolution by separating the respiratory and cardiac signals using a low-pass filter (LPF) and high-pass filter (HPF).

➉

Remove residual phase errors using a phase-adjustment technique for each of the separated signals.

⑪

Optionally apply zero padding and then use the MUSIC algorithm to obtain the super-resolution spectrum of the separated signals.

⑫

Estimate fr and fc by identifying the maximum peaks in each super-resolution spectrum.

3.2. Main Idea of the Proposed Method

3.2.1. Range Alignment

First, because ${\tau }_p\left(t_p\right)$ in (7) is heavily dependent on $v_0$ and $a_0$, the scatterer migrates into neighboring range bins. Consequently, $s_{cp}\left(t_p\right)$ cannot be formed by clipping $s_B(r,t_p)$ at a fixed range $r=ct/2$. Therefore, the RPs should be aligned such that the scatterer is located in the same range bin for all $t_p$. While exact estimation of motion parameters $v_0$ and $a_0$ is impossible without errors, we used a widely adopted non-parametric method. Among various methods, we employed entropy minimization, which is commonly used in inverse synthetic aperture radar imaging [30].

The entropy between two RPs ${\mathbf{rp}}_1\left(r\right)$ and ${\mathbf{rp}}_2(r-\epsilon )$ is defined by

$$H_{r{p}_1,r{p}_2}\left(\epsilon \right)=-{\int }_0^{r_{max}}{\mathbf{rp}}_{av}\left(r\right)ln{\mathbf{rp}}_{av}\left(r\right)dr,$$

(26)

where

$${\mathbf{rp}}_{av}\left(r\right)={\displaystyle \frac{|{\mathbf{rp}}_1\left(r\right)|+|{\mathbf{rp}}_2(r-\epsilon )|}{{\int }_0^{r_{max}}\left(\right|{\mathbf{rp}}_1\left(r\right)|+|{\mathbf{rp}}_2(r-\epsilon )\left|\right)dr}},$$

(27)

and $r_{max}=T_{PRI}·c/2$ is the maximum range. According to this criterion, the $\epsilon$ value that minimizes the 1D entropy is the shift that aligns well with ${\mathbf{rp}}_2\left(t\right)$. Consequently, the second RP of $s_B(r,t_p)$ and the relative shift $\epsilon$ minimizing (26) along with the average of the aligned RPs are determined. The RP is aligned with a shift by $\epsilon$.

3.2.2. Extraction of Vital Signals and Denosing Using PCA

Because the RPs of the body can be contaminated by noise or surrounding clutter, the phase of the RP can be distorted. In addition, determining the exact location of the chest is difficult, because the reflected signal is the sum of several body parts, mainly the torso, which can be buried within the torso signal. Therefore, a method that can denoise and separate vital signals is required. In this study, we applied PCA to fulfill these two requirements.

Assuming n RPs, let $\mathbf{A}$ be an $m\times n$ matrix constructed by sampling $\pm m/2$ range bins around the range bin with the maximum energy as follow:

$$\mathbf{A}={[{\mathbf{o}}_0,{\mathbf{o}}_1,...,{\mathbf{o}}_{m-1}]}^T,$$

(28)

where

$${\mathbf{o}}_k={[{\mathbf{rp}}_k\left(r_0\right),{\mathbf{rp}}_k\left(r_1\right),...,{\mathbf{rp}}_k\left(r_{n-1}\right)]}^T,$$

(29)

the covariance matrix is calculated by

$$\mathbf{Cov}={\displaystyle \frac{1}{m-1}}\sum _{k=0}^{m-1}{\mathbf{o}}_k{\mathbf{o}}_k^H,$$

(30)

where H is the conjugate transpose of vector. Subsequently, suing the eigenvalue decomposition of $\mathbf{Cov}$, given by

$$\mathbf{Cov}=\mathbf{E}\Lambda {\mathbf{E}}^T,$$

(31)

According to PCA, the orthonormal eigenvectors corresponding to large eigenvalues span the signal subspace of the $\mathbf{Cov}$, while those corresponding to small eigenvalues span the noise subspace. Therefore, by projecting $\mathbf{A}$ onto eigenvectors with large eigenvalues, the signal can be denoised. Additionally, because the large-valued torso signal corresponds to $\mathbf{A}$ projected onto the eigenvector $e_1$ with the largest eigenvalue, the small-valued respiratory and cardiac signals can be separated by projections onto eigenvectors with smaller eigenvalues. In this study, we use two eigenvectors, $e_2$ and $e_3$, corresponding to the second and third largest eigenvalues, respectively. The matrix $\mathbf{A}$ is projected onto ${\mathbf{E}}_{prj}=\left[e2e3\right]$ as

$${\mathbf{A}}_{prj}=[{\mathbf{a}}_1,{\mathbf{a}}_2]={\mathbf{AE}}_{prj},$$

(32)

Finally, the denoised sum of the respiratory and the cardiac signal is

$${\mathbf{s}}_{sum}={\mathbf{a}}_1+{\mathbf{a}}_2,$$

(33)

and it is used instead of the $s_{cp}\left(t_p\right)$ in (13). Note that $s_{sum}$ is a function of $t_p$.

3.2.3. Estimation and Removal of the Rigid-Body Doppler

The Doppler frequency $f_D$ corresponding to the rigid body velocity $v_0$ is estimated by

$${\widehat{f}}_D=\underset{f}{\mathrm{argmax}}\left\{\left|FFT\right\{s_{sum}\left\}\right|\right\},$$

(34)

where $FFT·$ denotes the FFT function. The rigid-body phase error is removed using

$$s_{cp}\left(t_p\right)={\mathbf{s}}_{sum}\left(t_p\right)exp\left(j2\pi {\widehat{f}}_D{t}_p\right).$$

(35)

3.2.4. MOCOM $\#1$

Although the rigid-body Doppler is removed using (35), the spectrum can be blurred by the effect of the acceleration $a_0$ of the rigid body. In addition, because the FT of the acceleration component yields a Doppler bandwidth with an irregular passband (Figure 1), the estimation of v0 may not be accurate. Consequently, Equation (35) may yield an estimation error, resulting in a shift in the frequency components of $S_{cp}\left(f\right)$.

This problem can be solved using only unwrapped $\theta \left(t_p\right)$. $\theta \left(t_p\right)$ includes the linear sum of ${\theta }_w\left(t_p\right)$, ${\theta }_r\left(t_p\right)$, and ${\theta }_c\left(t_p\right)$, and $(v_0,a_0)$ can be addressed separately. Because the spectrum is severely blurred in the frequency domain of the phase, extracting this component by applying signal-separation methods is difficult. In addition, because the shape of the unwrapped phase includes the effects of motion, removing the motion component from the time domain is more profitable.

The error can be perfectly removed once the motion component is accurately estimated using a proper polynomial, such as $P\left(t_p\right)=a_0+a_1{t}_p^1+a_2{t}_p^2+...+a_q{t}_p^q$ to describe $\theta \left(t_p\right)$, and a curve-fitting algorithm to find the parameter $P\left(t_p\right)$. However, the curve-fitting algorithm is time consuming, and numerous phase errors can occur if $P\left(t_p\right)$ does not match $\theta \left(t_p\right)$ accurately. Therefore, we coarsely compensated for the motion using a smoothing kernel function to smooth $\theta \left(t_p\right)$. Among the various methods, we used the following Gaussian kernel:

$$K\left(t_p\right)={\displaystyle \frac{d{t}_p}{\sqrt{2\pi }\sigma }}exp\left(-{\displaystyle \frac{{(t_p-{\mu }_{t_p})}^2}{2{\sigma }^2}}\right),$$

(36)

where $d{t}_p$ is the sampling interval of $t_p$, $\sigma$ is the standard deviation of the kernel, and ${\mu }_{t_p}$ is the mean value of $t_p$.

The smoothed phase for the MOCOM is given by

$${\theta }_s\left(t_p\right)=\theta \left(t_p\right)\otimes K\left(t_p\right).$$

(37)

Equation (37) is time consuming; therefore, it is converted in the frequency domain as follows:

$${\theta }_s\left(t_p\right)=IFT\left\{FT\left\{\theta \left(t_p\right)\right\}\times FT\left\{K\left(t_p\right)\right\}\right\},$$

(38)

where $IFT·$ is the inverse FT. Finally, by directly subtracting ${\theta }_s\left(t_p\right)$ from $\theta \left(t_p\right)$ as follows:

$${\theta }_{M1}\left(t_p\right)=\theta \left(t_p\right)-{\theta }_s\left(t_p\right),$$

(39)

the effect of ${\theta }_w\left(t_p\right)$ is eliminated. This method is direct and relatively simple (Figure 5). However, ${\theta }_{M1}\left(t_p\right)$ includes errors because various errors caused by phase unwrapping and noise cause a mismatch between ${\theta }_s\left(t_p\right)$ and the actual phase trajectory. Therefore, additional MOCOM is required. Therefore, we define this MOCOM as MOCOM $\#1$ to distinguish it from the additional MOCOM, MOCOM $\#2$.

The MOCOM $\#2$ proposed in this study reduces the residual error, as shown in Figure 6. As the exact value of the residual error is not known a priori, using the characteristics of the respiratory and cardiac signals is more effective. Owing to the sinusoidal nature of the two signals, the sum of the motion-compensated signals in the time domain fluctuates within constant upper and lower limits (Figure 5). Therefore, the ideal MOCOM flattens the phase by estimating the unknown residual error, as shown in Figure 3. A parametric method is ineffective for estimating that curve. Instead, we successively use three phase envelopes after MOCOM $\#1$: the upper, lower, and mean envelopes (Figure 7).

3.2.5. MOCOM $\#2$

Various methods for detecting the envelope of a signal in a communication area exist, such as squaring, low-pass filtering, and the Hilbert transform [31,32]. In this study, we use a simple detection method by using the maximum (or minimum) within a window with length L. By moving a window from left to right with a given time step $▵t_w$, the maximum value within the window is stored as ${\theta }_{max}(k▵t_w)$, and the minimum value is stored as ${\theta }_{min}(k▵t_w)$, where $k=0,1,2,...,round(T_f/▵t_w)$. $T_f$ is the length of the time frame, and $round\left(\right)$ is the function that determines the nearest integer (Figure 8). Zeros are added when the number of signals is smaller than L near the start and end of ${\theta }_{M1}\left(t_p\right)$.

The upper envelope $E_{upper}\left(t_p\right)$ and lower envelope $E_{lower}\left(t_p\right)$ are then obtained by curve-fitting ${\theta }_{max}(k▵t_w)$ and ${\theta }_{min}(k▵t_w)$. Among various methods, we use cubic spline interpolation method [33]. Finally, the mean envelope is obtained as

$$E_{mean}\left(t_p\right)=(E_{upper}\left(t_p\right)+E_{lower}\left(t_p\right))/2,$$

(40)

and the error is reduced by

$${\theta }_{M2}\left(t_p\right)={\theta }_{M1}\left(t_p\right)-E_{mean}\left(t_p\right).$$

(41)

However, depending on the length of the window L and time step $▵t_w$, the envelope may not describe the overall shape of ${\theta }_{M1}\left(t_p\right)$. In particular, near $t_p=0$ and $t_p=T_p$, the estimated envelope may increase or decrease significantly, owing to the zeros added to the window. In addition, the maximum and minimum values at a certain time may be the same, yielding a large difference from the near maximum or minimum value. In this case, ${\theta }_{M2}\left(t_p\right)$ may not be flat but curved near that point. This problem can be solved by applying MOCOM $\#2$ successively until the variance of $E_{mean}\left(t_p\right)$ is close to 0; that is,

$$var\left(E_{mean}\left(t_p\right)\right)\approx 0.$$

(42)

3.2.6. Reconstruction of Signals and Separation of Vital Signals Using LPF and HPF

To further remove the residual phase errors and obtain super-resolution for $f_r$ and $f_c$, the complex signal is first reconstructed as follows:

$$s_{M2}\left(t_p\right)=\left|s_{cp}\left(t_p\right)\right|exp\left(j{\theta }_{M2}\left(t_2\right)\right).$$

(43)

Subsequently, (43) is separated into respiratory and cardiac signals, that is, $s_R\left(t_p\right)$ and $s_C\left(t_p\right)$, respectively. Because the amplitudes of the cardiac spectrum are very small and can be influenced by the ripples of the passband, we used a maximally flat Butterworth filter: LPF, for the respiratory signal and HPF for the cardiac signal, in the following form:

$$HPF\left(\omega \right)={\displaystyle \frac{b_1+b_2{e}^{-j\omega }+b_3{e}^{-j2\omega }+...+b_{N_B}{e}^{-j(N_B-1)\omega }}{a_1+a_2{e}^{-j\omega }+a_3{e}^{-j2\omega }+...+a_{N_B}{e}^{-j(N_B-1)\omega }}},$$

(44)

where $N_B$ is the filter order.

3.2.7. Phase Adjustment

Despite its maximum flatness, ${\theta }_{M2}\left(t_p\right)$ contains only CR and RR frequencies. The residual phase errors may occur owing to estimation error, noise, and clutter. Consequently, residual errors should be removed using an appropriate method. The residual error is significantly reduced by using MOCOM $\#1$ and $\#2$; therefore, we used the 1D entropy of the following signal spectrum:

$$s_{fil}\left(t_p\right)=\left|s_{fil}\left(t_p\right)\right|exp\left(j{\theta }_{fil}\left(t_p\right)\right),$$

(45)

where $s_{fil}\left(t_p\right)$ is the filtered signal, representing either $s_R\left(t_p\right)$ and $s_C\left(t_p\right)$. Subsequently, the phase-error $\widehat{\theta }\left(t_p\right)$ is determined by minimizing the entropy, defined by

$$H\left(s_{PA}\left(t_p\right)\right)=-{\int }_{-\infty }^{\infty }\left|S_{PA}\left(f\right)\right|ln|S_{PA}\left(f\right)|df,$$

(46)

where

$$s_{PA}\left(t_p\right)=s_{fil}\left(t_p\right)exp\left(\widehat{\theta }\left(t_p\right)\right),$$

(47)

and $S_{PA}\left(f\right)=FT\left\{s_{PA}\left(t_p\right)\right\}$. This method is a major topic of interest and is known as phase adjustment in radar imaging. In this study, we applied the widely used minimum entropy phase adjustment (MEPA) [34]. The MEPA is fast and can be conducted in real time using the entropy gradient.

3.2.8. Zero-Padding and Application of MUSIC

While the spectrum of the vital signal can be significantly focused via phase adjustment, the vital signal can still be distorted by the side lobes of neighboring frequency components; this is particularly true for the small value at $f_c$, which is vulnerable to this effect. Additionally, because the frequency resolution is inversely proportional to the observation time $T_0$, a large $T_0$ is required, and the radar signal during this extended $T_0$ may contain significant movement. Consequently, the estimation accuracy of $f_r$ and $f_c$ may be degraded.

To improve estimation accuracy, the phase-adjusted signal is zero-padded to upsample the spectrum. However, this step is optional, as increasing the number of sampling points in the frequency domain does not increase the frequency resolution. Subsequently, MUSIC is applied to the zero-padded (optional) signal to obtain a super-resolution spectrum. MUSIC uses the characteristic that the direction vector containing the signal is orthogonal to the noise subspace consisting of noise eigenvectors.

Using, the vector $\mathbf{y}$ containing N samples of the phase ${\theta }_{PA}\left(t_p\right)$ of the radar signal after the phase adjustment, which given by

$$\begin{array}{cc}\hfill \mathbf{y}& =[y_0,y_1,...,y_{N-1}]\hfill \\ \hfill & =[{\theta }_{PA}\left(0\right),{\theta }_{PA}\left(T_{PRI}\right),...,{\theta }_{PA}\left((N-1)T_{PRI}\right)],\hfill \end{array}$$

(48)

the covariance matrix ${\mathbf{R}}_{yy}=E\left[{\mathbf{yy}}^H\right]$, where $E[·]$ is the ensemble average, is estimated by applying the modified spatial smoothing preprocessing to decorrelate signals from various frequency components as follows:

$${\mathbf{R}}_{yy}={\displaystyle \frac{1}{2M}}\sum _{k=0}^{M_{sub}-1}({\mathbf{R}}_k+{\mathbf{JR}}_k^{*}\mathbf{J}),$$

(49)

where

$${\mathbf{R}}_k={\mathbf{y}}_k{\mathbf{y}}_k^H,$$

(50)

$${\mathbf{y}}_k={[y_k,y_{k+1},...,y_{k+N_{sub}}-1]}^T,$$

(51)

and

$$\mathbf{J}=\left[\begin{array}{cccc}0& 0& ...& 1\\ 0& ...& 1& 0\\ ⋮& \ddots & ⋮& ⋮\\ 1& 0& ...& 0\end{array}\right].$$

(52)

In (49), $M_{sub}$ is the number of subarrays and $N_{sub}$ is the subarray dimension. (47) is the average ${\mathbf{R}}_k$ obtained by ${\mathbf{y}}_k$ clipped by two windows of size $N_{sub}$ moving forward and backward [26].

Assuming the number of frequency components L, the noise subspace matrix ${\mathbf{S}}_{noise}$ consisting $N_{sub}-L$ eigenvectors of ${\mathbf{R}}_{yy}$ corresponding to the smallest $N_{sub}-L$ eigenvalues is obtained, and the super-resolution spectrum is obtained as

$$sp\left(f\right)={\displaystyle \frac{1}{{\mathbf{v}}^H\left(f\right){\mathbf{S}}_{noise}{\mathbf{S}}_{noise}^H\mathbf{v}\left(f\right)}},$$

(53)

where the direction vector $\mathbf{v}\left(f\right)$ at a frequency f is defined by

$$\mathbf{v}\left(f\right)\left[n\right]=exp\left(j2\pi fn{T}_{PRI}\right),n=0,1,2,...,N_{sub}-1,$$

(54)

In (54), f can be set with infinite resolution, which is a major advantage of the MUSIC algorithm. In this study, we divided the signal into respiratory and cardiac signals and set f around $f_r$ and $f_c$ with a very high resolution, which may require very long observations if the FT is employed.

4. Experimental Results

4.1. Experimental Condition

To demonstrate the efficiency of the proposed method for VSE in the presence of body movement using measurements from general low-cost short-range radar, we used both IR-UWB and FMCW radars. For the IR-UWB radar, we employed the parameters of the X4M03 produced by Novelda, and for the FMCW radar, we used those of the Distance2GoL by Infineon (see Table 1). The radar signals of both radars were modeled, and $f_r$ and $f_c$ were estimated using the proposed method. A time frame of 40 s was used to obtain the spectrum with a 0.025 Hz resolution, and the time frame was moved with a step of 1 s to ensure overlapping neighboring frames.

Measurements were conducted using both the X4M03 and Distance2GoL radars. The subject was positioned approximately 2 m away from the radar sensor. Using the same frame time and frame step as in the experiment, the $f_r$s and $f_c$s of the received signal were estimated. To evaluate the accuracy of these estimations, actual $f_r$s and $f_c$s were measured using a wearable metabolic system, the K5 produced by COSMED, and compared with the estimated values (see Figure 9).

To facilitate computation, we used MATLAB for the experiment and measurement processes. We employed the functions $unwrap(·)$, $eig(·)$, and $butter(·)$ for phase unwrapping, eigenvalue decomposition, and the Butterworth filter, respectively. For envelope detection, the parameters were set as $L=150$ and $▵t_w=$ 6 s for IR-UWB and $L=300$ and $▵t_w=$ 6 s for FMCW. The Butterworth filter parameters were defined as $10/60\le f_1\le 30/60$ and $55/60\le f_2\le 100/60$ for the ranges of RR and CR, respectively. For $K\left(t_p\right)$, $\sigma$ was set to 0.7. The MUSIC algorithm was applied to the zero-padded signal, with zeros equal to the length of the signal added; $M_{sub}$ was set to 801 for IR-UWB and 1601 for FMCW. For each radar, the frequency was sampled at 0.0033 Hz for RR and 0.005 Hz for CR.

In this study, we compared the estimation accuracy of the proposed method with that of four different previous methods cited in [16,20,25], and [35]. The methods are as follows:

• The method in [16] compensated for phase error by detecting constant Doppler shift due to random body movement.
• The conventional method in [20] used a fuzzy rule to mitigate the effects of random movement.
• The method in [25] employed empirical mode decomposition (EMD).
• The method in [35] used VMD to detect vital signs.

4.2. Estimation Accuracy and Robustness of the Proposed Scheme

The experimental results confirm the efficiency of the proposed method (Figure 10 and Figure 11). Entropy minimization successfully aligned the RPs of the IR-UWB and FMCW radars (Figure 10b and Figure 11b). However, the phase of the signal post-PCA was considerably distorted owing to motion, obscuring the sinusoidal patterns of the RR and CR (Figure 10c and Figure 11c).

The proposed MOCOM $\#1$ effectively removed motion artifacts. As the movement was accurately modeled by (38), the phase error was substantially reduced, revealing a clear sinusoidal pattern in the phase components’ envelope (Figure 10d and Figure 11d). Additionally, MOCOM $\#2$ further eliminated residual errors, resulting in an even more pronounced sinusoidal pattern (Figure 10e and Figure 11e).

After applying LPF and HPF, the phase data clearly represented RR and CR (Figure 10f,g, Figure 11f,g). Compared with the combined signal of RR and CR, the sinusoidal pattern of CR was distinctly observable. Ultimately, the MUSIC algorithm accurately estimated both RR and CR (Figure 10h,i, Figure 11h,i).

The robust VSE performance of the proposed method was validated by comparing the estimated performance of the proposed scheme with the four existing techniques, as shown in Figure 12 and Figure 13. Furthermore, to evaluate the accuracy of the proposed scheme quantitatively, we computed the root-mean-square error (RMSE) between the RR and CR measured by the wearable metabolic system and estimated by the two radar systems. The RMSE is calculated as

$$RMSE=\sqrt{{\displaystyle \frac{1}{N_m}}\sum _{i=1}^{N_m}{({\overline{f}}_c\left(i\right)-{\widehat{f}}_c\left(i\right))}^2}$$

(55)

where ${\overline{f}}_c$ and ${\widehat{f}}_c$ are the measured and estimated values, respectively, and $N_m=75$ denotes the number of frames.

The results of the experiments conducted using the two radar systems are summarized in Table 2. Both the proposed and existing methods described in [35] demonstrated accurate estimation of RRs with an RMSE of less than 0.1 Hz in scenarios involving IR-UWB and FMCW radar systems. In contrast, conventional methods outlined in [16,20], and [25] exhibited significantly higher RMSEs when estimating RRs using the IR-UWB radar system. The performance of the proposed method in CR estimation was also significantly superior to that of the conventional methods.

In the Fuzzy rule-based method presented in [20], the initially estimated RRs and CRs were used as threshold values for the membership function in fuzzy logic rules. This method depends solely on fuzzy rules without incorporating any compensation techniques for VSE in subjects with random body movements. Consequently, the accuracy of VSE using this approach is highly dependent on the selected fuzzy logic thresholds. Determining appropriate thresholds using initial estimates is particularly challenging in this experiment, as vital signs were estimated for a subject with random body movements, even in the initial phases of the experiment.

Conventional methods referenced in [16,25,35] identified the desired RR and CR by selecting the first maximum peak in each spectrum. These peaks were isolated using bandpass filters in [16], EMD in [25], and VMD in [35]. However, because these methods do not account for motion compensation, their VSE accuracy is compromised by phase distortions caused by random body movements. Consequently, these conventional methods are inadequate for detecting vital sign in individuals with significant body movement.

5. Conclusions

In this study, we present a novel algorithm designed to enhance VSE accuracy in the presence of body movement by incorporating a motion compensation technique. The proposed scheme follows a four-step process: 1. Radar signal denoising using PCA: This step reduces noise and enhances signal clarity. 2. Phase distortion compensation with two MOCOM methods: These techniques correct distortions caused by body movements. 3. Separation of respiratory and cardiac signals using a bandpass filter: This step isolates the vital signs for a more accurate analysis. 4. VSE using the MUSIC method: This algorithm step ensures a precise estimation of vital signs.

Experimental results with IR-UWB and FMCW radar systems demonstrated that our method provides more accurate and robust VSE than conventional methods, especially in scenarios involving random body movements.

Future studies will focus on the biomechanical analysis of walking humans. We plan to extend our algorithm to estimate not only vital signs but also walking parameters. The presence of various micro-movements related to respiration, heartbeat, and limb movements leads to mutual interference and phase ambiguity in the received signals. To address these challenges, we aim to employ a MIMO radar system.