Heart interval estimation method based on FMCW radar

Abstract

The invention relates to the technical field of biological radar signal processing, in particular to an FMCW radar-based heart beat interval estimation method, which comprises the following steps of: S1, converting a chest vibration echo phase time sequence obtained by irradiating a chest area of a monitored object by an FMCW radar into an acceleration time sequence by using a second-order time derivative; s2, dynamically mapping the center frequency and the wavelet order in a preset frequency analysis interval, constructing a self-adaptive wavelet dictionary, then executing multi-order wavelet time domain convolution operation on the acceleration time sequence by using the dictionary, and aggregating time-frequency energy distribution to generate a time-frequency energy diagram; and S3, performing a deconvolution operation reconstruction strategy based on an energy-guided multi-stage time-frequency feature screening technology and wavelet function conjugation, generating an approximate time domain signal of heart beat vibration, and completing heart beat information inversion. According to the method, the accuracy and robustness of IBI extraction can be effectively improved under the condition of low signal-to-noise ratio, so that stable monitoring of light and moderate HRV (heart rate variability) is supported.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of biological radar signal processing, in particular to a heart beat interval estimation method based on FMCW radar. BACKGROUND

[0002] HRV (heart rate variability) is a key indicator for assessing cardiovascular health, mental stress and sleep quality, and its importance drives the preventive shift in health management paradigm. Despite the urgent demand, mainstream solutions such as smart watches cause discomfort due to contact wear and high night removal rate, resulting in data loss; emerging non-contact technologies (such as remote PPG) are limited by lighting conditions and privacy concerns.

[0003] As a leading technology for non-contact monitoring, millimeter wave radar can capture subtle movements in the chest cavity with high sensitivity, providing a new path for HRV monitoring. However, existing HRV extraction methods based on millimeter wave radar (including spectral analysis, deep learning, etc.) generally have limitations: either the error increases dramatically under noise and long distance, or it is difficult to be practical due to the limitation of computing resources, especially the real-time, low complexity and anti-interference ability. The modal aliasing problem of heartbeat and breathing signals further increases the difficulty of high-precision HRV feature extraction. Therefore, developing an HRV monitoring system that can operate remotely, with low error and real-time in actual scenarios and protect privacy has become a key bottleneck to promote the development of active health management. SUMMARY

[0004] (I) Technical problems solved

[0005] In view of the deficiencies of the prior art, the present application provides a heart beat interval estimation method based on FMCW radar, which solves the problems proposed in the background art.

[0006] (II) Technical solutions

[0007] In order to achieve the above purpose, the present application provides a heart beat interval estimation method based on FMCW radar, which combines an energy-guided signal reconstruction mechanism to realize non-contact and stable IBI estimation. The method mainly includes the following steps:

[0008] S1: First, the thoracic region of the monitored object is irradiated by a FMCW (Frequency Modulated Continuous Wave) radar, and radar echo signals containing heart and lung activity information are obtained. The echo signals are processed along the slow time dimension using coherent accumulation technology to improve the SNR (Signal-to-Noise Ratio) of the target signal. Then, the processed signals are subjected to zero padding and windowing operations, and the distance dimension vector corresponding to each receiving antenna channel is obtained by one-dimensional FFT (Fast Fourier Transform). After channel coherent accumulation of the distance FFT vectors of multiple receiving antenna channels, the target echo in the array normal direction can be further enhanced. In the enhanced data, the complex phase of the target distance unit is extracted, and a continuous phase time sequence representing the micro-motion of the thoracic cavity is recovered by a phase unwrapping algorithm.

[0009] In order to enhance the time domain feature description of the heart beat micro-vibration signal, the phase sequence is further processed by second-order time difference to obtain an acceleration time sequence. The sequence reflects the change trend of the thoracic cavity acceleration driven by the heart beat, and is the input signal for subsequent time-frequency analysis.

[0010] S2: In the time-frequency analysis stage, the present application proposes a decomposition method based on an adaptive wavelet dictionary. Specifically, first, a sequence of center frequencies {f i} is generated by interpolation in the preset frequency analysis interval [f min , f max ]. For each center frequency f i in the set, a dynamically mapped wavelet order n i is obtained. Thus, a set of frequency-order matched wavelet functions is constructed.

[0011] Specifically, the dynamic mapping method is:

[0012]

[0013] where f i is the i-th center frequency (Hz), f min , f max are the lower / upper limits of the analysis band, n i is the wavelet order matched with f i (dimensionless), n min , n max are the minimum / maximum order (dimensionless), and [·] represents the upward rounding.

[0014] Specifically, the Morlet type wavelet function is selected, and its expression is:

[0015]

[0016] where t represents time, f i is the center frequency, and j is the imaginary unit. is the scale factor related to the order n i

[0017] The acceleration time series is convolved with the wavelet function in the dictionary to calculate the convolution coefficients of each order wavelet. The results of all orders at the same center frequency are aggregated using geometric mean to improve the concentration and robustness of energy distribution:

[0018]

[0019] where C(f i ,t) is the aggregated wavelet coefficient, x(t) is the acceleration time series, and * represents convolution operation, is the wavelet function with the center frequency f i and the order n i , and K i is the number of wavelet orders participating in aggregation at f i .

[0020] The energy distribution at time t and frequency f i can be calculated by performing modulus square processing on the aggregated coefficient:

[0021] E(f i ,t) = |C(f i ,t)| 2

[0022] where |·| represents the modulus of a complex number. After traversing the entire frequency range, the complete time-frequency energy map E(f i ,t) can be obtained.

[0023] S3: After obtaining the time-frequency energy map, the present application further performs approximate reconstruction of the target heart beat signal based on the energy-guided multi-level time-frequency feature screening and reconstruction strategy. First, ridge line monitoring is performed on the time-frequency energy map, specifically including: determining the continuity of the energy distribution of each frequency component in the time-frequency energy map over the entire time range, and if a frequency component meets the continuity condition over the entire time interval, it is taken as a candidate ridge line component, wherein the continuity condition can be jointly defined by a coverage threshold and a maximum discontinuity length threshold. The center frequency position of the ridge line is determined according to the determination result.

[0024] Based on the weighted energy curve, searching is performed in the upward and downward frequency directions of the center frequency position until the energy value is below a preset proportion threshold, to determine the neighborhood frequency range of the ridge line. The neighborhood frequency range and the entire time range form a rectangular frequency band region, which serves as the basis for constructing the subsequent binary mask M(f,t).

[0025] ​In the rectangular frequency band region, the energy value of each time-frequency point is compared with a preset energy threshold, time-frequency points with energy value greater than the threshold are retained and assigned a value of 1, and the rest are assigned a value of 0, to generate a binary mask M(f, t), wherein f represents frequency and t represents time.

[0026] After extracting the wavelet convolution coefficient in the binary mask M(f, t), the conjugate of the corresponding wavelet function is deconvoluted to obtain an approximate reconstruction signal of each frequency component:

[0027]

[0028] wherein ψ * represents the complex conjugate of the wavelet function, and Re(·) represents the real part.

[0029] After energy weighting and superimposition of the reconstruction signals of all frequency components, a complete time-domain approximate heart beat signal is obtained:

[0030]

[0031] wherein, is the final reconstructed heart beat time-domain signal, f is the frequency index, is a narrow-band time-domain component obtained by deconvolution of the frequency channel f, ω f is the energy weight of the full time period of frequency f, E(f i is the energy of time t, frequency f i , ∑ t E(t,f) is the energy sum of frequency f in the entire time sequence, ∑ f,t E(f i is the energy sum of all frequencies and the full time range.

[0032] S4: Finally, the approximate time-domain signal of the generated heart beat vibration is subjected to peak detection, and the time interval between adjacent heart beat peaks is extracted to form a complete IBI sequence.

[0033] (Three) beneficial effects

[0034] Compared with the existing FMCW radar physiological signal processing method, the IBI estimation method based on adaptive wavelet dictionary time-frequency energy decomposition FMCW radar proposed in the present application has the following beneficial effects in the IBI estimation task:

[0035] The present application aims at the problem of modal aliasing and instability caused by the proximity of heart beat and respiratory modal frequency and large energy difference, proposes a dynamic matching mechanism of center frequency-wavelet order, adaptively maps the wavelet order according to the center frequency within the preset analysis band, constructs a wavelet dictionary matched in frequency and order, and enhances the representation ability of non-stationary and non-strict periodic heart beat signal from the source and relieves the energy leakage and modal separation difficulty. In the energy extraction stage, the geometric mean aggregation is used for the multi-order wavelet convolution coefficient of the same center frequency, and the time-frequency energy graph is calculated, so as to improve the concentration and robustness of energy distribution and provide a stable energy spectrum basis for subsequent feature extraction. In the approximate reconstruction stage, an energy-guided multi-level time-frequency feature screening and reconstruction strategy is proposed, which can focus on the heart beat main modal and reduce the risk of peak detection failure or drift when the heart rate drifts for a short time and the respiratory interference coexists. BRIEF DESCRIPTION OF DRAWINGS

[0036] Figure 1 The signal processing flow block diagram of the FMCW radar heart beat interval estimation method based on the adaptive wavelet dictionary of the present application is shown in the figure.

[0037] Figure 2 The center frequency and adaptive wavelet order mapping relationship diagram of the embodiment of the present application is shown in the figure.

[0038] Figure 3 The Morlet wavelet time domain waveform comparison diagram corresponding to part of different center frequencies in the embodiment of the present application is shown in the figure.

[0039] Figure 4 The time-frequency energy distribution diagram obtained by adaptive time-frequency energy decomposition of the thoracic acceleration signal obtained from the second derivative of the phase time sequence in the embodiment of the present application is shown in the figure.

[0040] Figure 5 The time domain waveform and peak detection result diagram of the subject heart beat signal collected by the PPG device in the embodiment of the present application is shown in the figure.

[0041] Figure 6 The time sequence comparison diagram of the IBI data extracted by the method of the present application and the IBI data measured by the PPG device in the embodiment of the present application is shown in the figure. DETAILED DESCRIPTION

[0042] The technical solutions in the embodiments of the present application will be described clearly and completely in combination with the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the protection scope of the present application.

[0043] EMBODIMENT

[0044] AsFigures 1-6 As shown, an embodiment of the present invention proposes a method for estimating the cardiac interval based on FMCW radar, which includes the following steps:

[0045] S1: First, the chest area of ​​the monitored object is illuminated by the FMCW radar to obtain radar echo signals containing cardiopulmonary activity information.

[0046] In this embodiment, an IWR6843 millimeter-wave radar front-end and a DCA1000 high-speed data acquisition board are used to collect vibration signals from the chest region of the monitored object. The radar is configured as a three-transmitter, four-receiver antenna array, with an operating frequency range of 60 GHz to 64 GHz, a fast time dimension sampling rate of 5209 kS / s, 256 fast time dimension sampling points, and a slow time dimension sampling rate of 50 Hz.

[0047] During the data acquisition process, only a single person is monitored within the radiation range of the radar antenna. The person is lying flat on the bed with their chest cavity located within the effective radiation coverage area of ​​the radar antenna, and the distance between the chest cavity and the radar antenna is 20cm to 30cm.

[0048] Coherent accumulation technology is used to process the echo signal along the slow time dimension to improve the SNR of the target's surface heartbeat vibration signal. Subsequently, zero-padding and windowing operations are performed on the processed signal, and the range vector corresponding to each receiving antenna channel is obtained through one-dimensional FFT. By superimposing the FFT results of multiple channels point by point, the target echo in the array normal direction can be further enhanced. From the enhanced data, the complex phase of the target range cell is extracted, and a phase unwrapping algorithm is used to recover the continuous phase time series representing the minute movements of the chest cavity.

[0049] To enhance the characterization of the time-domain features of the cardiac micro-vibration signal, the phase sequence was further processed using second-order time-difference to obtain an acceleration time sequence. This sequence reflects the changing trend of thoracic cavity acceleration driven by cardiac contraction and serves as the input signal for subsequent time-frequency analysis.

[0050] S2: In the time-frequency analysis stage, this invention proposes a decomposition method based on an adaptive wavelet dictionary. Specifically, firstly, within the preset frequency analysis interval [f min ,f max The center frequency sequence {f} is generated through interpolation. i For each center frequency f in the set... i Dynamically map out the appropriate wavelet order n i This allows for the construction of a frequency-order matched set of wavelet functions. In this embodiment, the preset frequency analysis interval is [1,2].

[0051] Specifically, the dynamic mapping method is as follows:

[0052]

[0053] f i is the i-th center frequency (Hz), f min is the i-th center frequency (Hz), f max is the lower / upper limit of the analysis band, n i is the i-th center frequency (Hz), f i is the matching wavelet order (dimensionless), n min is the i-th center frequency (Hz), f max is the minimum / maximum order (dimensionless), respectively, and [·] represents the upward rounding. In this embodiment, the maximum value of the order is 5 and the minimum value is 3, and the mapping relationship between the center frequency and the adaptive wavelet order is shown in FIG. 2. Figure 2

[0054] Specifically, the Morlet type wavelet function is selected, and its expression is as follows:

[0055]

[0056] where t represents time, f i is the center frequency, j is the imaginary unit, and is the scale factor related to the order n i . The time domain waveforms of the Morlet wavelet corresponding to different center frequencies are shown in FIG. 3. Figure 3

[0057] The wavelet function in the above dictionary is used to perform convolution operation on the acceleration time series, and the convolution coefficients of each order wavelet are calculated. For the results of multiple orders under the same center frequency, the geometric mean is used for aggregation to improve the concentration and robustness of the energy distribution:

[0058]

[0059] where C(f i ,t) is the aggregated wavelet coefficient, x(t) is the acceleration time series, * represents convolution operation, is the wavelet function with the center frequency f i and the order n i , and K i is the number of wavelet orders participating in aggregation at f i .

[0060] The energy distribution of the frequency f i at the time t can be calculated by performing modulus square processing on the aggregated coefficient:

[0061] E(f i ,t) = |C(f i ,t)| 2 ​​

[0062] Where |·| represents the complex modulus. After traversing the entire frequency range, the complete time-frequency energy map E(f) can be obtained. i The resulting time-frequency energy distribution diagram is shown in Figure 1. Figure 4 As shown.

[0063] S3: After obtaining the time-frequency energy map, this invention further uses an energy-guided multi-level time-frequency feature screening and reconstruction strategy to approximate the reconstruction of the target heartbeat signal, such as... Figure 4 As shown. Specifically, on the time-frequency energy map E(f,t), let the discrete frequency indices be i = 1, 2, ... N. f The discrete-time index is m = 1, 2, ..., N t , where N f For the number of frequency sampling points, N t Let f be the number of time sampling points, corresponding to the frequency sampling points. i At time t m The energy is:

[0064] e i (t m )=E(f i ,t m )

[0065] For frequency sampling point f i Determine its corresponding energy sequence Whether the continuity condition is met throughout the entire time span. Preferably, the continuity determination can take the following form:

[0066]

[0067] Where 1(·) is the indicator function, η is the lower limit of energy monitoring, ρ is the continuity ratio threshold, and N t The set of frequency sampling points that satisfy the condition, where the number of time sampling points is denoted as the ridge set, is called the ridge set. And its center position is defined as the ridge center frequency f. ridge With f ridge Expand k upwards and downwards from the center. up With k down With each frequency sampling step size, the neighborhood frequency range is obtained.

[0068]

[0069] The aforementioned neighborhood frequency range and the entire time range constitute a rectangular frequency band region R.

[0070] Within a rectangular frequency band R, all energy values ​​are sorted in descending order of amplitude, and their cumulative energy percentage is calculated. The amplitude corresponding to a percentage of 99% is taken as the energy threshold θ. Based on this, a binary mask M(f,t) is constructed.

[0071] After extracting the wavelet convolution coefficients in the binary mask M(f, t), the approximate reconstruction signal of each frequency component is obtained by deconvolution with the conjugate of the corresponding wavelet function:

[0072]

[0073] where ψ * denotes the complex conjugate of the wavelet function, and Re(·) denotes the real part.

[0074] After energy weighting and superimposition of the reconstruction signals of all frequency components, the complete time-domain approximate heart beat signal is obtained:

[0075]

[0076] where, is the final reconstructed heart beat time-domain signal, f is the frequency index, is the narrow-band time-domain component obtained by deconvolution of the frequency channel f, ω f is the full-time energy weight of frequency f, E(t, f) is the energy of time t and frequency f, ∑ t E(t, f) is the total energy of frequency f over the entire time sequence, ∑ f,t E(f i , t) is the total energy of all frequencies and the full time range.

[0077] S4: Finally, the approximate time-domain signal of the generated heart beat vibration is subjected to peak detection, and the time interval between adjacent heart beat peaks is extracted to form a complete IBI sequence.

[0078] The time-domain waveform of the subject's heart beat signal collected by the PPG device in this embodiment and the peak detection result are shown in FIG. 3. Figure 5 The PPG device is a HKG-07A infrared pulse sensor, and peak detection and IBI calculation are realized based on the minimum interval constraint principle.

[0079] The time sequence comparison chart of the IBI data extracted based on the radar sensor and the IBI data measured by the PPG device in this embodiment is shown in FIG. 4. Figure 6

[0080] The monitoring error calculated based on the mean absolute error in this embodiment is 15.28 ms.

[0081] ​Finally, it should be noted that the above only describes the preferred embodiments of the present application and is not intended to limit the present application. Although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art will appreciate that the technical solutions described in the foregoing embodiments can be modified or some technical features thereof can be replaced by equivalent ones. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A method for estimating cardiac interval based on FMCW radar, characterized in that, Includes the following steps: S1: The phase time series of thoracic vibration echoes obtained by FMCW radar illuminating the thoracic region of the monitored object is transformed into an acceleration time series using the second time derivative. S2: Perform adaptive time-frequency energy decomposition within the preset frequency analysis interval [f min ,f max The center frequency sequence {f} is generated using interpolation. i }, for each center frequency f in the set i Dynamically map out the appropriate wavelet order n i To construct a wavelet dictionary, multi-order wavelet time-domain convolution operations are performed on the acceleration time series obtained from S1 using the wavelet dictionary. Aggregation operations are performed on the wavelet convolution results of each order at the same center frequency, and the time-frequency energy distribution is calculated to generate a time-frequency energy map. S3: Perform energy-guided multi-level time-frequency feature filtering and reconstruction on the time-frequency energy map. Determine the continuity of the frequency components across the entire time range to identify the center frequency position of the ridge. Based on the weighted energy distribution, determine the neighborhood frequency range of the ridge. Construct a rectangular frequency band within this neighborhood frequency range. Generate a binary mask M(f,t) by filtering feature points according to a preset threshold. Extract the aggregated wavelet convolution coefficients within M(f,t). Perform deconvolution operations on the wavelet convolution coefficients and the wavelet function of the corresponding center frequency to obtain the reconstructed signals for each frequency component. A weighted superposition operation is performed on the reconstructed signals of each frequency component to generate an approximate time-domain signal of the heartbeat vibration. S4: Approximate time-domain signal for generating heartbeat vibrations Peak detection is performed to extract the time interval between adjacent heartbeat peaks, thus obtaining heartbeat interval data.

2. The method for estimating cardiac interval based on FMCW radar according to claim 1, characterized in that, The acquisition of the phase time sequence of the thoracic cavity vibration echo in S1 includes the following steps: Ⅰ: Coherently accumulate the FMCW radar echo signal along the slow time dimension; II: Perform zero-padding and windowing operations on the accumulated results, and obtain the distance FFT vector of each receiving antenna channel through one-dimensional fast Fourier transform; III: Perform channel coherent accumulation of the distance FFT vectors of multiple receiving antenna channels to enhance the target echo signal in the array normal direction; IV: Extract the complex phase sequence of the target distance cell from the accumulated result, and perform phase unwinding processing on it to obtain the time-domain phase sequence of the thoracic vibration.

3. The method for estimating cardiac interval based on FMCW radar according to claim 1, characterized in that, The dynamic mapping method for constructing the wavelet dictionary by dynamically mapping the center frequency value to the wavelet function order in S2 is as follows: Among them, f i For the i-th center frequency (Hz), f min f max To analyze the lower / upper frequency band, n i To be with f i The wavelet order (dimensionless) of the matching, n min n max These are the minimum and maximum orders (dimensionless), respectively, and [·] indicates rounding up; The aggregation operation method performed on the wavelet convolution results of different orders at the same center frequency is as follows: Among them, C(f) i ,t) represents the aggregated wavelet coefficients, x(t) represents the acceleration time series, and * indicates convolution operation. The center frequency is f i The order is n i wavelet function, K i For in f i The number of wavelet orders participating in the aggregation; The method for calculating energy distribution is: E(f i ,t)=|C(f i ,t)| 2 ; Where |·| represents the complex modulus, and after traversing the entire frequency range, the complete time-frequency energy map E(f,t) can be obtained.

4. The method for estimating cardiac interval based on FMCW radar according to claim 1, characterized in that, The full-time range continuity determination in S3 includes: calculating the coverage and maximum hole length based on the comparison results of the energy values ​​of the frequency components at each time point with the time-adaptive threshold, and determining the center frequency position of the ridge line based on the coverage threshold and the maximum hole length threshold.

5. The method for estimating cardiac interval based on FMCW radar according to claim 1, characterized in that, The determination of the S3 neighborhood frequency range includes: searching in the upper and lower frequency directions according to the weighted energy value of the ridge center frequency position to the position where the energy value is lower than the preset ratio threshold, so as to determine the upper and lower boundaries of the neighborhood.

6. The method for estimating cardiac interval based on FMCW radar according to claim 3, characterized in that, The wavelet function is the Morlet wavelet, and its expression is: Where t represents time, f i The center frequency is j, where j is the imaginary unit. For order n i Related scaling factors.

Images