A radar respiratory heart rate detection method, device and storage medium based on multi-unit detection and harmonic fusion statistics

Abstract

This invention discloses a radar breathing and heart rate detection method based on multi-unit detection and harmonic fusion statistics. The method includes the following steps: selecting an area of ​​interest (AOI) in a buffer region for subsequent processing; performing phase calculations on all range cells within the AOI to obtain an unwound phase matrix and a phase difference matrix; selecting suitable range cells for breathing and heart rate detection based on the phase matrix and phase difference matrix, and saving them to a breathing index array and a heart rate index array, respectively; performing multi-unit breathing detection on the selected breathing detection range cells and heart rate detection range cells, calculating the confidence level of each breathing detection result and heart rate detection result, and caching the breathing detection results and heart rate detection results of multiple frames and multiple units, respectively; performing statistical analysis on the cached breathing detection results, performing harmonic fusion statistical analysis on the cached heart rate detection results, and smoothing the results of both.

Description

Technical Field

[0001] The present invention relates to the technical field of non-contact human feature detection radar, and in particular to a radar respiratory heart rate detection method, device and storage medium based on multi-unit detection and harmonic fusion statistics. Background Art

[0002] Against the backdrop of an aging population and the continued improvement of the public's health requirements, radar-based non-contact respiratory and heart rate detection technology has received increasing attention. Compared with traditional contact detection technology, non-contact detection is more human-friendly, simple to operate, easy to deploy, and suitable for long-term health monitoring. The main principle of radar-based non-contact detection technology is to sense the rise and fall of the human chest through radar echoes, and further extract the respiratory and heart rate information therein. However, this technology still faces great challenges. Compared with respiratory detection, heart rate detection is much more complicated. There are two main reasons. First, the chest rise and fall caused by heart rate is much lower than that of breathing, and the reflection in signal characteristics is inherently weak. Second, the frequency bands of the various harmonics of breathing and heart rate overlap and the intensity is not weaker than that of heart rate, which seriously interferes with heart rate detection. How to improve the detection rate and accuracy of heart rate is the main problem facing current non-contact detection technology.

[0003] The existing technology mainly has the following disadvantages:

[0004] (1) The existing technology usually only tracks the energy peak position and detects respiration and heart rate based on the echo characteristics of the peak position. However, the energy peak position is not necessarily the position where the respiration and heart rate characteristics are most obvious. In fact, most positions of the human body contain both respiration and heart rate characteristics, but the proportions of respiration and heart rate characteristics are different. The ideal approach is to select the position with the most prominent respiration characteristics to detect respiration, and select the position with the most prominent heart rate characteristics to calculate the heart rate.

[0005] (2) Existing technologies often only focus on the fundamental frequency and treat the harmonic frequency as interference. Although the existence of harmonics does greatly increase the complexity of heart rate detection, harmonics also contain additional information of valid signals. When harmonics cannot be avoided, if this additional information can be fully utilized, the heart rate detection effect can also be improved. Summary of the invention

[0006] Technical purpose: In view of the shortcomings of existing non-contact radar detection of breathing and heart rate, the present invention discloses a radar breathing and heart rate detection method, equipment and storage medium based on multi-unit detection and harmonic fusion statistics, which can improve the radar breathing and heart rate detection effect, and focus on improving the detection rate and accuracy of heart rate detection.

[0007] Technical solution: To achieve the above technical objectives, the present invention adopts the following technical solution:

[0008] A radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics, the method comprising the following steps:

[0009] Use radar to collect target vital signs information and obtain the returned radar echo signal;

[0010] Receive ADC data frame by frame from the radar echo signal and perform FFT operation, and store the FFT operation results in a buffer of a specified length. The buffer contains a number of distance units, and each distance unit contains a number of frames;

[0011] In the buffer area, an area of ​​interest (AOI) is selected based on energy characteristics for subsequent processing, and an AOI contains multiple distance units;

[0012] The phase of all distance units in the AOI is calculated in the direction of frame number growth to obtain the unwrapping phase matrix and phase difference matrix;

[0013] According to the phase matrix and the phase difference matrix, the distance units suitable for detecting breathing and heart rate are respectively selected and saved in the breathing index array and the heart rate index array respectively;

[0014] Perform multi-unit breathing detection on the screened breathing detection distance units, calculate the confidence of each breathing detection result, and cache the breathing detection results of multiple frames and multiple units;

[0015] Perform statistical analysis on the cached breathing detection results and perform result smoothing;

[0016] Perform multi-unit heart rate detection on the screened heart rate detection distance units, calculate the confidence of each heart rate detection result, extract the heart rate and its harmonic information, and cache the multi-frame and multi-unit detection results;

[0017] Perform harmonic fusion statistical analysis on the cached heart rate detection results and perform result smoothing;

[0018] The final breathing and heart rate detection results are output to complete the entire detection process.

[0019] An electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, a radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics as described above is implemented.

[0020] A computer-readable storage medium stores computer-executable instructions for executing the radar respiration and heart rate detection method based on multi-unit detection and harmonic fusion statistics as described above.

[0021] Beneficial effects: The radar respiratory heart rate detection method, device and storage medium based on multi-unit detection and harmonic fusion statistics provided by the present invention have the following beneficial effects:

[0022] (1) The present invention is not limited to single-point detection at the energy peak, but performs multi-unit detection within the AOI, and selects appropriate detection units for respiration and heart rate according to the differences in different location characteristics, and measures respiration and heart rate separately. This method improves the detection effect of respiration and heart rate, effectively reduces the interference of respiratory harmonics on heart rate detection, improves the heart rate detection rate and accuracy, and enhances the stability of heart rate detection.

[0023] (2) Compared with the prior art that completely treats harmonics as interference, the present invention makes full use of the information in the harmonics. By fusing the information of multiple frames, multiple units, and multiple harmonics and adopting a weighted statistical strategy, the heart rate harmonics are folded into the fundamental wave interval to enhance the heart rate characteristics, thereby further improving the heart rate detection effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art are briefly introduced below.

[0025] Figure 1 It is a flow chart of the radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics of the present invention;

[0026] Figure 2 It is a structural schematic diagram of the buffer area of ​​the present invention;

[0027] Figure 3 This is a schematic diagram of AOI selection detection of the present invention;

[0028] Figure 4 A schematic diagram of the corresponding relationship between the detection window and the target of the sliding window weighted statistics of the present invention;

[0029] Figure 5 Schematic diagram of heart rate detection points in an embodiment of the present invention;

[0030] Figure 6 Schematic diagram of respiratory monitoring points in an embodiment of the present invention;

[0031] Figure 7 This is a schematic diagram of the heart rate output result of a single round of sliding window histogram statistics in an embodiment of the present invention;

[0032] Figure 8 Schematic diagram of heart rate output after harmonic fusion using sliding window histogram statistics in an embodiment of the present invention;

[0033] Fig. 9This is a schematic diagram of heart rate output without harmonic fusion and using sliding window histogram statistics;

[0034] Fig.10 This is a schematic diagram of a heart rate result finally output after smoothing in an embodiment of the present invention;

[0035] Fig.11 Schematic diagram of breathing output using sliding window histogram statistics in an embodiment of the present invention;

[0036] Fig.12 It is a schematic diagram of the breathing result finally output after smoothing in an embodiment of the present invention. DETAILED DESCRIPTION

[0037] The present invention will be described more clearly and completely below by way of a preferred embodiment in conjunction with the accompanying drawings, but the present invention is not limited to the scope of the embodiment.

[0038] In one embodiment, Figure 1 As shown, a radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics is provided. In this embodiment, a set of data of about 4 and a half minutes with an average heart rate of about 77 is cited. The method includes:

[0039] S101. Collect target vital signs information using radar and obtain returned radar echo signals.

[0040] The target's vital signs information includes breathing, heart rate and other information, and the returned radar echo signal is an analog signal carrying this information.

[0041] In one embodiment, a receiver in a radar system receives a radar echo signal and converts the analog signal into a digital signal through an analog-to-digital converter to generate ADC data.

[0042] In one embodiment, the frame rate f of the radar system st =20Hz, frame rate is the number of signal frames collected by the radar system in a unit of time (usually one second). A frame rate of 20Hz means collecting 20 frames of data per second. The higher the frame rate, the higher the frequency of data collection, which can capture more detailed dynamic changes. For example, for a frame rate of 20Hz, data is collected 20 times per second, which can better reflect changes in breathing and heart rate.

[0043] S102, receiving ADC data frame by frame from the radar echo signal and performing FFT operation, and storing the FFT operation result in a buffer of a specified length, the buffer containing a number of distance units, each distance unit containing a number of frames.

[0044] Since the frequency of breathing and heart rate is much lower than the radar frame rate, it is necessary to cache data frames of sufficient length to obtain good frequency resolution and accuracy. Generally, longer data lengths can often achieve higher signal-to-noise ratio gains. Of course, the longer the data, the better. The longer the data, the higher the cache delay.

[0045] like Figure 2 As shown, dbox is a buffer with a length of DBOX_LEN and a width of DBOX_WIDTH. DBOX_LEN corresponds to the maximum number of cached frames, and DBOX_WIDTH corresponds to the number of distance units of a single frame. Both are given constants, and the actual values ​​are related to the specific system design. In this embodiment, the maximum length DBOS_LEN of the buffer dbox is 512, and the width DBOX_WIDTH of the buffer dbox is 48, that is, there are 48 distance units in total.

[0046] In one embodiment, step S102 specifically includes the following steps:

[0047] S201 , receiving ADC data of several points in each frame.

[0048] In one embodiment, N0 points of raw ADC data are received, where N0 is 200.

[0049] S202: Perform windowing processing on the number of ADC data points in each frame.

[0050] Before performing FFT, windowing is performed on the ADC data to reduce spectrum leakage. Commonly used windows include Hanning window, Hamming window, etc. In one embodiment, a Hanning window is added to the ADC data.

[0051] S203 , performing one-dimensional FFT on the windowed ADC data, and extracting distance unit data corresponding to the distance range of target detection from the FFT operation result.

[0052] Perform N1-point FFT operation on the windowed ADC data to get the result OUT 1d , N1 is 256, so the FFT operation result OUT 1d is an array of complex numbers with a length of 256.

[0053] S204: Store the extracted distance unit data in a buffer.

[0054] In one embodiment, to save memory, there is no need to convert OUT 1d All units in the buffer are cached. Usually, the target detection range of respiratory and heart rate applications is sufficient in [0.3, 2] meters, so only OUT 1dThe unit corresponding to the interval can be saved in the buffer dbox. In this embodiment, the distance units corresponding to the distance range are the 8th to 55th units, that is, the number of units to be cached N r is 48, that is, DBOX_WIDTH is 48.

[0055] In one embodiment, the buffer caches at most DBOX_LEN frames of FFT results, and when the number of frames exceeds DBOX_LEN, the earlier cached data is deleted. dbox Record the current cache length of the buffer dbox, N dbox The initial value is zero. Each time a set of data is stored, N dbox Increase by 1 until the length DBOX_LEN is reached and no longer increases. After that, whenever a new set of data arrives, first delete the oldest set of data in the buffer, and then append the new data to the end of the buffer.

[0056] S103 , selecting an AOI in the buffer area based on energy characteristics for subsequent processing, where one AOI includes multiple distance units.

[0057] When the data buffer length reaches DBOX_LEN, the first algorithm calculation will be triggered, and then the algorithm calculation will be triggered with ΔT as the period. In this embodiment, ΔT is 1 second. Since the scheme is a detection algorithm based on the AOI area, the first step of the algorithm is to select the AOI area. There are many methods for selecting AOI, which can be based on peak selection, that is, selecting a continuous interval from a peak unit to a unit lower than the peak energy ΔE0dB, or based on background noise selection, that is, selecting a continuous interval higher than the noise level ΔE1dB and including the peak, or other methods that can obtain high-energy features and continuous areas. The present invention selects the AOI area based on energy.

[0058] In one embodiment, step S103 specifically includes the following steps:

[0059] S301 , calculating the RMS value of each distance unit in the buffer in the direction of increasing frame number to obtain an energy feature array.

[0060] In one embodiment, the energy characteristic may be root mean square (RMS), signal average amplitude, or amplitude mean square, etc.

[0061] For each distance unit in dbox, select a row of data along the frame number direction, such as Figure 2 For the dark part, calculate its RMS value. After all distance units are calculated, we get an E with a length of DBOX_WIDTH. list Array, E list The array represents the energy feature array, and each element of the array represents the signal energy level of the corresponding distance unit.

[0062] S302, searching for the maximum peak value in the energy feature array, and obtaining a peak unit index and a peak unit energy index.

[0063] This embodiment selects the AOI area based on the peak value, and selects the continuous interval between the peak unit and the unit below the peak energy ΔE0dB as the AOI. list Find the maximum peak in the array and get the peak unit index idx max and peak unit energy rms max .

[0064] S303, calculating the lower limit of the energy threshold, and finding the first distance unit below the lower limit of the energy threshold from the peak unit index low distance unit direction, and determining the starting distance unit of the AOI as the first distance unit below the lower limit of the energy threshold plus 1.

[0065] In one embodiment, the energy level lower than the peak energy ΔE0dB is set as the energy threshold lower limit th rms , th rms =(rms max -ΔE0), in this embodiment, the value of ΔE0 is 10dB.

[0066] In one embodiment, from idx max Low distance unit direction finds the first one below th rms The unit index s, then the starting distance unit index idx of AOI s The value is s+1; if the energy in the low-distance unit direction is higher than th rms , then the starting distance unit index idx s Set to starting index 1.

[0067] S304 , searching for the first distance unit below the lower limit of the energy threshold from the peak unit index toward the high distance unit direction, and determining the end distance unit of the AOI as the first distance unit below the lower limit of the energy threshold minus 1.

[0068] From idx max Start looking towards the high distance unit for the first cell below th rms The unit index e is then the end distance unit index idx of the AOI e The value is e-1; if the energy in the low-distance unit direction is higher than th rms , then the end distance unit index idx e Set to the end index DBOX_WIDTH.

[0069] S305 , if the number of distance units included in the obtained AOI exceeds a preset value, a unit clipping operation is performed, and the starting distance unit and the ending distance unit are adjusted according to the energy level until the distance units included in the AOI meet the preset distance unit requirements.

[0070] In one embodiment, if the number of distance units contained in the AOI obtained in the above step exceeds W aoi In this embodiment, W aoi If the setting value is 15, then unit clipping needs to be performed. The basic method of unit clipping is to compare idx s and idx e The energy of the position, if idx s Position energy is less than or equal to idx e Position energy, then idx s Increase by 1, otherwise idx e Decrease by 1, and the loop iterates until idx e and idx s The difference is W aoi -1.

[0071] AOI inspection diagram Figure 3 As shown in the figure, the red line part in the upper sub-figure represents the minimum energy threshold, the selected distance unit is circled in the box, and the box in the lower sub-figure circles the complete AOI area, corresponding to all the detection data of the selected distance unit in dbox.

[0072] S104, performing phase calculation on all distance units in the AOI in the direction of frame number growth to obtain an unwrapping phase matrix and a phase difference matrix.

[0073] In one embodiment, the heart rate is calculated based on the phase difference matrix dpbox, and the breathing rate is calculated based on the phase matrix pbox. At the same time, the phase matrix pbox is also the main basis for subsequent unit screening. The goal of this step is to calculate the phase difference matrix dpbox and the phase matrix pbox.

[0074] In one embodiment, the specific steps of step S104 are as follows:

[0075] S401, calculating the phase of each distance unit in the AOI frame by frame, and obtaining a phase matrix with winding.

[0076] The phase of the frame data in each distance unit of AOI is calculated using the four-quadrant inverse tangent method to obtain the phase matrix with winding abox. The size of the matrix is ​​W aoi ×DBOX_LENGTH, the phase value range in abox is (-π,π].

[0077] S402 , calculating the phase difference for each distance unit in the band-wrapped phase matrix frame by frame to obtain a band-wrapped phase difference matrix.

[0078] For each distance unit, starting from the second frame, according to the direction of frame number growth, the phase of the current frame is subtracted from the phase of the previous frame, and the phase difference is calculated frame by frame to obtain the winding phase difference matrix dabox. The calculation method of the phase difference dabox(i,j) of the jth frame of the i-th distance unit is:

[0079]

[0080] Where i represents the range unit number, j represents the frame number, dabox(i,j) represents the phase difference of the jth frame of the ith range unit, abox(i,j) represents the phase of the jth frame of the ith range unit, and abox(i,j-1) represents the phase of the j-1th frame of the ith range unit.

[0081] S403 , performing unwrapping processing on the phase difference in the wrapped phase difference matrix to obtain an unwrapped phase difference matrix.

[0082] The value range of the elements in abox is (-π,π], from which it can be inferred that the value range of the elements in dabox is (-2π,2π). At a sampling rate of 20Hz, the phase difference caused by the respiratory heart rate is usually much lower than π. Therefore, it is necessary to unwrap the elements in dabox whose absolute values ​​exceed π to obtain the unwrapped phase difference matrix dpbox. The unwrapped phase difference dpbox(i,j) of the j-th frame element in the i-th distance unit is calculated as follows:

[0083]

[0084] Where i represents the range unit number, j represents the frame number, dabox(i,j) represents the phase difference of the jth frame element in the ith range unit, and dpbox(i,j) represents the unwrapped phase difference of the jth frame in the ith range unit.

[0085] S404 , integrating the phase difference in the phase difference matrix frame by frame to obtain a dewrapped phase matrix.

[0086] For each distance unit in the phase difference matrix dpbox, the integral operation is performed in the direction of the frame number growth to obtain the phase matrix pbox. The calculation formula for the j-th frame element in the i-th distance unit in the phase matrix pbox is:

[0087]

[0088] Among them, k represents the frame number, dpbox(i,k) represents the k-th frame of the i-th distance unit in the phase difference matrix dpbox, and pbox(i,j) represents the j-th frame element of the i-th distance unit in the phase matrix pbox.

[0089] S105 . Filter out distance units suitable for detecting respiration and heart rate according to the phase matrix and the phase difference matrix, and save them in the respiration index array and the heart rate index array, respectively.

[0090] This step selects appropriate distance units for heart rate calculation and respiration calculation respectively. In one embodiment, step S105 specifically includes the following steps:

[0091] S501, calculating the phase mean and phase standard deviation of each distance unit in the phase matrix to obtain a phase mean array and a phase standard deviation array.

[0092] Calculate the phase mean for each distance unit of the phase matrix pbox in turn to obtain the phase mean array mP list , calculate the phase standard deviation for each distance unit of the phase matrix pbox in turn, and get the phase standard deviation array sdP list .

[0093] S502: Filter out distance units with mean values ​​lower than a threshold value based on the phase mean array as stable distance units and save them in a stable unit array.

[0094] Based on mPlist, select the values ​​with a mean lower than G mean The distance unit is saved in the array RBIN stable If no lower than G is found, mean If the distance unit is , then select the 5 smallest values ​​in mPlist and save their corresponding distance units to RBIN stable middle.

[0095] S503, selecting distance units suitable for breathing detection from the stable distance unit array based on the phase standard deviation array, and saving them in the breathing index array;

[0096] In RBIN stable In the corresponding distance unit, based on sdP list Further, the phase standard deviation is higher than the threshold G. std0 The distance unit is selected and the final filtered distance unit is saved to the breathing index array B idx If there is no distance unit that meets the conditions, then B idx The maximum peak cell index is saved in .

[0097] S504 , selecting distance units suitable for heart rate detection from the stable distance unit array based on the phase standard deviation array, and saving them in the heart rate index array.

[0098] In RBIN stable In the corresponding distance unit, based on sdP list Further select the phase standard deviation within the threshold interval [G std1 ,G std2 ] and save the final filtered distance units to the heart rate index array H idx If there is no distance unit that meets the conditions, then H idx The maximum peak cell index is saved in .

[0099] Among them G mean , G std0 represents the pre-set threshold, [G std1 ,G std2 ] indicates a preset threshold range. The threshold and threshold range are selected based on the actual application situation, and the selection is different in different situations; RBIN stable Represents a stable cell array.

[0100] S106, performing multi-unit breathing detection on the screened breathing detection distance units, calculating the confidence of each breathing detection result, and caching the breathing detection results of multiple frames and multiple units; performing statistical analysis on the cached breathing detection results, and performing result smoothing processing.

[0101] Since the breathing characteristics are usually more obvious, that is, the signal-to-noise ratio is higher, the peak counting method or the frequency domain analysis method can be used for breathing detection. The peak counting method has higher real-time performance, but also has higher requirements for the signal-to-noise ratio, while the frequency domain analysis method has a certain delay, but can provide a signal-to-noise ratio gain. The heart rate fluctuation characteristics are usually lower, so the frequency domain analysis method is generally used. In one embodiment, the steps of calculating breathing using the frequency analysis method are as follows:

[0102] S601. Design a bandpass filter according to the target respiratory frequency range, and then use the bandpass filter to filter the data corresponding to each distance unit in the respiratory index array to obtain the filtered data of each distance unit and save it in the respiratory data matrix.

[0103] Bandpass filter F for breathing detection br The design of the bandpass filter F needs to consider the target respiratory frequency range, and the filter parameters are set according to actual needs. br B idx A set of data corresponding to each distance unit in is passed through a bandpass filter F br Perform filtering.

[0104] S602 , performing windowing processing on the filtered data of each distance unit in the respiratory data matrix, and then performing FFT to obtain the amplitude spectrum of each distance unit and save it in the respiratory amplitude spectrum matrix.

[0105] The filtered data is subjected to windowing processing, such as using a Hanning window, and then an FFT operation is performed to obtain an amplitude spectrum matrix.

[0106] S603, calculating the noise array of the amplitude spectrum corresponding to each distance unit in the breathing amplitude spectrum matrix, and the noise arrays of all distance units constitute a breathing noise matrix.

[0107] The noise is calculated for the amplitude spectrum corresponding to each distance unit in the amplitude spectrum matrix in turn to obtain a noise matrix.

[0108] S604, searching for the peak position in the amplitude spectrum corresponding to each distance unit of the breathing amplitude spectrum matrix, calculating the signal-to-noise ratio of each peak based on the breathing noise matrix, and saving several peaks with the highest signal-to-noise ratios and their signal-to-noise ratios in the breathing detection result buffer.

[0109] The detection results of each distance unit include peak frequency and signal-to-noise ratio, which are saved as a set of data in B info in the buffer.

[0110] The peak value of the amplitude spectrum corresponding to each distance unit in the amplitude spectrum matrix is ​​searched in turn, and the peak position is interpolated to improve the accuracy of peak detection.

[0111] Calculate the signal-to-noise ratio for all detected peaks and sort them in descending order of signal-to-noise ratio. br The peak values ​​and their signal-to-noise ratios are saved in B info In the buffer zone, due to the relatively obvious breathing characteristics, usually N br Set it to 1. When all distance units have completed the calculation, B info The length of B increases by 1, that is, info Each unit of contains the respiratory peak detection results of all distance units in one calculation.

[0112] S605 , performing sliding window histogram statistical analysis on the data in the breathing detection result buffer area, finding the frequency area with the most concentrated number and taking the mean of all breathing results in the area as the breathing statistical result.

[0113] Set the window size and sliding step size, and use the sliding window to info The detection results in the buffer are counted, and the number of detection results in each window is counted to form a sliding window histogram.

[0114] The detection results are weighted by quantity. The results with high signal-to-noise ratio have larger weights to improve the accuracy of the statistical results. The weighted results are statistically analyzed by sliding window histogram. The weighted detection results are statistically analyzed by sliding window to find the frequency area with the largest number, as shown in the attached figure. Fig.11 The figure shows a schematic diagram of respiratory output using sliding window histogram statistics. Since the respiratory characteristics are much stronger than the heart rate, the statistical results are relatively smooth and there are basically no abnormal points.

[0115] Different from the conventional histogram, each unit does not overlap, but there is overlap between the units of the sliding window histogram. The sliding window histogram is counted according to the sliding window, and each unit corresponds to a window. Because there is an overlap between the windows, there is also an overlap between the corresponding histogram units. Figure 4 The correspondence between the sliding window statistical detection window and the target is shown.

[0116] S606: Smoothing the statistical results to obtain the final breathing detection results.

[0117] Combine historical data to smooth the test results and reduce fluctuations in a short period of time. Fig.12 The figure shows a schematic diagram of the breathing detection result after smoothing. It can be seen that once the initial breathing tracking is completed, a stable breathing value can be output continuously in the subsequent step.

[0118] In one embodiment, weighted averaging or other smoothing algorithms are used to obtain smoothed detection results.

[0119] Find the frequency region with the largest number in the frequency histogram after smoothing, calculate the mean of all results in the frequency region, combine the historical detection data, further smooth the results, and obtain the final breathing detection result.

[0120] S107, perform multi-unit heart rate detection on the screened heart rate detection distance units, calculate the confidence of each heart rate detection result, extract the heart rate and its harmonic information, and cache the multi-frame and multi-unit detection results; perform harmonic fusion statistical analysis on the cached heart rate detection results, and perform result smoothing.

[0121] The steps of calculating heart rate using frequency analysis are similar to those of calculating respiration, except that different filters are used. Because the present invention requires the use of heart rate harmonics to improve detection, the filter design needs to consider the passability of its heart rate harmonics. The specific steps are as follows:

[0122] S701. Design a bandpass filter according to the target heart rate range, and then use the bandpass filter to filter the data corresponding to each distance unit in the heart rate index array to obtain the filtered data of each distance unit and save it in the data matrix.

[0123] According to the passband range [BPSD min ,3*BPSD max ] Design filter F hr , where BPSD min Is the minimum heart rate to pay attention to, BPSD max is the maximum value of the heart rate. After the filter design is completed, H idx A set of data corresponding to each distance unit in is passed through a bandpass filter F hr The filtered data matrix is ​​recorded as SBOX, and each column of SBOX represents the filtered data of one distance unit.

[0124] S702 , performing windowing processing on the filtered data of each distance unit in the data matrix, and then performing FFT to obtain a heart rate amplitude spectrum matrix.

[0125] Perform FFT operations on each column in SBOX in turn. Before calculating the frequency, you need to apply a Hanning window to the data. Here, the number of FFT points is set to 1024. After the FFT operation is completed, the frequency interval [BPSD min ,3*BPSD max ] The corresponding frequency points are saved as the amplitude spectrum matrix ASPECT, where each column of ASPECT corresponds to the amplitude spectrum of a distance unit.

[0126] S703, calculating the noise array of the amplitude spectrum corresponding to each distance unit in the heart rate amplitude spectrum matrix, and the noise arrays of all distance units constitute the heart rate noise matrix.

[0127] The noise array is calculated column by column for the amplitude spectrum corresponding to each distance unit in ASPECT, and the noise arrays of all distance units constitute the noise matrix NoiseMatrix. Assuming that the amplitude spectrum of a distance unit is the array Spect, then the noise Noise(i) of the ii-th unit in the noise array Noise is calculated as:

[0128]

[0129] Among them, Noise l (ii) represents the noise level of the ii-th unit in the left window, Noise r(ii) represents the noise level of the ii-th unit in the right window, the GetOrder(l,r) operator represents the output of the value indexed as round(len*r) after arranging the sequence l with a sequence length of len from small to large, Spect(ii-guard-wlen:ii-guard) represents a new array composed of the elements from index ii-guard-wlen to index ii-guard in the array Spect, and the meaning of Spect(ii+guard:ii+guard+wlen) is similar to it; ratio represents the ratio, which defines the ratio used to select the value in the window as the noise value, and is usually used to obtain the median or other statistical values ​​in the window, wlen represents the window length, which defines the window size used to calculate the noise, guard represents the guard unit, which defines the number of frequency points skipped from the current ii unit to the left and right sides to avoid the detection result being interfered by the frequency near the current unit; in this embodiment, guard=1, wlen=20, and ratio=0.4.

[0130] That is, for the current unit ii, define the left and right window ranges as ii-guard-wlen to ii-guard and ii+guard to ii+guard+wlen respectively, then in the left window, sort all the values ​​in the window in ascending order, and select the value at the ratio position as Noise l (ii) Calculate Noise in the same way in the right window r (ii), finally take Noise l (ii) and Noise r The larger value in (ii) is used as the noise Noise(ii) of the current unit ii.

[0131] S704, searching for the peak position in the amplitude spectrum corresponding to each distance unit of the heart rate amplitude spectrum matrix, calculating the signal-to-noise ratio of each peak based on the heart rate noise matrix, and saving several peaks with the highest signal-to-noise ratios and their signal-to-noise ratios in the heart rate detection result buffer.

[0132] Perform peak detection on the amplitude spectrum Spect corresponding to each distance unit in ASPECT column by column, and calculate the signal-to-noise ratio of each peak value in combination with the noise matrix calculated in the previous step. Finally, select the N peaks with the highest signal-to-noise ratio. hr The peak values ​​and their signal-to-noise ratios are saved in H info in the buffer.

[0133] In this case, N hr The value is 6, Figure 5This is the complete heart rate detection graph of the data set in this case. The horizontal axis in the figure represents the number of rounds of heart rate statistics. Since this case is calculated once per second, the unit of the horizontal axis can also be understood as seconds. The vertical axis is the heart rate result detected directly from the spectrum. The actual heart rate is between 75 and 85. Obviously, there are a large number of interference points in the direct detection graph, such as respiratory harmonics and the second and third harmonics of heart rate.

[0134] S705. Perform harmonic fusion statistical analysis on the heart rate detection results cached in the heart rate detection result buffer area, fold the harmonic frequency into the fundamental frequency interval, and find the frequency area with the most concentrated number through sliding window weighted histogram statistics, and use the mean of all heart rate results in the area as the heart rate statistical result.

[0135] In this embodiment, each time the peak detection system in step S704 is executed once, a set of results will be cached in H info In the process, when the number of cached results reaches N sts In this embodiment, N sts When the value is 30, fusion statistics are started. The so-called fusion statistics is to fuse the harmonic and fundamental information of multiple frames and multiple units to calculate the heart rate result. The specific steps are as follows:

[0136] Remove H info The most recent N sts The group detection results are weighted and the quantity of each detection point is weighted. The weighting method is as follows. The results are saved in the weighted detection array D w middle.

[0137]

[0138] Wherein, Cnt represents the detection point, and SNR represents the signal-to-noise ratio.

[0139] Right D w The median frequency is higher than 1.9 times the minimum heart rate BPSD min and less than 2.1 times BPSD max The points are folded twice to obtain the secondary weighted folding matrix D wf2 , corresponding to D w Higher than 2.9 times BPSD min and less than 3.1 times BPSD max The points are folded three times to obtain the three-fold weighted folding matrix D wf3 , and D w , D wf2 and D wf3 Combine to form the complete weighted folding matrix D wf .

[0140] The quadratic weighted folding matrix is ​​calculated as follows:

[0141] D wf2· bps = D w· bps / 2,D w· bps>1.9*BPSD min ∩D w· bps<2.1*BPSD max

[0142] The cubic weighted folding matrix is ​​calculated as follows:

[0143] D wf3· bps = D w· bps / 3,D w· bps>2.9*BPSD min ∩D w· bps<3.1*BPSD max

[0144] In the fold weight matrix D wf Before performing quantitative statistics, we first remove points with a signal-to-noise ratio less than 2.5, and then perform quantitative statistics on the remaining data. The unit interval length of the histogram statistics is selected as 8, and the difference between each unit is 1, that is, if the heart rate interval corresponding to the first unit is [60,67], the interval corresponding to the second unit is [61,68], the interval corresponding to the third unit is [62,69], and so on. After completing the weighted histogram statistics, the heart rate mean corresponding to the number peak of each round of statistics is output as the heart rate statistics result.

[0145] Attached Figure 7 This is one of the test results of the data in this embodiment. It can be seen from the statistical chart that the maximum peak in the weighted statistical chart can correctly reflect the heart rate result, and the true value is about 77, while the second largest peak corresponds to the respiratory harmonics. It can also be seen from the figure that respiratory harmonics are the main factor interfering with heart rate detection.

[0146] To compare the improvement brought by harmonic fusion, Figure 8 and attached Fig. 9 The weighted histogram heart rate output results using harmonic fusion and the weighted histogram heart rate output results without harmonic fusion are given in the attached Fig. 9 It can be seen that the weighted histogram statistics without harmonic fusion has enhanced respiratory harmonic interference, and about 30% of the data are erroneous due to respiratory harmonic interference. It can be seen that the harmonic fusion solution can significantly improve the accuracy of heart rate detection.

[0147] S706: Smoothing the heart rate statistical results to obtain a final heart rate detection result.

[0148] Since interference is unavoidable and the human body cannot be completely still during the test, all direct test results will inevitably be wrong. At this time, the final test results can be improved through smoothing strategies. The smoothing process is mainly divided into three parts:

[0149] a) Starting point selection

[0150] The results of multiple consecutive heart rate tests are counted. If the number of points in the interval [-σ, +σ] reaches N start , then determine the starting heart rate value fhr, where σ represents the heart rate value.

[0151] b) Data retention

[0152] If the current heart rate detection result is frh cur Relative to the tracking value fhr last If the jump exceeds the threshold, the current detection results will not be fused immediately, and the final result will remain fhr last , the maximum number of steps to be maintained is N keep , where the tracking value fhr last Refers to the last heart rate detection result. The jump threshold is set according to the actual measurement. The maximum number of steps to be maintained is N keep It means if the current heart rate detection result is frh cur The number of jumps exceeding the threshold is greater than the maximum number of hold steps N keep After that, the heart rate detection results are allowed to be updated.

[0153] c) Alpha filtering

[0154] fhr=α*frh cur +(1-α)*fhr last

[0155] Wherein, α represents the filtering factor. In this embodiment, N start and N keep Set them to 12, and get the final heart rate smoothing output as shown in the attached figure. Fig.10 As shown, it can be seen that once the initial heart rate tracking is completed, a stable heart rate value can be output continuously thereafter.

[0156] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics, characterized in that: The method comprises the following steps: Use radar to collect target vital signs information and obtain the returned radar echo signal; Receive ADC data frame by frame from the radar echo signal and perform FFT operation, and store the FFT operation results in a buffer of a specified length. The buffer contains a number of distance units, and each distance unit contains a number of frames; In the buffer area, AOI is selected based on energy characteristics for subsequent processing, and one AOI contains multiple distance units; The phase of all distance units in the AOI is calculated in the direction of frame number growth to obtain the unwrapping phase matrix and phase difference matrix; According to the phase matrix and the phase difference matrix, the distance units suitable for detecting breathing and heart rate are respectively selected and saved in the breathing index array and the heart rate index array respectively; Perform multi-unit breathing detection on the screened breathing detection distance units, calculate the confidence of each breathing detection result, and cache the breathing detection results of multiple frames and multiple units; Perform statistical analysis on the cached breathing detection results and perform result smoothing; Perform multi-unit heart rate detection on the screened heart rate detection distance units, calculate the confidence of each heart rate detection result, extract the heart rate and its harmonic information, and cache the multi-frame and multi-unit detection results; Perform harmonic fusion statistical analysis on the cached heart rate detection results and perform result smoothing; The final breathing and heart rate detection results are output to complete the entire detection process.

2. According to claim 1, a radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics is characterized in that: The receiving ADC data frame by frame from the radar echo signal and performing FFT operation, and storing the FFT operation result in a buffer of a specified length includes: Receive ADC data of several points per frame; Perform windowing on the ADC data of each frame; Perform one-dimensional FFT on the windowed ADC data, and extract the distance unit data corresponding to the distance range of target detection from the FFT operation result; Store the extracted distance cell data into a buffer.

3. The radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics according to claim 1 is characterized in that: The selecting of AOI based on energy characteristics for subsequent processing in the buffer area includes: For each distance unit in the buffer, the RMS value is calculated in the direction of increasing frame number to obtain the energy feature array; Find the maximum peak in the energy feature array and obtain the peak unit index and peak unit energy index; Calculate the lower limit of the energy threshold, and find the first distance unit below the lower limit of the energy threshold from the peak unit index low distance unit direction, and determine the starting distance unit of the AOI as the first distance unit below the lower limit of the energy threshold plus 1; Starting from the peak unit index, search for the first distance unit below the lower limit of the energy threshold in the direction of the high distance unit, and determine that the end distance unit of the AOI is the first distance unit below the lower limit of the energy threshold minus 1; If the number of distance units contained in the obtained AOI exceeds the preset value, a unit clipping operation is performed, and the starting distance unit and the ending distance unit are adjusted according to the energy level until the distance units contained in the AOI meet the preset distance unit requirements.

4. The radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics according to claim 1 is characterized in that: The phase of all distance units in the AOI is calculated in the direction of frame number growth, and the unwrapped phase matrix and phase difference matrix are obtained, including: Calculate the phase of each distance unit in the AOI frame by frame to obtain a phase matrix with wrapping; The phase difference is calculated frame by frame for each distance unit in the band-wrapped phase matrix to obtain a band-wrapped phase difference matrix; Performing an unwrapping process on the phase difference in the wrapped phase difference matrix to obtain an unwrapped phase difference matrix; The phase difference in the phase difference matrix is ​​integrated frame by frame to obtain an unwrapped phase matrix.

5. The radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics according to claim 1 is characterized in that: According to the phase matrix and the phase difference matrix, the distance units suitable for detecting breathing and heart rate are selected respectively, and saved in the breathing index array and the heart rate index array respectively, including: Calculate the phase mean and phase standard deviation of each distance unit in the phase matrix to obtain a phase mean array and a phase standard deviation array; Based on the phase mean array, the distance units with mean values ​​lower than the threshold are selected as stable distance units and saved in the stable unit array; From the stable distance unit array, the distance unit suitable for breathing detection is selected based on the phase standard deviation array and saved in the breathing index array; From the stable distance unit array, the distance units suitable for heart rate detection are screened out based on the phase standard deviation array and saved in the heart rate index array.

6. The radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics according to claim 1 is characterized in that: Perform multi-unit breathing detection on the screened breathing detection distance units, calculate the confidence of each breathing detection result, and cache the breathing detection results of multiple frames and multiple units; Performing statistical analysis on the cached breathing detection results and smoothing the results includes: A bandpass filter is designed according to the target respiratory frequency range, and then the bandpass filter is used to filter the data corresponding to each distance unit in the respiratory index array, and the filtered data of each distance unit is obtained and saved in the respiratory data matrix; The filtered data of each distance unit in the respiratory data matrix is ​​windowed, and then FFT is performed to obtain the amplitude spectrum of each distance unit and save it in the respiratory amplitude spectrum matrix; Calculate the noise array of the amplitude spectrum corresponding to each distance unit in the breathing amplitude spectrum matrix, and the noise arrays of all distance units constitute the breathing noise matrix; Searching for the peak position in the amplitude spectrum corresponding to each distance unit of the breathing amplitude spectrum matrix, calculating the signal-to-noise ratio of each peak based on the breathing noise matrix, and saving several peaks with the highest signal-to-noise ratio and their signal-to-noise ratios in the breathing detection result buffer; Perform sliding window histogram statistical analysis on the data in the breathing detection result buffer area to find the frequency area with the most concentrated number and use the mean of all breathing results in the area as the breathing statistical result; The statistical results are smoothed to obtain the final breathing detection results.

7. The radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics according to claim 1 is characterized in that: Perform multi-unit heart rate detection on the screened heart rate detection distance units, calculate the confidence of each heart rate detection result, extract the heart rate and its harmonic information, and cache the multi-frame and multi-unit detection results; Performing harmonic fusion statistical analysis on the cached heart rate detection results and smoothing the results includes: Design a bandpass filter according to the target heart rate range, and then use the bandpass filter to filter the data corresponding to each distance unit in the heart rate index array, obtain the filtered data of each distance unit and save it in the data matrix; The filtered data of each distance unit in the data matrix is ​​windowed and then FFT is performed to obtain the heart rate amplitude spectrum matrix; Calculate the noise array of the amplitude spectrum corresponding to each distance unit in the heart rate amplitude spectrum matrix, and the noise arrays of all distance units constitute the heart rate noise matrix; Searching for the peak position in the amplitude spectrum corresponding to each distance unit in the heart rate amplitude spectrum matrix, calculating the signal-to-noise ratio of each peak based on the heart rate noise matrix, and saving several peaks with the highest signal-to-noise ratio and their signal-to-noise ratios in the heart rate detection result buffer; Perform harmonic fusion statistical analysis on the heart rate detection results cached in the heart rate detection result buffer area, fold the harmonic frequency into the fundamental frequency interval, and find the frequency area with the most concentrated number through sliding window weighted histogram statistics, and use the mean of all heart rate results in the area as the heart rate statistical result; The heart rate statistical results are smoothed to obtain the final heart rate detection results.

8. An electronic device, characterized in that: It includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements a radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics as described in any one of claims 1 to 7.

9. A computer-readable storage medium, characterized in that: Computer executable instructions are stored, and the computer executable instructions are used to execute a radar respiratory heart rate detection method based on multi-unit detection and harmonic fusion statistics as described in any one of claims 1-7.

Images