An Ultra-Wideband Radar Vital Sign Signal Denoising Algorithm Based on PE and TVF-EMD

Abstract

The present invention discloses an ultra-wideband radar vital sign signal denoising algorithm based on PE and TVF-EMD, belonging to the technical field of radar life detection. This method first preprocesses clutter and noise signals in the original radar echo signal by means of methods such as channel signal subtraction, subtraction averaging, linear trend suppression, automatic control gain, fast-time filtering, and slow-time filtering; then, by analyzing statistical features such as kurtosis, standard deviation, and variance of data in the slow-time direction, the PE algorithm is used to detect the position of the human target. This optimized algorithm can conveniently and accurately locate the moment when the signal mutates and can amplify the subtle changes in the signal. The TVF-EMD algorithm is used to decompose the extracted vital sign signal, calculate the respective energy ratios of each IMF component after decomposition within the respiration frequency band and the heartbeat frequency band, select appropriate signals for reconstruction, and finally perform spectral analysis on the reconstructed signal using FFT transformation. Using the radar vital sign denoising method proposed by the present invention, the target can be accurately located and vital information can be extracted, and the denoising result is very impressive.

Description

Technical Field

[0001] The present invention relates to the technical field of ultra-wideband radar life detection, and particularly relates to an ultra-wideband radar vital sign signal denoising algorithm based on PE and TVF-EMD. Background Technique

[0002] As an emerging vital sign detection means, ultra-wideband radar detection technology essentially emits electromagnetic waves, so it has extremely strong penetration ability, can penetrate non-metallic obstacles, and is reflected on the human body surface. At the same time, it also has advantages such as high range resolution, strong penetration ability, and strong anti-interference ability, and is not easily affected by external environmental factors such as weather, temperature, and light. It is a very ideal vital sign detection method; however, in the actual detection environment, there are not only the observed objects, but also other interferences, and it is necessary to perform denoising optimization processing on the radar raw echo to improve the signal-to-noise ratio. Therefore, the research on the ultra-wideband radar vital sign signal denoising algorithm is of great significance for improving the rescue efficiency of radar detecting vital signs and maintaining social stability. Summary of the Invention

[0003] Aiming at the deficiencies of the prior art, the present invention proposes an ultra-wideband radar vital sign signal denoising algorithm based on PE and TVF-EMD.

[0004] The technical solution adopted by the present invention is an ultra-wideband radar vital sign signal denoising algorithm based on PE and TVF-EMD, and the overall process is as Figure 1 shown, including the following 8 steps.

[0005] Step 1: Establish a mathematical model of the vital sign signal, and the specific steps are as follows.

[0006] Step 1.1: Assume that d0 is the distance from the surface of the human chest to the antenna, Δd is the periodic change of the chest cavity caused by human breathing and heartbeat, d r is the periodic change of the chest cavity caused by breathing, d h is the periodic change of the chest cavity caused by heartbeat, then the expression of the instantaneous distance from the radar antenna to the human chest is.

[0007] d(t) = d0 + Δd = d0 + d r + d h = d0 + A r cos(2πf r t) + A h cos(2πf h t) (1)

[0008] Wherein, t is the slow time, A r and A hThe amplitudes generated by respiratory motion and heartbeat motion respectively, f r and f h are the frequency of respiratory motion and the frequency of heartbeat motion respectively.

[0009] Step 1.2: Assume that in the detection scenario, except for the respiratory motion and heartbeat motion of the human body, other objects remain stationary, then the impulse response of the radar signal is.

[0010]

[0011] where τ is the fast time, a v δ(τ - τ v (t)) is the human target response, a v represents the amplitude of the human micro-motion echo signal, τ v (t) is the time delay of the human micro-motion echo in the fast time direction, is the sum of the responses of surrounding stationary targets, a i represents the amplitude of the echo signal of surrounding stationary objects, τ i is the time delay of the echo of surrounding stationary objects in the fast time direction, then τ v (t) can be expressed as.

[0012]

[0013] where v is the propagation speed of electromagnetic waves, τ0 is the fixed time delay between the radar antenna and the human body, τ r is the time delay of human respiratory motion, τ h is the time delay of heartbeat motion.

[0014] Step 1.3: Assume that the signal emitted by the transmitting antenna of the ultra-wideband pulsed radar is p(τ), then the echo signal received by the receiving antenna is.

[0015]

[0016] Step 1.4: Discretize the radar echo signal to obtain the mathematical model of the vital sign signal.

[0017]

[0018] where δ T is the fast time sampling interval, τ = mδ T m = 0, 1, 2…M - 1 is the discretization of the fast time, M is the number of fast time sampling points; T s is the slow time sampling interval, t = nT s n = 0, 1, 2…N - 1 is the discretization of the slow time, N is the number of slow time sampling points.

[0019] Step 2: Denoise and preprocess the original radar echo to improve the signal-to-noise ratio. The specific steps are as follows.

[0020] Step 2.1: Trace signal subtraction. Subtract each row of the radar echo signal from the previous row of the radar echo signal, and the quantities that remain constant during this process will be eliminated. For the radar echo matrix S(m,n), m is the number of sampling points of the radar echo delay time, and n is the number of sampling points of the scanning time. Assume the fast time series is Y j , j = 1, 2, 3…, N, where N is the number of sampling points. Then the result Y j ' of the trace signal subtraction is

[0021] Y j ' = Y j+1 - Y j (6)

[0022] Step 2.2: Subtraction of the average value. Take the average value of the radar echo data along the scanning time direction, that is, the DC component. The DC component is generated by stationary objects. Since the distance between the stationary object and the detection device remains fixed, while the breathing signal and heartbeat signal of the human body change periodically, the DC component can be removed by subtracting this average value from the original data. Then the result of eliminating the DC component is

[0023]

[0024] where can be approximately regarded as the DC component in the background clutter.

[0025] Step 2.3: Linear trend suppression. Estimate the potential linear trend in the background clutter obtained according to the linear least squares fitting in the scanning time direction, and then subtract this linear trend from the radar echo data to achieve the purpose of removing the background clutter. The calculation formula for the removal result is as follows.

[0026] W = Ω T - X(X T X) -1 X T Ω T (8)

[0027] where X = [x1, x2], x1 = [1, 2,…, N] T , x2 = [1, 1,…, 1] T .

[0028] Step 2.4: Automatic control gain. First, select a time window of 2λ + 1, and calculate the corresponding gain coefficient according to the energy therein to achieve the purpose of adaptive control. Let the n1th frame echo signal be r(τ, n1), then

[0029]

[0030] r E g(t,τ) = m g(t,τ) × r(t,τ) (10)

[0031] where g m (t,τ i ) is the gain coefficient, and r E (t,τ) is the processed human vital sign signal.

[0032] Step 2.5: Fast-time filtering. Wavelet threshold denoising based on wavelet transform is a common signal denoising method. How to select an appropriate threshold is the key to radar echo signal denoising. The generalized cross-validation function is used to determine the threshold. This method has greater advantages compared with the traditional D-J threshold as it does not rely on input-output data. The GCV formula is expressed as follows.

[0033]

[0034] where N represents the number of wavelet coefficients, N0 represents the number of wavelet coefficients equal to 0, w represents the echo signal with noise, and w t represents the echo signal after denoising.

[0035] Step 2.6: Slow-time filtering. A Butterworth low-pass filter is used to filter out high-frequency noise signals. Its squared magnitude function is defined as.

[0036]

[0037] where N is the order of the filter, ω c is the 3dB cut-off frequency, and ε is a parameter controlling the passband ripple amplitude. Since the human breathing frequency range is 0.1 - 0.5Hz and the heart rate frequency range is 0.8 - 3.0Hz, its frequency band is selected as 0.1 - 3.0Hz.

[0038] Step 3: The signal amplitudes in the target area range and non-target area range in the slow-time sequence of the radar echo matrix are different. The permutation entropy (PE) algorithm is used to detect the position of the human target. The specific steps are as follows.

[0039] Assume that the fast time is m1. For the slow-time sequence X[m1,i] = {y(1), y(2), …, y(i)} with length N, i = 1, 2, …, N, Y ∈ W M×N , phase space reconstruction can be obtained as follows.

[0040]

[0041] where τ d is the delay time, q is the embedding dimension, and K = N - (q - 1)τ d ; each row of this matrix is a reconstructed component, and there are K reconstructed components in total. By arranging the elements in the j-th reconstructed component in ascending order, a new reconstructed component can be obtained.

[0042] x(i + (j1 - 1)τ d ) ≤ x(i + (j2 - 1)τ d ) ≤ … ≤ x(i + (j q - 1)τ d ) (14)

[0043] If there are equal components among them, then they are arranged according to the size of the subscript i in j i , that is, the arrangement order of j1 and j2. For any component in the reconstructed matrix obtained by space reconstruction of an arbitrary time series, a set of sequences can be obtained.

[0044] S(l) = (j1, j2, …, j q ) (15)

[0045] where l = 1, 2, …, k, k ≤ q!, and the m-dimensional phase space mapping can generate m! kinds of S(l) orderings.

[0046] Calculate the occurrence probability P k of different permutation orders of each sequence S(l). The sum of the probabilities is 1. Then, in the form of Shannon entropy, the permutation entropy of the k different symbol sequences of the time series X[m1, i] can be defined as.

[0047]

[0048] When , that is, when the probabilities of each permutation order exist and are equal, the value of the permutation entropy is maximized at ln(m!), so the permutation entropy H P is normalized as follows.

[0049]

[0050] The value of H P represents the complexity of the slow time series X[m1, i]. The smaller the value of H P , the more regular the time series changes, and vice versa, the more random it is; in the radar echo matrix, the data change of the target position range in the slow time direction is relatively regular, and the corresponding PE value of the human target is smaller. The position P pos of the human target can be judged by finding the lowest point of the PE value result. The formula is as follows.

[0051] range = (v × P pos × Tf ) / 2 (18)

[0052] where \(v = 3\times10\) 8 m|s, \(T\) f is the sampling interval in the fast time direction.

[0053] Assume the lateral distance of the human chest cavity is \(d\), and the number of points occupied by the chest cavity distance in the signal is calculated as.

[0054]

[0055] Furthermore, the expression of the vital sign signal matrix \(\delta\) can be obtained as.

[0056]

[0057] Step 4: Range estimation, select the signal at the human body position and recombine the signals.

[0058] Step 5: Perform band-pass filtering on the combined signals.

[0059] Step 6: Perform TVF-EMD decomposition on the preprocessed signals. The specific operation steps are as follows.

[0060] Step 6.1: Local cut-off frequency

[0061] The cut-off frequency of the B-spline fitting filter changes with time. It mainly approximates the input signal by constructing polynomial splines and can be expressed in the following form.

[0062]

[0063] where \(\beta\) n (t) is the B-spline function, \(c(k)\) is the B-spline correlation coefficient, is the approximation error. When determining the B-spline correlation coefficient using the B-spline order \(n\) and nodes \(m\), the approximation error can be minimized, and the formula is expressed as.

[0064]

[0065] where [] ↑m is the upsampling operation, that is, inserting \(m\) sampling points between every two sampling points.

[0066] According to formula (21), \(c(k)\) can be expressed as.

[0067]

[0068] where [] ↓m is the downsampling operation, that is, sampling every \(m\) points; It is a pre-filter, and its formula is as follows.

[0069]

[0070] In summary, the B-spline fitter can be expressed in the following form.

[0071]

[0072] According to the above formula, the B-spline approximation can be regarded as a special form of low-pass filter. Generally, the node m is used to determine the local cut-off frequency of the B-spline filter, and the time for the signal to perform time-varying filtering is determined by the cut-off frequency. The specific implementation steps of time-varying filtering empirical mode decomposition are as follows.

[0073] Step 6.1.1: Calculate the instantaneous amplitude A(t) and instantaneous frequency of the input signal x(t) using the Hilbert transform

[0074]

[0075]

[0076] where X(t) is the Hilbert transform of the input signal x(t).

[0077] Step 6.1.2: Determine the local maximum value sequence A({t max}) and local minimum value sequence A({t min}) of the instantaneous amplitude A(t).

[0078] z(t) as the analytic signal of the multi-component signal can be expressed as.

[0079]

[0080] Thus, we can obtain.

[0081]

[0082]

[0083] where represents the instantaneous phase of the i-th component, and a i (t) represents the amplitude of the i-th component.

[0084] Assuming that a local minimum is obtained at t min , we can obtain.

[0085]

[0086] Substituting formula (31) into (29) and formula (30) gives:

[0087] A(t min ) = |a1(t min ) - a2(t min )| (32)

[0088]

[0089] According to the derivative operation, we have:

[0090] a′1(t min ) - a'2(t min ) = 0 (34)

[0091] Similarly, the following formula can be obtained:

[0092]

[0093] A(t max ) = |a1(t max ) + a2(t max )| (36)

[0094]

[0095] a′1(t max ) + a'2(t max ) = 0 (38)

[0096] Step 6.1.3: Perform interpolation fitting on A({t max}) and A({t min}) respectively to obtain the local extreme value fitting functions β1(t) and β2(t), and calculate the instantaneous mean function a1(t) and the instantaneous envelope function a2(t).

[0097]

[0098]

[0099] Step 6.1.4: Calculate the instantaneous frequency components and

[0100]

[0101] where η1(t) and η2(t) are the results after interpolation of respectively.

[0102] Step 6.1.5: Calculate the local cut-off frequency

[0103]

[0104] Step 6.1.6: Readjust the local cut-off frequency to solve the intermittent problem.

[0105] The local cut-off frequency may be affected by noise. Therefore, reconstruction can be performed to obtain a new signal.

[0106]

[0107] Step 6.2: Local mean function

[0108] Filter the input signal using a time-varying filter to obtain the local mean. Taking the extreme points of f(t) as the nodes for constructing the time-varying filter can make the cut-off frequency of the filter consistent with the local cut-off frequency. B-spline interpolation approximation of the input signal x(t) can obtain the approximation result m(t).

[0109] Step 6.3: Stop criterion

[0110] The judgment criterion is as follows.

[0111]

[0112]

[0113]

[0114] where B L (t) is the instantaneous bandwidth, is the weighted mean instantaneous frequency, ξ is the given bandwidth threshold. If the formula (45) is satisfied, the signal is a locally narrowband signal; according to the human breathing frequency of 0.1Hz - 0.5Hz and the heart rate of 0.8Hz - 3.0Hz, the echo signal can be divided into n IMF components and the percentage of breathing and heart rate energy for each component can be calculated in the frequency domain. The calculation formula is as follows.

[0115]

[0116]

[0117] where Q(j) is the frequency domain energy of the j-th IMF, Q r (j) is the breathing frequency domain energy of the j-th IMF, Q h (j) is the heart rate frequency domain energy of the j-th IMF, δ r and δ h are the judgment thresholds for breathing and heart rate respectively.

[0118] Step 7: Reconstruct the IMF respiration component and heartbeat component that meet the requirements, and the formula is as follows.

[0119] S r (t) = ∑IMF(t) (50)

[0120] S h (t) = ΣIMF(t) (51)

[0121] Step 8: Perform spectral analysis on the reconstructed signal using the Fast Fourier Transform (FFT).

[0122] The beneficial effects of adopting the above technical solutions are as follows: The present invention provides an ultra-wideband radar vital sign signal denoising algorithm based on PE and TVF-EMD. This method estimates the distance between the human target and the radar by calculating the PE value of the radar received pulse along the slow time direction. The combined radar echo signal is adaptively decomposed into IMFs using the TVF-EMD algorithm. According to the energy percentage of each IMF within the respiration and heartbeat frequency bands, the respiration and heartbeat signals are reconstructed and FFT is performed on them to obtain the frequencies of respiration and heartbeat. Description of the Drawings

[0123] Figure 1 It is a flowchart of an ultra-wideband radar vital sign signal denoising algorithm based on PE and TVF-EMD of the present invention.

[0124] Figure 2 It is a diagram of the vital sign signal model of the present invention.

[0125] Figure 3 It is a performance graph of the preprocessing algorithm of the present invention.

[0126] Among them, (a) Original radar echo matrix diagram; (b) RPS processing result diagram; (c) MS processing result diagram; (d) AGC processing result diagram; (e) Linear trend suppression result comparison diagram; (f) Fast time filtering result comparison diagram; (g) Slow time filtering result comparison diagram.

[0127] Figure 4 It is a schematic diagram showing the human body position with PE values at different positions of the present invention.

[0128] Among them, (a) PE value distance curve at 0.5 m; (b) PE value distance curve at 1.0 m; (c) PE value distance curve at 1.5 m.

[0129] Figure 5 It is a decomposition diagram of other comparison algorithms of the present invention.

[0130] Among them, (a) EMD algorithm decomposition diagram; (b) EEMD algorithm decomposition diagram; (c) CEEMD algorithm decomposition diagram.

[0131] Figure 6 This is the TVF-EMD decomposition diagram of the vital sign signals of the present invention.

[0132] Figure 7 This is the diagram of the proportion of respiratory and heartbeat energy of each IMF component of the present invention.

[0133] Figure 8 These are the reconstructed respiratory signal and heartbeat signal of the present invention.

[0134] Figure 9 These are the spectrograms of the respiratory signal and heartbeat signal of the present invention.

[0135] Among them, (a) is the spectrogram of the respiratory signal; (b) is the spectrogram of the heartbeat signal. Detailed implementation manners

[0136] The following further describes in detail the specific implementation manners of the present invention with reference to the accompanying drawings.

[0137] In the actual detection process of the radar, the radar echo signal mainly includes the vital sign signals of the target human body, noise interference in the surrounding environment, and some other clutter signals. The main purpose of life detection is to extract the vital sign information of the human body from the complex echo signal according to the characteristics of human respiratory, heartbeat and other life signals, and adopt a suitable denoising algorithm to filter and reduce the noise and clutter interference generated in this process to improve the signal-to-noise ratio. The overall algorithm flow is as Figure 1 shown; first, in the preprocessing stage, the background clutter is simply eliminated by the channel signal subtraction method; the DC component in the background clutter is eliminated by the direct subtraction average method; the linear trend in the radar echo data is eliminated by the linear trend suppression method; then the multipath effect is reduced by automatic control gain, and the weak respiratory signal and heartbeat signal are enhanced in the slow time direction to improve the signal-to-noise ratio of the human vital signs; the generalized cross-validation is used to determine the irrelevance between the threshold and the real data and noise energy; finally, the high-frequency noise in the echo signal is filtered by the Butterworth low-pass filter; then, by analyzing the statistical characteristics such as kurtosis, standard deviation, and variance of the data in the slow time direction, the PE algorithm is used to detect the position of the human target; and the TVF-EMD algorithm is used to decompose the extracted vital sign signals, calculate the respective energy proportions of each IMF component in the respiratory frequency band and the heartbeat frequency band after decomposition, select the appropriate signal for reconstruction, and finally perform spectral analysis on the reconstructed signal by using the FFT transform.

[0138] The vital sign signal model of this embodiment is as Figure 2 shown, and the specific steps are as follows.

[0139] Assume that d0 is the distance from the surface of the human chest to the antenna, and Δd is the periodic change of the chest cavity caused by human respiration and heartbeat (dr is the periodic change of the chest cavity caused by breathing, d h is the periodic change of the chest cavity caused by the heartbeat), then the expression for the instantaneous distance from the radar antenna to the human chest is as follows.

[0140] d(t) = d0 + Δd = d0 + d r + d h = d0 + A r cos(2πf r t) + A h cos(2πf h t) (1)

[0141] where t is the slow time, A r and A h are the amplitudes caused by the breathing motion and the heartbeat motion respectively, f r and f h are the frequencies of the breathing motion and the heartbeat motion respectively.

[0142] Assume that in the detection scenario, except for the breathing motion and the heartbeat motion of the human body, other objects remain stationary, then the impulse response of the radar signal is as follows.

[0143]

[0144] where τ is the fast time, a v δ(τ - τ v (t)) is the human target response, a v represents the amplitude of the human micro-motion echo signal, τ v (t) is the time delay of the human micro-motion echo in the fast time direction, is the sum of the responses of surrounding stationary targets, a i represents the amplitude of the echo signal of surrounding stationary objects, τ i is the time delay of the echo of surrounding stationary objects in the fast time direction, then τ v (t) can be expressed as follows.

[0145]

[0146] where v is the propagation speed of the electromagnetic wave, τ0 is the fixed time delay between the radar antenna and the human body, τ r is the time delay of the human breathing motion, τ h is the time delay of the heartbeat motion.

[0147] Assume that the signal emitted by the transmitting antenna of the ultra-wideband pulsed radar is p(τ), then the echo signal received by the receiving antenna is as follows.

[0148]

[0149] A mathematical model for obtaining vital sign signals by discretely processing radar echo signals.

[0150]

[0151] Among them, δ T is the fast-time sampling interval, τ = mδ T m = 0, 1, 2…M - 1 is the discretization of fast time, and M is the number of fast-time sampling points; T s is the slow-time sampling interval, t = nT s n = 0, 1, 2…N - 1 is the discretization of slow time, and N is the number of slow-time sampling points.

[0152] The performance diagram of the preprocessing algorithm of this embodiment is as Figure 3 shown, and the specific steps are as follows.

[0153] Channel signal subtraction. The simplest and most direct method to eliminate the background clutter in the radar echo signal is to subtract each row of the radar echo signal from the previous row of the radar echo signal. In this process, the quantities that remain constant will be eliminated. For the radar echo matrix S(m, n), m is the number of sampling points of the radar echo delay time, and n is the number of sampling points of the scanning time. Assuming the fast-time sequence is Y j , j = 1, 2, 3…, N, where N is the number of sampling points, then the result Y j ' of the channel signal subtraction is.

[0154] Y j ' = Y j+1 - Y j (6)

[0155] Subtraction-averaging method. Take the average value of the radar echo data along the scanning time direction, that is, the DC component; the DC component is generated by stationary objects. Since the distance between the stationary object and the detection device remains fixed, and the breathing signal and heartbeat signal of the human body change periodically, the DC component can be removed by subtracting this average value from the original data. Then the result of eliminating the DC component is.

[0156]

[0157] Among them, can be approximately regarded as the DC component in the background clutter.

[0158] Linear trend suppression. For the background clutter mixed with stationary signals and linear trends, the method of linear trend suppression can be adopted. The specific steps are as follows: The potential linear trend in the obtained background clutter can be estimated according to the linear least squares fitting in the scanning time direction, and then this linear trend is subtracted from the radar echo data to achieve the purpose of removing the background clutter. The calculation formula for the removal result is as follows.

[0159] W = Ω T -X(X T X) -1 X T Ω T (8)

[0160] Where X = [x1, x2], x1 = [1, 2,..., N] T , x2 = [1, 1,..., 1] T .

[0161] Automatic control gain. This method mainly plays the role of increasing the energy value of the vital sign signal. It will be automatically adjusted according to the feedback result of the energy value in the echo signal matrix, so as to strengthen the vital sign signal and ignore other useless signals to achieve the purpose of increasing the signal-to-noise ratio. The specific steps are as follows: First, select a time window of 2λ + 1, and calculate the corresponding gain coefficient according to the energy therein to complete the purpose of adaptive control; Let the echo signal of the n1th frame be r(τ, n1), then.

[0162]

[0163] r E (t, τ) = g m (t, τ) × r(t, τ) (10)

[0164] Where g m (t, τ i ) is the gain coefficient, and r E (t, τ) is the processed human vital sign signal.

[0165] Fast time filtering. Wavelet threshold denoising based on wavelet transform is a common signal denoising method. How to select an appropriate threshold is the key to radar echo signal denoising; The generalized cross-validation function is used to determine the threshold. This method has greater advantages compared with the traditional D-J threshold as it does not rely on input-output data. The GCV formula expression is as follows.

[0166]

[0167] Where N represents the number of wavelet coefficients, N0 represents the number of wavelet coefficients that are 0, w represents the echo signal with noise, w tRepresents the echo signal after noise reduction.

[0168] Slow-time filtering, using a Butterworth low-pass filter to filter out high-frequency noise signals, and its squared magnitude function is defined as.

[0169]

[0170] Where N is the order of the filter, ω c Is the 3dB cut-off frequency, and ε is a parameter that controls the band-pass ripple amplitude. Since the human breathing frequency range is 0.1 - 0.5Hz and the heartbeat frequency range is 0.8 - 3.0Hz, its frequency band is selected as 0.1 - 3.0Hz.

[0171] The schematic diagram of the human body position shown by the PE values at different positions in this embodiment is as Figure 4 Shown. The signal amplitudes in the target area range and the non-target area range in the slow-time sequence of the radar echo matrix are different. The permutation entropy algorithm is used to detect the position of the human target. Permutation entropy is an index that can be used to represent the complexity of the signal and measure the uncertainty of information. By comparing the values of adjacent time series to detect the dynamic changes of the time series, the more regular the time series, the smaller its corresponding permutation entropy; on the contrary, the more complex the time series, the larger its corresponding permutation entropy. Since the vital sign signals of the human body are relatively simple compared to the noise clutter in the detection environment, the PE value corresponding to the human target will be smaller. This characteristic can be used as the main basis for detecting the position of the human target. The specific steps are as follows.

[0172] Assume that the fast time is m1, for the slow-time sequence X[m1,i] = {y(1), y(2), …, y(i)}, i = 1, 2, …, N, Y ∈ W M×N , phase space reconstruction can be obtained.

[0173]

[0174] Where τ d Is the delay time, q is the embedding dimension, K = N - (q - 1)τ d , each row of this matrix is a reconstructed component, and there are K reconstructed components in total. By arranging the elements in the jth reconstructed component in ascending order, a new reconstructed component can be obtained.

[0175] x(i+(j1 - 1)τ d ) ≤ x(i+(j2 - 1)τ d ) ≤ … ≤ x(i+(j q - 1)τ d ) (14)

[0176] If there are equal components among them, then according to j iArrange in ascending order of the subscript i of the lower and upper indices, i.e., the arrangement order of j1 and j2. For any component in the reconstructed matrix obtained by spatial reconstruction of any time series, a set of sequences can be obtained.

[0177] S(l) = (j1, j2, …, j q ) (15)

[0178] where l = 1, 2, …, k, k ≤ q!, and the m-dimensional phase space mapping can generate m! kinds of S(l) orderings.

[0179] Calculate the occurrence probability P of different permutation orders of each sequence S(l) k , and the sum of probabilities is 1. Then, in the form of Shannon entropy, the permutation entropy of the k different symbol sequences of the time series X[m1, i] can be defined as.

[0180]

[0181] When that is, when the probabilities of each permutation order exist and are equal, the value of the permutation entropy is maximized at ln(m!), so the permutation entropy H P is normalized to obtain.

[0182]

[0183] H P represents the complexity of the slow time series X[m1, i]. The smaller the value of H P , the more regular the time series changes, and vice versa, the more random it is; in the radar echo matrix, the data change in the slow time direction within the target position range is relatively regular, and the corresponding PE value of the human target is smaller. The position P of the human target can be judged by finding the lowest point of the PE value result pos , and the formula is as follows.

[0184] range = (v × P pos × T f ) / 2 (18)

[0185] where v = 3 × 10 8 m / s, and T f is the sampling interval in the fast time direction.

[0186] Assume that the lateral distance of the human chest is d, and the number of points occupied by the chest distance in the signal is calculated as.

[0187]

[0188] Furthermore, the expression of the vital sign signal matrix δ can be obtained as.

[0189]

[0190] The decomposition diagrams of other comparison algorithms in this embodiment are as Figure 5 shown, and the specific steps are as follows.

[0191] Use the EMD, EEMD, and CEEMD algorithms to decompose the simulated mixed signal of the formula h(t) = sin(20πf1t) + 0.5sin(2πf2t), where f1 = 10, f2 = 100, and the sampling frequency is 500 Hz, and then compare and explain the advantages and disadvantages of each algorithm.

[0192] Empirical Mode Decomposition (EMD) is a commonly used method for local spectral analysis of non-stationary signals whose frequency changes with time. It mainly decomposes the signal into several Intrinsic Mode Functions (IMFs) using the extreme point information of the signal. From Figure 5 the EMD algorithm decomposition diagram in (a), it can be seen that this algorithm has the problem of mode mixing, mainly including two phenomena: a single IMF signal contains different time scales or the same time scale appears in different IMFs.

[0193] To avoid the problem of mode mixing, the concept of Ensemble Empirical Mode Decomposition (EEMD) is proposed, that is, adding white noise with a normal distribution to the original signal and regarding it as a whole, and then performing EMD decomposition on it. From Figure 5 the EEMD algorithm decomposition diagram in (b), it can be seen that its disadvantage is that although the added white noise is basically cancelled after ensemble averaging, there is still residue, and this residual noise cannot be ignored, and the number of ensemble averaging times is generally more than several hundred times, which is very time-consuming.

[0194] In response to the above problems, Complete Ensemble Empirical Mode Decomposition (CEEMD) is proposed. This algorithm adds a group of white noise and then subtracts a group of white noise. In this way, the noise will decrease as the number of ensemble averaging times increases, but at the same time, there is also the problem of too many iteration times. The decomposition diagram of the CEEMD algorithm is as Figure 5 shown in (c).

[0195] The TVF-EMD decomposition diagram of the vital sign signal in this embodiment is as Figure 6As shown, the signal is adaptively decomposed into a series of intrinsic mode functions, arranged in order from low order to high order and from high frequency to low frequency. Each intrinsic mode function represents the oscillation characteristics of the signal within different frequency ranges, and the order decreases as the IMF frequency increases. That is, the higher the IMF frequency, the lower the order; conversely, the lower the frequency, the higher the order. The specific steps are as follows.

[0196] Local cut-off frequency: The cut-off frequency of the B-spline fitting filter changes with time. It mainly approximates the input signal by constructing polynomial splines and can be expressed in the following form.

[0197]

[0198] Where β n (t) is the B-spline function, c(k) is the B-spline correlation coefficient, is the approximation error. When determining the B-spline correlation coefficient using the B-spline order n and nodes m, the approximation error can be minimized, and the formula is expressed as.

[0199]

[0200] Where [] ↑m is the upsampling operation, that is, inserting m sampling points between every two sampling points.

[0201] According to formula (21), c(k) can be expressed as.

[0202]

[0203] Where [] ↓m is the downsampling operation, that is, sampling every m points; is the pre-filter, and the formula is expressed as follows.

[0204]

[0205] In summary, the B-spline fitting device can be expressed in the following form.

[0206]

[0207] According to the above formula, the B-spline approximation can be regarded as a special form of low-pass filter. Usually, the local cut-off frequency of the B-spline filter is determined by the nodes m, and the time for the signal to perform time-varying filtering is determined by the cut-off frequency. The specific implementation steps for time-varying filtering empirical mode decomposition are as follows.

[0208] 1. Use the Hilbert transform to calculate the instantaneous amplitude A(t) and instantaneous frequency of the input signal x(t)

[0209]

[0210] where X(t) is the Hilbert transform of the input signal x(t).

[0211] 2. Determine the local maximum sequence A({t max}) and the local minimum sequence A({t min}) of the instantaneous amplitude A(t).

[0212] The analytic signal z(t) of the multi-component signal can be expressed as.

[0213]

[0214] Thus, it can be obtained that.

[0215]

[0216]

[0217] where represents the instantaneous phase of the i-th component, and a i (t) represents the amplitude of the i-th component.

[0218] Assuming that a local minimum is obtained at t min it can be obtained that.

[0219]

[0220] Substituting formula (31) into (29) and formula (30), it can be obtained that.

[0221] A(t min ) = |a1(t min ) - a2(t min )| (32)

[0222]

[0223] According to the derivative operation, it can be obtained that.

[0224] a′1(t min ) - a'2(t min ) = 0 (34)

[0225] Similarly, the following formula can be obtained.

[0226]

[0227] A(t max ) = |a1(t max) + a2(t max ) | (36)

[0228]

[0229] 3. Interpolate and fit A({t max}) and A({t min}) respectively to obtain the local extreme value fitting functions β1(t) and β2(t), and calculate the instantaneous mean function a1(t) and the instantaneous envelope function a2(t).

[0230]

[0231]

[0232] 4. Calculate the instantaneous frequency components and

[0233]

[0234]

[0235] where η1(t) and η2(t) are the results after interpolation of respectively.

[0236] 5. Calculate the local cut-off frequency

[0237]

[0238] 6. Readjust the local cut-off frequency to solve the intermittent problem.

[0239] The local cut-off frequency may be affected by noise, so reconstructing can obtain a new signal.

[0240]

[0241] Local mean function: Filter the input signal using a time-varying filter to obtain the local mean, and use the extreme points of f(t) as the nodes for constructing the time-varying filter, which can make the cut-off frequency of the filter consistent with the local cut-off frequency. Approximating the input signal x(t) by B-spline interpolation can obtain the approximation result m(t).

[0242] The stopping criterion is as follows, where B L (t) is the instantaneous bandwidth, is the weighted mean instantaneous frequency, and ξ is the given bandwidth threshold.

[0243]

[0244]

[0245]

[0246] The proportion of the breathing and heartbeat energy of each IMF component in this embodiment is as Figure 7 shown. The horizontal axis represents the serial number of the IMF, and the vertical axis represents the energy percentage of each IMF component. For the heartbeat signal sequence, the energy proportions of IMF5 and IMF6 both reach more than 90%. The energy proportions of the remaining IMFs in the heartbeat frequency band are relatively small. Therefore, the above two intrinsic mode functions can be selected for reconstructing the heartbeat signal. For the breathing signal sequence, the energy proportions of IMF8, IMF9, IMF10, IMF11, and IMF12 all reach 90%. The energy proportions of the remaining IMFs in the breathing frequency band are relatively small. Therefore, the above five intrinsic mode functions can be selected to reconstruct the breathing signal. The specific steps are as follows.

[0247] Divide the echo signal into n IMF components and calculate the breathing and heartbeat energy percentages for each component in the frequency domain. The calculation formula is as follows.

[0248]

[0249]

[0250] where Q(j) is the frequency domain energy of the jth IMF, Q r (j) is the breathing frequency domain energy of the jth IMF, Q h (j) is the heartbeat frequency domain energy of the jth IMF, δ r and δ h are the judgment thresholds for breathing and heartbeat respectively.

[0251] The reconstructed breathing signal and heartbeat signal in this embodiment are as Figure 8 shown. The specific steps are as follows.

[0252] Reconstruct the IMF breathing components and heartbeat components that meet the requirements. The formula is as follows.

[0253] S r (t) = ∑IMF(t) (50)

[0254] S h (t) = ∑IMF(t) (51)

[0255] After reconstructing the breathing and heartbeat signals, use the FFT method for frequency domain analysis respectively. The spectra of the breathing signal and heartbeat signal in this embodiment are as Figure 9 shown.

[0256] Although the specific embodiments of the present invention have been described above, those skilled in the art should understand that these are only examples, and various changes or modifications can be made to these embodiments without departing from the principles and essence of the present invention. The scope of the present invention is only defined by the appended claims.

Claims

1. A denoising algorithm for ultra-wideband radar vital sign signals based on PE and TVF-EMD, characterized in that, The method includes the following eight steps: Step 1: Establish a mathematical model of vital sign signals; The mathematical model of vital sign signals in Step 1 refers to a mathematical formula obtained by superimposing electromagnetic wave signals of ultra-wideband radar detection technology on human respiratory and heartbeat signals during rescue using ultra-wideband radar. Step 2: Perform denoising preprocessing on the original radar echo to improve the signal-to-noise ratio; Step 3: In the slow time series of the radar echo matrix, the signal amplitudes in the target area range and non-target area range are different. Use permutation entropy, i.e., the PE parameter, to detect the position of the human target; Step 4: Range estimation, select the signals at the human position and recombine the signals; Step 5: Perform band-pass filtering on the combined signals; Step 6: Perform TVF-EMD decomposition on the preprocessed signals; Step 7: Reconstruct the IMF respiratory component and heartbeat component that meet the requirements; Step 8: Perform spectral analysis on the reconstructed signals using FFT transformation.

2. A denoising algorithm for ultra-wideband radar vital sign signals based on PE and TVF-EMD according to claim 1, characterized in that, The process of Step 1 is as follows: Step 1.1: Assume that d0 is the distance from the surface of the human chest to the antenna, and Δd is the periodic change of the chest cavity caused by human respiration and heartbeat. Then the expression for the instantaneous distance from the radar antenna to the human chest is: d(t) = d0 + Δd = d0 + d r + d h = d0 + A r cos(2πft r )+ A h cos(2πft h t) (1) where t is the slow time, A r and A h are the amplitudes generated by respiratory motion and heartbeat motion respectively, and f r and f h are the frequencies of respiratory motion and heartbeat motion respectively; Step 1.2: Assume that in the detection scenario, except for the respiratory and heartbeat movements of the human body, other objects remain stationary. Then the impulse response of the radar signal is: where τ is the fast time, a v δ(τ - τ v (t)) is the human target response, a v represents the amplitude of the human micro - motion echo signal, τ v (t) is the time delay of the human micro - motion echo in the fast - time direction, is the sum of the responses of surrounding stationary targets, a i represents the amplitude of the echo signal of surrounding stationary objects, τ i is the time delay of the echo of surrounding stationary objects in the fast - time direction, then τ v (t) can be expressed as: where \(v\) is the propagation speed of the electromagnetic wave, \(\tau_0\) is the fixed time delay between the radar antenna and the human body, \(\tau\) r is the time delay of the human respiratory motion, and \(\tau\) h is the time delay of the heartbeat motion; Step 1.3: Assume that the signal emitted by the transmitting antenna of the ultra-wideband pulsed radar is p(τ). Then the echo signal received by the receiving antenna is: Step 1.4: Discretize the radar echo signal to obtain the mathematical model of vital sign signals: Among them, δ T is the fast time sampling interval, τ = mδ T where m = 0, 1, 2…M - 1 is the discretization of the fast time, and M is the number of fast time sampling points; T s is the slow time sampling interval, t = nT s where n = 0, 1, 2…N - 1 is the discretization of the slow time, and N is the number of slow time sampling points.

3. A denoising algorithm for ultra-wideband radar vital sign signals based on PE and TVF-EMD according to claim 1, characterized in that, The process of Step 2 is as follows: Step 2.1: Trace signal subtraction. The simplest and most direct method for eliminating background clutter in radar echo signals is to subtract each row of radar echo signals from the previous row. In this process, the quantities that remain constant will be eliminated. For the radar echo matrix S(m,n), m is the number of sampling points of the radar echo delay time, that is, the number of fast time points; n is the number of sampling points of the scanning time, also called the number of slow time points; assuming the fast time series is Y j , j = 1, 2, 3…, N, where N is the number of sampling points, then the result Y j ' is: Y’ j = Y j+1 - Y j (6) Step 2.2: Subtraction averaging method. Take the average value of the radar echo data along the scanning time direction, i.e., the DC component. The DC component is generated by stationary objects. Since the distance between stationary objects and the detection device remains fixed, while the respiratory and heartbeat signals of the human body are periodic, the DC component can be removed by subtracting this average value from the original data. The result of removing the DC component is: Among them, can be approximately regarded as the DC component in the background clutter; Step 2.3: Linear trend suppression. Estimate the potential linear trend in the background clutter obtained according to the linear least squares fitting in the scanning time direction, and then subtract this linear trend from the radar echo data to achieve the purpose of removing background clutter. The calculation formula for the removal result is as follows: W = Ω T -X(X T X) -1 X T Ω T (8) where X = [x1, x2], x1 = [1, 2, …, N] T , x2 = [1, 1, …, 1] T ; Step 2.4: Automatic control gain. Select a time window of 2λ + 1 and calculate the corresponding gain coefficient according to the energy therein to complete the purpose of adaptive control. Let the n1-th frame echo signal be r(τ, n1), then: r E (t, τ) = g m (t, τ) × r(t, τ) (10) where g m (t,τ i ) is the gain coefficient, and r E (t,τ) is the processed human vital sign signal; Step 2.5: Fast time filtering. Wavelet threshold denoising based on wavelet transform is a common signal denoising method. How to select the appropriate threshold is the key to radar echo signal denoising. Use the generalized cross-validation function to determine the threshold. This method has greater advantages compared with the traditional D-J threshold as it does not rely on input and output data. The GCV formula expression is as follows: Where N represents the number of wavelet coefficients, N0 represents the number of wavelet coefficients equal to 0, w represents the echo signal with noise, and w t represents the echo signal after noise reduction; Step 2.6: Slow-time filtering. A Butterworth low-pass filter is used to filter out high-frequency noise signals. Its squared magnitude function is defined as: where N is the order of the filter, ω c is the 3dB cut-off frequency, and ε is the parameter controlling the amplitude of the band-pass ripple. Since the human body's breathing frequency range is 0.1 - 0.5Hz and the heart rate frequency range is 0.8 - 3.0Hz, its frequency band is selected as 0.1 - 3.0Hz.

4. A denoising algorithm for ultra-wideband radar vital sign signals based on PE and TVF-EMD according to claim 1, characterized in that, The process of Step 3 is as follows: Assume that the fast time is m1. For the slow time series X[m1, i] = {y(1), y(2), …, y(i)} with length N, where i = 1, 2, …, N and Y ∈ W M×N , phase space reconstruction can be obtained as follows: where τ d is the delay time, q is the embedding dimension, and K = N - (q - 1)τ d ; each row of this matrix is a reconstructed component, and there are K reconstructed components in total. By arranging the elements in the j-th reconstructed component in ascending order, a new reconstructed component can be obtained: x(i+(j1-1)τ d )≤x(i+(j2-1)τ d )≤…≤x(i+(j q -1)τ d ) (14) If there are equal components among them, then arrange them according to the magnitude of the subscript i in j i in the arrangement order of j1 and j2. For any component in the reconstructed matrix obtained after the spatial reconstruction of any time series, a set of sequences can be obtained: S(l) = (j1, j2, …, j q ) (15) where l = 1, 2, …, k, k ≤ q!, and the m-dimensional phase space mapping can generate m! kinds of S(l) orderings; Calculate the occurrence probability P of different permutation orders of each sequence S(l). k , if the sum of probabilities is 1, then in the form of Shannon entropy, the permutation entropy of k different symbol sequences of the time series X[m1, i] can be defined as: When the value of the permutation entropy is maximized at ln(m!) when the probability of each permutation order exists and is equal. Therefore, for the permutation entropy H P the following normalization can be obtained: H P The value represents the complexity of the slow time series X[m1, i]. The smaller the value of H P , the more regular the time series changes. Conversely, the more random it is. In the radar echo matrix, the data changes in the slow time direction within the target position range are relatively regular, and the corresponding PE value of the human target is smaller. The position P of the human target can be judged by finding the lowest point of the PE value result pos , and the formula is as follows: range=(v×P pos ×T f ) / 2 (18) where v = 3×10 8 m / s, T f is the sampling interval in the fast time direction; Assume that the lateral distance of the human chest is d, and the number of points occupied by the chest distance in the signal is calculated as: Furthermore, the expression of the vital sign signal matrix δ can be obtained as:

5. A denoising algorithm for ultra-wideband radar vital sign signals based on PE and TVF-EMD according to claim 1, characterized in that, The process of Step 6 is as follows: Step 6.1: Local cut-off frequency The cut-off frequency of the B-spline fitting filter changes with time. It mainly approximates the input signal by constructing polynomial splines and can be expressed in the following form: where β n (t) is a B-spline function, and c(k) is the B-spline correlation coefficient. is the approximation error. When determining the B-spline correlation coefficient using the B-spline order n and the knots m, the approximation error can be minimized, and the formula is expressed as: Among them [] ↑m is an upsampling operation, that is, m sampling points are inserted between every two sampling points; According to formula (21), c(k) can be expressed as: where[] ↓m is a downsampling operation, that is, sampling every m points; is a pre-filter, and the formula is as follows: In summary, the B-spline fitting device can be expressed in the following form: According to the above formula, the B-spline approximation can be regarded as a special form of low-pass filter. Generally, the node m is used to determine the local cut-off frequency of the B-spline filter, and the time for the signal to perform time-varying filtering is determined by the cut-off frequency. The specific implementation steps of time-varying filtering empirical mode decomposition are as follows: Step 6.1.1: Calculate the instantaneous amplitude A(t) and the instantaneous frequency of the input signal x(t) by using the Hilbert transform where X(t) is the Hilbert transform of the input signal x(t); Step 6.1.2: Determine the local maximum value sequence A({t max}) and the local minimum value sequence A({t min}) of the instantaneous amplitude A(t). The analytical signal of the multi-component signal can be expressed as: From this, we can obtain: where represents the instantaneous phase of the i-th component, and a i (t) represents the amplitude of the i-th component; Assume that a local minimum is obtained at t min to obtain: Substituting formula (31) into (29) and formula (30), we can get: A(t min ) = |a1(t min ) - a2(t min )| (32) According to the derivative operation, we can obtain: a’1(t min ) - a'2(t min ) = 0 (34) Similarly, the following formula can be obtained: A(t max ) = |a1(t max ) + a2(t max )| (36) a’1(t max )+a'2(t max )=0 (38) Step 6.1.3: Interpolate and fit A({t max}) and A({t min}) respectively to obtain local extreme value fitting functions β1(t) and β2(t), and calculate the instantaneous mean function a1(t) and the instantaneous envelope function a2(t) Step 6.1.4: Calculate the instantaneous frequency component and where η1(t) and η2(t) are respectively the results after interpolating ; Step 6.1.5: Calculate the local cut-off frequency Step 6.1.6: Readjust the local cut-off frequency to solve the intermittent problem The local cut-off frequency may be affected by noise, so reconstructing it can obtain a new signal: Step 6.2: Local mean function Use the time-varying filter to filter the input signal to obtain the local mean, and use the extreme points of f(t) as the nodes for constructing the time-varying filter, which can make the cut-off frequency of the filter consistent with the local cut-off frequency. By performing B-spline interpolation approximation on the input signal x(t), the approximation result m(t) can be obtained; Step 6.3: Stop criterion The judgment criterion is as follows: Among them, B L (t) is the instantaneous bandwidth, is the weighted mean instantaneous frequency, ξ is the given bandwidth threshold, and if the formula (45) is satisfied, the signal is a locally narrowband signal; according to the human breathing frequency of 0.1Hz - 0.5Hz and the heartbeat frequency of 0.8Hz - 3.0Hz, the echo signal can be divided into n IMF components, and the percentage of breathing and heartbeat energy of each component can be calculated in the frequency domain. The calculation formula is as follows: where Q(j) is the frequency-domain energy of the j-th IMF, Q r (j) is the respiratory frequency-domain energy of the j-th IMF, Q h (j) is the heartbeat frequency-domain energy of the j-th IMF, δ r and δ h are the judgment thresholds for respiration and heartbeat, respectively.

6. A denoising algorithm for ultra-wideband radar vital sign signals based on PE and TVF-EMD according to claim 1, characterized in that, The process of Step 7 is as follows: Reconstruct the IMF respiration component and heartbeat component that meet the requirements. The formula is as follows: S r (t) = ∑IMF(t) (50) S h (t) = ∑IMF(t) (51).

Images