A Novel Real-Time Human Heart Rate Estimation Method for Noncontact Vital Sign Radar Detection

This paper proposes a novel real-time human heart rate (HR) estimation method for noncontact vital sign radar detection. The proposed method combines the Hough transform based respiratory harmonics suppression algorithm and linear predictive coding (LPC) based HR estimation algorithm. Since respiration signals can cause serious interference to HR estimation, the Hough transform based respiratory harmonics suppression algorithm is first applied to successively identify and filter the respiration signals and the higher order harmonics by their distributions on the 2-D spectrum. After suppression, HR can be derived in real time based on the linear predictor coefficients of the filtered signals. Compared to the traditional method, the proposed method can estimate HR more precisely by reducing the interference of respiratory harmonics and avoiding the problem of insufficient frequency resolution caused by short-time windows. Experimental results are presented to demonstrate the performance of the proposed method.


I. INTRODUCTION
In recent years, the fast development of remote sensing techniques has promoted the wide application of noncontact vital sign measurement in numerous fields, such as the daily monitoring of vital signs [1]- [5], human and animal health care [6]- [9], and disaster relief [10]- [13]. Heart rate (HR) is considered as an important vital sign. The most popular noncontact HR measurement method is using the Doppler radar, which can be realized by detecting the micro motion caused by the human physiological movement through phase modulation effect. However, the use of the Doppler radar causes two challenges in detecting human HR. The first is how to deal with the harmonics interference caused by respiration signals. When heartbeat signals and respiratory harmonics are close or overlapping in the frequency domain, the respiratory harmonics may be mistaken as the heartbeat signals. The second is how to achieve HR extraction in real time so that it can be more applicable and practical.
The associate editor coordinating the review of this manuscript and approving it for publication was Seung-Hyun Kong .
This requires HR to be extracted in short time windows with low computation complexity.
To overcome these obstacles, various studies have been conducted on HR acquisition using the Doppler radar. A novel signal processing technique for Doppler radar cardiopulmonary sensing is applied, which uses minimum mean square error to estimate respiratory harmonics thus reducing its interference [14]. A respiratory harmonics cancellation method [15] is applied to avoid the effect of respiratory harmonics in complex signal demodulation [16]. A continuouswavelet filter and ensemble empirical mode decomposition (EEMD) based algorithm is applied for cardiopulmonary signal recovery and separation so that more accurate HR can be obtained in the time domain [17]. A leakage cancellation technique is used for the long-range detection of human HR, using arctangent demodulation without dc offset compensation [18]. The algorithms of the above studies well deal with the interference of respiratory harmonics, but all require a long time window to achieve HR acquisition. A time-window-variation technique, which can estimate HR by using 2-5 s short-period time windows, is applied to fast HR acquisition [19]. However, the measurement accuracy VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ of this algorithm decreases when the heartbeat signals and respiratory harmonics are close or overlapping in the frequency domain. Based on [19], a wavelet-transform-based data-length-variation technique is explored, which can distinguish heartbeat signals from respiratory harmonics and estimate HR with 3-5 s data length [20]. Nevertheless, this method is computationally intensive since the high-resolution wavelet frequency spectrum requires multiple wavelet transforms. The goal of this paper is to explore a novel real-time HR estimation method for long-time HR monitoring, using Hough transform and linear predictive coding (LPC) algorithm, which can both suppress respiratory harmonics and achieve HR estimation in a shorter time window with a less intensive computation. When the breathing of subjects is excessively strong, the heartbeat signals will be covered by respiratory harmonics in the frequency domain, which leads to the inaccurate estimation of HR. Therefore, a respiratory harmonics suppression method using the 2-D spectrum is presented to distinguish respiratory harmonics and the heartbeat signal in the frequency domain. In the 2-D spectrum, the respiratory harmonics and heartbeat signal, which are close or overlapping in the frequency domain, can be distinctly separated by their different rates of change of frequency. In the conventional spectrum analysis methods, the frequency resolution is proportional to the length of time window. The shorter length of time window may lead to lower frequency resolution. However, most real-time applications require HR acquisition in a very short time window, which means the results of the traditional methods may not be practical. To solve this problem, the LPC based HR derivation algorithm is applied to estimate HR with high frequency resolution in each 3 s time window. This method uses the M -order linear prediction coefficients to derive the parameters of the M components of the signals, thereby estimating the real-time HR. The experimental data are collected from the Doppler radar vital sign detection system. Experimental results are presented to demonstrate the performance of the proposed method.
The paper is organized as follows: In Section II, the architecture of the Doppler radar vital sign detection system and Doppler radar signal model are introduced. Section III proposes a respiratory harmonics suppression method and uses LPC to estimate human heart rate. In Section IV, experimental results illustrate the performance of the proposed method and the conclusion is made in Section V.

II. DOPPLER RADAR
The Doppler radar vital sign detection system acquires breathing and heartbeat information by detecting chest-wall movement. Fig. 1 shows the architecture of the Doppler radar vital sign detection system. It transmits a single frequency sine wave continuously towards a subject and obtains information of chest-wall movement through the phase modulation of the reflected wave. The single frequency transmitted signal can be expressed as where A is the signal amplitude, f the radar carrier frequency and ϕ 0 the initial phase. After being reflected by the human body, the phase of the transmitted signal is modulated by the chest-wall movement w(t) and the received signal is expressed as where k is the scattering coefficient, d 0 the initial distance between antenna and human body, λ the wavelength of the signal, c the velocity of light, and w(t) can be expressed as [20] w(t) = m h sin(ω h t + ϕ h ) + m r sin(ω r t + ϕ r ) In this formula, m h , ω h , and ϕ h are the amplitude, the angular frequency and the initial phase of the heartbeat movement, respectively. m r , ω r , and ϕ r are the amplitude, the angular frequency and the initial phase of the respiratory movement, respectively. After demodulation by using the complex signal demodulation method [16], the baseband signal can be expressed as Process the baseband signal as follows where d is the direct current component. The respiration and heartbeat information is included in the signal s(t).

A. RESPIRATORY HARMONICS SUPPRESSION METHOD
Most of the Doppler radar vital sign detection methods estimate HR in the frequency domain. But there are many challenges in real-time HR acquisition by using the traditional time-frequency analysis methods. According to (3), chestwall movement is related to breathing rate, breathing amplitude, heart rate, and heartbeat amplitude. The amplitude of breathing signals ranges from 4 to 12 mm and the normal breathing rate is 5-25 beats per minute(bpm) [21]. The amplitude of heartbeat signals ranges from 0.2 to 0.5 mm and the normal heart rate is 50-120 bpm [22]. Since the breathing amplitude is much larger than the heartbeat amplitude, respiratory harmonics can interfere with the identification of the heartbeat signals. It's hard to distinguish between respiratory harmonics and heartbeat signals when their frequencies are close. In addition, respiration and heartbeat signals are not standard sinusoidal signals and their frequencies are variable. Therefore, respiration and heartbeat signals have broad bandwidths in the frequency domain, which may cause the heartbeat signals and respiratory harmonics to overlap in the frequency domain. When breathing is strong, harmonics can even cover the heartbeat signals, which can reduce the accuracy of HR estimation. Since the bases of the Fourier transform (FT) are sine waves with a fixed frequency, it is not suitable for processing respiration and heartbeat signals with varying frequency. Therefore, we use Hough transform to process the signals. The capability of the Hough transform to accumulate signal energy along a preset frequency trajectory provides an effective approach to identifying and extracting components needed [23], [24]. Compared with the Fourier transform, the Hough transform can make energy more concentrated. Demodulate signal s(t) with the frequency model F(α, t) and the Fourier transform is then applied as where the frequency model F(α, t) can be expressed as As shown in Fig. 2, α is the angle between the linear frequency model and the x axis, and t end is duration of signal s(t). As angle α increases, the linear frequency model F(α, t) rotates counterclockwise around the midpoint (t end /2, 0). According to (6), Hough transform accumulates energy along the linear frequency model F(α, t). Assume the frequencies of the respiration signal and the second harmonic are changing as shown by the curves in Fig. 2, and their linear fitting results are shown as dashed lines F 1 and F 2 , respectively. The ordinates of the midpoints of F 1 and F 2 are f 1 and f 2 , respectively. When α = α 1 , the frequency model F(α 1 , t) is parallel to F 1 , which means that the respiration signal can generate maximum energy accumulation along frequency model F(α 1 , t). Since F 1 can be obtained by translating  F(α 1 , t) by distance f 1 , the peak point of energy accumulation of the respiration signal is at S(α 1 , f 1 ). Similarly, when α = α 2 , the second harmonic will peak at the point S(α 2 , f 2 ). According to the properties of harmonics [15], the peak point S(α i , f i ) of i-th harmonic should approximately satisfy the relationship as The simulation results in Fig. 3 illustrate this feature. In order to make respiratory harmonics more obvious, the simulated chest-wall movement caused by breathing is modeled by [25]. The respiration rate increases from 15 to 21 bpm. The simulated chest-wall movement caused by heartbeat is modeled by [26]. The heart rate decreases from 66 to 54 bpm and the amplitude is one tenth of the respiration signal. The length VOLUME 8, 2020 of the simulation data is 20 s. The Fourier transform of the signal is shown in Fig. 3(a). Although the peak of respiration signal is easy to identify, it is difficult to distinguish the second harmonic, the third harmonic, and the heartbeat signal which have multiple peaks. The Hough transform is shown in Fig. 3(b) when α = 0.005. The main peak of respiration signal (f = 0.3 Hz) is more concentrated than in Fig. 3(a), which has a higher amplitude and a narrower bandwidth. The Hough transform is shown in Fig. 3(c) when α = 0.01. It is easy to see that the second harmonic has a distinct peak (f = 0.58 Hz). The Hough transform is shown in Fig. 3(d) when α = 0.015. It is also easy to see that the energy of third harmonic is obviously more concentrated when f = 0.92 Hz.
The simulation results show that when the rotation angle of the frequency model F(α, t) is α i , the peak of i-th harmonic dose appear near S(α i , f i ). Since the frequency of respiration signals and heartbeat signals are both varying parameters and their rates of change of frequency are different, their peaks will appear in different positions. We can use this feature to identify and filter respiration signals and heartbeat signals.
According to the position (α, f ) of peak, a band stop filter can be applied to suppress the corresponding signal components as where s (t) is the filtered signal and f B (f ) is the filter frequency response. To suppress the thermal noise interference of Doppler radar, Gaussian filter is applied in this study to suppress components which can be expressed as where f is the frequency of peak and σ is a scale coefficient which can be used to adjust the filter band-width. In order to illustrate this feature, the Hough transform of the simulation signal is calculated multiple times with α as the independent variable. Then, with α as the x axis and f as the y axis, a 2-D spectrum S(α, f ) is obtained, which is shown in Fig. 4(a). The most obvious peak in Fig. 4(a), located at (0.005 rad, 0.3 Hz), is the respiration signal. The second obvious peak, located at (0.011 rad, 0.6 Hz), is the second respiratory harmonic. Since the third respiratory harmonic and heartbeat signal are too weak, it is hard to identify them in Fig. 4(a). Fig. 4(b), (c), and (d) are 2-D spectrums after filtering the respiration signal, second harmonic, and third harmonic by (9), respectively. As shown in Fig. 4(c), the peak of the third harmonic can be found at (0.014 rad, 0.9 Hz). As shown in Fig. 4(d), the α, f coordinates of peak point are (-0.01 rad, 1 Hz), which match the frequency of the heartbeat signal we set. According to the definition of α, the peak of the increasing frequency signal is on the right side of the 2-D spectrum. In contrast, the peak of the decreasing frequency signal is on the left side. Since the rate of change of frequency of each signal component is different, their peaks will appear in different positions, which can suppress respiratory harmonics more accurately. Since α is small, (8) can be simplified to α i ≈ i · α 1 and f i ≈ i · f 1 . Therefore, the peaks of the respiration signal and higher-order harmonics are distributed at nearly equal intervals with (0 rad, 0 Hz) as the origin. However, the actual chest wall movement is not a standard sine wave. Moreover, because of the influence of insufficient data length and noise interference, the peak distribution does not strictly follow this rule. Therefore, the prediction peak point (α i , f i ) is first obtained by (8), and then, points (α i , f i ), namely the maximum value near (α i , f i ), are selected as the 88692 VOLUME 8, 2020  (6) is based on Fourier transform, the frequency resolution f is inversely proportional to the length of the signal. The frequency resolution can be improved by using the zeros padding method before Fourier transform. In this work, f = 1/60 Hz (1 bpm). The 2-D search is performed only in one of the gray areas at a time instead of being conducted globally, which can greatly reduce the computational cost and ensure the obtaining of results in real time. However, when multiple components are extracted simultaneously, the energy leakage of the strong component will interfere with the parameter estimation of the weak component. In order to improve the performance of the proposed method, the filtering algorithm is improved by combining CLEAN algorithm [27], which can reduce the effect of sidelobes of strong signals on weak signals. The respiration signal and higher-order respiratory harmonics are extracted sequentially rather than simultaneously. After filtering the respiration signal, the second harmonic is then extracted. Next, filter the second harmonic from the original signal, re-estimate the respiration signal and correct its peak point. After filtering the new respiration signal from the original signal, re-estimate the second harmonic and correct its peak point. By analogy, from strong to weak, complete the suppression of each component. The flowchart of the respiratory harmonics suppression method is shown in Fig. 6, where s 1 (t) is the respiration signal, s i (t) is the ithorder respiratory harmonic, s h (t) is the heartbeat signal, and T is the threshold used to control the number of harmonics to be suppressed. When the modulus of the peak point (α i , f i ) is less than the threshold T (|S(α i , f i )| < T ), the suppression of the harmonics is stopped. The threshold T is a priori parameter. The number of suppressed harmonics is usually controlled at 3 or 4.

B. HEART RATE REAL-TIME ESTIMATION
In order to use the above method to identify and suppress respiratory harmonics more accurately, the data length of the signal s(t) needs to be about 20 s. The reasons are shown as follows. First, the resolution of the Hough transform is inversely proportional to the length of time window. When the length of time window is short, the resolution is not enough to recognize respiratory harmonics. Second, when the length of time window is small, each time window contains only few respiration cycles and the rate of change of frequency of the respiration signal is small, which will cause respiratory harmonics to fail to have obvious energy accumulation on the 2-D spectrum. However, real-time HR detection requires that HR should be estimated at least every 3 s. Thus, the main steps of real-time HR detection at t 0 can be summarized as follows: 1. Use the respiratory harmonics suppression method (introduced in Section III-A) to process the signal s(t) (t = t 0 -20 to t 0 ) and then obtain the signal s h (t).
2. Intercept the last 3 s of signal s h (t) to get a new heartbeat signal s h (t) (t = t 0 -17 to t 0 ).

Calculate the frequency of signal s h (t) as the real-time HR at t 0 .
Since the length of time window is too small, the traditional spectral analysis method is limited by the insufficient frequency resolution to accurately obtain the frequency spectrum of signal s h (t). Therefore, the linear predictive coding (LPC) of the signal is used to directly derive the frequency information of the signal, which is not disturbed by the frequency resolution. In addition to the heartbeat signal, s h (t) also contains some unsuppressed respiratory harmonic components and noise, so data sequence s h (n) can be expressed as where A i is amplitude, f i the frequency, T s the sampling period, θ i the initial phase, and M the number of components. According to the information prediction theory, the forward and backward linear prediction model of the data sequence s h (n) can be expressed as [28], [29] Forward : where M is the predictor order, a M (i) (i = 1, . . . , M ) are the M -order forward predictor coefficients, and a * M (i) are the conjugate of a M (i). Among a number of different LPC algorithms, the Burg-based LPC algorithm, which determines the optimum predictor parameters based on forward and backward average error power, is one of the most costeffective algorithms [30], [31]. According to the Levinson recursive relationship [32], which is shown in (13), the M -order forward and backward prediction errors can be expressed as where E f ,M and E b,M are the M -order forward and backward prediction error, respectively. The predictor reflection coefficient K M can be expressed as (15) where N is the length of data sequence s h (n). Substitute (15) in (14) and the M -order average predictor error power P M can be expressed as From the above formula, in order to solve the M -order predictor coefficients a M (i), first set the initial value as follows Then, according to (15), (13), and (16), estimate the predictor reflection coefficient K 1 , the predictor coefficients a 1 , and the average predictor error power P 1 , respectively. Finally, update the forward and backward prediction error E f ,1 and E b,1 according to (14). Repeat the above steps until a M (i) is obtained. Generally, a higher prediction order would lead to higher prediction accuracy. However, when prediction order is too large, an overfitting issue may occur. Thus, we set a preset threshold τ , which controls the predictor accuracy. Iterative calculation will be stopped when τ satisfies the following relationship The most proper τ usually ranges from 0.1 to 0.3 [31]. The smaller τ is, the more the number of iterations is. The threshold τ in this work is set as 0.3. Burg-based LPC algorithm starts from 0-order and solves a M (i) by iterative calculation.
Since the heartbeat signal s h (n) contains fewer signal components, Burg-based LPC algorithm is suitable for solving its linear predictor coefficients. The HR acquisition method using linear predictor coefficients is introduced as follows.

19) Since s h (n) is obtained by linearly superposing each component, each component satisfies the following relationship
After dividing A k e j[2πf k (n−M −1)T s +θ k ] on both sides and transposition, (20) can be expressed as Then (21) can be expressed as It can be known from (23) that X k (k = 1, . . . , M ) is the M roots of the following unary M -degree equation Therefore, after obtaining the M -order predictor coefficients a M (i) by using Burg-based LPC algorithm, a unary M -degree equation can be constructed based on (24) and its roots are X k . Then, according to (22), f k can be obtained as follows If we substitute (22) in (11), a system of equations can be obtained as follows where Y k = A k e jθ k , k = 1, . . . , M (27) and N is the length of data sequence s h (n). The least squares solution Y LS of (26) can be obtained by the following normal equation where X H is the Hermitian conjugate of matrix X and X + is the Moore-Penrose pseudoinverse matrix of X. Then, according to (27), A k and θ k can be obtained as follows Therefore, the parameters of each component of the signal s h (n) can be derived by the M -order predictor coefficients a M (i). It can be known from (24) that the M -order linear predictor coefficients can solve M signal components. The number of signal components is determined by (18). Finally, the frequency f k of the signal component, whose amplitude A k is the largest, will be taken as real-time HR at t 0 .
The above process presents that the proposed method is not limited by the frequency resolution and can obtain HR in real time with low computational cost.

IV. EXPERIMENTAL RESULTS
To evaluate the performance of the proposed method, a Doppler radar vital sign detection system is constructed in our laboratory. The architecture of Doppler radar vital sign detection system is shown in Fig. 1. Considering the size, portability, and price of the Doppler radar vital sign detection system, it is designed with a single transmit antenna and a single receive antenna. The radar parameters are listed in Table 1. Fig. 7(a) shows the transmitter module and Fig. 7(b) shows the receiver module. The radar receiver samples the received signals, transmits the data to the computer, and then processes the signal on the computer. Fig. 7(c) shows the experimental environment. The subject sat still in front of the radar in a distance of about 0.5-1.5 m. Meanwhile, an ECG sensor was tied under the subject's clothes to collect HR as a reference standard value. Fig. 8 shows one set of measured experimental results. The subject breathed calmly and the influence of respiratory harmonics on HR estimation was not severe. Fig. 8(a) shows the 2-D spectrum of the original signal. x, y, z, and { are respectively the peaks of the respiration signal, the second, the third harmonic, and the heartbeat signal. The third harmonic and heartbeat signal, which have similar frequencies, are clearly separated in the 2-D spectrum. Thus, the respiratory harmonics can be suppressed according to the flowchart shown in Fig. 6. And then, the real-time HR is calculated according to the three steps described in the first paragraph of Section III-B. The final frequency derivation results are shown by the red lines in Fig. 8(b), and the frequency with the largest amplitude (60 bpm) is selected as HR. The black triangle presents the reference HR, and it is consistent with the estimated results of the proposed method. The blue line shows Fourier transform of the original signal with a 3 s time window. On the line, it is confusing to decide the HR between the peak | and peak }. According to the reference HR, the peak | is the combined result of the heartbeat and the respiratory harmonic. It is not accurate due to the influence of the respiratory harmonic and the overly short time window. Moreover, the peak } can be identified as the respiratory harmonic whose amplitude can be reduced by the proposed method in this study. Fig. 9 shows HR estimation results versus time with different estimation methods. The total length of the data is 60 s and HR is estimated each 3 s. The black line shows the reference HR. The blue line presents the results of the proposed method, which is close to the reference HR. The magenta line shows the experimental results of the time-window-variation method [19], and the error is also small. It should be noted that the subject's breathing was unstable between 18 s to 30 s. This makes the results of time-window-variation method fluctuate. The green line shows the results of the Wavelet transform (WT). In order to improve the estimation accuracy, the WT algorithm compared in this study uses a length of previous data, which, however, causes computational complexity. The red line shows the results the traditional Fourier transform and the error is large. Table 2 shows the experimental errors of different methods in the case of calm breathing with more measured data. The mean relative error (MRE) and rootmean-squared error (RMSE) are used to quantify the accuracy of different methods, which is expressed as (32) where N is the number of HR data, HR meas the estimated HR, and HR ref the reference HR. Compared with the other three methods, the proposed method reduces the mean relative error to 2.1%. In addition, Table 2 shows the average time required for HR estimation by each method in a 3 s time window. The processing time is obtained by using MATLAB R2014b. As shown in Table 2, WT requires a long processing time, which cannot meet the real-time requirements. Although the processing time of the proposed method is not as short as the FT and time-window-variation method, its computational complexity still satisfies the realtime requirements.
In order to reflect the performance of the proposed method more clearly, the subjects were required to take strong  breathes deliberately so that the respiratory harmonics could cause severe interference. Fig. 10 shows one set of experimental results. Fig. 10(a) shows the 2-D spectrum of the original signal. x, y, z, and { are respectively the peaks of the respiration signal, the second, the third, and the fourth  harmonic. Their position distribution conforms to (8). Since HR changes greatly, it has no obvious peak. The final frequency derivation results are shown in Fig. 10(b). It can be seen from the blue line that since the breathing is too strong, the heartbeat signal is covered by respiratory harmonics, and only two peaks | and } can be identified. However, after suppression, the proposed method can more accurately estimate the HR. Fig. 11 shows the results of different HR estimation methods. Compared with the other three methods, the proposed method presents a more accurate result. Table 3 shows the experimental errors of different methods in the case of strong breathing with more measured data. Compared with the other three methods, the proposed method reduces the mean relative error to 3.5%.

V. CONCLUSION
In this paper, a real-time human HR estimation method is proposed for noncontact vital sign radar detection. Hough transform is used to suppress respiratory harmonics. The proposed 2-D spectrum method can significantly distinguish between the heartbeat signal and higher-order harmonics. Based on this, a respiratory harmonics suppression process combining CLEAN algorithm is proposed. Then the LPC algorithm is applied to estimate HR in real time. The proposed derivation process can use the linear predictor coefficients of the filtered signal to obtain the HR, which can avoid the problem of insufficient resolution caused by short-time windows. Combining the results of two experiments, the proposed method can accurately estimate the heart rate both in the cases of calm and strong breathing. Compared with other methods, the proposed method has less error, which demonstrates its good performance. Overall, the superiority of the proposed HR estimation method indicates considerable potential in real-time HR detection applications. From 2012 to 2015, he was a Lecturer with Central South University, Changsha, where he is currently an Associate Professor. His current research interests include target detecting, tracking, and identification methods for radar applications.
XIALI YU was born in Shenzhen, China, in 1995. He received the B.S. degree in optoelectronics information science and engineering from Central South University, Changsha, China, where he is currently pursuing the master's degree with the School of Physics and Electronics, with a focuses on the research of radar signal processing and target detection.