Night respiration and heartbeat monitoring signal processing method based on FMCW radar

Abstract

The invention belongs to the technical field of intelligent health monitoring, and particularly relates to an FMCW radar-based night respiration and heartbeat monitoring signal processing method, which comprises the following steps of S1, acquiring a chest vibration displacement signal of a human body by using an FMCW radar, and completing the acquisition of a radar signal; according to the recovery separation method based on the two-parameter LMS filter, monitoring and target positioning of multiple life signals can be achieved by the FMCW radar under the condition of low signal-to-noise ratio, and compared with a band-pass filtering (BPF) method and a wavelet transform (WT) method, more accurate breathing and heartbeat signals can be obtained by the method based on the two-parameter LMS filter, and the method has the advantages of being simple in structure, convenient to operate and high in practicability. Particularly, when the breathing signal has multiple harmonic waves and the frequency changes, the method is more accurate compared with other methods, breathing and heartbeat signals can be rapidly converged by setting a proper step length, and small steady-state deviation is kept when the signals are received next time.

Description

Technical Field

[0001] The present invention relates to the technical field of intelligent health monitoring, and particularly to a method for processing signals of night respiration and heartbeat monitoring based on FMCW radar. Background Art

[0002] Respiration and heart rate are important physiological parameters for clinical diagnosis and disease prevention; they directly reflect the physiological states of the human respiratory system and circulatory system. Nowadays, the vast majority of devices for measuring respiration and heartbeat are contact devices. Compared with traditional methods, non-contact micro-power radar is more convenient for patients to monitor respiration and heart rate, and does not increase the workload of nurses. These advantages make non-contact micro-power radar a favorable choice for hospital monitoring.

[0003] At the present stage, there are usually the following solutions for signal acquisition:

[0004] 1. Band-pass filtering is mostly used to extract the respiration and heartbeat signals measured by FMCW radar. However, the respiration signals generated by periodic chest wall movements often have many harmonics and variable instantaneous frequencies, which will interfere with the reliable detection of weaker heartbeat signals. Therefore, the band-pass filtering method cannot effectively separate these two signals and is not suitable for short-time signal processing.

[0005] 2. There is also a method that uses wavelet transform to separate respiration and heartbeat signals from mixed physiological signals. Compared with the band-pass filtering (BPF) method, the wavelet transform method can retain the boundaries and energies of the reconstructed signals in signals with short durations. However, when the respiration signal contains harmonics or the instantaneous frequency of the respiration signal changes severely, the reconstructed signals will be interfered.

[0006] 3. There is also a method that realizes adaptive cancellation of respiration harmonics and enhances heartbeat signals by minimizing the mean square error cost function. However, since the residual signal power related to the heartbeat signal in the error cost function is minimized, the heartbeat signals cannot be successfully retained.

[0007] 4. In order to suppress unnecessary periodic fluctuation components, some scholars have proposed a projection matrix method to extract heart rate. However, for correlation matrices that usually have large dimensions, eigenvalue decomposition is very time-consuming.

[0008] 5. There is also a method that uses continuous wavelet filters and an algorithm based on ensemble empirical mode decomposition (EEMD) to be applied to the separation of cardiopulmonary signals; this algorithm obtains accurate beat intervals in the time domain. However, when there are multiple harmonics in the respiration signal, since the EEMD method may not be able to completely solve the Gaussian white noise in the case of multiple harmonics, a high accuracy cannot be achieved. Therefore, we propose a method for processing signals of night respiration and heartbeat monitoring based on FMCW radar to solve the above problems. Summary of the Invention

[0009] (1) Technical problem to be solved

[0010] In view of the deficiencies of the prior art, the present invention provides a method for processing signals of night respiration and heartbeat monitoring based on an FMCW radar, and solves the problems raised in the above-mentioned background art.

[0011] (2) Technical solution

[0012] In order to achieve the above object, the present invention specifically adopts the following technical solutions:

[0013] A method for processing signals of night respiration and heartbeat monitoring based on an FMCW radar includes the following steps:

[0014] S1: Use an FMCW radar to obtain the thoracic vibration displacement signal of the human body and complete the acquisition of radar signals;

[0015] S2: Perform de-ramping processing on the signals collected by the radar, and design a low-pass filter to suppress the high-frequency components in the de-ramping processing to obtain the actual beat frequency signal;

[0016] S3: Orthogonalize the actual beat frequency signal into a complex beat frequency signal;

[0017] S4: Perform static clutter suppression on the complex beat frequency signal;

[0018] S5: Use a Hanning window to reduce the sidelobes of the signal after static clutter suppression, and use the fast Fourier transform to obtain the spectrum;

[0019] S6: Perform CFAR constant false alarm detection on the spectrum obtained after the fast Fourier transform;

[0020] S7: Estimate the angle of arrival of the signal after CFAR constant false alarm detection through the MUSIC algorithm;

[0021] S8: For the signal after estimating the angle of arrival through the MUSIC algorithm, use the PCA algorithm to extract the respiration signal and heartbeat signal of multiple targets;

[0022] S9: Use the two-parameter LMS filtering method to separate the respiration signal and the heartbeat signal.

[0023] Further, the specific content of S1 is as follows: Use the IWR6843 radar + DCA1000 high-speed data acquisition board, turn on 2 transmit antennas and 4 receive antennas, control the FMCW radar to transmit a linearly frequency-modulated signal with a linearly increasing frequency from its transmit antenna to the chest area of the target to be measured, with the frequency change range set at 60 GHz - 64 GHz. Set the sampling time to 50 us, and uniformly sample 256 points within this time, so that each sampling of the echo signal can collect subtle breathing signals. And set the radar sampling frequency to 100 Hz, so that the radar can extract the subtle changes in the chest vibration of the target person. Finally, extract the phase signal of the radar receive antenna.

[0024] Further, the specific steps of S2 are as follows:

[0025] S2.1: De-ramp the signal collected by the radar, and use linear regression to remove the trend component from the signal;

[0026] S2.2: Design a Butterworth low-pass filter to remove high-frequency components, select a cut-off frequency of 40 Hz and a sampling frequency of 100 Hz, and retain the low-frequency signal;

[0027] S2.3: Conduct spectrum analysis, and finally obtain the actual beat frequency signal.

[0028] Further, the specific steps of S3 are as follows:

[0029] S3.1: Extract the sine and cosine components from the actual beat frequency signal. The original signal contains multiple beat frequencies;

[0030] S3.2: Separate the sine wave and cosine wave from it by analyzing the frequency components of the signal;

[0031] S3.3: Generate a complex-form beat frequency signal. Combine the sine and cosine components, and a complex-form signal can be obtained. The real part of this complex signal represents the cosine component of the signal, and the imaginary part represents the sine component.

[0032] Further, in S4, static clutter suppression is achieved by subtracting the average value of all waveforms from each waveform to remove the DC component in the data; the formula is:

[0033]

[0034] where Y(t,τ) is the radar map data composed of the received radar signals, with a size of N s ×N c , N s is the number of samples of the fast sampling point d, and N c is the number of slow-time sampling points composed of the chirp frequencies.

[0035] Further, the specific steps of S5 are as follows:

[0036] S5.1: When obtaining the spectrum of the signal using the fast Fourier transform, in order to reduce the sidelobe effect in the spectrum and retain the vibration phase information of the signal, the signal is processed using a Hanning window.

[0037] S5.2: Calculate the weighting value of the Hanning window, which gradually changes with time from the start to the end of the signal, making the two ends of the signal smoother.

[0038] S5.3: Apply the Hanning window to the original signal.

[0039] S5.4: Multiply each point of the signal by the corresponding window function value to obtain the windowed signal.

[0040] S5.5: The amplitude distribution of the windowed signal changes under the influence of the window function, removing the abrupt edge effect.

[0041] S5.6: Perform FFT transformation on the windowed signal. Through FFT, the spectrum of the signal is obtained. This spectrum contains the frequency components of the signal. The spectrum of the signal shows the amplitude distribution at different frequencies, and the amplitudes of the heartbeat and respiration at different frequencies can be distinguished from the spectrogram.

[0042] Further, the specific steps of S6 are as follows:

[0043] S6.1: After the input signal, set a detection unit and a neighborhood unit; Detection unit: Set a time - space position as i, that is, x i (t) is the unit to be detected, and t is a moment of the radar signal; Neighborhood unit: Select a group of neighborhood units on both sides of the detection unit for estimating the power of the background noise. Select symmetric neighborhood units, set as M units on the left side of the detection unit and M units on the right side of the detection unit, forming a total neighborhood window of 2M. The neighborhood units should not contain the target signal, ensuring that they are only used for background noise estimation.

[0044] S6.2: Background noise estimation: In CA - CFAR, the estimation of the background noise is completed by averaging the signal powers of the neighborhood units; Set the units of the left neighborhood and the right neighborhood as: x1(t), x2(t),...x M (t), x i+1 (t), x i+2 (t),...x i+M (t), and calculate the average power of these 2M neighborhood units.

[0045]

[0046] Obtained is the background noise estimation of the i-th detection unit;

[0047] S6.3: Determine the estimation of the detection threshold background noise The detection threshold used to determine the target is calculated by setting a fixed threshold factor γ, and the threshold value is calculated as follows:

[0048] S6.4: Target detection: Compare the detection threshold with the power value of the detection unit; if the signal power x of the detection unit i (t) is greater than the threshold λ after background noise estimation threshold , then it is determined that this unit is a target unit, indicating that the target signal is detected. If x i (t) is less than the threshold, then it is determined that this unit is background noise and not a target. The judgment condition: If x i (t) ≥ λ threshold then the target is detected.

[0049] Furthermore, the specific content of S7 is: Let the signal vector received by the radar be where In the formula, Y(t) is the signal intensity received by the radar antenna; is the matrix composed of the steering vector a(θ), representing the gain of the signal at different arrival angles; X(t) is the intensity of the signal reaching the receiving array, and N(t) is the random error;

[0050] According to the principle of MUSIC algorithm, the signal arrival angle is determined based on the orthogonality characteristic of the vector and the noise subspace. The determination method is as follows:

[0051] a T (θ)U N = 0, where θ = {θ1, θ2,... θ M} In the formula, a T (θ) is the transpose of the steering vector a(θ); the matrix U N is the noise subspace separated from Y(t) using the MUSIC algorithm;

[0052] Since the noise subspace cannot fully satisfy orthogonality, a spectral peak optimization search is performed on the AOA estimation, and the expression is:

[0053] where P MUSIC (θ) is the probability spectral peak value with the given angle θ as the signal arrival angle.

[0054] Furthermore, the specific steps of S8 are as follows:

[0055] S8.1: Data preprocessing: In this step, first, it is necessary to determine the window size W for collecting data around the peak point. Then, select an appropriate number of reference points N, choosing several points before and after the peak point. Next, based on the number of selected reference points and the time length T, sample data of each reference point are collected on a slow time scale to form a data matrix of dimension N×T. To eliminate the static path component in the signal, the constant offset of each feature in the input signal is removed. This constant offset is obtained by calculating the average amplitude of each stream within the time window over a time scale of T seconds;

[0056] S8.2: Correlation estimation: After the data preprocessing is completed, the next step is to calculate the covariance matrix R, that is, to perform a correlation analysis on the various features of the input signal. The form of the covariance matrix is R = H×H T , where H is the data matrix after preprocessing, and H T is its transpose; the dimension of the covariance matrix is N×N, where N is the number of reference points. The covariance matrix can reflect the linear relationship or dependence between different features. Therefore, by calculating the covariance matrix, it can be known which features have a strong correlation. In the application of this system, we choose N = 7, that is, 3 points are selected before and after the reference point respectively. Such a span is 30 cm, which can effectively cover the width of the human target and ensure that the data can fully reflect the information of the thoracic region. By performing eigenvalue decomposition on the covariance matrix, the principal components with the largest variance and the most representativeness in the data can be extracted, thus achieving the goal of dimensionality reduction and feature selection;

[0057] S8.3: Eigenvalue decomposition: In the previous step, the covariance matrix R is obtained. Perform eigenvalue decomposition on this covariance matrix, that is, solve the eigenvalues and eigenvectors of matrix R. Specifically, first calculate the eigenvalues λ and the corresponding eigenvectors ν of the covariance matrix R. These eigenvectors and eigenvalues describe the main directions and variance magnitudes of the data; then, sort them in descending order according to the magnitudes of the eigenvalues, and arrange the eigenvectors corresponding to these sorted eigenvalues in order to form a new matrix P. Each column in this matrix is an eigenvector, and each eigenvector is arranged as a row vector in the matrix, representing the most representative and important directions in the data. Through eigenvalue decomposition, not only are various eigenvectors obtained, but also the importance of each principal component can be understood according to the magnitudes of the eigenvalues. The larger the eigenvalue, the greater the contribution of the corresponding eigenvector to the variance of the data, which also means that the direction corresponding to this eigenvector contains more effective information;

[0058] S8.4: Motion signal reconstruction: The data is transformed into the new space constructed by the new feature vectors, that is, Y = PH, and the corresponding y i = q i H constructs the principal components, qi and y i are the i-th eigenvector and the i-th principal component respectively; applying the PCA algorithm to process the received data, finally N principal components are obtained, {Y} = {y1, y2,... y n}, and then select the principal component y with the most information i for subsequent processing.

[0059] Further, the specific steps of S9 are as follows:

[0060] S9.1: Calculate the reference signal: For the physiological signal processed by the PCA algorithm, its peak value can be extracted by the peak detection method. Since the heartbeat signal is weaker than the respiration signal, the extracted reference signal can be regarded as the instantaneous respiration signal. The interval between consecutive measurements usually corresponds to an oscillation period. Without loss of generality, assume that the frequency of the respiration signal remains constant within a single period, but may change between different periods. At the sampling frequency of f s the instantaneous frequency and phase of the n-th target respiration signal can be calculated;

[0061] S9.2: Rewrite the respiration signal model: Ensure the effective harmonic number M of the respiration signal by checking whether there are peaks in the adjacent frequency intervals of multiples of the fundamental frequency on the spectrum of the n-th target respiration signal. Each harmonic and are time-varying amplitudes and phases, equivalent to the in-phase component and the quadrature component of the trigonometric function transformation and the amplitudes and phases and Represent the in-phase component and the quadrature component of the M-th harmonic as row vectors, and then arrange the amplitudes and phases of the in-phase component and the quadrature component into vectors;

[0062] The respiration signal model containing only the fundamental frequency of the n-th target can be expressed as:

[0063]

[0064] S9.3: Allocate variable step sizes: The adaptive filtering step size decreases as the harmonic number decreases. All step sizes are grouped as follows in the diagonal matrix:

[0065]

[0066] S9.4. Estimate the respiration and heartbeat signals;

[0067] The initial settings of the two parameters of the LMS filter are 0. Only the respiration signal of the fundamental frequency is calculated by the following formula:

[0068]

[0069] The remaining signals, including harmonics and heartbeat signals, are expressed as:

[0070]

[0071] Respiratory signals and their harmonics Calculated as:

[0072]

[0073] The estimated heartbeat signal is calculated as:

[0074]

[0075] The updated LMS filter coefficient equation is:

[0076]

[0077] If a suitable initial step size is selected, through two-parameter LMS filtering, the respiratory and heartbeat signals can converge quickly and maintain a small steady-state offset in the next step, so that the respiratory and heartbeat signals can be accurately estimated and the precise frequency components can be obtained;

[0078] The LMS filter is equivalent to a cascade of M notch filters, and these notch filters are centered on the harmonics of the fundamental frequency;

[0079] The 3dB bandwidth of each notch filter is Therefore, if the frequency difference between some harmonics of the respiratory signal and the heartbeat signal is less than then the amplitude attenuation of the heartbeat signal obtained by our method is greater than 3dB. If all the frequency components of the heartbeat signal remain unchanged, it will not affect the heart rate variability analysis. However, if the frequency difference between some harmonics of the respiratory signal and the heartbeat signal is much greater than then the amplitude of the heartbeat signal obtained by our method will not be severely suppressed either.

[0080] (III) Beneficial effects

[0081] Compared with the prior art, the present invention provides a method for processing signals for nighttime respiratory and heartbeat monitoring based on FMCW radar, and has the following beneficial effects:

[0082] The recovery and separation method of the present invention is based on a two-parameter least mean square filter, enabling the FMCW radar to monitor multi-life signals and target positioning under low signal-to-noise ratio conditions. Compared with the band-pass filtering (BPF) method and the wavelet transform (WT) method, the method based on the two-parameter LMS filter can obtain more accurate respiration and heartbeat signals. Especially when the respiration signal has multiple harmonics and its frequency changes, this method is more accurate than other methods. By setting an appropriate step size, the respiration and heartbeat signals can converge quickly and maintain a small steady-state offset when the next signal is received. Compared with the least mean square error cost function method, the two-parameter LMS method can well retain the respiration and heartbeat signals. Compared with the projection matrix method, the matrix dimension of LMS is small and it is easy to perform eigenvalue decomposition. Compared with the empirical mode decomposition method, the two-parameter LMS method can solve white noise in the case of multiple harmonics and can accurately separate the respiration and heartbeat signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] Figure 1 is the overall flowchart of the algorithm of the present invention;

[0084] Figure 2 is the flowchart of the two-parameter LMS algorithm in the algorithm of the present invention;

[0085] Figure 3 is the schematic diagram of the two-parameter LMS algorithm principle. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0086] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0087] Embodiment

[0088] As Figures 1-3 shown, a method for processing signals for monitoring respiration and heartbeat at night based on an FMCW radar proposed in an embodiment of the present invention includes the following steps: S1: Use the FMCW radar to obtain the thoracic vibration displacement signal of the human body and complete the acquisition of radar signals;

[0089] In this example, an IWR6843 radar + DCA1000 high-speed data acquisition board is used to collect thoracic vibration signals. Two transmitting antennas and four receiving antennas are enabled. The FMCW radar is controlled to transmit a linearly frequency-modulated signal with a linearly increasing frequency from its transmitting antenna to the thoracic region of the target to be measured. The frequency change range is set to 60 GHz - 64 GHz. The sampling time is set to 50 μs, and 256 points are evenly sampled within this time, so that subtle breathing signals can be collected for each sampling of the echo signal. Moreover, the radar sampling frequency is set to 100 Hz, so that the radar can extract the subtle changes in the thoracic vibration of the target person. Finally, the phase signal of the radar receiving antenna is extracted.

[0090] S2: Perform de-ramping processing on the signals collected by the radar, and design a low-pass filter to suppress the high-frequency components in the de-ramping processing to obtain the actual frequency signal. The specific steps are as follows:

[0091] S2.1: First, perform de-ramping on the signals, and use linear regression to remove the trend components from the signals. The radar signals contain trend components, which will affect subsequent signal processing. By removing the trend components from the original signals and retaining their true changing parts, subsequent spectral analysis and Fourier transform can be made easier.

[0092] S2.2: Design a Butterworth low-pass filter. Since it has a flat passband response and a smooth attenuation characteristic at the cut-off frequency, it can well suppress the high-frequency part in the radar signals. The order of the filter determines the steepness of the filter and the change of the frequency response. The higher the order, the steeper the response of the filter, but the computational complexity will increase. Since the cut-off frequency is 40 Hz and the sampling frequency is 100 Hz, set the passband frequency A p to 1 dB, and the stopband frequency A s to 40 dB (to make the stopband have a higher attenuation).

[0093] Calculate the required order n, and through it is calculated that n = 6.

[0094] Design the poles of the Butterworth filter: The poles of the Butterworth low-pass filter are evenly distributed on the unit circle. For an n = 6-order filter, its pole positions are: The specific values of the poles are:

[0095] The complex representation forms of these poles are:

[0096] Calculate the transfer function of the filter: For the Butterworth filter, the transfer function H(s) is:

[0097] For each pole p kSubstitute the above formula to obtain the transfer function of the filter.

[0098] Bilinear transformation: Convert the analog transfer function to a digital transfer function through bilinear transformation. The transformation formula is: where T is the sampling period, T = 1 / f s , z -1 is the z-domain delay operator. Substitute the obtained analog transfer function into the conversion formula to obtain the digital filter transfer function. Finally, the normalized cut-off frequency is w c = 2.513 rad / s;

[0099] S2.3: Perform spectral analysis on the signal after removing the high-frequency components, and finally obtain the actual frequency signal.

[0100] S3: Orthogonalize the actual frequency signal into a complex frequency signal; the specific steps are as follows:

[0101] S3.1: Extract the sine and cosine components from the actual frequency signal. The original signal contains multiple beat frequencies.

[0102] S3.2: Separate the sine wave and cosine wave by analyzing the frequency components of the analytic signal.

[0103] S3.3: Generate a complex form of the beat frequency signal: Combine the sine and cosine components to obtain a complex form of the signal. The real part of this complex signal represents the cosine component of the signal, and the imaginary part represents the sine component. In this way, the phase information of the signal is retained, and the signal becomes a complex form (complex beat frequency signal).

[0104] S4: Perform static clutter suppression on the complex beat frequency signal;

[0105] In a nighttime sleep environment, to accurately extract vital sign signals through radar, background clutter needs to be suppressed. Vital sign signals are weak and periodic signals. If the background clutter cannot be filtered out, it will interfere with subsequent data processing. Static clutter suppression removes the DC component in the data by subtracting the average value of all waveforms from each waveform. The formula is:

[0106]

[0107] where Y(t,τ) is the radar map data composed of received radar signals, with a size of N s ×N c . N s is the number of samples of the fast sampling point d, and N c is the number of slow-time sampling points composed of chirp frequencies.

[0108] S5: Apply a Hanning window to the signal after static clutter suppression to reduce the sidelobes, and use the Fast Fourier Transform (FFT) to obtain the spectrum. The specific process is as follows:

[0109] S5.1: When using the FFT to obtain the spectrum of the signal, to reduce the sidelobe effect in the spectrum and retain the vibration phase information of the signal, use the Hanning window to process the signal.

[0110] S5.2: Calculate the weighting values of the Hanning window, which change gradually over time from the start to the end of the signal, making the two ends of the signal smoother.

[0111] S5.3: Apply the Hanning window to the original signal.

[0112] S5.4: Multiply each point of the signal by the corresponding window function value to obtain the windowed signal.

[0113] S5.5: The amplitude distribution of the windowed signal changes under the influence of the window function, removing the sudden edge effect.

[0114] S5.6: Perform FFT transformation on the windowed signal. Through FFT, obtain the spectrum of the signal. This spectrum contains the frequency components of the signal. The spectrum of the signal shows the amplitude distribution at different frequencies, and the amplitudes of heartbeat and respiration at different frequencies can be distinguished from the spectrogram.

[0115] S6: Perform CFAR (Constant False Alarm Rate) detection on the spectrum obtained after the Fast Fourier Transform. The specific process is as follows:

[0116] S6.1: After the input signal, set a switching unit and a neighborhood unit: Switching unit (Cell Under Test, CUT): Set a time - space position as i, that is, x i (t) is the switching unit to be detected, and t is a moment of the radar signal. Neighborhood unit (Training Cells): Select a group of neighborhood units on both sides of the switching unit to estimate the power of the background noise. Select symmetric neighborhood units, set as M units on the left side of the switching unit and M units on the right side of the switching unit, forming a total neighborhood window of 2M. The neighborhood units should not contain the target signal, ensuring that they are only used for background noise estimation;

[0117] S6.2: Background noise estimation: In CA - CFAR, the estimation of the background noise is completed by averaging the signal powers of the neighborhood units. Set the units of the left neighborhood and the right neighborhood as: x1(t), x2(t),...x M (t) (left neighborhood), x i+1 (t), x i+2 (t),...xi+M (t) (right neighborhood), calculate the average power of these 2M neighborhood cells

[0118]

[0119] Obtain which is the background noise estimate of the i-th switching unit;

[0120] S6.3: Determine the estimate of the detection threshold background noise which will be used to determine the detection threshold of the target. By setting a fixed threshold factor γ, calculate this threshold value:

[0121] S6.4: Target detection: Compare the detection threshold with the power value of the switching unit. If the signal power x i (t) of the switching unit is greater than the threshold λ after background noise estimation threshold , then determine that this unit is a target unit, indicating that a target signal has been detected. If x i (t) is less than the threshold, then determine that this unit is background noise and not a target. Judgment condition: If x i (t) ≥ λ threshold then a target is detected.

[0122] S7: Estimate the angle of arrival of the signal after CFAR constant false alarm detection through the MUSIC algorithm; MUSIC is a spectrum-based algorithm and is often used for the estimation of the angle of arrival. By matching the received signal vector with a preset direction vector, the angle of arrival can be effectively identified. Specifically, for each received signal vector, compare the direction vector a(θ) corresponding to its angle, and a likelihood function can be obtained, indicating the probability of the actual arrival position of this angle. Let the signal vector received by the radar be where in the formula, Y(t) is the signal intensity received by the radar antenna; is the matrix composed of steering vectors a(θ), representing the gain of the signal for different arrival angles; X(t) is the intensity of the signal arriving at the receiving array, and N(t) is the random error.

[0123] According to the principle of the MUSIC algorithm, determine the signal arrival angle based on the orthogonality characteristic of the vector and the noise subspace. The determination method is as follows:

[0124] a T (θ)U N = 0, where θ = {θ1, θ2,... θ M}}, in the formula, a T (θ) is the transpose of the steering vector a(θ); the matrix U Nis the noise subspace separated from Y(t) using the MUSIC algorithm.

[0125] Since the noise subspace does not fully satisfy orthogonality, a spectral peak optimization search is performed for AOA estimation, and the expression is:

[0126] where P MUSIC (θ) is the probability spectral peak value with the given angle θ as the signal arrival angle.

[0127] S8: For the signal whose arrival angle is estimated by the MUSIC algorithm, the PCA algorithm is used to extract the respiratory signal and heartbeat signal of multiple targets; the aim is to find the point with the strongest energy containing thoracic cavity information from multiple target points. PCA is a classic algorithm widely used in data dimensionality reduction and feature extraction. Its core idea is to map high-dimensional feature data to a new low-dimensional space. In the new space, the dimension of the data is reduced while retaining as much information as possible. Specifically, PCA projects the data from the N-dimensional space to the K-dimensional space through the dimensionality reduction process. The features of the new K-dimensional data are called principal components, which are orthogonal to each other and can retain the variance or energy in the original data to the greatest extent. PCA is implemented through the following steps:

[0128] S8.1: Data preprocessing: In this step, first, it is necessary to determine the window size W for collecting data around the peak point, and then select an appropriate number of reference points N, usually several points before and after the peak point. Then, based on the selected number of reference points and the time length T, and collect the sample data of each reference point on the slow time scale, so as to form a data matrix with the dimension of N×T. To eliminate the static path component in the signal, the constant offset of each feature in the input signal is removed. This constant offset is obtained by calculating the average amplitude of each stream within the time window, usually on the time scale of T seconds. The purpose of removing the offset is to eliminate the interference from static environmental factors and make the extracted signal better reflect the dynamic changing features.

[0129] S8.2: Correlation estimation: After the data preprocessing is completed, the next step is to calculate the covariance matrix R, that is, to perform a correlation analysis on each feature of the input signal. The form of the covariance matrix is R = H×H T , where H is the data matrix after preprocessing, H Tis its transpose. The covariance matrix has dimensions N×N, where N is the number of reference points. The covariance matrix can reflect the linear relationships or dependencies between different features. Thus, by calculating the covariance matrix, we can understand which features have strong correlations. In the application of this system, we choose N = 7, that is, 3 points are selected before and after the reference point respectively. Such a span is 30 cm, which can effectively cover the width of the human target and ensure that the data can fully reflect the information of the thoracic region. By performing eigenvalue decomposition on the covariance matrix, the principal components with the largest variance and most representativeness in the data can be extracted, thus achieving the goals of dimensionality reduction and feature selection.

[0130] S8.3: Eigenvalue Decomposition: The covariance matrix R was obtained in the previous step. Perform eigenvalue decomposition on this covariance matrix, that is, solve the eigenvalues and eigenvectors of matrix R. Specifically, first calculate the eigenvalues λ and the corresponding eigenvectors ν of the covariance matrix R. These eigenvectors and eigenvalues describe the main directions and variance magnitudes of the data. Then, sort them in descending order according to the magnitudes of the eigenvalues, and arrange the eigenvectors corresponding to these sorted eigenvalues in order to form a new matrix P. Each column in this matrix is an eigenvector. Each eigenvector is arranged as a row vector in the matrix, representing the most representative and important directions in the data. Through eigenvalue decomposition, not only are various eigenvectors obtained, but also the importance of each principal component can be understood based on the magnitudes of the eigenvalues. The larger the eigenvalue, the greater the contribution of the corresponding eigenvector to the variance of the data, which also means that the direction corresponding to this eigenvector contains more valid information. Therefore, in this way, we can project the high-dimensional data into a new feature space, thereby extracting the most important information components. This process of eigenvalue decomposition is one of the core steps of the PCA algorithm. It maps the original data to a new orthogonal space, so that each new dimension contains the maximum variance information of the data, thus achieving effective dimensionality reduction and information compression.

[0131] S8.4: Motion Signal Reconstruction: The data is transformed into the new space constructed by the new eigenvectors, that is, Y = PH, and the corresponding y i = q i H constructs the principal components, q i and y i are the i-th eigenvector and the i-th principal component respectively.

[0132] Apply the PCA algorithm to process the received data, and finally obtain N principal components, {Y} = {y1, y2,... y n}, and then select the principal component y i with the most information for subsequent processing.

[0133] S9: Separate the respiratory signal and the heartbeat signal using the two-parameter LMS filtering method; the specific process is as follows:

[0134] S9.1: Calculate the reference signal: The peak of the physiological signal processed by the PCA algorithm is detected by the peak detection method. Since the heartbeat signal is much weaker than the respiratory signal, the reference signal can be considered as the instantaneous respiratory signal. The interval between consecutive instantaneous measurements corresponds to one oscillation period. Without loss of generality, assume that the frequency of the respiratory signal is constant within one period, but may vary between two periods. At the sampling rate of f s , calculate the instantaneous frequency and phase of the nth target respiratory signal.

[0135] S9.2: Reconstruct the respiratory signal model: Ensure the effective harmonic number M of the respiratory signal by checking whether there are peaks in the adjacent frequency intervals of the multiples of the fundamental frequency on the spectrum of the nth target respiratory signal. Each harmonic and are time-varying amplitudes and phases, equivalent to the in-phase and quadrature components of the trigonometric function transformation and of the amplitudes and phases and Represent the in-phase and quadrature components of the Mth harmonic as row vectors, and then arrange the amplitudes and phases of the in-phase and quadrature components into vectors.

[0136] The respiratory signal model containing only the fundamental frequency of the nth target can be expressed as:

[0137]

[0138] S9.3: Allocate variable step sizes: The adaptive filtering step size decreases as the harmonic number decreases, and all step sizes are grouped as follows in a diagonal matrix:

[0139]

[0140] S9.4: Estimate the respiratory and heartbeat signals

[0141] The initial settings of the two parameters of the LMS filter are 0, and only the respiratory signal of the fundamental frequency is calculated by the following formula:

[0142]

[0143] The remaining signals, including harmonics and heartbeat signals, are expressed as:

[0144]

[0145] The respiratory signal and its harmonics are calculated as:

[0146]

[0147] The estimated heart rate signal is calculated as:

[0148]

[0149] The updated LMS filter coefficient equation is:

[0150]

[0151] If an appropriate initial step size is selected, through two-parameter LMS filtering, the respiration and heart rate signals can converge quickly and maintain a small steady-state offset in the next step. In this way, the respiration and heart rate signals can be accurately estimated and the precise frequency components can be obtained.

[0152] The LMS filter is equivalent to a cascade of M notch filters, and these notch filters are centered on the harmonics of the fundamental frequency.

[0153] The 3dB bandwidth of each notch filter is Therefore, if the frequency difference between some harmonics of the respiration signal and the heart rate signal is less than then the amplitude attenuation of the heart rate signal obtained by our method is greater than 3dB. If all the frequency components of the heart rate signal remain unchanged, it will not affect the heart rate variability analysis. However, if the frequency difference between some harmonics of the respiration signal and the heart rate signal is much greater than then the amplitude of the heart rate signal obtained by our method will not be severely suppressed either.

[0154] Finally, it should be noted that the above are only the preferred embodiments of the present invention and are not used to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements for some of the technical features. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for processing signals of monitoring respiration and heartbeat at night based on FMCW radar, characterized in that: It includes the following steps: S1: Use an FMCW radar to obtain the thoracic vibration displacement signal of the human body and complete the acquisition of radar signals; S2: Perform de-ramping processing on the signals collected by the radar, and design a low-pass filter to suppress the high-frequency components in the de-ramping processing to obtain the actual beat frequency signal; S3: Orthogonalize the actual beat frequency signal into a complex beat frequency signal; S4: Perform static clutter suppression on the complex beat frequency signal; S5: Apply a Hanning window to the signal after static clutter suppression to reduce sidelobes, and use the fast Fourier transform (FFT) to obtain the spectrum; S6: Perform CFAR constant false alarm detection on the spectrum obtained after the fast Fourier transform; S7: Estimate the angle of arrival of the signal after CFAR constant false alarm detection through the MUSIC algorithm; S8: For the signal after estimating the angle of arrival through the MUSIC algorithm, use the PCA algorithm to extract the respiratory signal and heartbeat signal of multiple targets; S9: Use the two-parameter LMS filtering method to separate the respiratory signal and the heartbeat signal.

2. The method for processing signals of nocturnal respiration and heartbeat monitoring based on FMCW radar according to claim 1, wherein: The specific content of S1 is as follows: Use an IWR6843 radar + DCA1000 high-speed data acquisition board, turn on 2 transmitting antennas and 4 receiving antennas, control the FMCW radar to transmit a linearly frequency-modulated signal with a linearly increasing frequency from its transmitting antenna to the thoracic region of the target to be measured, with the frequency change range set to 60 GHz - 64 GHz, set the sampling time to 50 us, and uniformly sample 256 points within this time, so that each sampling of the echo signal can collect subtle respiratory signals, and set the radar sampling frequency to 100 Hz, so that the radar can extract the subtle changes in the thoracic vibration of the target person to be measured, and finally extract the phase signal of the radar receiving antenna.

3. A method for processing signals of monitoring respiration and heartbeat at night based on FMCW radar according to claim 1, characterized in that: The specific steps of S2 are as follows: S2.1: Perform de-ramping on the signals collected by the radar and use linear regression to remove the trend component from the signals; S2.2: Design a Butterworth low-pass filter to remove high-frequency components, select the cut-off frequency to be 40 Hz, the sampling frequency to be 100 Hz, and retain the low-frequency signals; S2.3: Perform spectrum analysis to finally obtain the actual beat frequency signal.

4. A method for processing signals of night respiration and heartbeat monitoring based on FMCW radar according to claim 1, characterized in that: The specific steps of S3 are as follows: S3.1: Extract the sine and cosine components from the actual beat frequency signal, and the original signal contains multiple beat frequencies; S3.2: Separate the sine wave and cosine wave by analyzing the frequency components of the signal; S3.3: Generate a complex beat frequency signal. Combine the sine and cosine components to obtain a complex signal. The real part of this complex signal represents the cosine component of the signal, and the imaginary part represents the sine component.

5. A method for processing signals for monitoring night-time respiration and heartbeat based on an FMCW radar according to claim 1, characterized in that: In S4, static clutter suppression is achieved by subtracting the average value of all waveforms from each waveform to remove the DC component in the data; the formula is: Among them, Y(t,τ) is the radar map data composed of received radar signals, with a size of N s ×N c , N s is the number of samples of the fast sampling point d, and N c is the number of slow-time sampling points composed of chirp frequencies.

6. A method for processing signals of nocturnal respiration and heartbeat monitoring based on FMCW radar according to claim 1, characterized in that: The specific steps of S5 are as follows: S5.1: When using the fast Fourier transform to obtain the spectrum of the signal, in order to reduce the sidelobe effect in the spectrum and retain the vibration phase information of the signal, use a Hanning window to process the signal; S5.2: Calculate the weighting value of the Hanning window, and this weighting value gradually changes with time from the start to the end of the signal, making the two ends of the signal smoother; S5.3: Apply the Hanning window to the original signal; S5.4: Multiply each point of the signal by the corresponding window function value to obtain the windowed signal; S5.5: The amplitude distribution of the windowed signal changes under the influence of the window function, removing the sudden edge effect; S5.6: Perform FFT transformation on the windowed signal. The spectrum of the signal is obtained through FFT. This spectrum contains the frequency components of the signal. The spectrum of the signal shows the amplitude distribution at different frequencies. The amplitudes of the heartbeat and respiration at different frequencies can be distinguished from the spectrogram.

7. A method for processing signals for monitoring night-time respiration and heartbeat based on an FMCW radar according to claim 1, characterized in that: The specific steps of S6 are as follows: S6.1: After the input signal, set the detection unit and the neighborhood unit; Detection unit: Set a spatio-temporal position as i, that is, x i (t) is the unit to be detected, and t is a moment of the radar signal; Neighborhood unit: Select a group of neighborhood units on both sides of the detection unit for estimating the power of the background noise. Select symmetric neighborhood units, set as M units on the left side of the detection unit and M units on the right side of the detection unit, forming a total neighborhood window of 2M. The neighborhood units should not contain the target signal, ensuring that they are only used for background noise estimation; S6.2: Background noise estimation: In CA-CFAR, the background noise is estimated by averaging the signal powers of the neighboring cells; the cells in the left neighborhood and the right neighborhood are set as: x1(t), x2(t),... x M (t), x i+1 (t), x i+2 (t),... x i+M (t), and calculate the average power of these 2M neighboring cells Obtain which is the background noise estimate of the i-th detection unit; S6.3: Determine the estimation of the detection threshold background noise The detection threshold for determining the target will be calculated by setting a fixed threshold factor γ to calculate this threshold value: S6.4: Target detection: Compare the detection threshold with the power value of the detection unit; if the signal power x i (t) of the detection unit is greater than the threshold λ after background noise estimation threshold , then determine that this unit is a target unit, indicating that a target signal has been detected. If x i (t) is less than the threshold, then determine that this unit is background noise and not a target. The judgment condition: If x i (t) ≥ λ threshold then a target is detected.

8. A method for processing signals of monitoring night respiration and heartbeat based on FMCW radar according to claim 1, characterized in that: The specific content of S7 is as follows: Let the signal vector received by the radar be where In the formula, Y(t) is the signal intensity received by the radar antenna; is the matrix composed of the steering vector a(θ), representing the gain of the signal for different arrival angles; X(t) is the intensity of the signal arriving at the receiving array, and N(t) is the random error; According to the MUSIC algorithm principle, determine the signal arrival angle based on the characteristic that the vector is orthogonal to the noise subspace. The determination method is as follows: a T (θ)U N = 0, where θ = {θ1, θ2,... θ M}, in the formula, a T (θ) is the transpose of the steering vector a(θ); the matrix U N is the noise subspace separated from Y(t) using the MUSIC algorithm; Since the noise subspace does not fully satisfy orthogonality, perform spectral peak optimization search on the AOA estimation. The expression is: Among them, P MUSIC (θ) is the peak value of the probability spectrum with the given angle θ as the angle of arrival of the signal.

9. A method for processing signals for monitoring respiration and heartbeat at night based on FMCW radar according to claim 1, characterized in that: The specific steps of S8 are as follows: S8.1: Data preprocessing: In this step, first determine the window size W for collecting data around the peak point, then select an appropriate number of reference points N, select several points before and after the peak point. Then, based on the selected number of reference points and the time length T, collect the sample data of each reference point on the slow time scale to form a data matrix with a dimension of N×T. To eliminate the static path component in the signal, remove the constant offset of each feature in the input signal. This constant offset is obtained by calculating the average amplitude of each stream within the time window and is performed on the time scale of T seconds; S8.2: Correlation Estimation: After the data preprocessing is completed, the next step is to calculate the covariance matrix R, that is, to perform a correlation analysis on each feature of the input signal. The form of the covariance matrix is R = H × H T , where H is the data matrix after preprocessing, and H T is its transpose; the dimension of the covariance matrix is N × N, where N is the number of reference points. The covariance matrix can reflect the linear relationship or dependence between different features. Therefore, by calculating the covariance matrix, we can find out which features have a strong correlation. In the application of this system, we choose N = 7, that is, 3 points are selected before and after the reference point respectively. Such a span is 30 cm, which can effectively cover the width of the human target and ensure that the data can fully reflect the information of the thoracic region. By performing eigenvalue decomposition on the covariance matrix, the principal components with the largest variance and the most representativeness in the data can be extracted, thus achieving the goals of dimensionality reduction and feature selection; S8.3: Eigenvalue decomposition: The covariance matrix R obtained in the previous step is subjected to eigenvalue decomposition, that is, solve the eigenvalues and eigenvectors of the matrix R. Specifically, first calculate the eigenvalues λ and the corresponding eigenvectors ν of the covariance matrix R. These eigenvectors and eigenvalues describe the main directions and variance magnitudes of the data. Then, sort them in descending order according to the magnitudes of the eigenvalues, and arrange the eigenvectors corresponding to these sorted eigenvalues in order to form a new matrix P. Each column in this matrix is an eigenvector, and each eigenvector is arranged as a row vector in the matrix, representing the most representative and important directions in the data. Through eigenvalue decomposition, not only can each eigenvector be obtained, but also the importance of each principal component can be understood according to the magnitudes of the eigenvalues. The larger the eigenvalue, the greater the contribution of the corresponding eigenvector to the variance of the data, which also means that the direction corresponding to this eigenvector contains more effective information; S8.4: Motion signal reconstruction: The data is transformed into a new space constructed by a new feature vector, i.e., Y = PH, corresponding to y i = q i H constructs the principal components, q i and y i are the i-th feature vector and the i-th principal component respectively; Apply the PCA algorithm to process the received data, and finally obtain N principal components, {Y} = {y1, y2,... y n}, and then select the principal component y i with the most information for subsequent processing.

10. A method for processing signals of monitoring respiration and heartbeat at night based on FMCW radar according to claim 1, characterized in that: The specific steps of S9 are as follows: S9.1: Calculate the reference signal: For the physiological signal processed by the PCA algorithm, its peak value can be extracted by the peak detection method. Since the heartbeat signal is weaker than the respiration signal, the extracted reference signal can be considered as the instantaneous respiration signal. The interval between consecutive measurements usually corresponds to one oscillation period. Without loss of generality, it is assumed that the frequency of the respiration signal remains constant within a single period, but may vary between different periods. In the case of a sampling frequency of f s , the instantaneous frequency and phase of the nth target respiration signal can be calculated; S9.2: Rewrite the respiratory signal model: Ensure the effective harmonic number M of the respiratory signal by checking whether there are peaks in the adjacent frequency intervals of the multiples of the fundamental frequency on the spectrum of the nth target respiratory signal. Each harmonic and has time-varying amplitude and phase, which is equivalent to the in-phase component and the quadrature component of the trigonometric function transformation and in terms of amplitude and phase and Represent the in-phase components and quadrature components of the M harmonics as row vectors, and then arrange the amplitudes and phases of the in-phase components and quadrature components into vectors; The fundamental frequency respiration signal model containing only the nth target can be expressed as: S9.3: Allocate variable step sizes: The adaptive filtering step size decreases as the number of harmonics decreases. All step sizes are grouped according to the following diagonal matrix: S9.

4. Estimate the respiration and heartbeat signals; The initial settings of the two parameters of the LMS filter are 0. Only the fundamental frequency respiration signal is calculated by the following formula: The remaining signals, including harmonics and heartbeat signals, are expressed as: Respiratory signal and its harmonics Calculated as: The estimated heartbeat signal is calculated as: The updated LMS filter coefficient equation is: If a suitable initial step size is selected, through two-parameter LMS filtering, the respiration and heartbeat signals can converge quickly and maintain a small steady-state offset in the next step, so that the respiration and heartbeat signals can be accurately estimated and the precise frequency components can be obtained; The LMS filter is equivalent to a cascade of M notch filters, and these notch filters are centered on the harmonics of the fundamental frequency; The 3dB bandwidth of each notch filter is Therefore, if the frequency difference between some harmonics of the respiration signal and the heartbeat signal is less than then the amplitude attenuation of the heartbeat signal obtained by our method is greater than 3dB. If all frequency components of the heartbeat signal remain unchanged, it will not affect the heart rate variability analysis. However, if the frequency difference between some harmonics of the respiration signal and the heartbeat signal is much greater than then the amplitude of the heartbeat signal obtained by our method will not be severely suppressed either.

Images