A method for estimating inter-beat intervals based on FMCW radar

Abstract

The present application relates to biological radar signal processing technical field, especially for a kind of based on FMCW radar heartbeat interval estimation method, including the following steps: S1: the chest vibration echo phase time series obtained by using second-order time derivative transformation into acceleration time series using FMCW radar irradiation monitoring object chest region;S2: in preset frequency analysis interval, dynamic mapping center frequency and wavelet order, construct adaptive wavelet dictionary, then using dictionary executes multi-order wavelet time domain convolution operation to acceleration time series and aggregates time-frequency energy distribution to generate time-frequency energy graph;S3: based on energy-guided multistage time-frequency feature screening technology and wavelet function conjugate are carried out inverse convolution operation reconstruction strategy, generate the approximate time domain signal of heart beat vibration, complete heart beat information inversion.The present application can effectively improve the accuracy and robustness of IBI extraction under low signal-to-noise ratio condition, to support the stable monitoring of light and moderate HRV (heart rate variability).

Description

Technical Field

[0001] This invention relates to the field of bio-radar signal processing technology, specifically to a method for estimating the cardiac interval based on FMCW radar. Background Technology

[0002] HRV (Heart Rate Variability) is a key indicator for assessing cardiovascular health, mental stress, and sleep quality, and its importance has driven a shift in health management paradigms towards prevention. Despite the urgent need, mainstream solutions such as smartwatches suffer from discomfort due to contact wearing and high rates of removal at night, resulting in data gaps; emerging contactless technologies (such as remote PPG) are limited by lighting conditions and privacy concerns.

[0003] Millimeter-wave radar, as a cutting-edge technology for non-contact monitoring, offers a new approach to HRV monitoring due to its high sensitivity in capturing minute chest cavity movements. However, existing millimeter-wave radar-based HRV extraction methods (including spectral analysis and deep learning) generally have limitations: errors increase dramatically under noise and long distances, or they are difficult to implement due to computational resource constraints, especially in balancing real-time performance, low complexity, and anti-interference capabilities. The modal aliasing problem between heartbeat and respiratory signals further increases the difficulty of high-precision HRV feature extraction. Therefore, developing an HRV monitoring system that can operate remotely, with low error, in real-time, and while protecting privacy in real-world scenarios has become a key bottleneck in promoting the development of proactive health management. Summary of the Invention

[0004] (a) Technical problems to be solved

[0005] To address the shortcomings of existing technologies, this invention provides a method for estimating cardiac intervals based on FMCW radar, which solves the problems mentioned in the background section.

[0006] (II) Technical Solution

[0007] To achieve the above objectives, this invention proposes a method for estimating the intercardiac interval (IBI) based on FMCW radar, which, combined with an energy-guided signal reconstruction mechanism, achieves non-contact and stable IBI estimation. The method mainly includes the following steps:

[0008] S1: First, the chest region of the monitored object is illuminated by FMCW (Frequency Modulated Continuous Wave) radar to acquire radar echo signals containing cardiopulmonary activity information. Coherent accumulation techniques are used to process the echo signals along the slow time dimension to improve the SNR (Signal-to-Noise Ratio) of the target signal. Subsequently, zero-padding and windowing operations are performed on the processed signal, and a range vector corresponding to each receiving antenna channel is obtained through a one-dimensional FFT (Fast Fourier Transform). Channel coherent accumulation of the range FFT vectors from multiple receiving antenna channels further enhances the target echo in the array normal direction. 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 minute chest cavity movements.

[0009] 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.

[0010] 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.

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

[0012]

[0013] 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.

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

[0015]

[0016] Where t represents time, f i The center frequency is j, where j is the imaginary unit. For order n i Related scaling factors.

[0017] Using the wavelet functions from the aforementioned dictionary, convolution operations are performed on the acceleration time series to calculate the convolution coefficients for each order of wavelet. For all orders at the same center frequency, geometric averaging is used for aggregation to improve the concentration and robustness of the energy distribution.

[0018]

[0019] 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.

[0020] The frequency f can be calculated by performing modulo-square processing on the aggregated coefficients. i Energy distribution at time t:

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

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

[0023] 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. First, ridge 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; if a frequency component satisfies the continuity condition throughout the entire time interval, it is considered a candidate ridge component, where the continuity condition can be jointly limited by a coverage threshold and a maximum discontinuity length threshold. The center frequency position of the ridge is determined based on the determination result.

[0024] Based on the weighted energy curve, the frequency range of the ridge is determined by searching in the upper and lower frequency directions above and below the center frequency position until the energy value is lower than a preset proportional threshold. This neighborhood frequency range is then combined with the entire time interval to form a rectangular frequency band, which serves as the basis for constructing the subsequent binary mask M(f,t).

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

[0026] After extracting the wavelet convolution coefficients within the binary mask M(f,t), deconvolve them with the conjugate of the corresponding wavelet function to obtain the approximate reconstructed signal for each frequency component:

[0027]

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

[0029] After energy-weighting and superimposing the reconstructed signals of all frequency components, a complete time-domain approximate heartbeat signal is obtained:

[0030]

[0031] in, For the final reconstructed time-domain heartbeat signal, f is the frequency index. ω represents the narrowband time-domain component obtained by deconvolution of the frequency channel f. f E(f) is the energy weight for the entire time period at frequency f. i (,t) represents time t and frequency f i Energy, ∑ t E(t,f) is the total energy of frequency f over the entire time series, ∑ f,t E(f i ,t) represents the total energy across all frequencies and the entire time range.

[0032] S4: Finally, the approximate time-domain signal that generates the heartbeat vibration. Peak detection is performed to extract the time interval between adjacent heartbeat peaks, forming a complete IBI sequence.

[0033] (III) Beneficial Effects

[0034] Compared with existing FMCW radar physiological signal processing methods, the time-frequency energy decomposition method for estimating the intercardia of FMCW radar based on adaptive wavelet dictionary proposed in this invention has the following beneficial effects in IBI estimation tasks:

[0035] This invention addresses the challenges of modal aliasing and instability caused by the proximity of cardiac and respiratory modal frequencies and large energy differences. It proposes a dynamic matching mechanism between center frequency and wavelet order, adaptively mapping the wavelet order according to the center frequency within a preset analysis band to construct a frequency-order matched wavelet dictionary. This enhances the characterization ability of non-steady-state and non-strictly periodic cardiac signals from the source and alleviates the difficulties of energy leakage and modal separation. In the energy extraction stage, multi-order wavelet convolution coefficients at the same center frequency are geometrically averaged and used to calculate the time-frequency energy map, thereby improving the concentration and robustness of energy distribution and providing 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. This strategy can focus on the main cardiac modality and reduce the risk of peak detection failure or drift when short-term heart rate drift and respiratory interference coexist. Attached Figure Description

[0036] Figure 1 This is a flowchart of the signal processing of the FMCW radar intercardia estimation method based on an adaptive wavelet dictionary according to the present invention.

[0037] Figure 2 This is a schematic diagram illustrating the mapping relationship between the center frequency and the adaptive wavelet order in an embodiment of the present invention;

[0038] Figure 3 This is a schematic diagram comparing the time-domain waveforms of Morlet wavelets corresponding to different center frequencies in some embodiments of the present invention;

[0039] Figure 4 This is a time-frequency energy distribution map obtained by adaptive time-frequency energy decomposition of the thoracic acceleration signal obtained from the second derivative of the phase time series in an embodiment of the present invention.

[0040] Figure 5 This is a schematic diagram of the time-domain waveform and peak detection results of the subject's heartbeat signal collected by the PPG device in an embodiment of the present invention;

[0041] Figure 6 This is a time series comparison chart of the IBI data extracted by the method of the present invention and the IBI data measured by the PPG device in an embodiment of the present invention. Detailed Implementation

[0042] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0043] Example

[0044] like Figure 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 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 represent the minimum and maximum orders (dimensionless), respectively, with [·] indicating rounding up. In this embodiment, the maximum order is 5 and the minimum is 3. The mapping relationship between the center frequency and the adaptive wavelet order is as follows: Figure 2 As shown.

[0054] Specifically, the wavelet function is of type Morlet, and its expression is:

[0055]

[0056] Where t represents time, f i The center frequency is j, where j is the imaginary unit. For order n i Related scaling factors. Examples of Morlet wavelet time-domain waveforms corresponding to different center frequencies, such as... Figure 3 As shown.

[0057] Using the wavelet functions from the aforementioned dictionary, convolution operations are performed on the acceleration time series to calculate the convolution coefficients for each order of wavelet. For multiple orders of results at the same center frequency, geometric averaging is used for aggregation to improve the concentration and robustness of the energy distribution.

[0058]

[0059] 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.

[0060] The frequency f can be calculated by performing modulo-square processing on the aggregated coefficients. i Energy distribution at time t:

[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 Let be the number of time sampling points, and let the set of frequency sampling points that satisfy the condition be denoted as 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 within the binary mask M(f,t), deconvolve them with the conjugate of the corresponding wavelet function to obtain the approximate reconstructed signal for each frequency component:

[0072]

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

[0074] After energy-weighting and superimposing the reconstructed signals of all frequency components, a complete time-domain approximate heartbeat signal is obtained:

[0075]

[0076] in, For the final reconstructed time-domain heartbeat signal, f is the frequency index. ω represents the narrowband time-domain component obtained by deconvolution of the frequency channel f. f Let E(t,f) be the energy weight for frequency f over the entire time period, and E(t,f) be the energy at time t and frequency f. ∑ t E(t,f) is the total energy of frequency f over the entire time series, ∑ f,t E(f i ,t) represents the total energy across all frequencies and the entire time range.

[0077] S4: Finally, the approximate time-domain signal that generates the heartbeat vibration. Peak detection is performed to extract the time interval between adjacent heartbeat peaks, forming a complete IBI sequence.

[0078] In this embodiment, the time-domain waveform and peak value detection results of the subject's heartbeat signal acquired by the PPG device are as follows: Figure 5 As shown, the PPG device model is HKG-07A infrared pulse sensor, which realizes peak detection and IBI calculation based on the minimum interval constraint principle.

[0079] The time series comparison chart of IBI data extracted from radar sensors and IBI data measured by PPG equipment in this embodiment is shown below. Figure 6 As shown.

[0080] In this embodiment, the monitoring error was calculated to be 15.28 ms based on the mean absolute error method.

[0081] Finally, it should be noted that the above descriptions are merely preferred embodiments of the present invention and are not intended 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 described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

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