CN114052709B - Robust millimeter wave radar vital sign measurement method - Google Patents
Robust millimeter wave radar vital sign measurement method Download PDFInfo
- Publication number
- CN114052709B CN114052709B CN202111460385.6A CN202111460385A CN114052709B CN 114052709 B CN114052709 B CN 114052709B CN 202111460385 A CN202111460385 A CN 202111460385A CN 114052709 B CN114052709 B CN 114052709B
- Authority
- CN
- China
- Prior art keywords
- signal
- respiratory
- heartbeat
- period
- signals
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 39
- 238000009528 vital sign measurement Methods 0.000 title claims abstract description 10
- 230000029058 respiratory gaseous exchange Effects 0.000 claims abstract description 52
- 238000001914 filtration Methods 0.000 claims abstract description 20
- 238000005259 measurement Methods 0.000 claims abstract description 10
- 210000000115 thoracic cavity Anatomy 0.000 claims abstract description 5
- 230000000241 respiratory effect Effects 0.000 claims description 73
- 238000005311 autocorrelation function Methods 0.000 claims description 33
- 238000012545 processing Methods 0.000 claims description 14
- 238000001228 spectrum Methods 0.000 claims description 9
- 230000001121 heart beat frequency Effects 0.000 claims description 8
- 230000036391 respiratory frequency Effects 0.000 claims description 8
- 230000001427 coherent effect Effects 0.000 claims description 6
- 230000036387 respiratory rate Effects 0.000 claims description 6
- 238000009825 accumulation Methods 0.000 claims description 5
- 238000009531 respiratory rate measurement Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 230000003068 static effect Effects 0.000 claims description 3
- 230000035565 breathing frequency Effects 0.000 claims description 2
- 238000004364 calculation method Methods 0.000 claims description 2
- 238000002592 echocardiography Methods 0.000 claims description 2
- 238000012544 monitoring process Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 7
- 230000008859 change Effects 0.000 description 4
- 238000001514 detection method Methods 0.000 description 4
- 238000010009 beating Methods 0.000 description 2
- 210000000038 chest Anatomy 0.000 description 2
- 235000011888 snacks Nutrition 0.000 description 2
- 208000035473 Communicable disease Diseases 0.000 description 1
- 230000002612 cardiopulmonary effect Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000009532 heart rate measurement Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000474 nursing effect Effects 0.000 description 1
- 238000013441 quality evaluation Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/08—Detecting, measuring or recording devices for evaluating the respiratory organs
- A61B5/0816—Measuring devices for examining respiratory frequency
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
- A61B5/024—Detecting, measuring or recording pulse rate or heart rate
- A61B5/0245—Detecting, measuring or recording pulse rate or heart rate by using sensing means generating electric signals, i.e. ECG signals
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7221—Determining signal validity, reliability or quality
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/7253—Details of waveform analysis characterised by using transforms
- A61B5/7257—Details of waveform analysis characterised by using transforms using Fourier transforms
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Pathology (AREA)
- Physiology (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Signal Processing (AREA)
- Biophysics (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Psychiatry (AREA)
- Pulmonology (AREA)
- Cardiology (AREA)
- Mathematical Physics (AREA)
- Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
Abstract
The invention provides a millimeter wave radar vital sign measurement method, which uses millimeter wave radar to measure thoracic skin fluctuation near the heart to obtain radar original echo signals, and comprises the following steps: the radar original echo signal is processed, and the initial position of the detected human body target is obtained; obtaining slow time dimension signals of an initial position and surrounding adjacent units, carrying out low-pass filtering and band-pass filtering on phase signals of each position of the adjacent units to obtain respiration signals and heartbeat signals of each position, evaluating the quality of the respiration signals and the heartbeat signals of N positions, taking the position with the best quality of the respiration signals and the heartbeat signals as a measurement point of the respiration heart rate, and outputting the respiration heart rate of the position.
Description
Technical Field
The invention belongs to the technical field of medical equipment, and particularly relates to a robust millimeter wave radar vital sign measurement method.
Background
Vital signs are the basis for human health assessment, wherein cardiopulmonary parameters are the core physiological indicators that measure the basic state of vital signs. Vital sign monitoring techniques can be broadly divided into two categories, contact and non-contact, in terms of the form of monitoring. For some special monitoring scenarios, such as: the non-contact physical sign monitoring mode has more advantages in nursing homes, patients with damaged skin, infectious diseases and the like. Common non-contact sign monitoring methods include infrared, video, acoustic, electromagnetic, and the like. The invention researches the non-contact vital sign detection of the focusing millimeter wave radar (electromagnetic wave).
Briefly, millimeter wave radar performs spatial environment sensing by transmitting millimeter wave signals to a detection space and receiving spatial echoes. If an object is present in the detection space, the object will generate a reflected signal and be received by the radar. The distance between the target and the radar causes a corresponding time delay between the radar transmit signal and the target echo signal. In practice, by measuring this time delay, distance information between the target and the radar can be deduced. When a human body target breathes normally, the chest cavity can generate corresponding fluctuation along with the breathing process, and the beating of the heart can also bring about the vibration near the outer skin of the heart. When the radar detects a human body target, fluctuation of the chest cavity and beating of the heart of the human body can cause weak change of the distance between the human body target and the radar. Such a weak change in distance causes a corresponding phase change in the received radar signal. By measuring the phase change, the respiratory state and the heartbeat state of the human body can be estimated.
Chinese patent publication CN112674738A discloses a method and apparatus for detecting respiratory heartbeat signals, and chinese patent publication CN112168152a discloses a method, apparatus and computer-readable storage medium for detecting respiratory heartbeat. Both patents determine the target position by directly selecting the highest energy position from radar range dimension signals as the echo signal of the human body target. Then, phase information is extracted from the target echo of the human body, fourier transformation is carried out, and then the highest frequency point is searched from the frequency spectrum distribution to be used as frequency estimation of respiratory and heartbeat signals.
When the radar is used for vital sign monitoring, the selection of the measuring points directly influences the estimated quality of the respiration and heartbeat signals. In the two patents, the selection mode of the monitoring point is that the echo signal amplitude of the target position of the human body is the largest, namely: the maximum amplitude of the echo signal is considered as the chest position of the human body target. However, in practice, the position where the radar echo amplitude is maximum is not necessarily the chest position due to the limitation of radar resolution, the difference between individuals, and the like, and thus, respiratory and heartbeat signal estimation is inaccurate or invalid.
In addition, in the above two patents, the estimation of the respiratory and heartbeat frequency of the human body is obtained by performing fourier transform on the corresponding signals. Due to the characteristic of Fourier transform, if the observation time is too short, the precision of the Fourier transform is poor, and the frequency estimation is inaccurate; if the observation time is too long, the measured human body needs to be kept still for a long time, and the measurement difficulty is increased.
Disclosure of Invention
In order to solve the technical problems, the invention provides a millimeter wave radar vital sign measurement method, which uses millimeter wave radar to measure thoracic skin fluctuation near the heart to obtain radar original echo signals, and comprises the following steps of
Step 1, processing an original radar echo signal to obtain an initial position of a detected human body target;
step 2, taking the initial position as a center, and acquiring slow time dimension signals of the initial position and surrounding adjacent units;
step 3, carrying out low-pass filtering on the phase signals of each position of the adjacent units to obtain breathing signals of each position, calculating an autocorrelation function of the breathing signals, and estimating the period of the breathing signals and the intensity of the breathing signals on the autocorrelation function; i.e. breathing cycle and breathing cycle intensity;
step 4, according to the period estimation result in the step 3, if the breathing period of the nth position is not zero, converting the breathing period into breathing frequency; calculating the signal-to-noise ratio of the respiratory signal according to the signal-to-noise ratio;
step 5, according to the calculation results of the step 3 and the step 4, if the respiration period of the nth position is not zero, evaluating the respiration signal quality of the nth position;
sequencing all positions with the estimated respiratory signal quality according to the respiratory signal quality QBREATn, taking the position with the best respiratory quality as a respiratory measurement point, and outputting the respiratory rate of the position as a final respiratory measurement value;
step 6, carrying out band-pass filtering on the phase signals PHn of each position of the adjacent units, further obtaining heartbeat signals, calculating an autocorrelation function of the heartbeat signals, and estimating the heartbeat signal period and the heartbeat signal intensity on the heartbeat autocorrelation function; namely, the heart cycle and the heart cycle intensity;
step 7, according to the period estimation result in the step 6, if the heartbeat period of the nth position is not equal to 0, converting the heartbeat period into the heartbeat frequency; calculating the signal-to-noise ratio of the heartbeat signal;
step 8, judging the estimation results of the step 6 and the step 7, and if the heartbeat period of the nth position is not equal to 0, evaluating the quality of the snack jump signal;
all the positions with the assessed heart beat signal quality are ranked according to the heart beat signal quality QHeartn, the position with the best heart beat signal quality is used as a heart rate measuring point, and the heart rate of the position is output as a final heart rate measuring value.
Further, the processing the radar original echo in step 1 includes: and carrying out wave beam formation, distance-azimuth dimension Fourier transformation, removing space static clutter, and carrying out coherent accumulation or non-coherent accumulation on the signals, thereby carrying out constant false alarm detection on the original echo signals.
Further, step 2 comprises the sub-steps of:
step 2.1, determining the adjacent units, wherein the adjacent units comprise initial positions and are N positions, and the slow time signal vector of the nth position is recorded as:
Sn＝[sn(1),sn(2),…sn(T*Fs)],
wherein, n=1, 2, … N, sn has a signal length of t×fs, where T is a slow time processing time, and Fs is a slow time sampling rate;
step 2.2, extracting phase information from each Sn according to the time sequence, and performing unwrapping operation on the phase signals, wherein the result is recorded as PHn;
PHn＝unwrap(phase(Sn))；
wherein phase (·) represents the phase signal; unwrap (·) represents unwrapping of the phase signal.
Further, step 3 includes the sub-steps of:
step 3.1, after low-pass filtering, the respiration signal is expressed as:
Breathn＝filter(PHn,LP)，
wherein filter (·) represents the filtering operation on the signal, LP is a low-pass filter parameter;
step 3.2, obtaining an autocorrelation function of the respiratory signal to obtain an autocorrelation curve acf_bn of the respiratory signal:
ACF_Bn＝autocorr(Breathn)，
wherein, autocorr (·) is the calculated autocorrelation function;
step 3.3, estimating the period of the respiration signal, i.e. the respiration signal period, on the autocorrelation function curve acf_bn, while the intensity at the peak point of the respiration period represents the intensity of the respiration period.
Further, step 4 includes the sub-steps of:
step 4.1, the reciprocal of the respiratory cycle is recorded as the respiratory frequency BreathFreqn;
step 4.2, calculating the spectrum BreathFFTn of the respiratory signal Breathn at the nth position,
BreathFFTn＝FFT(Breathn)，
wherein: FFT (& gt) is to perform Fourier transform on the signal;
step 4.3, according to the respiratory rate BreathFreqn, finding the amplitude of the corresponding frequency point on BreathFFTn, and taking the amplitude as the strength of the respiratory signal, and marking as BreathPSn:
BreathPSn＝BreathFFTn(BreathFreqn)，
in BreathFFTn, the noise is noted as BreathPNn,
BreathPNn＝SUM(BreathFFTn)–BreathPSn；
wherein SUM (·) is a SUM operation;
step 4.4, calculating the signal-to-noise ratio BreathSNRn of the respiratory signal at the nth position:
BreathSNRn＝10*log10(BreathPSn/BreathPNn)。
further, assessing the respiratory signal quality at the nth location in step 5 includes: the respiratory quality is expressed as:
QBreathn＝A*BreathStengthn+B*BreathSNRn,
wherein BreathStengthn represents the respiratory cycle intensity, breathSNRn represents the signal-to-noise ratio of the respiratory signal, A and B are weighting coefficients of the respiratory cycle intensity and the respiratory signal-to-noise ratio,
further, step 6 includes the sub-steps of:
step 6.1, after filtering by a band-pass filter, the heartbeat signal is expressed as:
Heartn＝filter(PHn,BP)，
wherein, filter (·) represents filtering operation on the signal, BP is a band-pass filter parameter;
step 6.2, obtaining an autocorrelation function of the heartbeat signal to obtain an autocorrelation curve ACF_Hn of the heartbeat signal:
ACF_Hn＝autocorr(Heartn)；
and 6.3, taking the peak point interval of the ACF_Hn curve as the period of the heartbeat signal, and simultaneously, representing the intensity of the period of the heartbeat signal by the intensity at the peak point.
Further, the low-pass filter uses a low-pass filter with a cutoff frequency of 0.8 Hz; the band pass filter uses a band pass filter with a passband frequency bp=0.8-2.5 Hz.
Further, step 7 includes the sub-steps of:
step 7.1, recording the reciprocal of the heartbeat period as heart frequency HeartFreqn;
step 7.2, calculating a spectrum HeartFFTn of the heartbeat signal Heartn at the nth position,
HeartFFTn＝FFT(Heartn)，
wherein, FFT (·) is to make Fourier transform for the signal;
step 7.3, according to the heart beat frequency, finding the amplitude of the corresponding frequency point on the HeartFFTn, and taking the amplitude as the intensity of the heart beat signal, and recording as HeartPSn:
HeartPSn＝HeartFFTn(HeartFreqn)；
in the HeartFFTn, the noise is noted as HeartPNn,
HeartPNn＝SUM(HeartFFTn)–HeartPSn；
wherein SUM (·) is a SUM operation;
step 7.4, calculating the signal-to-noise ratio HeartSNRn of the respiratory signal at the nth position;
HeartSNRn＝10*log10(HeartPSn/HeartPNn)。
further, step 8 includes: evaluating an nth location heartbeat signal quality, the heartbeat signal quality expressed as:
QHeartn＝A*HeartStengthn+B*HeartSNRn，
wherein HeartStengthn represents the heart cycle intensity, heartSNRn represents the signal-to-noise ratio of the heart signal, and A and B are weighting coefficients of the heart cycle intensity and the signal-to-noise ratio of the heart signal.
By adopting the method, the Fourier transform result is not used as the sole basis of respiratory/heartbeat frequency estimation, so that the respiratory/heartbeat frequency estimation accuracy and robustness are improved; the respiratory/heartbeat signal quality of each spatial position is comprehensively evaluated, so that the problem of inaccurate or wrong respiratory/heart rate measurement caused by the selection deviation of the measurement position is solved.
Drawings
FIG. 1 is a flow chart of the present invention;
fig. 2 is a graph of a typical respiratory signal and its corresponding autocorrelation function and spectrum.
Detailed Description
The invention provides a millimeter wave radar vital sign measurement method for signal time-frequency domain combined processing, which can be used for stably estimating respiratory frequency and heart rate in signals and overcoming inherent problems of a Fourier transform frequency estimation method; and secondly, the signal quality of each spatial position can be evaluated, and the problem of inaccurate or wrong measurement of breathing and heartbeat caused by the selection deviation of the measurement position is solved.
The following detailed description of specific embodiments of the invention refers to the accompanying drawings.
The main functions of the invention can be summarized as follows: on the basis of the known target position P, the quality of vital sign signals of the target position P and the positions nearby the target position P is evaluated, and the optimal vital sign measurement position P is selected from the target position P, so that the robustness and the accuracy of vital sign measurement are improved. The process can be divided into 8 steps, and the flow is shown in figure 1.
Step S101 is implemented by the following method: in the present invention, it is assumed that basic radar signal processing has been completed, and an initial position of a measured human body target (hereinafter referred to as target position P) has been obtained. Basic radar signal processing, including but not limited to: the method comprises the steps of carrying out wave beam formation (or distance-azimuth dimension Fourier transform) on an original echo signal, removing space static clutter, carrying out coherent (or non-coherent) accumulation on the signal, detecting constant false alarm and the like.
Step S102 is implemented by the following method: after the initial position P of the target is obtained, a slow time dimension signal (complex number) of the target position P and neighboring units (including P, total N positions) around the target position P is obtained with the target position P as a center, and the slow time signal of the nth position is denoted as sn= [ Sn (1), sn (2), … Sn (t×fs) ], and n=1, 2, … N. The signal length of Sn is t×fs, where T is the slow time processing time and Fs is the slow time sampling rate.
The phase information is extracted for each Sn in time sequence and the phase signal is subjected to an unwrapping operation, the result being denoted PHn.
PHn＝unwrap(phase(Sn))；
Wherein phase (·) represents the phase signal; unwrap (·) represents unwrapping of the phase signal.
Step S103 is implemented by the following method: the phase signal PHn for each position is low pass filtered to obtain a respiratory signal, then an autocorrelation function of the respiratory signal is calculated, and respiratory signal period and intensity are estimated on the autocorrelation function.
After low pass filtering, the respiratory signal can be expressed as:
Breathn＝filter(PHn,LP)；
where filter (·) represents the filtering operation on the signal, LP is a low-pass filter parameter with cut-off frequency=0.8 Hz.
Then, the respiratory signal is subjected to autocorrelation function to obtain a respiratory signal autocorrelation curve, which is marked as ACF_Bn:
ACF_Bn＝autocorr(Breathn)。
wherein, autocorr (·) is the calculated autocorrelation function.
The respiratory signal period and intensity process is estimated on the autocorrelation function curve acf_bn as: the period of normal respiration of human body is 2 s-10 s (corresponding to frequency 0.1 Hz-0.5 Hz), so on the autocorrelation function curve of normal respiration, at least one peak point (named PosBreath) should exist in the range of 2 s-10 s. In addition, there should also be a maximum negative correlation point before the peak point of ACF_Bn (at about PosBreath/2). In summary, if the position of the first peak point of the-acf_bn curve is exactly 1/2 of the position of the peak point of the acf_bn curve, the respiratory signal Breathn is considered to have a respiratory cycle, and the peak point (PosBreathn) of the acf_bn curve is the respiratory signal cycle, and the intensity at the peak point represents the intensity of the respiratory cycle.
[ValBreathn,PosBrethn]＝findpeak(ACF_Bn)
[ValBreathn_neg,PosBreathn_neg]＝findpeak(-ACF_Bn)
If 1.9<(PosBreathn/PosBreathn_neg)<2.1
BreathCirclen＝PosBreathn
BreathStengthn＝ValBreathn
Else
BreathCirclen＝0
BreathStengthn＝0
End
findpeak (·) represents the search peak position, breathCirclen and breathstenthn represent the respiratory signal period and respiratory signal period intensity at the nth position, n=1, 2,3 … N, respectively.
Step S104 is implemented by the following method: from the result in step S103, if the respiratory cycle breathcircle of the nth position is not equal to 0, it is indicated that there is a respiratory signal at the position, while the respiratory cycle breathcircle is converted into the respiratory frequency BreathFreqn.
BreathFreqn＝1/BreathCirclen。
The spectrum of the respiratory signal Breathn at the nth position is calculated, denoted BreathFFTn, breathFFTn =fft (Breathn), FFT (·) being the fourier transform of the signal.
According to the respiratory frequency BreathFreqn, the amplitude of the corresponding frequency point is found on BreathFFTn, and the amplitude is taken as the strength of the respiratory signal and is recorded as BreathPSn. Breathpsn= BreathFFTn (BreathFreqn). Meanwhile, in BreathFFTn, signals other than BreathPSn are regarded as noise, and are denoted as BreathPNn.
BreathPNn＝SUM(BreathFFTn)–BreathPSn；
Wherein SUM (·) is the SUM operation.
Next, the signal-to-noise ratio BreathSNRn of the respiratory signal at the nth position is calculated.
BreathSNRn＝10*log10(BreathPSn/BreathPNn)。
Step S105 is implemented by the following method: according to the results of steps S103 and S104, if the respiratory cycle breathcircle of the nth position is not equal to 0, the respiratory signal quality of the point is evaluated, and the respiratory quality qbreathn=a×breathstenthn+b×breathsnrn. Wherein BreathStengthn, breathSNRn represents the respiration cycle intensity and the signal-to-noise ratio of the respiration signal, and A and B are weighting coefficients of the respiration cycle intensity and the signal-to-noise ratio of the respiration signal, which can be empirically set.
Finally, the breathing mass QBreathn is ordered for N positions, n=1, 2, … N. The position with the best breathing quality is taken as a measuring point of breathing, and the breathing rate of the position is output as a final breathing measuring value.
To facilitate an understanding of the above process, a typical respiration signal, and its corresponding autocorrelation curves and spectra are shown in the left and right of fig. 2, respectively.
Step S106 is implemented by the following method: the processing of the heartbeat signal is similar to the processing of the respiration signal. And carrying out band-pass filtering on the phase signal PHn of each position to obtain a heartbeat signal, then calculating an autocorrelation function of the heartbeat signal, and estimating the period and the intensity of the heartbeat signal on the heartbeat autocorrelation function.
After filtering, the heartbeat signal may be expressed as:
Heartn＝filter(PHn,BP)；
wherein, filter (·) represents a filtering operation on the signal, and BP is a band-pass filter parameter with passband frequency=0.8-2.5 Hz.
Then, the heart beat signal is subjected to autocorrelation function to obtain a heart beat signal autocorrelation curve which is marked as ACF_Hn:
ACF_Hn＝autocorr(Heartn)。
the heart beat signal period and intensity process is estimated on the autocorrelation function curve ACF_Hn as follows:
since the normal heart beat of human body ranges from 0.4s to 1.25s (0.8 Hz to 2.5 Hz), at least one peak point (marked as PosHeartn) should exist in the range from 0.4s to 1.25s on a normal heart beat autocorrelation function curve. In addition, there should also be a maximum negative correlation point before this peak point (at about PosHeartn/2). Therefore, if the position of the first peak point of the acf_hn curve is exactly 1/2 of the position of the peak point of the acf_hn curve, the heartbeat signal is considered to have a period, and the peak point (posheart) of the ACFHn curve is the period of the heartbeat signal. Meanwhile, the intensity at the peak point (posheart) represents the intensity of the heart beat signal period.
[ValHeartn,PosHeartn]＝findpeak(ACF_Hn)
[ValHeartn_neg,PosHeartn_neg]＝findpeak(-ACF_Hn)
If 1.9<(PosHeartn/PosHeartn_neg)<2.1
HeartCirclen＝PosHeartn
HeartStengthn＝ValHeartn
Else
HeartCirclen＝0
HeartStengthn＝0
End
Heart cycle and heart cycle intensities at the nth position, n=1, 2,3 … N, respectively.
Step S107 is implemented by the following method: according to the result in step S106, if the heartbeat cycle heart cycle at the nth position is not equal to 0, it is indicated that there is a heartbeat signal at the position, and the heartbeat cycle heart cycle is converted into the heartbeat frequency heart freqn.
HeartFreqn＝1/HeartCirclen。
The spectrum of the signal heart at the nth position is calculated, denoted HeartFFTn, heartFFTn =fft (heart), FFT (·) being the fourier transform of the signal.
According to heart rate HeartFreqn, finding the amplitude of the corresponding frequency point on HeartFFTn, and taking the amplitude as the power of heart rate signals, and recording as HeartPSn.
HeartPSn＝HeartFFTn(HeartFreqn)。
Meanwhile, in the HeartFFTn, signals other than HeartPSn are regarded as noise, denoted as HeartPNn.
HeartPNn＝SUM(HeartFFTn)–HeartPSn；
Wherein SUM (·) is the SUM operation.
Next, the signal-to-noise ratio heart snrn of the respiratory signal at the nth position is calculated.
HeartSNRn＝10*log10(HeartPSn/HeartPNn)。
Step S108 is implemented by the following method: according to the results of step S106 and step S107, if the heartbeat cycle heart cycle at the nth position is not equal to 0, the snack hop signal quality is evaluated, the heartbeat signal quality:
QHeartn＝A*HeartStengthn+B*HeartSNRn。
wherein HeartStengthn, heartSNRn represents the heart cycle intensity and the signal-to-noise ratio of the heart signal, and A and B are weighting coefficients of the heart cycle intensity and the signal-to-noise ratio of the heart signal, which can be set empirically.
Finally, the heart beat signal quality qheart N of the N positions is ordered, n=1, 2, … N, the position with the best heart beat signal quality is taken as the heart rate measuring point, and the heart rate of the position is output as the final heart rate measuring value.
Compared with the prior art, the invention firstly does not use the Fourier transformation result as the only basis for respiratory (heartbeat) frequency estimation, and improves the estimation accuracy and robustness of respiratory (heartbeat) rate;
secondly, the respiratory (heartbeat) signal quality of each spatial position is comprehensively evaluated, so that the problem of inaccurate or wrong measurement of the respiratory (heartbeat) rate caused by the selection deviation of the measurement position is solved.
Taking respiratory rate estimation as an example (step S103-step S105), in the process of estimating respiratory quality near a target, the invention provides a measurement method of joint processing (step S105) of a signal time domain (step S103) and a frequency domain (step S104) (the heart rate estimation process (step S106-step S108) is similar to respiratory rate estimation). The method firstly does not take the Fourier transform result as the only basis of respiratory (heart beat) frequency estimation, and improves the estimation accuracy and robustness of respiratory (heart) rate; secondly, the respiratory (heartbeat) signal quality of each spatial position is comprehensively evaluated, so that the problem of inaccurate or wrong measurement of the respiratory (heartbeat) rate caused by the selection deviation of the measurement position is solved.
1. The method for judging the signal period by the autocorrelation function comprises the following steps: if the position of the first peak point of the ACF curve is exactly 1/2 of the position of the peak point of the ACF curve, the signal is considered to have a period, and the peak point of the ACF curve is the period of the signal.
2. The mode of carrying out time-frequency domain joint processing on the vital sign signals is as follows: and (5) carrying out signal quality evaluation by combining the signal period intensity in the signal autocorrelation function with the signal to noise ratio.
3. In the form of a fine search of vital sign measurement locations in the vicinity of the target.
Finally, it should be noted that the above-mentioned embodiments are merely for illustrating the technical solution of the embodiment of the present invention, and not for limiting, and although the embodiment of the present invention has been described in detail with reference to the above-mentioned preferred embodiments, it should be understood by those skilled in the art that modifications and equivalent substitutions can be made to the technical solution of the embodiment of the present invention without departing from the spirit and scope of the technical solution of the embodiment of the present invention.
Claims (7)
1. A millimeter wave radar vital sign measurement method is characterized in that the method uses a millimeter wave radar to measure thoracic skin fluctuation near the heart to obtain a radar original echo signal, and the method comprises the following steps of
Step 1, processing an original radar echo signal to obtain an initial position of a detected human body target;
step 2, taking the initial position as the center, acquiring slow time dimension signals of the initial position and surrounding adjacent units, and specifically comprising the following substeps:
step 2.1, determining the adjacent units, wherein the adjacent units comprise initial positions and are N positions, and the slow time signal vector of the nth position is recorded as:
Sn＝[sn(1),sn(2),…sn(T*Fs)],
wherein, n=1, 2, … N, sn has a signal length of t×fs, where T is a slow time processing time, and Fs is a slow time sampling rate;
step 2.2, extracting phase information from each Sn according to the time sequence, and performing unwrapping operation on the phase signals, wherein the result is recorded as PHn;
PHn＝unwrap(phase(Sn))；
wherein phase (·) represents the phase signal; unwrap (·) represents unwrapping the phase signal;
step 3, carrying out low-pass filtering on the phase signals of each position of the adjacent units to obtain breathing signals of each position, calculating an autocorrelation function of the breathing signals, and estimating the period of the breathing signals and the intensity of the breathing signals on the autocorrelation function; i.e. breathing cycle and breathing cycle intensity;
step 4, according to the period estimation result in the step 3, if the breathing period of the nth position is not zero, converting the breathing period into breathing frequency; calculating the signal-to-noise ratio of the respiratory signal according to the signal-to-noise ratio;
step 5, according to the calculation results of the step 3 and the step 4, if the respiration period of the nth position is not zero, evaluating the respiration signal quality of the nth position;
sequencing all positions with the estimated respiratory signal quality according to the respiratory signal quality QBREATn, taking the position with the best respiratory quality as a measurement point of respiration, and outputting the respiratory rate of the position as a final respiratory measurement value, wherein the respiratory signal quality is expressed as:
QBreathn＝A*BreathStengthn+B*BreathSNRn,
wherein, breathStengthn represents the respiratory cycle intensity, breathSNRn represents the signal-to-noise ratio of the respiratory signal, and a and B are the weighted coefficients of the respiratory cycle intensity and the respiratory signal-to-noise ratio;
step 6, carrying out band-pass filtering on the phase signals PHn of each position of the adjacent units, further obtaining heartbeat signals, calculating an autocorrelation function of the heartbeat signals, and estimating the heartbeat signal period and the heartbeat signal intensity on the heartbeat autocorrelation function; namely, the heart cycle and the heart cycle intensity;
step 7, according to the period estimation result in the step 6, if the heartbeat period of the nth position is not equal to 0, converting the heartbeat period into the heartbeat frequency; calculating the signal-to-noise ratio of the heartbeat signal;
step 8, judging the estimation results of the step 6 and the step 7, and if the heartbeat period of the nth position is not equal to 0, evaluating the heartbeat signal quality of the position;
sequencing all positions with the assessed heart beat signal quality according to the heart beat signal quality QHeartn, taking the position with the best heart beat signal quality as a heart rate measuring point, and outputting the heart rate of the position as a final heart rate measuring value;
the heartbeat signal quality qheart is expressed as:
QHeartn＝A1*HeartStengthn+B1*HeartSNRn，
wherein HeartStengthn represents the heart cycle intensity, heartSNRn represents the signal-to-noise ratio of the heart beat signal, and A1 and B1 are the weighted coefficients of the heart beat signal cycle intensity and the heart beat signal-to-noise ratio.
2. The method of claim 1, wherein the processing of the radar raw echoes of step 1 comprises: and carrying out wave beam formation, distance-azimuth dimension Fourier transformation, removing space static clutter, carrying out coherent accumulation or non-coherent accumulation on the signals, and detecting constant false alarm.
3. The method of claim 1, wherein step 3 comprises the sub-steps of:
step 3.1, after low-pass filtering, the respiration signal is expressed as:
Breathn＝filter(PHn,LP)，
wherein filter (·) represents the filtering operation on the signal, LP is a low-pass filter parameter;
step 3.2, obtaining an autocorrelation function of the respiratory signal to obtain an autocorrelation curve acf_bn of the respiratory signal:
ACF_Bn＝autocorr(Breathn)，
wherein, autocorr (·) is the calculated autocorrelation function;
step 3.3, estimating the period of the respiration signal, i.e. the respiration signal period, on the autocorrelation function curve acf_bn, while the intensity at the peak point of the respiration period represents the intensity of the respiration period.
4. The method of claim 1, wherein step 4 comprises the sub-steps of:
step 4.1, the reciprocal of the respiratory cycle is recorded as the respiratory frequency BreathFreqn;
step 4.2, calculating the spectrum BreathFFTn of the respiratory signal Breathn at the nth position,
BreathFFTn＝FFT(Breathn)，
wherein: FFT (& gt) is to perform Fourier transform on the signal;
step 4.3, according to the respiratory rate BreathFreqn, finding the amplitude of the corresponding frequency point on BreathFFTn, and taking the amplitude as the strength of the respiratory signal, and marking as BreathPSn:
BreathPSn＝BreathFFTn(BreathFreqn)，
in BreathFFTn, the noise is noted as BreathPNn,
BreathPNn＝SUM(BreathFFTn)–BreathPSn；
wherein SUM (·) is a SUM operation;
step 4.4, calculating the signal-to-noise ratio BreathSNRn of the respiratory signal at the nth position:
BreathSNRn＝10*log10(BreathPSn/BreathPNn)。
5. the method of claim 1, wherein step 6 comprises the sub-steps of:
step 6.1, after filtering by a band-pass filter, the heartbeat signal is expressed as:
Heartn＝filter(PHn,BP)，
wherein, filter (·) represents filtering operation on the signal, BP is a band-pass filter parameter;
step 6.2, obtaining an autocorrelation function of the heartbeat signal to obtain an autocorrelation curve ACF_Hn of the heartbeat signal:
ACF_Hn＝autocorr(Heartn)；
and 6.3, taking the peak point interval of the ACF_Hn curve as the period of the heartbeat signal, and simultaneously, representing the intensity of the period of the heartbeat signal by the intensity at the peak point.
6. The method of claim 5, wherein the low pass filtering uses a low pass filter with a cutoff frequency of 0.8 Hz; the band pass filter uses a band pass filter with a passband frequency bp=0.8-2.5 Hz.
7. The method of claim 1, wherein step 7 comprises the sub-steps of:
step 7.1, recording the reciprocal of the heartbeat period as heart frequency HeartFreqn;
step 7.2, calculating a spectrum HeartFFTn of the heartbeat signal Heartn at the nth position,
HeartFFTn＝FFT(Heartn)，
wherein, FFT (·) is to make Fourier transform for the signal;
step 7.3, according to the heart beat frequency, finding the amplitude of the corresponding frequency point on the HeartFFTn, and taking the amplitude as the intensity of the heart beat signal, and recording as HeartPSn:
HeartPSn＝HeartFFTn(HeartFreqn)；
in the HeartFFTn, the noise is noted as HeartPNn,
HeartPNn＝SUM(HeartFFTn)–HeartPSn；
wherein SUM (·) is a SUM operation;
step 7.4, calculating the signal-to-noise ratio HeartSNRn of the respiratory signal at the nth position;
HeartSNRn＝10*log10(HeartPSn/HeartPNn)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111460385.6A CN114052709B (en) | 2021-12-02 | 2021-12-02 | Robust millimeter wave radar vital sign measurement method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111460385.6A CN114052709B (en) | 2021-12-02 | 2021-12-02 | Robust millimeter wave radar vital sign measurement method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114052709A CN114052709A (en) | 2022-02-18 |
| CN114052709B true CN114052709B (en) | 2024-03-15 |
Family
ID=80228320
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111460385.6A Active CN114052709B (en) | 2021-12-02 | 2021-12-02 | Robust millimeter wave radar vital sign measurement method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114052709B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114527463B (en) * | 2022-04-24 | 2022-08-12 | 南京楚航科技有限公司 | In-vehicle living body detection method and device by utilizing phase matching |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110115592A (en) * | 2018-02-07 | 2019-08-13 | 英飞凌科技股份有限公司 | The system and method for the participation level of people are determined using millimetre-wave radar sensor |
| CN110192850A (en) * | 2019-05-31 | 2019-09-03 | 湖南省顺鸿智能科技有限公司 | Extracting method and system based on heartbeat signal under radar return strong noise background |
| CN112401856A (en) * | 2020-11-16 | 2021-02-26 | 武汉烽火凯卓科技有限公司 | Nursing home monitoring method and system based on millimeter wave radar |
| CN112965060A (en) * | 2021-02-19 | 2021-06-15 | 加特兰微电子科技(上海)有限公司 | Detection method and device for vital sign parameters and method for detecting physical sign points |
| CN113534141A (en) * | 2021-07-01 | 2021-10-22 | 深圳晶华相控科技有限公司 | Remote vital sign detection method and device based on phased array radar technology |
-
2021
- 2021-12-02 CN CN202111460385.6A patent/CN114052709B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110115592A (en) * | 2018-02-07 | 2019-08-13 | 英飞凌科技股份有限公司 | The system and method for the participation level of people are determined using millimetre-wave radar sensor |
| CN110192850A (en) * | 2019-05-31 | 2019-09-03 | 湖南省顺鸿智能科技有限公司 | Extracting method and system based on heartbeat signal under radar return strong noise background |
| CN112401856A (en) * | 2020-11-16 | 2021-02-26 | 武汉烽火凯卓科技有限公司 | Nursing home monitoring method and system based on millimeter wave radar |
| CN112965060A (en) * | 2021-02-19 | 2021-06-15 | 加特兰微电子科技(上海)有限公司 | Detection method and device for vital sign parameters and method for detecting physical sign points |
| CN113534141A (en) * | 2021-07-01 | 2021-10-22 | 深圳晶华相控科技有限公司 | Remote vital sign detection method and device based on phased array radar technology |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114052709A (en) | 2022-02-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109965858B (en) | Ultra-wideband radar-based human body vital sign detection method and device | |
| EP3701873B1 (en) | System and method for indicating risk of coronary artery disease | |
| Jezewski et al. | A novel technique for fetal heart rate estimation from Doppler ultrasound signal | |
| KR102025328B1 (en) | Apparatus and method for generating ultrasonic vector doppler image using plane wave synthesis | |
| JP6574524B2 (en) | Imaging system and method for determining translational speed of a catheter | |
| CN111856455A (en) | Multi-target heart rate and respiration measuring method and system matched with different radar bandwidths | |
| US20150238169A1 (en) | Ultrasonic measurement apparatus and ultrasonic measurement method | |
| JP2001286474A (en) | Dynamic measurement of subject's parameter | |
| KR101827522B1 (en) | Apparatus and method for measuring heart rate through non-contact | |
| CN114052709B (en) | Robust millimeter wave radar vital sign measurement method | |
| CN107106125B (en) | System and method for measuring arterial parameters | |
| CN113786176B (en) | Accurate millimeter wave radar breath and heartbeat measurement method, system and storage medium | |
| CN113273978B (en) | Ultra-wideband radar-based human body respiration and heartbeat frequency detection method | |
| CN111685760B (en) | Human body respiratory frequency calculation method based on radar measurement | |
| JP2022185752A (en) | Biological information detection system, program, and biological information detection method | |
| CN109452954A (en) | Ultrasonic imaging method and device | |
| KR20130095160A (en) | Ultrasound apparatus and method for generating ultrasound image | |
| Lu et al. | Accurate heart beat detection with doppler radar using bidirectional GRU network | |
| CN115040091A (en) | VMD algorithm-based millimeter wave radar life signal extraction method | |
| Lee et al. | Towards higher accuracy and better noise-tolerance for fetal heart rate monitoring using Doppler ultrasound | |
| CN114931371A (en) | Method and device for measuring lung function index, electronic device and storage medium | |
| JP2007090003A (en) | Ultrasonic diagnostic apparatus and controlling method thereof | |
| JP2007268023A (en) | Method and apparatus for observing embolus and ultrasonic diagnostic equipment using the same | |
| JP4357260B2 (en) | Acceleration pulse wave measuring device | |
| CN112043256A (en) | Radar-based multi-target heart rate real-time measurement method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |