Heart rate estimation method based on millimeter wave radar

Abstract

The present invention provides a heart rate estimation method based on millimeter wave radar. The method comprises the following steps: sampling a human body by using a millimeter wave radar to obtain a plurality of frequency modulated continuous wave sweep signals; forming a distance-slow time matrix and locating the distance between the human body and the radar; extracting the human body phase by using an inverse tangent function, unfolding the phase whose phase change exceeds a set threshold by using a phase unwrapping function, obtaining the phase difference of adjacent frames, and smoothing the difference exceeding the threshold by using an interpolation method; removing the noise of the smoothed phase signal by using a wavelet transform, and obtaining a reconstructed phase signal; combining an EEMD-ICA method to separate a heartbeat signal, a respiratory signal, a respiratory harmonic and noise in the reconstructed phase signal, and obtaining a heartbeat waveform; constructing a spatial spectrum function by using a multi-signal classification algorithm, and estimating the heartbeat frequency by searching the spectrum peak; the method can accurately separate the heartbeat signal and greatly improve the accuracy of the heartbeat frequency estimation.

Description

Technical Field

[0001] The invention relates to a heart rate estimation method based on millimeter wave radar, belonging to the technical field of radar detection. Background Art

[0002] Heart rate is an important reference indicator and basis for medical and health care. Traditional heart rate estimation methods mainly use contact wearable sensors or adhesive electrodes to monitor heart rate, which are basically wired devices and have many limitations in use. In addition, contact sensors are usually complicated to operate and the user experience is not good. In response to these problems, researchers proposed a heart rate detection technology based on millimeter-wave radar, which can monitor the user's heart rate over a long distance and contactlessly, making the whole process more convenient and comfortable.

[0003] At present, the use of millimeter wave radar for heart rate estimation mainly includes target positioning, phase extraction, noise removal, signal separation and frequency estimation. However, related research still faces the following difficulties and challenges:

[0004] (1) The reflection of surrounding static objects and the multipath effect in the environment will affect the echo signal. The random movement of the human body will also bring a lot of noise to the data collection.

[0005] (2) The chest displacement caused by breathing is much larger than the heartbeat, and its harmonics mask the heartbeat waveform, which will affect the accuracy of heartbeat signal separation.

[0006] For noise processing, usually arc tangent demodulation and phase compensation are used in the phase extraction stage, and bandpass filters are used in the signal separation stage. However, they are not effective in suppressing noise near the heartbeat frequency band, and cannot solve the problem of the second and third harmonics of breathing masking the heartbeat signal.

[0007] To solve the problem of respiratory harmonics, the Empirical Mode Decomposition (EMD) method can be used to separate the heartbeat signal in the signal separation stage. However, its theoretical system is not yet mature, and in the process of generating the Intrinsic Mode Functions (IMF), adjacent IMF component waveforms often alias, and the method's anti-noise ability is also average.

[0008] The above problems are issues that should be considered and solved in the process of heart rate estimation based on millimeter wave radar. Summary of the invention

[0009] The purpose of the present invention is to provide a heart rate estimation method based on millimeter wave radar to solve the problem that noise and respiratory harmonics existing in the prior art affect heartbeat waveform extraction and heart rate estimation, and the accuracy of heart rate estimation needs to be improved.

[0010] The technical solution of the present invention is:

[0011] A heart rate estimation method based on millimeter wave radar comprises the following steps:

[0012] S1. Use millimeter wave radar to sample the human body and obtain M frequency modulated continuous wave sweep signals S(1), S(2), …, S(k), …, S(M), where the kth sweep signal S(k) = [s1(k), s2(k), …, s N (k)], 1≤k≤M, N is the number of sampling points of a frequency sweep signal;

[0013] S2. Apply the distance-fast Fourier transform to each sampling point of the M frequency modulated continuous wave sweep signals S(1), S(2), ..., S(k), ..., S(M), obtain the corresponding distance unit for each sampling point, form a distance-slow time matrix, and locate the distance from the human body to the radar;

[0014] S3. After determining the distance between the human body and the radar, use the inverse tangent function to extract the human body phase, use the phase unwrapping function to unwrap the phase whose phase change exceeds the set threshold, make a difference between the unwrapped phases to obtain the phase difference between adjacent frames, use the interpolation method to smooth the difference that exceeds the threshold, and obtain the phase difference between adjacent frames after smoothing;

[0015] S4, using wavelet transform to remove noise from the smoothed phase signal, obtaining denoised wavelet coefficients, performing inverse wavelet transform to reconstruct the signal, and obtaining a reconstructed phase signal;

[0016] S5, combining the ensemble empirical mode decomposition and independent component analysis method, namely, the EEMD-ICA method, to separate the heartbeat signal, the respiratory signal, the respiratory harmonics and the noise in the reconstructed phase signal to obtain the heartbeat waveform;

[0017] S6. For the obtained heartbeat waveform, a multi-signal classification algorithm is used to construct a spatial spectrum function, and the heartbeat frequency is estimated by spectrum peak search.

[0018] Further, in step S1, M frequency modulated continuous wave sweep signals are obtained, specifically,

[0019] S11, millimeter wave radar transmission signal X t (t) and the received signal X r (t) are:

[0020] X t (t) = sin(2πf min +πSt 2 ) (1)

[0021] X r(t) = sin[2πf min (t-τ)+πS(t-τ) 2 ] (2)

[0022] Among them, f min is the initial frequency of the millimeter wave radar, S is the frequency growth slope, and τ is the signal round trip time;

[0023] S12, the mixer performs I / Q mixing and low-pass filtering on the transmit signal and the receive signal, and outputs the intermediate frequency signal Y(t):

[0024]

[0025] Among them, f min is the initial frequency of the millimeter wave radar, S is the frequency growth slope, and τ is the signal round trip time;

[0026] S13. Sample the intermediate frequency signal at equal intervals to obtain swept frequency signals S(1), S(2), ..., S(k), ..., S(M).

[0027] Further, in step S2, a distance-fast Fourier transform is applied to each sampling point of the M frequency modulated continuous wave sweep signals S(1), S(2), ..., S(k), ..., S(M), and each sampling point obtains a corresponding distance unit to form a distance-slow time matrix, and locates the distance from the human body to the radar, specifically,

[0028] S21, applying distance-fast Fourier transform to each sampling point of the M frequency modulated continuous wave swept frequency signals S(1), S(2), ..., S(k), ..., S(M), and obtaining a corresponding distance unit for each sampling point. The distance unit is a complex signal, and the distance units are put into a matrix by column to form a distance-slow time matrix;

[0029] S22, for each row in the distance-slow time matrix, i.e., each distance unit, the influence of static clutter and DC component is eliminated by subtracting the average value, and the distance unit du after removing the average value is obtained. pq `:

[0030] du pq `=du pq –mean(p) (4)

[0031] Among them, pq is the complex signal in the p-th row and q-th column of the distance-slow time matrix, 1≤q≤M, mean(p) is the average value of the complex signal in the p-th row of the distance-slow time matrix;

[0032] S23, the distance unit du after removing the average value pq`Construct a new distance-slow time matrix, obtain the summation result of each row of the new distance-slow time matrix, and use the distance unit with the largest summation result as the distance from the human body to the radar.

[0033] Further, in step S3, after determining the distance from the human body to the radar, use the arctangent function to extract the human body phase, use the phase unwrapping function to unwrap the phase with a phase change exceeding the set threshold, take the difference of the unwrapped phases to obtain the adjacent frame phase difference, and use the interpolation method to smooth the difference exceeding the threshold to obtain the smoothed adjacent frame phase difference. Specifically,

[0034] S31. After determining the distance from the human body to the radar, the arctangent function can be used for the complex signal of the distance unit where the human body is located to extract the human body phase φ(n):

[0035]

[0036] where I(n) is the imaginary part of the nth frame complex signal in the distance unit where the human body is located, and R(n) is the real part of the nth frame complex signal in the distance unit where the human body is located;

[0037] S32. Use the phase unwrapping function to unwrap the phase with a phase change exceeding ±π:

[0038]

[0039] where φ(n + 1) is the phase of the (n + 1)th frame;

[0040] S33. Calculate the adjacent frame phase difference φ diff (n):

[0041] φ diff (n) = φ(n + 1) - φ(n) (7)

[0042] where 0 < n < M, φ(n) is the phase of the nth frame, and φ(n + 1) is the phase of the (n + 1)th frame;

[0043] S34. If the absolute value of the adjacent frame phase difference φ diff (n) exceeds the set threshold, use the interpolation method for smoothing to obtain the smoothed adjacent frame phase difference φ diff (n):

[0044]

[0045] where φ diff (n - 1) is the (n - 1)th phase difference, and φ diff (n + 1) is the (n + 1)th phase difference.

[0046] Further, in step S4, wavelet transform is used to remove the smoothed adjacent frame phase difference signal φ diff (n)`, obtain the denoised wavelet coefficients, perform inverse wavelet transform to reconstruct the signal, and obtain the reconstructed phase signal, specifically,

[0047] S41, using wavelet transform to smooth the adjacent frame phase difference φ from multiple scales diff (n)` is decomposed to obtain the multi-scale wavelet transform coefficient w;

[0048] S42, perform soft threshold denoising on the multi-scale wavelet transform coefficient w to obtain the denoised wavelet coefficient w λ :

[0049]

[0050] Among them, λ wt is a given threshold, sgn() is a sign function;

[0051] S43, the denoised wavelet coefficients w λ Perform inverse wavelet transform to reconstruct the signal and obtain the reconstructed phase signal φ re (n).

[0052] Further, in step S5, the heartbeat signal, respiratory signal, respiratory harmonics and noise in the reconstructed phase signal are separated by combining the ensemble empirical mode decomposition and independent component analysis method, namely, the EEMD-ICA method, to obtain a heartbeat waveform, specifically,

[0053] S51, use ensemble empirical mode decomposition EEMD to reconstruct the phase signal φ re (n) performing signal separation to obtain a plurality of intrinsic mode functions (IMF) components including heartbeat signals, respiratory signals, respiratory harmonics, and noise, and obtaining a component set;

[0054] S52. Use independent component analysis method ICA to analyze the obtained component set to obtain a heartbeat waveform.

[0055] Furthermore, in step S51, the reconstructed phase signal φ is decomposed using EEMD. re (n) performing signal separation to obtain multiple IMF components including heartbeat signals, respiratory signals, respiratory harmonics, and noise, specifically,

[0056] S511, adding white noise with standard normal distribution to the reconstructed phase signal φ for K times re (n), we get the new signal after adding white noise for the i-th time

[0057]

[0058] Among them, n i (n) represents the i-th added white noise sequence, i = 1, 2, ..., K;

[0059] S512, adding white noise to the new signal obtained each time Perform empirical mode decomposition (EMD) to decompose the signal into a finite number of IMF components to separate the noise caused by unconscious micro-movements of the human body from the heartbeat signal. Assuming that each EMD decomposition obtains P IMF components, K decompositions will obtain K groups of IMF component sets, and each group contains P IMF components.

[0060] S513. According to the principle that the statistical average of unrelated sequences is zero, corresponding components are selected from the K groups of IMF component sets and averaged respectively to obtain a group of P component sets IMFS(n).

[0061] Further, in step S52, the obtained component set IMFS(n) is analyzed using the independent component analysis method ICA to obtain the heartbeat waveform, specifically,

[0062] S521, establishing an ICA model for the component set obtained in step S51:

[0063] X(n)=AE(n) (11)

[0064] Among them, X(n) is the signal observation value, that is, the component set IMFS(n), E(n) is the source signal to be estimated, which is composed of the real heartbeat waveform, respiratory waveform, respiratory harmonics and noise signal, and A is the mixing matrix;

[0065] S522, reversely estimate the source signal E(n) through the observed signal, and establish an ICA solution model:

[0066] Eest=WX(n) (12)

[0067] Where Eest is the estimated value of the source signal E(n), which consists of multiple components, and W is the unmixing matrix;

[0068] S523: Calculate the estimated value Eest of the source signal E(n) by solving the unmixing matrix W, and then extract the heartbeat waveform therefrom.

[0069] Further, in step S523, the estimated value Eest of the source signal E(n) is calculated by solving the unmixing matrix W, and then the heartbeat waveform is extracted therefrom, specifically,

[0070] S5231, zero-meaning the component set IMFS(n), that is, subtracting the average value of each IMF component in the component set from its own, to obtain a zero-meaning component set;

[0071] S5232, whitening the component set subjected to zero mean processing to obtain a component set subjected to zero mean processing and whitening processing;

[0072] S5233, using the maximum likelihood method to estimate each row of the unmixing matrix W for the component set IMFS(n) after zero mean and whitening, and multiplying them with the component set IMFS(n) respectively, the product Gaussian uncorrelation is maximized, and it is put into the estimated value Eest of the source signal E(n), and finally a set of independent components is obtained;

[0073] S5234. Perform fast Fourier transform on each independent component in the estimated value Eest of the source signal E(n) to obtain the frequency spectrum of each independent component. Then calculate the total frequency band energy of each independent component and the heartbeat frequency band energy in the set frequency range in the frequency domain, and calculate the ratio of the heartbeat frequency band energy to the total frequency band energy respectively, and take the independent component with the largest ratio as the heartbeat waveform.

[0074] The beneficial effects of the present invention are as follows: the heart rate estimation method based on millimeter wave radar can accurately separate the heartbeat signal and greatly improve the accuracy of heartbeat frequency estimation. By adopting the EEMD-ICA method, the signal is adaptively decomposed into a finite number of IMF components according to the characteristics of the input signal, and the waveform aliasing problem of EMD is solved by adding white noise during decomposition. The noise can be effectively removed on the basis of ensuring the integrity of the signal in the time domain and frequency domain, and the heartbeat signal, breathing signal, breathing harmonics and noise signal can be accurately and effectively separated, and the influence of breathing harmonics and noise on heartbeat waveform extraction and heart rate estimation can be overcome, effectively improving the recognition accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 is a flowchart of a heart rate estimation method based on millimeter wave radar according to an embodiment of the present invention;

[0076] Figure 2 is a schematic diagram for explaining how to use wavelet transform to remove noise from a smoothed phase signal in an embodiment;

[0077] Figure 3 3 is a schematic diagram illustrating the use of the EEMD-ICA method to separate the reconstructed phase signal to obtain a heartbeat waveform in an embodiment. DETAILED DESCRIPTION

[0078] The preferred embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0079] Example

[0080] A heart rate estimation method based on millimeter wave radar, such as Figure 1 , including the following steps,

[0081] S1. Use millimeter wave radar to sample the human body and obtain M frequency modulated continuous wave sweep signals S(1), S(2), …, S(k), …, S(M), where the kth sweep signal S(k) = [s1(k), s2(k), …, s N (k)], 1≤k≤M, N is the number of sampling points of a frequency sweep signal.

[0082] In step S1, M frequency-modulated continuous wave sweep signals are obtained, specifically,

[0083] S11, millimeter wave radar transmission signal X t (t) and the received signal X r (t) are:

[0084] X t (t) = sin(2πf min +πSt 2 ) (1)

[0085] X r (t) = sin[2πf min (t-τ)+πS(t-τ) 2 ] (2)

[0086] Among them, f min is the initial frequency of the millimeter wave radar, S is the frequency growth slope, and τ is the signal round trip time;

[0087] S12, the mixer performs I / Q mixing and low-pass filtering on the transmit signal and the receive signal, and outputs the intermediate frequency signal Y(t):

[0088]

[0089] Among them, f min is the initial frequency of the millimeter wave radar, S is the frequency growth slope, and τ is the signal round trip time;

[0090] S13. Sample the intermediate frequency signal at equal intervals to obtain swept frequency signals S(1), S(2), ..., S(k), ..., S(M).

[0091] S2. Apply distance-fast Fourier transform to each sampling point of the M frequency modulated continuous wave sweep signals S(1), S(2), ..., S(k), ..., S(M), and obtain the corresponding distance unit for each sampling point to form a distance-slow time matrix, and locate the distance from the human body to the radar.

[0092] S21, applying distance-fast Fourier transform to each sampling point of the M frequency modulated continuous wave swept frequency signals S(1), S(2), ..., S(k), ..., S(M), and obtaining a corresponding distance unit for each sampling point. The distance unit is a complex signal, and the distance units are put into a matrix by column to form a distance-slow time matrix;

[0093] S22, for each row in the distance-slow time matrix, i.e., each distance unit, the influence of static clutter and DC component is eliminated by subtracting the average value, and the distance unit du after removing the average value is obtained. pq `:

[0094] du pq `=du pq –mean(p) (4)

[0095] Among them, pq is the complex signal in the p-th row and q-th column of the distance-slow time matrix, 1≤q≤M, mean(p) is the average value of the complex signal in the p-th row of the distance-slow time matrix;

[0096] In step S22, due to the reflection of static objects and the DC component generated after Fourier transform, a higher power reflection value will be generated on the distance-slow time matrix. Formula (4) is used to eliminate the influence of static clutter and DC.

[0097] S23, the distance unit du after removing the average value pq `Construct a new distance-slow time matrix, obtain the sum of each row of the new distance-slow time matrix, and use the largest distance unit in the sum as the distance from the human body to the radar.

[0098] In step S3, after determining the distance between the human body and the radar, the inverse tangent function is used to extract the human body phase, and the phase unwrapping function is used to unwrap the phase whose phase change exceeds the set threshold. The unwrapped phase is subtracted to obtain the phase difference of adjacent frames, and the interpolation method is used to smooth the difference exceeding the threshold to obtain the phase difference of adjacent frames after smoothing, which is specifically,

[0099] S31. After determining the distance from the human body to the radar, the inverse tangent function can be used to extract the phase φ(n) of the human body from the complex signal of the distance unit where the human body is located:

[0100]

[0101] Wherein, I(n) is the imaginary part of the complex signal of the nth frame in the distance unit where the human body is located, and R(n) is the real part of the complex signal of the nth frame in the distance unit where the human body is located;

[0102] S32. Use the phase unwrapping function to unwrap the phase of the phase change exceeding ±π:

[0103]

[0104] Among them, φ(n + 1) is the phase of the (n + 1)-th frame;

[0105] S33. Calculate the phase difference φ diff (n) between adjacent frames:

[0106] φ diff (n) = φ(n + 1) - φ(n) (7)

[0107] Among them, 0 < n < M, φ(n) is the phase of the n-th frame, and φ(n + 1) is the phase of the (n + 1)-th frame;

[0108] S34. If the absolute value of the phase difference φ diff (n) between adjacent frames exceeds the set threshold, then use the interpolation method for smoothing processing to obtain the smoothed phase difference φ diff (n)` between adjacent frames:

[0109]

[0110] Among them, φ diff (n - 1) is the (n - 1)-th phase difference value, and φ diff (n + 1) is the (n + 1)-th phase difference value.

[0111] S4. Use wavelet transform to remove the noise of the smoothed phase signal, obtain the denoised wavelet coefficients, and perform inverse wavelet transform to reconstruct the signal to obtain the reconstructed phase signal. Such as Figure 2 :

[0112] S41. Use wavelet transform to decompose the smoothed phase difference φ diff (n)` between adjacent frames at multiple scales to obtain the multi-scale wavelet transform coefficients w;

[0113] S42. Perform soft threshold denoising processing on the multi-scale wavelet transform coefficients w to obtain the denoised wavelet coefficients w λ :

[0114]

[0115] Among them, λ wt is the given threshold, and sgn() is the sign function;

[0116] In step S42, after wavelet decomposition, most of the wavelet coefficients with larger absolute values are useful signals, while the coefficients with smaller amplitudes are generally noise. Therefore, when the absolute value of the wavelet coefficient is less than λ wt , directly set it to zero, which can eliminate the influence of noise. The wavelet coefficients w λ obtained after soft threshold denoising have better overall continuity.

[0117] S43, the denoised wavelet coefficients w λ Perform inverse wavelet transform to reconstruct the signal and obtain the reconstructed phase signal φ re (n). Thus, the reconstructed phase signal φ re (n) It removes the noise and retains the original heartbeat component.

[0118] S5. Combining the ensemble empirical mode decomposition and independent component analysis method, namely the EEMD-ICA method, separates the heartbeat signal, respiratory signal, respiratory harmonics and noise in the reconstructed phase signal, suppresses the influence of respiratory harmonics, and obtains the heartbeat waveform. Figure 3 :

[0119] S51, use ensemble empirical mode decomposition EEMD to reconstruct the phase signal φ re (n) performing signal separation to obtain a plurality of intrinsic mode functions (IMF) components including heartbeat signals, respiratory signals, respiratory harmonics, and noise, and obtaining a component set;

[0120] S511, adding white noise with standard normal distribution to the reconstructed phase signal φ for K times re (n), we get the new signal after adding white noise for the i-th time

[0121]

[0122] Among them, n i (n) represents the i-th added white noise sequence, i = 1, 2, ..., K;

[0123] S512, adding white noise to the new signal obtained each time Perform empirical mode decomposition (EMD) to decompose the signal into a finite number of IMF components to separate the noise caused by unconscious micro-movements of the human body from the heartbeat signal. Assuming that each EMD decomposition obtains P IMF components, K decompositions will obtain a total of K groups of IMF component sets, and each group contains P IMF components.

[0124] In step S512, EMD can adaptively decompose the signal into a finite number of IMF components according to the characteristics of the signal itself, thereby separating the noise caused by unconscious micro-movements of the human body from the heartbeat signal.

[0125] S513. According to the principle that the statistical average of unrelated sequences is zero, corresponding components are selected from the K groups of IMF component sets and averaged respectively to obtain a group of P component sets IMFS(n).

[0126] The above process is the process of EEMD. Compared with the single EMD processing, EEMD avoids the waveform aliasing problem of adjacent IMF components in EMD.

[0127] S52. Use independent component analysis method ICA to analyze the obtained component set to obtain a heartbeat waveform.

[0128] S521, establishing an ICA model for the component set obtained in step S51:

[0129] X(n)=AE(n) (11)

[0130] Among them, X(n) is the signal observation value, that is, the component set IMFS(n), E(n) is the source signal to be estimated, which is composed of the real heartbeat waveform, respiratory waveform, respiratory harmonics and noise signal, and A is the mixing matrix;

[0131] S522, reversely estimate the source signal E(n) through the observed signal, and establish an ICA solution model:

[0132] Eest=WX(n) (12)

[0133] Where Eest is the estimated value of the source signal E(n), which consists of multiple components, and W is the unmixing matrix;

[0134] S523: Calculate the estimated value Eest of the source signal E(n) by solving the unmixing matrix W, and then extract the heartbeat waveform therefrom.

[0135] S5231, zero-meaning the component set IMFS(n), that is, subtracting the average value of each IMF component in the component set from its own, to obtain a zero-meaning component set;

[0136] S5232, whitening the component set after zero mean processing, reducing the correlation between the features of the data, and all the features have the same variance, to obtain a component set after zero mean processing and whitening processing;

[0137] S5233, using the maximum likelihood method to estimate each row of the unmixing matrix W for the component set IMFS(n) after zero mean and whitening, and multiplying them with the component set IMFS(n) respectively, the product Gaussian uncorrelation is maximized, and it is put into the estimated value Eest of the source signal E(n), and finally a set of independent components is obtained;

[0138] In step S5233, since the independence of the signal is equivalent to Gaussian irrelevance, to obtain an independent source signal, we can start from obtaining the estimated value Eest of the Gaussian irrelevance source signal E(n).

[0139] S5234. Perform FFT on each independent component in the estimated value Eest of the source signal E(n) to obtain the frequency spectrum of each independent component. Then calculate the total frequency band energy of each independent component and the heartbeat frequency band energy in a set frequency range such as 0.8HZ to 3HZ in the frequency domain, and calculate the ratio of the heartbeat frequency band energy to the total frequency band energy respectively, and take the independent component with the largest ratio as the heartbeat waveform.

[0140] In step S5, the EEMD method is introduced to adaptively decompose the signal into a finite number of IMF components according to the characteristics of the signal itself, and the waveform aliasing problem of EMD is solved by adding white noise during decomposition. After that, the ICA is used to effectively remove the noise, and the heartbeat signal can be effectively separated, and the integrity of the heartbeat signal can be maintained in both the frequency domain and the time domain.

[0141] S6. For the obtained heartbeat waveform, a multi-signal classification algorithm is used to construct a spatial spectrum function, and the heartbeat frequency is estimated by spectrum peak search.

[0142] This heart rate estimation method based on millimeter wave radar can accurately separate the heartbeat signal and greatly improve the accuracy of heart rate frequency estimation. By adopting the EEMD-ICA method, the signal is adaptively decomposed into a finite number of IMF components according to the characteristics of the input signal. By adding white noise during decomposition, the waveform aliasing problem of EMD is solved, and the noise can be effectively removed on the basis of ensuring the integrity of the signal in the time domain and frequency domain. The heartbeat signal, breathing signal, breathing harmonics and noise signal can be accurately and effectively separated, and the influence of breathing harmonics and noise on heartbeat waveform extraction and heart rate estimation can be overcome, effectively improving the recognition accuracy.

[0143] This kind of heart rate estimation method based on millimeter wave radar, in the target positioning stage, for the reflection of static objects, and the reflection of higher power generated by the DC component after Fourier transform, the present invention eliminates the influence by subtracting the average value at each distance, which can improve the accuracy of target positioning. In the phase extraction stage, the present invention uses arc tangent demodulation to extract the phase, uses differential phase to intuitively restore the signal of human chest vibration, and introduces interpolation method to smooth the excessive difference, so as to obtain a stable phase signal and improve the accuracy of subsequent respiratory heartbeat signal separation. In the noise removal stage, the present invention uses wavelet transform multi-scale decomposition of the phase signal, and on the basis of retaining the useful part of the signal, the noise in the phase can be removed, further improving the accuracy of subsequent heartbeat signal separation. In the signal separation stage, the present invention combines EEMD and ICA algorithms, firstly performs collective empirical mode decomposition on the processed phase signal, obtains multiple IMF components containing respiratory heartbeat components, respiratory harmonics, and noise, and then uses ICA algorithm to separate useful signals and noise therefrom, so as to improve the signal-to-noise ratio of the signal and improve the accuracy of the final result.

[0144] This heart rate estimation method based on millimeter wave radar introduces the Ensemble Empirical Mode Decomposition (EEMD) and Independent Component Correlation Algorithm (ICA) method, combines the two, and proposes the EEMD-ICA algorithm. By using the modal decomposition and independent component analysis methods to process radar echo signals, the heart rate information can be accurately estimated.

[0145] The results of the experimental verification of the heart rate estimation method based on millimeter wave radar in the embodiment are as follows:

[0146] In order to verify the performance of the heart rate estimation method based on millimeter wave radar in the embodiment, a comparative test was conducted between the embodiment method and a plurality of existing methods, and the comparative methods include:

[0147] (1) BPF+GI. This method uses a band pass filter (BPF) to separate signals of different frequency bands during signal separation. After calculating the frequency using the fast Fourier transform, the heart rate is obtained using the Gaussian interpolation (GI). (2) IWT, Improved Wavelet Transform. In the signal separation stage, this method resamples the signal, restores the heartbeat signal through wavelet decomposition and reconstruction, uses the fast Fourier transform, and retains the peak value greater than the threshold for calculating the frequency. (3) EMD, Empirical Mode Decomposition (EMD). This method is used in the signal separation stage. (4) EEMD, Ensemble Empirical Mode Decomposition (EEMD).

[0148] On the same data set, the experimental results of the heart rate estimation method based on millimeter wave radar in the embodiment and the existing method are shown in Table 1:

[0149] Table 1 Prediction accuracy of the embodiment method and the existing method

[0150] algorithm Heart rate prediction accuracy BPF+GI 85.08％ IWT 86.71％ EMD 87.17％ EEMD 87.33％ Example Methods 96.48％

[0151] It can be seen from the results in Table 1 that the heart rate prediction accuracy of the embodiment method is 96.48%, which is better than the existing method, verifying the effectiveness of the embodiment method.

[0152] This heart rate estimation method based on millimeter wave radar can overcome the multipath effect in the environment and the influence of static object clutter, and can overcome the influence of respiratory harmonics and noise on heartbeat waveform extraction and heart rate estimation, effectively improving the recognition accuracy.

[0153] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A heart rate estimation method based on millimeter wave radar, characterized in that: The following steps are included: S1. Use millimeter wave radar to sample the human body and obtain M frequency modulated continuous wave sweep signals S(1), S(2), …, S(k), …, S(M), where the kth sweep signal S(k) = [s1(k), s2(k), …, s N (k)], 1≤k≤M, N is the number of sampling points of a frequency sweep signal; S2. Apply the distance-fast Fourier transform to each sampling point of the M frequency modulated continuous wave sweep signals S(1), S(2), ..., S(k), ..., S(M), obtain the corresponding distance unit for each sampling point, form a distance-slow time matrix, and locate the distance from the human body to the radar; S3. After determining the distance between the human body and the radar, use the inverse tangent function to extract the human body phase, use the phase unwrapping signal function to unwrap the phase whose phase change exceeds the set threshold, make a difference between the unwrapped phases to obtain the phase difference between adjacent frames, use the interpolation method to smooth the difference that exceeds the threshold, and obtain the phase difference between adjacent frames after smoothing; S4, using wavelet transform to remove noise from the smoothed phase signal, obtaining denoised wavelet coefficients, performing inverse wavelet transform to reconstruct the signal, and obtaining a reconstructed phase signal; S5, combining the ensemble empirical mode decomposition and independent component analysis method, namely, the EEMD-ICA method, to separate the heartbeat signal, the respiratory signal, the respiratory harmonics and the noise in the reconstructed phase signal to obtain the heartbeat waveform; S51, use ensemble empirical mode decomposition EEMD to reconstruct the phase signal φ re (n) performing signal separation to obtain a plurality of intrinsic mode functions (IMF) components including heartbeat signals, respiratory signals, respiratory harmonics, and noise, and obtaining a component set; S511, adding white noise with standard normal distribution to the reconstructed phase signal φ for K times re (n), we get the new signal after adding white noise for the i-th time Among them, n i (n) represents the i-th added white noise sequence, i = 1, 2, ..., K; S512, adding white noise to the new signal obtained each time Perform empirical mode decomposition (EMD) to decompose the signal into a finite number of IMF components to separate the noise caused by unconscious micro-movements of the human body from the heartbeat signal. Assuming that each EMD decomposition obtains P IMF components, K decompositions will obtain a total of K groups of IMF component sets, and each group contains P IMF components. S513, according to the principle that the statistical average value of unrelated sequences is zero, select corresponding components from the K groups of IMF component sets and take the average respectively to obtain a group of P component sets IMFS(n); S52, using an independent component analysis method ICA to analyze the obtained component set to obtain a heartbeat waveform; S521, establishing an ICA model for the component set obtained in step S51: X(n)=AE(n) (11) Among them, X(n) is the signal observation value, that is, the component set IMFS(n), E(n) is the source signal to be estimated, which is composed of the real heartbeat waveform, respiratory waveform, respiratory harmonics and noise signal, and A is the mixing matrix; S522, reversely estimate the source signal E(n) through the observed signal, and establish an ICA solution model: Eest=WX(n) (12) Where Eest is the estimated value of the source signal E(n), which consists of multiple components, and W is the unmixing matrix; S523, calculating the estimated value Eest of the source signal E(n) by solving the unmixing matrix W, and then extracting the heartbeat waveform therefrom; S5231, zero-meaning the component set IMFS(n), that is, subtracting the average value of each IMF component in the component set from its own, to obtain a zero-meaning component set; S5232, whitening the component set subjected to zero mean processing to obtain a component set subjected to zero mean processing and whitening processing; S5233, using the maximum likelihood method to estimate each row of the unmixing matrix W for the component set IMFS(n) after zero mean and whitening, and multiplying them with the component set IMFS(n) respectively, the product Gaussian uncorrelation is maximized, and it is put into the estimated value Eest of the source signal E(n), and finally a set of independent components is obtained; S5234. Perform a fast Fourier transform on each independent component in the estimated value Eest of the source signal E(n) to obtain the spectra of the respective independent components. Then, calculate the total frequency band energy of each independent component and the heartbeat frequency band energy within the set frequency range in the frequency domain, and calculate the ratio of the heartbeat frequency band energy to the total frequency band energy respectively. Select the independent component with the largest ratio as the heartbeat waveform. S6. For the obtained heartbeat waveform, use the multiple signal classification algorithm to construct a spatial spectrum function and estimate the heartbeat frequency through spectrum peak search.

2. The heart rate estimation method based on millimeter wave radar as claimed in claim 1, characterized in that: In step S1, obtain the frequency-swept signals of M frequency-modulated continuous waves. Specifically, S11, millimeter wave radar transmission signal X t (t) and the received signal X r (t) are: X t (t)=sin(2πf min +πSt 2 ) (1) X r (t)=sin[2πf min (t-τ)+πS(t-τ) 2 ] (2) Among them, f min is the initial frequency of the millimeter-wave radar, S is the frequency growth slope, and τ is the signal round-trip time; S12. After the mixer performs I / Q mixing and low-pass filtering on the transmitted signal and the received signal, it outputs the intermediate frequency signal Y(t): Among them, f min is the initial frequency of the millimeter wave radar, S is the frequency growth slope, and τ is the signal round trip time; S13. Perform equally-spaced sampling on the intermediate frequency signal to obtain the frequency-swept signals S(1), S(2), …, S(k), …, S(M).

3. The heart rate estimation method based on millimeter wave radar as claimed in claim 1, characterized in that: In step S2, apply the range-fast Fourier transform to each sampling point of the M frequency-swept signals of the frequency-modulated continuous waves S(1), S(2), …, S(k), …, S(M). Each sampling point obtains the corresponding range cell, forming a range-slow time matrix, and locate the distance from the human body to the radar. Specifically, S21. Apply the range-fast Fourier transform to each sampling point of the M frequency-swept signals of the frequency-modulated continuous waves S(1), S(2), …, S(k), …, S(M). Each sampling point obtains the corresponding range cell, which is a complex signal. Place the range cells into the matrix by column to form a range-slow time matrix. S22, for each row in the distance-slow time matrix, i.e., each distance unit, the influence of static clutter and DC component is eliminated by subtracting the average value, and the distance unit du after removing the average value is obtained. pq `: Go pq `=do pq -mean(p) (4) Among them, pq is the complex signal in the p-th row and q-th column of the distance-slow time matrix, 1≤q≤M, mean(p) is the average value of the complex signal in the p-th row of the distance-slow time matrix; S23, the distance unit du after removing the average value pq `Construct a new distance-slow time matrix, obtain the sum of each row of the new distance-slow time matrix, and use the largest distance unit in the sum as the distance from the human body to the radar.

4. The heart rate estimation method based on millimeter wave radar as claimed in claim 3, characterized in that: In step S3, after determining the distance from the human body to the radar, use the arctangent function to extract the human body phase, use the phase unwrapping function to unwrap the phase where the phase change exceeds the set threshold, take the difference of the unwrapped phases to obtain the adjacent frame phase difference, and use the interpolation method to smooth the difference exceeding the threshold to obtain the smoothed adjacent frame phase difference. Specifically, S31. After determining the distance from the human body to the radar, the arctangent function can be used on the complex signal of the range cell where the human body is located to extract the phase φ(n) of the human body: where I(n) is the imaginary part of the nth frame complex signal in the range cell where the human body is located, and R(n) is the real part of the nth frame complex signal in the range cell where the human body is located; S32. Use the phase unwrapping function to unwrap the phase where the phase change exceeds ±π: where φ(n + 1) is the phase of the (n + 1)th frame; S33, find the phase difference φ between adjacent frames diff (n): f diff (n)=φ(n+1)-φ(n) (7) where 0 < n < M, φ(n) is the phase of the nth frame, and φ(n + 1) is the phase of the (n + 1)th frame; S34, if the phase difference between adjacent frames is φ diff (n) If the absolute value exceeds the set threshold, the interpolation method is used for smoothing to obtain the phase difference φ of adjacent frames after smoothing. diff (n)`: Among them, φ diff (n-1) is the n-1th phase difference, φ diff (n+1) is the n+1th phase difference value.

5. The heart rate estimation method based on millimeter wave radar according to any one of claims 1 to 4, characterized in that: In step S4, wavelet transform is used to remove the smoothed adjacent frame phase difference signal φ diff (n)`, obtain the denoised wavelet coefficients, perform inverse wavelet transform to reconstruct the signal, and obtain the reconstructed phase signal, specifically, S41, using wavelet transform to smooth the adjacent frame phase difference φ from multiple scales diff (n)` is decomposed to obtain the multi-scale wavelet transform coefficient w; S42, perform soft threshold denoising on the multi-scale wavelet transform coefficient w to obtain the denoised wavelet coefficient w λ : Among them, λ wt is a given threshold, sgn() is a sign function; S43, the denoised wavelet coefficients w λ Perform inverse wavelet transform to reconstruct the signal and obtain the reconstructed phase signal φ re (n).

Images