Heart rate estimation method and system based on millimeter wave radar stft image spectrum line extraction and harmonic inversion

Abstract

This invention discloses a heart rate estimation method and system based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion. This method acquires the millimeter-wave radar complex IQ slow-time signal of the target, and obtains the surface micro-motion phase signal after phase enhancement processing. STFT analysis is performed on the surface micro-motion phase signal to obtain a time-frequency energy map covering the respiratory, fundamental, and higher-order heartbeat harmonic frequency bands. The time-frequency energy map is then used as a two-dimensional image for background normalization, vertical interference suppression, multi-directional linear enhancement, binarization, skeletonization, and connected component analysis to extract candidate image spectral line objects. Spectral line objects are marked according to the respiratory dominant frequency trajectory and the predicted respiratory harmonic trajectory, and a heartbeat harmonic spectral line family satisfying the harmonic relationship and order normalization consistency is constructed. After vertical interference or spectral line breakage regions, the higher-order heartbeat harmonic trajectory is tracked based on the joint tracking of second- and third-order heartbeat harmonics. Finally, the heartbeat fundamental trajectory is recovered by inversion from the heartbeat harmonic spectral line family, and the actual visible fundamental point is checked back to output the heartbeat fundamental frequency, heart rate, fundamental state, and confidence level. This invention can stably obtain the fundamental frequency of the heartbeat in scenarios where the fundamental frequency of the heartbeat is weakly visible, intermittent, or invisible.

Description

Technical Field

[0001] This invention relates to the fields of non-contact vital sign detection, millimeter-wave radar signal processing, and time-frequency image analysis, and particularly to a method and system for heart rate estimation based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion. This method is suitable for monitoring vital signs in small animals such as mice, and can also be extended to the detection of localized surface micro-movements in the human body, sleep monitoring, and non-contact physiological parameter observation. Background Technology

[0002] Millimeter-wave radar can obtain information about minute displacements on the surface of a target by transmitting and receiving electromagnetic waves. It features non-contact operation, no need for attached electrodes, and the ability to monitor continuously for extended periods, and has been used to detect vital signs such as respiration and heartbeat. However, for small laboratory animals such as mice, the phase changes related to the heartbeat in the radar echo are usually weak due to their small size, high heart rate, and limited amplitude of chest and abdominal micro-movements, making the detection of the heartbeat fundamental frequency quite challenging.

[0003] Current millimeter-wave radar methods for detecting vital signs typically extract the micro-motion phase signal from the body surface from the complex IQ slow-time signal, and then estimate the respiratory and heart rates through methods such as bandpass filtering, FFT, short-time Fourier transform, peak search, or ridge tracking. In scenarios where the respiratory signal is strong and the heart rate signal is weak, the fundamental heart rate is often obscured by low-frequency background, higher-order harmonics of respiration, multipath reflection, or body motion interference, resulting in discontinuities or even obscurity of the fundamental heart rate spectrum.

[0004] In actual experiments, another phenomenon exists: the fundamental frequency of the heartbeat is not obvious in the STFT time-frequency energy diagram, while the second and third order harmonics of the heartbeat form relatively clear linear spectral lines. If heart rate is still determined solely by the maximum spectral peak or a single ridge within the fundamental frequency band of the heartbeat, it is easy to miss detections, make false detections, or cause trajectory jumps.

[0005] On the other hand, the heartbeat, respiration, respiration harmonics, vertical body motion interference, and noise spectral lines in the STFT time-frequency energy map all exhibit different image morphologies. If only frame-by-frame peak detection is used, it is difficult to utilize object-level information such as the duration, linewidth, smoothness, contrast, and slope of the spectral lines. Therefore, it is necessary to process the STFT time-frequency energy map as a two-dimensional image, extract continuous candidate image spectral line objects from the image and perform semantic discrimination, and use higher-order heartbeat harmonic spectral lines to invert and recover the heartbeat fundamental wave trajectory when the heartbeat fundamental wave is weakly visible or invisible. Summary of the Invention

[0006] In view of this, the purpose of the present invention is to provide a heart rate estimation method and system based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion, which solves the problem that existing methods are difficult to stably extract the heart rate fundamental frequency when the heart rate fundamental frequency is weakly visible, intermittent, or blocked by respiratory harmonics, low-frequency background and vertical interference.

[0007] To achieve the above objectives, the present invention provides the following technical solution: The heart rate estimation method based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion provided by this invention includes the following steps: Step 1: Obtain the millimeter-wave radar complex IQ slow-time signal corresponding to the distance gate or detection channel of the target under test, and perform outlier repair, stable segment truncation and phase enhancement processing on the complex IQ slow-time signal to obtain the surface micro-motion phase signal; Step 2: Perform short-time Fourier transform on the micro-motion phase signal on the body surface to obtain an STFT time-frequency energy map covering the respiratory frequency band, the fundamental frequency band of the heartbeat, and the higher-order heartbeat harmonic frequency band. Step 3: The STFT time-frequency energy map is preprocessed and line structure enhanced as a two-dimensional image, and candidate image spectral line objects are obtained through binarization, skeletonization and connected component analysis. Step 4: Extract the respiratory dominant frequency trajectory and generate the respiratory harmonic prediction trajectory. Based on the positional relationship between the candidate image spectral line object and the respiratory harmonic prediction trajectory, mark the candidate image spectral line object as a respiratory spectral line, respiratory harmonic spectral line, respiratory / heartbeat blurred spectral line, or suspected heartbeat spectral line. Step 5: Search for a family of heartbeat harmonic spectral lines in the candidate image spectral line object set that satisfies integer multiple frequency relationships, time overlap relationships, short interval succession relationships, and order normalized trajectory consistency. Step 6: When vertical interference, broadband transient interference, or spectral line breakage region is detected, joint supplementary tracking is performed in the STFT time-frequency energy map based on at least two higher-order heartbeat harmonic orders, and the higher-order heartbeat harmonic trajectory obtained by supplementary tracking is incorporated into the heartbeat harmonic spectral line family. Step 7: Perform robust smoothing on the higher-order heartbeat harmonic spectral lines in the heartbeat harmonic spectral line family, and invert and recover the heartbeat fundamental wave trajectory based on the order-normalized frequency trajectory of at least one higher-order heartbeat harmonic spectral line object in the heartbeat harmonic spectral line family. Step 8: Look back at the actual visible fundamental point in the STFT time-frequency energy map near the recovered heartbeat fundamental trajectory, and output the heartbeat fundamental frequency, heart rate, fundamental state and confidence level by combining the inversion residual, the proportion of actual visible fundamental waves and the spectral family score.

[0008] Furthermore, the phase enhancement processing includes at least one of robust IQ center compensation, sliding center compensation, mean-to-DC compensation, differential phase or DAM phase extraction methods; The robust IQ center compensation includes: Iteratively fit the center of the complex IQ trajectory, calculate the radius sequence from the sampling point to the fitted center, remove outliers based on the median of the radius and the median of the absolute deviation of the radius, and iteratively obtain the stable compensation center c*; calculate the compensated complex IQ slow-time signal according to the following formula: ; And calculate according to the following formula: ; Then, x[n] is detrended, mean-removed, or normalized. in, This represents a complex slow-time signal after IQ center compensation; Indicates the first The original millimeter-wave radar complex IQ slow-time signal at each sampling moment; This refers to the micro-motion phase signal of the body surface obtained by extracting and expanding the compensated complex signal, and the unit is usually radians; Indicates phase expansion; This represents the complex argument extraction operation, used to obtain the wrapper phase of a complex signal.

[0009] Furthermore, the STFT time-frequency energy map is calculated based on the body surface micro-motion phase signal x[n] according to the following formula: ; ; ; in, Represents the complex time-frequency matrix; w[n] is the short-time analysis window; The discrete sampling point index represents the phase signal of the micro-motion on the body surface; Indicates the time frame position or the center position of the analysis window in the STFT; This indicates the STFT analysis frequency or frequency point, in Hz. Indicates the sampling rate; Indicates time position ,frequency The time-frequency energy at that point is equal to the square of the STFT coefficient magnitude; Indicates that it is composed of all The complex time-frequency matrix formed; This represents the time-frequency energy after logarithmic transformation and maximum value normalization, in dB. Indicates that it is composed of all The constructed time-frequency energy matrix; A constant used to prevent overflow in logarithmic operations.

[0010] Furthermore, the preprocessing described in step three includes target frequency band interception, frequency-direction background normalization, time-direction background normalization, and vertical interference suppression; The vertical interference suppression includes statistically analyzing the proportion of bright pixels in each time series. When the proportion of bright pixels is greater than a preset threshold, the time series and its neighboring time series are marked as vertical interference regions. In subsequent image spectral line extraction, spectral family scoring, joint supplementary tracking, or fundamental wave trajectory recovery, the vertical interference regions are shielded, weighted, linearly crossed, or have endpoint constraints applied.

[0011] Furthermore, the linear structure enhancement described in step three includes employing at least one of the following: multi-directional linear filter, Gabor filter, morphological top-hat filter, Frangi-type linear structure filter, and skeletonization processing.

[0012] Furthermore, each candidate image spectral line object is composed of a set of time-frequency points in the STFT time-frequency energy map. The method represents and calculates at least one object-level feature, including duration, frequency range, average frequency, center frequency trajectory, linewidth, frequency smoothness, background contrast, time coverage, slope, discontinuity rate, and whether it crosses a vertical interference region; the object-level features are used for candidate spectral line semantic labeling, heartbeat harmonic spectral line family construction, and confidence scoring; wherein, Represents a time frame. Indicates frequency point, This represents the corresponding time-frequency energy.

[0013] Furthermore, in step four, the dominant respiratory frequency trajectory is identified within the respiratory frequency band. and based on Constructing respiratory harmonic prediction trajectories for candidate image spectral line objects The median frequency distance between the center frequency trajectory f_i(t) and the predicted trajectories of each order of respiratory harmonics is calculated using the following formula: ; Where k is a positive integer, Represents candidate image spectral line objects The time support range; when Less than the preset threshold At that time, the candidate image spectral line objects are marked as respiratory harmonic spectral lines; When the spectral line object is located in the heartbeat candidate frequency band or its object-level features are similar to the heartbeat spectral line, it is marked as a respiratory / heartbeat fuzzy spectral line and used as a fuzzy candidate in subsequent comprehensive judgment.

[0014] Furthermore, the joint pursuit described in step six is ​​implemented through dynamic programming, specifically including: A candidate grid for the heartbeat fundamental frequency is established. For each candidate fundamental frequency f0, the local spectral peak contrast at its second harmonic position 2f0 and third harmonic position 3f0 is examined simultaneously. The joint observation score is calculated according to the following formula:

[0015] Where C(P,t,f) represents the contrast of the local peak near frequency f in the time-frequency energy map relative to the local background; Introducing adjacent fundamental frequency hopping penalties and maximum hopping constraints in the time direction, according to the formula:

[0016] The trajectory of the fundamental wave of the continuous heartbeat is solved by dynamic programming, with the constraint that |f0 - f0'| ≤ Δf_max; Heart rate is calculated based on the fundamental wave trajectory of the heartbeat.

[0017] The heart rate estimation system based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion provided by this invention includes: The radar data acquisition module is configured to acquire the complex IQ slow-time signal of the target range gate or detection channel; The signal enhancement and phase extraction module is configured to perform phase enhancement processing on the complex IQ slow-time signal and extract the micro-motion phase signal of the body surface. The STFT time-frequency analysis module is configured to calculate the STFT time-frequency energy map of the micro-motion phase signal on the body surface; The STFT image preprocessing module is configured to perform background normalization, vertical interference suppression, and linear structure enhancement on the STFT time-frequency energy map. The image spectral line object extraction module is configured to extract candidate image spectral line objects and object-level features. The respiratory spectral marking module is configured to extract the dominant respiratory frequency trajectory and generate a respiratory harmonic prediction trajectory. The heartbeat harmonic spectral line family construction module is configured to construct the heartbeat harmonic spectral line family based on integer multiple frequency relationships, time continuity relationships, and order normalization trajectory consistency. The joint tracking module is configured to perform joint tracking of the fundamental wave trajectory of the heartbeat based on at least two higher-order heartbeat harmonics after vertical interference, broadband transient interference, or spectral line breakage regions. The spectral line family post-processing module is configured to perform robust smoothing on the higher-order heartbeat harmonic spectral lines in the heartbeat harmonic spectral line family. The heartbeat harmonic inversion and verification module is configured to invert and recover the fundamental wave trajectory of the heartbeat from the heartbeat harmonic spectral line family, and to check the actual visible fundamental wave points. The results output module is configured to output the heartbeat fundamental frequency, heart rate, fundamental frequency state, visible harmonic order, inversion residual, actual visible fundamental frequency ratio, and confidence level.

[0018] Furthermore, the joint tracking module implements tracking through dynamic programming, specifically configured as follows: A candidate grid for the heartbeat fundamental frequency is established. For each candidate fundamental frequency f0, the local spectral peak contrast at its second harmonic position 2f0 and third harmonic position 3f0 is examined simultaneously. The joint observation score is calculated according to the following formula:

[0019] Where C(P,t,f) represents the contrast of the local peak near frequency f in the time-frequency energy map relative to the local background; Introducing adjacent fundamental frequency hopping penalties and maximum hopping constraints in the time direction, according to the formula:

[0020] The trajectory of the fundamental wave of the continuous heartbeat is solved by dynamic programming, with the constraint that |f0 - f0'| ≤ Δf_max; Heart rate is calculated based on the fundamental wave trajectory of the heartbeat.

[0021] The beneficial effects of this invention are as follows: This invention discloses a heart rate estimation method and system based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion. This method acquires the millimeter-wave radar complex IQ slow-time signal of the target, and obtains the surface micro-motion phase signal after phase enhancement processing. STFT analysis is performed on the surface micro-motion phase signal to obtain a time-frequency energy map covering the respiratory, fundamental, and higher-order heartbeat harmonic frequency bands. The time-frequency energy map is then used as a two-dimensional image for background normalization, vertical interference suppression, multi-directional linear enhancement, binarization, skeletonization, and connected component analysis to extract candidate image spectral line objects. Spectral line objects are marked according to the respiratory dominant frequency trajectory and the predicted respiratory harmonic trajectory, and a heartbeat harmonic spectral line family satisfying the harmonic relationship and order normalization consistency is constructed. After vertical interference or spectral line breakage regions, the higher-order heartbeat harmonic trajectory is tracked based on the joint tracking of second- and third-order heartbeat harmonics. Finally, the heartbeat fundamental trajectory is recovered by inversion from the heartbeat harmonic spectral line family, and the actual visible fundamental point is checked back to output the heartbeat fundamental frequency, heart rate, fundamental state, and confidence level. This invention can stably obtain the fundamental heartbeat frequency in scenarios where the fundamental heartbeat is weakly visible, intermittent, or invisible. Compared with existing technologies, it has the following advantages: 1. This invention processes the STFT time-frequency energy map as a two-dimensional image, and replaces simple frame-by-frame peak search with image spectral line object extraction. It can utilize object-level features such as the temporal continuity, smoothness, linewidth and background contrast of the spectral lines to improve the reliability of heartbeat spectral line recognition in weak heartbeat scenarios.

[0022] 2. This invention improves the effects of complex IQ trajectory eccentricity, outliers, and slow drift on phase demodulation through phase enhancement processing, enabling the heartbeat fundamental and higher harmonics to form clearer linear spectral lines in the STFT image, providing stable input for subsequent object extraction.

[0023] 3. This invention reduces the interference of broadband vertical bright bars caused by body movement, attitude changes, or acquisition transients on spectral line extraction and fundamental frequency recovery by vertical interference suppression, short-interval continuation, and high-order harmonic compensation tracking mechanism after obstruction, thereby improving the continuity of the heartbeat spectral line object.

[0024] 4. This invention semantically labels candidate spectral lines by predicting the main respiratory trajectory and its harmonics, avoiding misjudging respiratory harmonics as heartbeats, while allowing ambiguous respiratory / heartbeat spectral lines to be comprehensively evaluated in subsequent spectral line family determination.

[0025] 5. When the fundamental frequency of the heartbeat is weakly visible, intermittent, or blocked, this invention does not rely on the peak value of the fundamental frequency. Instead, it recovers the trajectory of the fundamental frequency of the heartbeat by inverting the second-order, third-order, or higher-order heartbeat harmonic spectral lines, thereby improving the stability of heartbeat fundamental frequency detection in weak heartbeat scenarios.

[0026] 6. This invention can output the recovered fundamental wave trajectory, the actual visible fundamental wave point, the visible harmonic order, the inversion residual, the fundamental wave state, and the confidence level, making the detection results interpretable and facilitating subsequent medical or animal experimental analysis.

[0027] The above and other objects, advantages, and features of the present invention will be more fully set forth and demonstrated through the following detailed description of specific embodiments in conjunction with the accompanying drawings. Those skilled in the art, upon referring to the following detailed description and the accompanying drawings, will be able to better understand and realize the above advantages of the present invention. Other objects, features, and advantages of the present invention will become clearer after being described in detail in the detailed description section in conjunction with the accompanying drawings. Attached Figure Description

[0028] To make the objectives, technical solutions, and beneficial effects of this invention clearer, the following drawings are provided for illustration.

[0029] Figure 1 This is an overall flowchart of the heartbeat fundamental frequency detection method of the present invention; Figure 2 This is a schematic diagram of the heartbeat fundamental frequency detection system module of the present invention; Figure 3 This is a schematic diagram of the STFT image preprocessing and image spectral line object extraction process of the present invention; Figure 4 This is a schematic diagram illustrating the spectral line object features and semantic tags of candidate images in this invention; Figure 5This is a schematic diagram illustrating the relationship between the heartbeat harmonic spectral line family and the fundamental frequency inversion according to the present invention. Figure 6 This is a schematic diagram of the superposition of the high-order heartbeat harmonic spectral line family, joint tracking results, and recovered fundamental wave trajectory of the present invention; Figure 7 This is a schematic diagram comparing the fundamental wave state determination and result output of the present invention. Detailed Implementation

[0030] The present invention will be further described below with reference to the accompanying drawings and specific embodiments to enable those skilled in the art to better understand and implement the present invention. The described embodiments are only for explaining the present invention and are not intended to limit the scope of protection of the present invention.

[0031] Example 1

[0032] like Figure 1 As shown, this embodiment provides a heart rate estimation method and system based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion. Taking mouse millimeter-wave FMCW radar vital sign detection as an example, this embodiment acquires complex IQ slow-time signals from the chest, abdomen, or body surface region of the mouse under test. The output heart rate fundamental frequency is then recovered through phase enhancement, STFT image processing, image spectral line object extraction, respiratory spectral line labeling, heart rate harmonic spectral line family construction, joint tracking after occlusion, and heart rate fundamental frequency inversion.

[0033] Where I represents the in-phase component and Q represents the quadrature component; the complex IQ slow-time signal refers to the complex baseband echo signal acquired in accordance with the continuous chirping period or radar frame sequence on the selected range gate or detection channel, which can be expressed as z[n]=I[n]+jQ[n]; Step 1: Millimeter-wave radar echo signal preprocessing and phase enhancement Acquire the complex IQ slow-time signal z[n] = I[n] + jQ[n] corresponding to the target range gate or detection channel. Perform outlier detection on the acquired data. Outliers include sampling points whose real or imaginary parts are NaN, Inf, or significantly deviate from the local median. If outliers are found, use linear interpolation or robust interpolation between adjacent valid sampling points to repair them. Select stable data segments according to the experimental settings, removing data from the radar startup phase, target placement and adjustment phase, or phases with significant body movement.

[0034] In this embodiment, the micro-motion phase signal of the body surface is extracted, and the complex IQ slow-time signal is subjected to phase enhancement processing.

[0035] Preferably, robust IQ center compensation is used: In each iteration, the center of the circle is fitted based on the currently retained complex IQ sampling points. And calculate the distance from each sampling point to the center of the circle. radius:

[0036] Calculate the median of the radius sequence:

[0037] And calculate the median of the absolute deviation of the radius:

[0038] When the radius corresponding to the sampling point satisfies:

[0039] If the sampling point is identified as an outlier, it will be removed from the next center-fitting iteration. After repeated iterations, a stable compensation center c* is obtained, and the compensated complex signal is calculated.

[0040]

[0041]

[0042] Subsequently, x[n] is detrended, mean-removed, or normalized to obtain the micro-motion phase signal of the body surface.

[0043] in, This represents a complex slow-time signal after IQ center compensation; Indicates the first The original millimeter-wave radar complex IQ slow-time signal at each sampling moment; This refers to the micro-motion phase signal of the body surface obtained by extracting and expanding the compensated complex signal, and the unit is usually radians; Indicates phase expansion; This represents the complex argument extraction operation, used to obtain the wrapper phase of a complex signal.

[0044] This processing reduces the impact of IQ trajectory eccentricity, low-frequency drift, and outliers on phase expansion, resulting in more continuous and clearer line spectral lines in the subsequent STFT time-frequency energy plot.

[0045] Step 2: Generation of time-frequency energy map based on STFT A short-time Fourier transform is performed on the micro-motion phase signal on the body surface to obtain the complex time-frequency matrix S(f,t), and the time-frequency energy map is calculated.

[0046] ; ; ; in, Represents the complex time-frequency matrix; w[n] is the short-time analysis window; The discrete sampling point index represents the phase signal of the micro-motion on the body surface; Indicates the time frame position or the center position of the analysis window in the STFT; This indicates the STFT analysis frequency or frequency point, in Hz. Indicates the sampling rate; Indicates time position ,frequency The time-frequency energy at that point is equal to the square of the STFT coefficient magnitude; Indicates that it is composed of all The complex time-frequency matrix formed; This represents the time-frequency energy after logarithmic transformation and maximum value normalization, in dB. Indicates that it is composed of all The constructed time-frequency energy matrix; A constant to prevent overflow in logarithmic operations; A constant used to prevent overflow in logarithmic operations.

[0047] Preferably, in detection scenarios targeting small laboratory animals such as mice, the window length of the short-time analysis window can be selected to be 10 to 30 seconds, the overlap rate between adjacent windows can be selected to be 80% to 98%, and the frequency display range of the STFT time-frequency energy graph can cover 0Hz to 25Hz or be extended to a wider frequency band according to actual needs. It should be noted that the above-mentioned window length, overlap rate, and frequency range are all exemplary preferred parameters. Those skilled in the art can make adaptive adjustments according to the radar sampling rate, the heartbeat frequency range of the target under test, and the actual experimental conditions, so that the fundamental wave and higher harmonics of the heartbeat present clear linear spectral lines on the time-frequency graph.

[0048] Human Body: In human vital sign detection scenarios, millimeter-wave radar can be directed towards the chest, abdomen, back, or other body surface areas that can reflect cardiopulmonary micromovements of the subject, and extract complex IQ slow-time signals from the target's distance gate or detection channel. The fundamental frequency of human respiration is typically located in a lower frequency band, while the fundamental frequency of heartbeat and its higher harmonics are located in a relatively higher frequency band. Therefore, the corresponding STFT display frequency band can be selected based on the physiological frequency range of the subject.

[0049] In one implementation of human body detection, the window length of the STFT short-time analysis window can be selected from 5 seconds to 20 seconds, preferably from 8 seconds to 15 seconds; the overlap rate between adjacent analysis windows can be selected from 70% to 98%, preferably from 80% to 95%; the display frequency range of the STFT time-frequency energy map can be selected from 0Hz to 10Hz, or extended to higher frequencies according to the human body status, sampling rate and required higher-order heartbeat harmonics.

[0050] Step 3: STFT image preprocessing and candidate spectral line extraction like Figure 3 As shown, firstly, a target frequency band containing the fundamental heartbeat and higher-order heartbeat harmonics, such as 4Hz to 25Hz, is extracted; then, background normalization is performed along the frequency and time directions to reduce low-frequency background and global energy drift; next, vertical interference columns are detected. First, a bright pixel determination threshold is set on the preprocessed STFT time-frequency energy map. The bright pixel determination threshold can be determined based on the mean, standard deviation, median, absolute offset median, or a preset quantile of the time-frequency energy within the entire target frequency band. For example, the following can be used:

[0051] As the threshold for determining bright pixels, where... This represents the time-frequency energy within the target frequency band. This represents the threshold adjustment coefficient. For the... Each time column is used to calculate the proportion of frequency points in that time column that are above the bright pixel determination threshold.

[0052] in, This indicates the total number of frequency points within the target frequency band. Indicates the characteristic function. When the ratio of bright pixels... Greater than the preset ratio threshold At that time, the first The time series is initially marked as vertical interference. The proportional threshold... Desirable to Preferred selection to The determination can also be made by combining the total energy of the time series, the frequency coverage width, or the degree of energy surge relative to adjacent time series. To cover the edges and energy diffusion areas of vertical interference, the initially determined vertical interference series can be extended forward and backward by multiple time frames; or the vertical interference series mask can be subjected to dilation, closing operations, or short-interval merging processing. Continuous or adjacent vertical interference series and their extended time series together constitute the vertical interference region. For the vertical interference region, masking or weight reduction can be performed during the extraction of candidate image spectral lines; during spectral line continuation and joint tracking, linear crossing, smooth interpolation, endpoint constraints, or reduced observation scores can be performed based on the effective trajectories before and after the interference region to avoid broadband vertical bright bars being misidentified as heartbeat latent wave spectral lines. Corresponding symbol explanation: Wherein, This indicates the threshold for determining bright pixels; Indicates the first The proportion of bright pixels in each time column; Indicates the number of frequency points within the target frequency band; Indicates the first One frequency point; Indicates the first One time frame; Indicates the threshold for determining the vertical interference ratio; This represents the bright pixel threshold adjustment coefficient.

[0053] In this embodiment, spectral line objects are extracted from the image. The preprocessed time-frequency image is enhanced with multi-angle linear structures using linear convolution kernels in multiple directions, Gabor filters, Frangi-type filters, or morphological methods to enhance near-horizontal or oblique thin line structures. Subsequently, binarization, skeletonization, and connected component analysis are performed to obtain candidate image spectral line objects. Each candidate image spectral line object can be represented as a set of time-frequency points.

[0054]

[0055] in, Represents the i-th candidate image spectral line object; Represents a time frame. Indicates frequency point, This represents the corresponding time-frequency energy. This represents the number of time-frequency points contained in the i-th candidate image spectral line object.

[0056] The center frequency trajectory can be calculated based on this set of time and frequency points. Duration, frequency range, linewidth, background contrast, frequency smoothness, time coverage, slope, and discontinuity.

[0057] Step 4: Extraction of the dominant respiratory frequency trajectory and labeling of candidate spectral lines like Figure 4 As shown, the trajectory of the dominant respiratory frequency is identified within the respiratory frequency band. and based on Construct the respiratory harmonic prediction trajectory, where k is a positive integer. For candidate image spectral line objects... Calculate the median frequency distance between its center frequency trajectory f_i(t) and the predicted trajectories of each order of respiratory harmonics.

[0058]

[0059] in, Represents candidate image spectral line objects The time support range.

[0060] when Less than the preset threshold When this happens, the candidate image spectral line object is marked as a respiratory harmonic spectral line; When a spectral line object is located in the heartbeat candidate frequency band or its object-level features are similar to those of a heartbeat spectral line, it can be marked as a respiratory / heartbeat ambiguous spectral line instead of being directly excluded, so that it can be further determined by the heartbeat spectral line family relationship.

[0061] Step 5: Construct a family of heartbeat harmonic spectral lines based on integer multiple relationships and order normalization consistency. like Figure 5 As shown, different heartbeat harmonic orders are tried for each suspected heartbeat spectral line object, and the order-normalized trajectory is calculated.

[0062]

[0063] If at least two spectral objects have the same order-normalized trajectory within an overlapping time interval, a short-interval contiguous interval, or a contiguous interval separated by a vertical interference region, and their inversion residuals are less than a preset threshold, they are identified as belonging to the same heartbeat harmonic spectral line family. This spectral line family may include the heartbeat fundamental spectrum, or it may only include second-, third-, or higher-order heartbeat harmonic spectra.

[0064]

[0065] Where G represents the heartbeat harmonic spectral line family, TG represents the time support range or time union range of each spectral line object, wi represents the weight determined based on the background contrast, duration and smoothness of the spectral line object, f0(t) represents the fused heartbeat fundamental trajectory after order normalization, and R represents the inversion residual.

[0066] Step 6: Joint tracking based on higher-order heartbeat harmonics after occlusion When vertical interference, broadband transient interference, or spectral line breakage causes discontinuities in the second and third order heartbeat harmonic spectral lines in the later part of the image skeletonization result, the system no longer relies solely on skeleton connected domain fragments. Instead, it establishes a candidate grid of heartbeat fundamental frequencies in the original STFT time-frequency energy map and examines the local spectral peak contrast at the second harmonic position 2f0 and the third harmonic position 3f0 for each candidate fundamental frequency f0.

[0067] In one implementation, for time frames and candidate fundamental frequency The joint observation score is:

[0068] Where C(P,t,f) represents the contrast of the local peak near frequency f in the STFT energy map relative to the local background.

[0069] In this embodiment, the local background refers to the background in the current time frame. Within, based on the frequency to be evaluated Centered on the target frequency range, the time-frequency energy contained within the preset background frequency range after excluding the central spectral peak protection zone. Specifically, in the frequency range... A peak search range, a central protection zone, and a background statistics range are set up nearby. The local peak is the maximum time-frequency energy within the peak search range, and the local background is the median, mean, truncated mean, or preset quantile of the time-frequency energy at each frequency point outside the central protection zone within the background statistics range.

[0070] Indicates the first Time frames At the candidate heart rate The corresponding higher-order heartbeat harmonic joint observation score; In the STFT time-frequency energy diagram P, the first... Time frames and second-order heartbeat harmonic frequency Local peak contrast calculated in the vicinity; In the STFT time-frequency energy diagram P, at the th Time frames and third-order heartbeat harmonic frequency Local peak contrast calculated in the vicinity; Subsequently, adjacent fundamental frequency jump penalties and maximum jump constraints are introduced in the time direction, and the continuous heartbeat fundamental frequency trajectory is obtained through dynamic programming:

[0071] And can be constrained | - '|≤Δfmax. Δfmax can be adaptively determined based on the time interval between adjacent STFT time frames, the rate of change of the heart rate of the target being measured, and the STFT frequency resolution. In one embodiment, Δfmax is 1 to 10 frequency resolution units, or 0.05Hz to 0.50Hz, preferably 0.10Hz to 0.30Hz.

[0072] in, Indicates as of the date The timeframes and the current candidate heartbeat base frequency are The optimal cumulative path score at that time; Indicates the first Candidate heartbeat fundamental frequencies in -1 time frames; Indicates as of the date -1 time frame and candidate heartbeat base frequency is The optimal cumulative path score at that time; The smoothing penalty coefficient represents the change in candidate heartbeat fundamental frequency between adjacent time frames; Indicates the STFT time frame index; This represents the candidate heartbeat baseband frequencies in the current time frame; Δfmax represents the maximum allowed heartbeat baseband frequency jump between adjacent time frames. For time frames marked as vertical interference columns or preset occlusion time windows, the system does not directly utilize the energy score of that time column. Instead, it performs linear traversal or smooth continuation based on the effective trajectories before and after the occlusion area to avoid the vertical bright bars causing deviation of higher-order heartbeat harmonic trajectories. This obtains the heartbeat fundamental wave trajectory. After (t), generate the second harmonic trajectory f2(t)=2. (t) and the third harmonic trajectory f3(t)=3 (t), and incorporate it as an additional spectral line object into the heartbeat harmonic spectral line family.

[0073] Step 7: Robust smoothing of spectral lines and heartbeat fundamental wave inversion and recovery When the fundamental heartbeat wave is clearly visible, the system directly uses the fundamental wave spectral line as the trajectory of the main heartbeat fundamental wave, and uses the consistency of integer multiples of second-order, third-order, and other higher-order heartbeat harmonics as verification. When the fundamental heartbeat wave is weakly visible, intermittent, or obscured by low-frequency background or respiratory harmonics, the system reconstructs the fundamental wave trajectory based on at least one higher-order heartbeat spectral line. For a higher-order heartbeat spectral line of order m, the calculation... The weights are set according to the contrast, duration, smoothness, discontinuity rate and breathing harmonic distance of the spectral lines, and weighted fusion is performed within the time union range of each spectral line object.

[0074] Before or during the recovery of the fundamental wave trajectory, robust smoothing can be performed on the final heartbeat harmonic spectral line family. This robust smoothing does not re-search for STFT peaks or alter the order relationship of the spectral line objects. Instead, it performs local trend estimation, outlier replacement, maximum offset limitation, and maximum slope limitation on the determined higher-order heartbeat harmonic spectral lines to reduce the impact of skeleton f0 sawtooth, local cusps, and segmentation boundaries on the recovered fundamental wave trajectory.

[0075] Step 8: Fundamental wave regression, confidence assessment, and heart rate result output. like Figure 6 and Figure 7 As shown, the system recovers the fundamental wave trajectory. Perform a local peak re-examination in the vicinity. If in ±Δ If a peak exists within the range where the relative local background is higher than the threshold and the frequency deviation is less than the threshold, it is taken as the actual visible fundamental frequency point. Based on the proportion of actual visible fundamental frequency points, inversion residuals, number of heartbeat harmonic spectral line families, and object-level characteristics, the heartbeat fundamental frequency, heart rate, visible harmonic order, fundamental frequency state, and confidence level are output.

[0076] BPM= In one implementation, the confidence score can be determined by the number of spectral family orders, the proportion of actual visible fundamental frequencies, the spectral family inversion residuals, the breathing harmonic proximity penalty, and the quality of the spectral object.

[0077] Score in, This indicates the overall quality score of the spectral line family. This indicates the proportion of the actual visible fundamental frequency. Indicates the number of visible harmonic orders. Indicates the inversion residual. This represents the proximity penalty for respiratory harmonics, with η1 to η5 being preset weighting coefficients. If the score is lower than the preset threshold, a message "Result is unstable" or "Fundamental frequency cannot be reliably confirmed" will be output; if the score reaches the preset threshold, the restored fundamental heart rate and heart rate will be output.

[0078] Combination Figure 6 The results of superimposing the higher-order heartbeat harmonic spectral line family and the recovered fundamental wave trajectory shown can be seen that when the second and third-order heartbeat harmonics are clearer in the STFT time-frequency energy diagram and the heartbeat fundamental wave is weaker, the system can still recover the heartbeat fundamental wave trajectory based on the consistency of the heartbeat harmonic spectral line family and the joint tracking results after occlusion. Figure 7 An example of the output result is given, in which the recovered fundamental heart rate is about 5.78 Hz, corresponding to a heart rate of about 346.8 beats / minute. The visible harmonic orders are second and third, the inversion residual is about 0.01 Hz, and the actual visible proportion of the fundamental frequency is about 62.6%, thus demonstrating that this method can simultaneously provide heart rate results and result reliability indicators.

[0079] Example 2

[0080] like Figure 2 As shown, this embodiment provides a heart rate estimation system based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion. The system includes: The system includes a radar data acquisition module, a signal enhancement and phase extraction module, an STFT time-frequency analysis module, an STFT image preprocessing module, an image spectral line object extraction module, a respiratory spectral line marking module, a heartbeat harmonic spectral line family construction module, a joint tracking module, a spectral line family post-processing module, a heartbeat harmonic inversion and verification module, and a result output module.

[0081] The radar data acquisition module acquires the complex IQ slow-time signal from the millimeter-wave radar at the target range gate or detection channel. The signal enhancement and phase extraction module performs robust center compensation, sliding center compensation, or other phase enhancement processing on the complex IQ slow-time signal and extracts the micro-motion phase signal from the unfolded surface. The STFT time-frequency analysis module calculates the STFT time-frequency energy map.

[0082] The STFT image preprocessing module performs target frequency band truncation, background normalization, vertical interference detection and suppression, and linear structure enhancement. The image spectral line object extraction module binarizes, skeletonizes, and performs connected component analysis on the enhanced time-frequency image, outputting candidate image spectral line objects and their object-level features. The respiratory spectral line marking module extracts the main respiratory trajectory and generates respiratory harmonic prediction trajectories.

[0083] The heartbeat harmonic spectral line family construction module is used to determine the heartbeat harmonic spectral line family based on integer multiple frequency relationships, temporal overlap or continuity, and consistency of order-normalized trajectories among candidate image spectral line objects. The joint tracking module implements tracking through dynamic programming, specifically configured as follows: a candidate grid for the heartbeat fundamental frequency is established; for each candidate fundamental frequency f0, the local spectral peak contrast at its second harmonic position 2f0 and third harmonic position 3f0 is simultaneously examined; according to the formula:

[0084] Calculate the joint observation score, where C(P,t,f) represents the contrast of a local peak near frequency f in the time-frequency energy map relative to the local background; introduce adjacent fundamental frequency jump penalties and maximum jump constraints in the time direction, according to the formula:

[0085] The fundamental trajectory of the continuous heartbeat is solved by dynamic programming, with the constraint |f0 - f0'| ≤ Δf_max; the heart rate is calculated based on the fundamental trajectory of the heartbeat.

[0086] The spectral family post-processing module robustly smooths the spectral lines in the final heartbeat harmonic spectral family. The heartbeat harmonic inversion and verification module recovers the fundamental heartbeat frequency from the higher-order heartbeat harmonic spectral family when the fundamental frequency is weakly visible or invisible, and retrieves the actual visible fundamental frequency point. The results output module outputs the fundamental heartbeat frequency, heart rate, fundamental frequency state, visible harmonic order, inversion residual, actual visible fundamental frequency proportion, and confidence level.

[0087] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, combinations, or improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A heart rate estimation method based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion, characterized in that, Includes the following steps: Step 1: Obtain the millimeter-wave radar complex IQ slow-time signal corresponding to the distance gate or detection channel of the target under test, and perform outlier repair, stable segment truncation and phase enhancement processing on the complex IQ slow-time signal to obtain the surface micro-motion phase signal; Step 2: Perform short-time Fourier transform on the micro-motion phase signal on the body surface to obtain an STFT time-frequency energy map covering the respiratory frequency band, the fundamental frequency band of the heartbeat, and the higher-order heartbeat harmonic frequency band. Step 3: The STFT time-frequency energy map is preprocessed and line structure enhanced as a two-dimensional image, and candidate image spectral line objects are obtained through binarization, skeletonization and connected component analysis. Step 4: Extract the respiratory dominant frequency trajectory and generate the respiratory harmonic prediction trajectory. Based on the positional relationship between the candidate image spectral line object and the respiratory harmonic prediction trajectory, mark the candidate image spectral line object as a respiratory spectral line, respiratory harmonic spectral line, respiratory / heartbeat blurred spectral line, or suspected heartbeat spectral line. Step 5: Search for a family of heartbeat harmonic spectral lines in the candidate image spectral line object set that satisfies integer multiple frequency relationships, time overlap relationships, short interval succession relationships, and order normalized trajectory consistency. Step 6: When vertical interference, broadband transient interference, or spectral line breakage region is detected, joint supplementary tracking is performed in the STFT time-frequency energy map based on at least two higher-order heartbeat harmonic orders, and the higher-order heartbeat harmonic trajectory obtained by supplementary tracking is incorporated into the heartbeat harmonic spectral line family. Step 7: Perform robust smoothing on the higher-order heartbeat harmonic spectral lines in the heartbeat harmonic spectral line family, and invert and recover the heartbeat fundamental wave trajectory based on the order-normalized frequency trajectory of at least one higher-order heartbeat harmonic spectral line object in the heartbeat harmonic spectral line family. Step 8: Look back at the actual visible fundamental point in the STFT time-frequency energy map near the recovered heartbeat fundamental trajectory, and output the heartbeat fundamental frequency, heart rate, fundamental state and confidence level by combining the inversion residual, the proportion of actual visible fundamental waves and the spectral family score.

2. The heart rate estimation method based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion according to claim 1, characterized in that, The phase enhancement processing includes at least one of robust IQ center compensation, sliding center compensation, mean DC removal compensation, differential phase or DAM phase extraction methods; The robust IQ center compensation includes: Iteratively fit the center of the complex IQ trajectory, calculate the radius sequence from the sampling point to the fitted center, remove outliers based on the median of the radius and the median of the absolute deviation of the radius, and iteratively obtain the stable compensation center c*; calculate the compensated complex IQ slow-time signal according to the following formula: ； And calculate according to the following formula: ； Then, x[n] is detrended, mean-removed, or normalized. in, This represents a complex slow-time signal after IQ center compensation; Indicates the first The original millimeter-wave radar complex IQ slow-time signal at each sampling moment; This refers to the micro-motion phase signal of the body surface obtained by extracting and expanding the compensated complex signal, and the unit is usually radians; Indicates phase expansion; This represents the complex argument extraction operation, used to obtain the wrapper phase of a complex signal.

3. The heart rate estimation method based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion according to claim 1, characterized in that, The STFT time-frequency energy map is calculated based on the micro-motion phase signal x[n] on the body surface according to the following formula: ； ； ； in, Represents the complex time-frequency matrix; w[n] is the short-time analysis window; The discrete sampling point index represents the phase signal of the micro-motion on the body surface; Indicates the time frame position or the center position of the analysis window in the STFT; This indicates the STFT analysis frequency or frequency point, in Hz. Indicates the sampling rate; Indicates time position ,frequency The time-frequency energy at that point is equal to the square of the STFT coefficient magnitude; Indicates that it is composed of all The complex time-frequency matrix formed; This represents the time-frequency energy after logarithmic transformation and maximum value normalization, in dB. Indicates that it is composed of all The constructed time-frequency energy matrix; A constant used to prevent overflow in logarithmic operations.

4. The heart rate estimation method based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion according to claim 1, characterized in that, The preprocessing described in step three includes target frequency band interception, frequency-direction background normalization, time-direction background normalization, and vertical interference suppression. Vertical interference suppression involves statistically analyzing the proportion of bright pixels in each time series. When the proportion of bright pixels exceeds a preset threshold, the time series and its neighboring time series are marked as vertical interference regions. These regions are then masked, weighted, linearly crossed, or have endpoints constrained in subsequent image spectral line extraction, spectral family scoring, joint tracking, or fundamental wave trajectory recovery.

5. The heart rate estimation method based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion according to claim 1, characterized in that, The linear structure enhancement described in step three includes at least one of the following: multi-directional linear filter, Gabor filter, morphological top-hat filter, Frangi-type linear structure filter, and skeletonization processing.

6. The heart rate estimation method based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion according to claim 1, characterized in that, Each candidate image spectral line object is composed of a set of time-frequency points in the STFT time-frequency energy map. The method represents and calculates at least one object-level feature, including duration, frequency range, average frequency, center frequency trajectory, linewidth, frequency smoothness, background contrast, time coverage, slope, discontinuity rate, and whether it crosses a vertical interference region; the object-level features are used for candidate spectral line semantic labeling, heartbeat harmonic spectral line family construction, and confidence scoring; wherein, Represents a time frame. Indicates frequency point, This represents the corresponding time-frequency energy.

7. The heart rate estimation method based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion according to claim 1, characterized in that, Step four involves identifying the dominant respiratory frequency trajectory within the respiratory frequency band. and based on Constructing respiratory harmonic prediction trajectories for candidate image spectral line objects The median frequency distance between the center frequency trajectory f_i(t) and the predicted trajectories of each order of respiratory harmonics is calculated using the following formula: ； Where k is a positive integer, Represents candidate image spectral line objects The time support range; when Less than the preset threshold At that time, the candidate image spectral line objects are marked as respiratory harmonic spectral lines; When the spectral line object is located in the heartbeat candidate frequency band or its object-level features are similar to the heartbeat spectral line, it is marked as a respiratory / heartbeat fuzzy spectral line and used as a fuzzy candidate in subsequent comprehensive judgment.

8. The heart rate estimation method based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion according to claim 1, characterized in that, The joint pursuit described in step six is ​​implemented through dynamic programming, specifically including: A candidate grid for the heartbeat fundamental frequency is established. For each candidate fundamental frequency f0, the local spectral peak contrast at its second harmonic position 2f0 and third harmonic position 3f0 is examined simultaneously. The joint observation score is calculated according to the following formula: Where C(P,t,f) represents the contrast of the local peak near frequency f in the time-frequency energy map relative to the local background; Introducing adjacent fundamental frequency hopping penalties and maximum hopping constraints in the time direction, according to the formula: The trajectory of the fundamental wave of the continuous heartbeat is solved by dynamic programming, with the constraint that |f0 - f0'| ≤ Δf_max; Heart rate is calculated based on the fundamental wave trajectory of the heartbeat.

9. A heart rate estimation system based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion, characterized in that, include: The radar data acquisition module is configured to acquire the complex IQ slow-time signal of the target range gate or detection channel; The signal enhancement and phase extraction module is configured to perform phase enhancement processing on the complex IQ slow-time signal and extract the micro-motion phase signal of the body surface. The STFT time-frequency analysis module is configured to calculate the STFT time-frequency energy map of the micro-motion phase signal on the body surface; The STFT image preprocessing module is configured to perform background normalization, vertical interference suppression, and linear structure enhancement on the STFT time-frequency energy map. The image spectral line object extraction module is configured to extract candidate image spectral line objects and object-level features. The respiratory spectral marking module is configured to extract the dominant respiratory frequency trajectory and generate a respiratory harmonic prediction trajectory. The heartbeat harmonic spectral line family construction module is configured to construct the heartbeat harmonic spectral line family based on integer multiple frequency relationships, time continuity relationships, and order normalization trajectory consistency. The joint tracking module is configured to perform joint tracking of the fundamental wave trajectory of the heartbeat based on at least two higher-order heartbeat harmonics after vertical interference, broadband transient interference, or spectral line breakage regions. The spectral line family post-processing module is configured to perform robust smoothing on the higher-order heartbeat harmonic spectral lines in the heartbeat harmonic spectral line family. The heartbeat harmonic inversion and verification module is configured to invert and recover the fundamental wave trajectory of the heartbeat from the heartbeat harmonic spectral line family, and to check the actual visible fundamental wave points. The results output module is configured to output the heartbeat fundamental frequency, heart rate, fundamental frequency state, visible harmonic order, inversion residual, actual visible fundamental frequency ratio, and confidence level.

10. The heart rate estimation method system based on millimeter-wave radar STFT image spectral line extraction and harmonic inversion according to claim 9, characterized in that, The joint tracking module implements tracking through dynamic programming, and its specific configuration is as follows: Establish a candidate grid for the fundamental frequency of heartbeat, and examine the local spectral peak contrast at the second harmonic position 2f0 and the third harmonic position 3f0 for each candidate fundamental frequency f0. The joint observation score is calculated using the following formula: Where C(P,t,f) represents the contrast of the local peak near frequency f in the time-frequency energy map relative to the local background; Introducing adjacent fundamental frequency hopping penalties and maximum hopping constraints in the time direction, according to the formula: The trajectory of the fundamental wave of the continuous heartbeat is solved by dynamic programming, with the constraint that |f0 - f0'| ≤ Δf_max; Heart rate is calculated based on the fundamental wave trajectory of the heartbeat.

Images