CN113951856A - 一种基于多普勒雷达心跳检测的频谱估计方法 - Google Patents
一种基于多普勒雷达心跳检测的频谱估计方法 Download PDFInfo
- Publication number
- CN113951856A CN113951856A CN202111190185.3A CN202111190185A CN113951856A CN 113951856 A CN113951856 A CN 113951856A CN 202111190185 A CN202111190185 A CN 202111190185A CN 113951856 A CN113951856 A CN 113951856A
- Authority
- CN
- China
- Prior art keywords
- signal
- heartbeat
- frequency
- chest wall
- algorithm
- 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.)
- Pending
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 44
- 238000001228 spectrum Methods 0.000 title claims abstract description 42
- 238000000034 method Methods 0.000 title claims description 44
- 210000000779 thoracic wall Anatomy 0.000 claims abstract description 59
- 238000006073 displacement reaction Methods 0.000 claims abstract description 54
- 238000004364 calculation method Methods 0.000 claims abstract description 11
- 238000001914 filtration Methods 0.000 claims abstract description 6
- 230000029058 respiratory gaseous exchange Effects 0.000 claims description 26
- 238000004458 analytical method Methods 0.000 claims description 23
- 238000005516 engineering process Methods 0.000 claims description 17
- 230000006870 function Effects 0.000 claims description 13
- 238000007781 pre-processing Methods 0.000 claims description 11
- 238000012545 processing Methods 0.000 claims description 9
- 238000000926 separation method Methods 0.000 claims description 9
- 230000009466 transformation Effects 0.000 claims description 8
- 238000010183 spectrum analysis Methods 0.000 claims description 7
- 230000000241 respiratory effect Effects 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- 230000008859 change Effects 0.000 claims description 4
- 230000008569 process Effects 0.000 claims description 4
- 230000004069 differentiation Effects 0.000 claims description 3
- 230000010363 phase shift Effects 0.000 claims description 3
- 238000002360 preparation method Methods 0.000 claims description 3
- 230000009467 reduction Effects 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims description 3
- 230000001629 suppression Effects 0.000 claims description 3
- 230000003111 delayed effect Effects 0.000 claims description 2
- 230000001121 heart beat frequency Effects 0.000 abstract description 6
- 230000033001 locomotion Effects 0.000 abstract description 6
- 238000004804 winding Methods 0.000 abstract description 3
- 238000010009 beating Methods 0.000 abstract description 2
- 210000000038 chest Anatomy 0.000 abstract description 2
- 230000036541 health Effects 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 238000000354 decomposition reaction Methods 0.000 description 4
- 238000003745 diagnosis Methods 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 4
- 238000000605 extraction Methods 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- 230000000737 periodic effect Effects 0.000 description 4
- 230000003595 spectral effect Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 238000009825 accumulation Methods 0.000 description 2
- 230000002612 cardiopulmonary effect Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 208000019622 heart disease Diseases 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000011084 recovery Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 238000010989 Bland-Altman Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 238000013186 photoplethysmography Methods 0.000 description 1
- 238000011897 real-time detection Methods 0.000 description 1
- 230000036391 respiratory frequency Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000033764 rhythmic process Effects 0.000 description 1
- 238000010187 selection method Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 210000000115 thoracic cavity Anatomy 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
- 210000000707 wrist Anatomy 0.000 description 1
Images
Classifications
-
- 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
-
- 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/02444—Details of sensor
-
- 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/726—Details of waveform analysis characterised by using transforms using Wavelet transforms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Cardiology (AREA)
- Public Health (AREA)
- Animal Behavior & Ethology (AREA)
- Remote Sensing (AREA)
- Physiology (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Radar, Positioning & Navigation (AREA)
- General Health & Medical Sciences (AREA)
- Veterinary Medicine (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Electromagnetism (AREA)
- Signal Processing (AREA)
- Psychiatry (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Artificial Intelligence (AREA)
- Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)
Abstract
本发明提出了一种基于多普勒雷达的非接触式精确心跳检测方法。首先，采用连续波多普勒雷达体征检测的原理捕捉人体心脏跳动引起的微弱胸腔位移引起的频移，即雷达向受试者通过发射端发射T(t)。随后T(t)信号的相位被受试者胸壁运动x(t)所调制，反射信号被雷达接收端所捕获；其次，反射信号经过前端的低噪声放大器后，得到基带信号，采用正交混频器解调基带信号，再经过增益放大器后，由反正切解调和接卷绕计算可得到胸壁信号；最后，胸壁模拟信号经过模数转换器得到数字信号x(t)，在数字信号处理后在PC端经过马勒特算法和匹配滤波算法分离提取心脏信号，再经过FTPR‑TWV算法后得到心跳频谱结果。
Description
技术领域
本发明涉及心率检测领域，特别是涉及一种在频域基于多普勒雷达的非接触式精确心跳检测方法。
背景技术
随着生活质量的提高，人们对生命健康越来越重视，而心率检测对于日常健康监测和临床医疗诊断有着重大作用。一方面，心跳的频率能够反映基本的健康状况；另一方面，心跳的节律和强度在很大程度上有利于心脏病的诊断。与传统的检测心跳的接触方式比如心电图、光电容积脉搏波、腕带式脉搏血氧仪相比，基于多普勒雷达的非接触式心跳检测在一些不能接触电极的独特场合具有优越性，例如抗震救灾时的生命检测、婴儿监护和疲劳检测，这些接触式心跳检测都非常不方便甚至不切实际。
调频连续波雷达作为一种雷达传感器，近几十年来被广泛应用于远程监测呼吸、心跳等生命体征。目前为止，由于具有I/Q通道解调的正交接收器架构的广泛应用，许多正交解调方法相继出现。小角度近似仅通过预设目标位置就可以很好地执行。复信号解调引入谐波干扰，严重降低心率检测精度。此外，反正切解调对目标位置和谐波干扰不敏感，但仍然具有很强的非线性和不连续性。扩展的微分和交叉乘法算法克服了相位缠绕和谐波干扰的缺点，其存在的问题是需要解决由于离散积分导致的误差累积问题。
雷达心跳检测的挑战之一是在呼吸及其谐波强干扰下准确估计心率。这是因为心跳引起的胸壁位移(0.1～0.5mm)远小于呼吸引起的胸壁位移(1～12mm)。连续小波变换、集合经验模态分解、奇异谱分析和变分模态分解等一系列代表性方法已用于心肺信号分离和提取以消除呼吸干扰。虽然这些方法对于呼吸和心跳的分离和恢复已经证明是有效的，但仍然存在一些不足。上述方法对信噪比非常敏感，而且其计算成本大大增加。最近，基于匹配滤波器可以最大化输出信噪比的特性来恢复心跳信号的能力和鲁棒性已经在存在大规模随机身体运动的情况下得到证实。此外，提出的多分辨率分析已广泛应用于轴承故障诊断、雷达系统和图像处理。相比之下，使用多普勒雷达进行生命体征研究的方面还有所欠缺。
另一个挑战是短期时间窗口中频谱分辨率的不足。使用超过10秒的长周期时间窗口可以实现足够的频谱分辨率。但是时间窗口越长，心跳细节的损失就会越多。因此，对于频域分析中的HRV检测，需要在短周期时间窗口内处理雷达数据以实现快速HR采集。
发明内容
基于以上两个技术难题，本发明提出一种基于多普勒雷达心跳准确检测的频谱估计方法。
一种基于多普勒雷达心跳准确检测的频谱估计方法，其步骤如下：
步骤一：采用多普勒雷达发射雷达信号获取回波信号，并用信号解调技术进行正交信号解调；
步骤二：通过前端的低通滤波器进行杂波抑制和胸壁位移微分得到预处理后的胸壁位移信号；
步骤三：胸壁信号经过模数转换器得到数字信号；
步骤四：利用多频率分析算法(MRA)实现呼吸信号与心跳信号的分离；
步骤五：利用模板匹配算法(TMF)恢复隐藏在呼吸信号中的心跳信号；
步骤六：设计并采用频率时间相位回归技术(FTPR)与时间窗口变化技术(TWV)相结合的算法估计心跳频谱；
所述的获取回波信号的方法是：通过芯片上的雷达模块发射端发射频率为77GHz的毫米波T(t)到目标特征体处，雷达信号可表示为
由频域中相移的变化来捕捉由呼吸和心跳引起的胸壁位移信号，反射信号R(t)被雷达接收端所捕获，反射信号可表示为：
其中c为无线电的波速，λ＝c/f表示雷达波长；
把R(t)近似看作经过2d0/c时间延迟并相位调制后的T(t)，反射信号R(t)再经过前端的低噪声放大器后，得到基带信号B(t)，
其中，
为常数相位偏移，与接收机本身参数以及d0相关；残余相位噪声简化为
然后，由正交混频器(I/Q Mixer)解调基带信号B(t)，同相和正交信号分别表示为：
再经过增益放大器(Gain)后，由呼吸和心跳引起的胸壁信号x(t)由反正切解调和解卷绕技术计算可得，表达式如下：
最后，模拟信号x(t)经过模数转换器(ADC)得到数学信号x(t)，用于数字信号处理(DSP)以及PC端心跳检测算法的开发；
对得到的雷达所捕获的胸壁位移信号进行预处理，预处理步骤包括去除背景噪声和多项式趋势、带通滤波，
多普勒雷达的前端采用了正交接收机结构，首先，通过正交信号解调技术获得相位，然后，根据与相位成正比的关系计算胸壁位移，将反正弦解调算法引入到正交信号的解调中，解调后的胸壁位移信号表示为：
其中N是雷达原始信号的采样数，λ是多普勒雷达的波长，Φ[0]是初始相位，I[n]是基带的同相信号，Q[n]是正交信号，x[n]是胸壁位移信号；
通过从原始胸壁位移中减去最佳拟合线来消除多项式趋势，预处理后的胸壁位移信号可表示为：
多分辨率分析(MRA)建立了小波变换和数字滤波器之间的联系，基于多分辨率分析理论，采用基于最大重叠离散算法(MODWT)的MRA方法对多频段预处理胸壁位移信号进行分解，实现心跳和呼吸信号的初步分离；
小波变换本质上是通过母小波ψ(t)的尺度和以为操作实现的多尺度分析，对于任意函数f(t)∈L2(R)，母小波的中心频率记为f0，则
因此小波变换的时域形式可表示为：
其中a是与频率信息对应的比例因子，b是与时空信息相关的位移因子，通过设置尺度参数a来改变窗口的形状，为心跳等较快的时间提供良好的时间分辨率，为呼吸较慢的事件提供出色的频率分辨率，这是小波变换的多分辨率特性；
然而，连续小波变换基函数ψa，b(t)对于快速心率检测具有很大的计算冗余，通过对参数a和b进行离散化，转化为离散小波变换，计算效率大大提高，所以让
得到离散小波的表达式为：
使a0＝2，b0＝1，上述表达式可以简化为二次型小波，其表达式为：
FTPR算法检测心跳信号的具体步骤如下：
1)将心跳信号的时域表达式xh[n]经过转化变为xh win[n]，即在时域中增加时间窗口可以有效的减少信号在频域中的流失，
2)对xh win[n]进行傅里叶变换对应的频谱，并确定频谱峰值，保留峰值对应的频率及其相邻频率，舍弃剩余的频谱后获得重构频谱，
3)对重构频谱进行反傅里叶变换来计算复杂信号xh rec[n]，然后对xh rec[n]使用解调技术获得相位
4)最后通过公式
得到所需要的频率，
基于频率时间相位回归技术与时间窗口变化技术相结合的算法(FTPR-TWV)的心跳频谱估计方法具体流程是：
1)首先，将区间为T的时间窗内心跳信号分别2fs份，
2)对2fs个时间窗分别进行频谱分析，
3)再把2fs份频谱分析结果组合起来，
4)取最优时间窗，
5)最后用FTPR算法进行心跳信号的估计。
本发明的优点在于：
1)使用多普勒雷达，提高检测灵敏度。
2)MRA和TMF算法实现了心肺信号的分离和隐藏心跳信号的恢复。
3)将频谱细化算法引入心跳频谱估计，该方法实现了快速、精确的心率估计，并解决了短时窗下频谱分辨率不足的难题。
4)本发明的算法在短时间内也可以精确且快速地提取心率特征，并具有很好的有效性和鲁棒性。
附图说明
图1是本发明基于多普勒雷达心跳检测流程图。
图2a是本发明基于马勒特算法的Coiflet5的标度函数图。
图2b是本发明基于马勒特算法的Coiflet5的小波函数图。
图3是本发明解调前，预处理后、MRA算法处理后以及TMF算法处理后的胸壁位移信号及其频谱图。
图4是本发明在具体实验环境下的各项实验参数图。
图5是本发明FTPR-TWV算法的流程图。
图6a是本发明在T＝10s时HRV分析结果图。
图6b是本发明在T＝3s时HRV分析结果图。
具体实施方式
下面结合附图对本发明作进一步详细的阐述。
本发明提出了一种基于多普勒雷达的非接触式精确心跳检测方法。首先，采用连续波多普勒雷达体征检测的原理捕捉人体心脏跳动引起的微弱胸腔位移引起的频移，即雷达向受试者通过发射端发射T(t)。随后T(t)信号的相位被受试者胸壁运动x(t)所调制，反射信号被雷达接收端所捕获；其次，反射信号经过前端的低噪声放大器后，得到基带信号，采用正交混频器解调基带信号，再经过增益放大器后，由反正切解调和接卷绕计算可得到胸壁信号；最后，胸壁模拟信号经过模数转换器得到数字信号x(t)，在数字信号处理后在PC端经过马勒特算法和匹配滤波算法分离提取心脏信号，再经过FTPR-TWV算法后得到心跳频谱结果。
如图1所示：一种基于FTPR和TWV组合的多普勒雷达心跳准确检测的频谱估计
方法，其步骤如下：
步骤一：采用多普勒雷达发射雷达信号获取回波信号，并用信号解调技术进行正交信号解调；
步骤二：通过前端的低通滤波器进行杂波抑制和胸壁位移微分得到预处理后的胸壁位移信号；
步骤三：胸壁信号经过模数转换器得到数字信号；
步骤四：利用多频率分析算法(MRA)实现呼吸信号与心跳信号的分离；
步骤五：利用模板匹配算法(TMF)恢复隐藏在呼吸信号中的心跳信号；
步骤六：设计并采用频率时间相位回归技术(FTPR)与时间窗口变化技术(TWV)相结合的算法估计心跳频谱；
所述的获取回波信号的方法是：主控芯片采用IWR1443芯片，运算速度快，拥有丰富的外设结构。通过IWR1443芯片上的雷达模块发射端发射毫米波T(t)到目标特征体处(雷达信号可表示为
中，f为雷达工作频率，
为初始相位。由频域中相移的变化来捕捉由呼吸和心跳引起的胸壁位移信号，反射信号R(t)被雷达接收端所捕获，反射信号可表示为：
其中c为无线电的波速，λ＝c/f表示雷达波长；
把R(t)近似看作经过2d0/c时间延迟并相位调制后的T(t)，反射信号R(t)再经过前端的低噪声放大器后，得到基带信号B(t)，
其中，
为常数相位偏移，与接收机本身参数以及d0相关；残余相位噪声简化为
然后，由正交混频器(I/Q Mixer)解调基带信号B(t)，同相和正交信号分别表示为：
再经过增益放大器(Gain)后，由呼吸和心跳引起的胸壁信号x(t)由反正切解调和解卷绕技术计算可得，表达式如下：
最后，模拟信号x(t)经过模数转换器(ADC)得到数学信号x(t)，用于数字信号处理(DSP)以及PC端心跳检测算法的开发；
对得到的雷达所捕获的胸壁位移信号进行预处理，预处理步骤包括去除背景噪声和多项式趋势、带通滤波。由于雷达测量范围内除受试者胸壁运动外还存在大量静止物体，通过减去原始雷达回波数据的平均值来去除静止背景噪声。另外，由于热噪声和时间漂移引起的发射单元幅值的不稳定导致回波数据呈现严重的多项式趋势，通过减去最小二乘法拟合的曲线来去除多项式趋势。此外，基于呼吸频率典型范围为0.1～0.3Hz、HR范围为1～3Hz这一事实，带通滤波器选择巴特沃斯滤波器，其带通范围设置为0.8～4Hz，并应用于去噪后的胸壁位移信号，就得到了预处理后的胸壁位移信号。
使用多普勒雷达检测生命体征是基于这样一个事实，即由于多普勒效应，由呼吸和心跳引起的胸壁位移对传输的微波信号进行相位调制。因此，多普勒雷达的前端采用了正交接收机结构。首先，通过正交信号解调技术获得相位。然后，可以根据与相位成正比的关系计算胸壁位移。我们将反正弦解调算法引入到正交信号的解调中。由于摒弃了传统的反正切方法，低采样率的多普勒雷达系统可以避免相位缠绕和误差累积问题。因此，解调后的胸壁位移信号可表示为：
其中N是雷达原始信号的采样数，λ是多普勒雷达的波长，Φ[0]是初始相位，I[n]是基带的同相信号，Q[n]是正交信号，x[n]是胸壁位移信号；
信号预处理包括去除背景噪声和消除多项式趋势信号去除环境杂波，增强生命体征三个部分。本发明假设雷达测量范围内的传播环境是静态的，也就是说需要测量生命体征的生命体是静止不动的。除主体的胸壁运动外，可以通过平均胸壁位移数据来捕获静止散射体的恒定特征，然后可以通过减去平均值来去除背景杂波。另外由热噪声和时间漂移引起的发射单元的幅度不稳定性会导致胸壁位移信号的多项式趋势。严重的多项式趋势会造成低频分量影响甚至掩盖主频。因此，有必要通过从原始胸壁位移中减去最佳拟合线来消除多项式趋势。预处理后的胸壁位移信号可表示为：
采用心跳检测算法获取心跳频谱的方法是：为实现复杂环境下多普勒雷达对心跳信号的精确估计，本发明提出了一种雷达体征数据处理的具体实现方案，主要包括三部分：雷达信号预处理、心跳信号提取、心跳频谱估计。雷达捕获的受试者胸壁位移信号x(t)作为本发明算法的输入。
多分辨率分析(MRA)建立了小波变换和数字滤波器之间的联系，基于多分辨率分析理论，采用基于最大重叠离散算法(MODWT)的MRA方法对多频段预处理胸壁位移信号进行分解，实现心跳和呼吸信号的初步分离；
小波变换本质上是通过母小波ψ(t)的尺度和以为操作实现的多尺度分析，对于任意函数f(t)∈L2(R)，母小波的中心频率记为f0，则
因此小波变换的时域形式可表示为：
其中a是与频率信息对应的比例因子，b是与时空信息相关的位移因子，通过设置尺度参数a来改变窗口的形状，为心跳等较快的时间提供良好的时间分辨率，为呼吸较慢的事件提供出色的频率分辨率，这是小波变换的多分辨率特性；
然而，连续小波变换基函数ψa，b(t)对于快速心率检测具有很大的计算冗余，通过对参数a和b进行离散化，转化为离散小波变换，计算效率大大提高，所以让
得到离散小波的表达式为：
使a0＝2，b0＝1，上述表达式可以简化为二次型小波，其表达式为：
与传统的傅里叶变换不同，小波变换具有小波基函数，如Haar小波、Daubechies小波、Morlet小波、Symlet小波和Coiflet小波。本发明根据心跳信号的特点，选择Coiflet小波作为母小波。原因是Coiflet小波的小波函数和尺度函数具有良好的对称性，可以避免相位失真。其次，Coiflet小波函数具有合适的紧支撑长度，并且在有效支撑区域外迅速衰减。最后，Coiflet小波的形状类似于心跳信号的形状。Coiflet5小波的标度函数和小波函数分别如图2a，图2b所示。
选取Coiflet小波作为离散小波变换的母小波，选取消失矩阶为5阶，分解层数为6层。经过MODWT预处理后得到胸壁位移信号
计算出的小波变换矩阵表示为WT(L+1)×N。其中L使分解层数，N是位移信号的样本数。然后对小波变换矩阵进行MRA算法变换，计算出的MRA矩阵为WTmra，其不同行对应不同尺度的MRA结果。
由于心跳频谱通常在0.8Hz到2Hz之间变化，其频率范围对应于矩阵WTmra的第k行到第l行之间。此外，小波系数的能量分布特性与信号的频率特性密切相关。因此，基于能量比的不同频段的加权重建可以突出心跳的基频分量。每个尺度的能量可以表示为：
其中i表示小波分解的第i层，i＝1，2，...，L+1。
表示第i行第n个小波系数，E(i)表示第i层的能量。基于马勒特算法，能量加权的重构心跳信号可以表示为：
xmra[n]表示使用MRA算法进行初步分离和重建后的胸壁位移信号。
模板过滤匹配器：在信号处理领域，匹配滤波器是一种最优滤波器，其最优准则是由输出信号的最大信噪比来决定的。对于雷达系统，使用发射信号作为模板信号的匹配滤波器已普遍应用于脉冲雷达。尽管如此，匹配滤波器很少用于多普勒雷达的生命体征检测领域。模板匹配滤波器通过输入信号与模板信号之间的卷积运算，可以最大限度地恢复与隐藏在原始输入信号中的模板信号相似的信号。模板匹配滤波器的关键是模板信号的选择，直接选择低速段信号作为模板信号。然后通过多项式拟合去除模板信号中的身体运动，得到心跳模板信号。该方法的难点在于多项式拟合阶数的选择和提取精度低。因此，本发明提出了一种模板选择方法来避免呼吸干扰。在实验过程中，要求受试者保持呼吸几秒钟，已导出仅由心跳引起的胸壁位移信号作为模板信号。随后，能量重构的位移信号与模板信号在时域进行卷积运算。然后，通过模板匹配滤波器后的信号可以表达为：
xh[n]＝xmra[n]*h*[-n] (3f)
其中h*[-n]是模板信号的取反共轭。为了直观地说明我么们提出的分离和提取算法，通过检查胸壁位移波形的可视化和相应的频谱来研究差异，如图3所示，图3左侧图所示为不同处理阶段的胸壁位移波形，图3右侧图所示为对应的谱图。将图3a和3b与图3c和3d进行比较，可以很容易地看出，预处理有效地抑制了背景杂波，呼吸和呼吸引起的胸壁位移信号得到了有效抑制，心跳信号明显增强。如图3e和图3f所示，位移幅度减小了约5倍，呼吸频率基频分量被消除，表明应用MRA算法实现了呼吸和心跳的初步分离。此外，图3g和图3h表明位移信号具有很强的周期性，其主频谱峰值非常显著，表明TMF可以恢复隐藏在呼吸信号中的心跳信号。
使用多普勒雷达进行心跳检测本质上一个频率估计问题，频率分辨率决定了心跳估计的精度。对于小于5秒的短时间窗口，传统的傅里叶算法缺乏频率分辨率，难以快速准确地估计心率。
基于当信号存在主频率相位随时间线性变化的客观事实，FTPR算法将传统的频域峰值检测转换为时域中相位和时间的线性回归拟合，相位的斜率与信号的主频率成正比。以时间为10秒的时间窗口为例，传统傅里叶算法的频率分辨率是0.1Hz，而FTPR算法的频率分辨率为0.01Hz。由于频率分辨率的提高，心跳信号检测精度大大提高。
FTPR算法检测心跳信号的具体步骤如下：
1)将心跳信号的时域表达式xh[n]经过转化变为xh win[n]，即在时域中增加时间窗口可以有效的减少信号在频域中的流失，
2)对xh win[n]进行傅里叶变换对应的频谱，并确定频谱峰值，保留峰值对应的频率及其相邻频率，舍弃剩余的频谱后获得重构频谱，
3)对重构频谱进行反傅里叶变换来计算复杂信号xh rec[n]，然后对xh rec[n]使用解调技术获得相位
4)最后通过公式
得到所需要的频率。
然而，在短时窗内，心跳频谱的主峰会严重偏离真实的心跳频率。原因在于首先提取的心跳信号中仍有残留的谐波噪声；其次，由于频率分辨率不够，真实的心跳频率不是频率分辨率的整数倍。针对这些问题，本发明提出了一种结合FTPR和TWV的频谱估计新方法。
对于时间为T秒的时间窗，频率分辨率
为了准确地测量心跳频率，必须满足以下两个条件或者满足两个条件中的一个：(1)Δf＜＜HR，(2)HR＝kΔf，其中k为整数。显然，增加时间窗的长度可以满足检测精度，但实时检测的性能却大大降低。另外，对于周期信号，如果期望频率与频率分辨率的关系不是整数倍，则会产生频谱泄露，降低基频能力。该比值越接近整数倍，基频能力泄露越小。基于这一特性，可以确定最佳时间窗口以提高心率检测精度。
为显示TWV组合频谱的直观优势，以T＝3s和fs＝32Hz的时间窗为例，生成了一系列不同的时间窗。如图4所示，3s时间窗的频谱峰值为0.215，对应的频率为1.333Hz。然后，组合频谱中的最大峰值为0.223，对应于3.179秒的时间窗口和1.255Hz的频率，所以3.179秒的时间窗口是最合适的。基于上述分析，引入TWV算法可以避免严重的主频估计误差。
基于频率时间相位回归技术与时间窗口变化技术相结合的算法(FTPR-TWV)的心跳频谱估计方法具体流程是：
1)首先，将区间为T的时间窗内心跳信号分别2fs份，
2)对2fs个时间窗分别进行频谱分析，
3)再把2fs份频谱分析结果组合起来，
4)取最优时间窗，
5)最后用FTPR算法进行心跳信号的估计
图5为FTPR-TWV算法的流程图。
实验参数介绍
本发明体征检测实验在真实的办公环境下实施，实验的参考心跳数据由传统的接触式ECG模块测量，ECG电极贴于受试者心脏附近，采样频率设为500Hz。另外，实验采用德州仪器公司设计的IWR1443Boost毫米波雷达平台来测量体征信息，其相关实验参数如表1所示。实验过程中，受试者静坐与雷达正前方0.5m处，胸腔与雷达基本处于同一水平线，保持正常呼吸，记录两分钟实验数据。
此外，在实际应用中，需要选择合适的时间分辨率与频率分辨率。一般而言，时间窗口越短，时间分辨率越好，呈现更多心跳细节，但同时频率分辨率越低，检测精度越差。为验证本文算法的普适性，选择3s和10s两种不同时间窗口长度用于实时心率检测。
表1
性能分析：
为验证本发明提出的心跳检测算法的有效性与鲁棒性，我们在不同的实验条件下测量了多组实验数据，并采用多个评价指标对本文方法的性能进行评估。
此外，在实际应用中，需要选择合适的时间分辨率与频率分辨率。一般而言，时间窗口越短，时间分辨率越好，呈现更多心跳细节，但同时频率分辨率越低，检测精度越差。本发明基于多普勒雷达非接触式检测心跳方法具有普适性。
心率(HR)和心率变异性(HRV)是人们在心脏跳动方面衡量身体健康的两个重要指标。
心率(HR)是指每分钟的平均心跳次数，常被用做健康检测的基本指标。心率准确度(HRA)定义为多普勒雷达检测到的心率与ECG检测到的参考心率偏差在±2％以内的时间百分比，我们也将其作为心率指标。
心率变异性(HRV)是指连续心跳之间的时间间隔，揭示了每次心跳的规律性。心率变异性分析对心脏疾病的诊断和治疗具有重要意义。不同于在时域中通过检测相邻心跳的时间间隔来确定心率，本发明采用在频域中搜索心跳频率的心率检测方法，其中连续心跳时间间隔BBIs可用BPM来表示：
另外，平均相对误差(MRE)能被表示为：
HRV分析评价指标：RR间期标准差(SDNN)、连续RR间期的均方根差(RMSSD)，其中，NBBI为连续心跳间隔数，可分别表示为：
此外，Bland-Altman分析也常用于医学临床研究，以评估不同参数的一致性。我们也采用这种方法来评估雷达和心电图BBI的一致性。其表示为：
LoAL＝Bias-1.96SD (8)
LoAU＝Bias+1.96SD (9)
表2a和表2b为FTPR-TWV算法与FFT算法和FFT-TWV算法的HR分析结果。其中，三种算法的HR准确率均在92％以上，尤其是FTPR-TWV的平均HR准确率为99.7％，与参考值吻合度较高。FFT-TWV、FTPR和FTPR-TWV的平均HR准确度分别为82.204％、42.938％、92.091％。FTPR-TWV的HR精度明显高于其他两种。此外，FTPR-TWV的平均RMSE为0.901BPM，低于FFT-TWV，进一步验证了本发明算法在短时窗下也可以呈现更多的心跳细节。
表2a
表2b
HRV分析结果如表3a、表3b所示。从表中可知，FTPR-TWV算法的MRE范围为0.358％～0.636％，显示出较高的一致性。此外，小的RMSSD表明BBI的波动很小，意味着HR细节的损失很小。FFT-TWV、FTPR和FTPR-TWV算法的MRE平均值分别为1.168BPM、2.229BPM和0.910BPM。这些算法和ECG之间HRV评估指标的平均差异分别为SDNNs(1.856ms，3.688ms，0.910ms)，RMSSDs(4.408ms，4.610ms，2.839ms)。与FTPR的性能相比，由于将TWV引入FTPR，FTPR，FTPR-TWV的上述错误会大大减少。最后，FTPR-TWV的所有偏差幅度都低于0.8ms，进一步验证了我们提出的方法在短期时间窗中的有效性和鲁棒性。
表3a
表3b
为了直观地展示我们提出的算法在不同时段窗口中的性能优势，图6a和图6b展示了不同算法在不同时间窗口中的HR估计结果，ECG作为参考。如图6a所示，三种算法在偏差范围内的HR准确率约为98％。但是，由于主频发散，FTPR的HR曲线有两个严重的偏差点还有待发现，心跳细节的损失也很严重。如图6b所示，FRPR算法的HR准确度仅为43.220％，这是由于短周期时间窗中存在严重的主频估计偏差。因此，在FTPR算法中引入TWV算法来避免该问题，并确定最佳时间窗口。从图中可以发现，FTPR-TWV的HR准确度高达94.068％。此外，FTPR-TWV算法在跟踪心跳大突变点和呈现心跳多样性方面可以实现高一致性，具有强大的有效性和鲁棒性。
Claims (1)
1.一种基于多普勒雷达心跳准确检测的频谱估计方法，其特征在于，其步骤如下：
步骤一：采用多普勒雷达发射雷达信号获取回波信号，并用信号解调技术进行正交信号解调；
步骤二：通过前端的低通滤波器进行杂波抑制和胸壁位移微分得到预处理后的胸壁位移信号；
步骤三：胸壁信号经过模数转换器得到数字信号；
步骤四：利用多频率分析算法实现呼吸信号与心跳信号的分离；
步骤五：利用模板匹配算法恢复隐藏在呼吸信号中的心跳信号；
步骤六：设计并采用频率时间相位回归技术与时间窗口变化技术相结合的算法估计心跳频谱；
所述的获取回波信号的方法是：通过芯片上的雷达模块发射端发射频率为77GHz的毫米波T(t)到目标特征体处，雷达信号可表示为
由频域中相移的变化来捕捉由呼吸和心跳引起的胸壁位移信号，反射信号R(t)被雷达接收端所捕获,反射信号可表示为：
其中c为无线电的波速，λ＝c/f表示雷达波长；
把R(t)近似看作经过2d0/c时间延迟并相位调制后的T(t)，反射信号R(t)再经过前端的低噪声放大器后，得到基带信号B(t)，
其中，
为常数相位偏移，与接收机本身参数以及d0相关；残余相位噪声简化为
然后，由正交混频器(I/Q Mixer)解调基带信号B(t)，同相和正交信号分别表示为：
再经过增益放大器后，由呼吸和心跳引起的胸壁信号x(t)由反正切解调和解卷绕技术计算可得，表达式如下：
最后，模拟信号x(t)经过模数转换器得到数学信号x(t)，用于数字信号处理以及PC端心跳检测算法的开发；
对得到的雷达所捕获的胸壁位移信号进行预处理，预处理步骤包括去除背景噪声和多项式趋势、带通滤波，
多普勒雷达的前端采用了正交接收机结构，首先，通过正交信号解调技术获得相位，然后，根据与相位成正比的关系计算胸壁位移，将反正弦解调算法引入到正交信号的解调中，解调后的胸壁位移信号表示为：
其中N是雷达原始信号的采样数，λ是多普勒雷达的波长，Φ[0]是初始相位，I[n]是基带的同相信号，Q[n]是正交信号，x[n]是胸壁位移信号；
通过从原始胸壁位移中减去最佳拟合线来消除多项式趋势，预处理后的胸壁位移信号可表示为：
多分辨率分析(MRA)建立了小波变换和数字滤波器之间的联系，基于多分辨率分析理论，采用基于最大重叠离散算法的MRA方法对多频段预处理胸壁位移信号进行分解，实现心跳和呼吸信号的初步分离；
小波变换本质上是通过母小波ψ(t)的尺度和以为操作实现的多尺度分析，对于任意函数f(t)∈L2(R)，母小波的中心频率记为f0，则
因此小波变换的时域形式可表示为：
其中a是与频率信息对应的比例因子，b是与时空信息相关的位移因子，通过设置尺度参数a来改变窗口的形状，为心跳等较快的时间提供良好的时间分辨率，为呼吸较慢的事件提供出色的频率分辨率，这是小波变换的多分辨率特性；
然而，连续小波变换基函数ψa，b(t)对于快速心率检测具有很大的计算冗余，通过对参数a和b进行离散化，转化为离散小波变换，计算效率大大提高，所以让
得到离散小波的表达式为：
使a0＝2，b0＝1，上述表达式可以简化为二次型小波，其表达式为：
FTPR算法检测心跳信号的具体步骤如下：
1)将心跳信号的时域表达式xh[n]经过转化变为xh win[n]，即在时域中增加时间窗口可以有效的减少信号在频域中的流失，
2)对xh win[n]进行傅里叶变换对应的频谱，并确定频谱峰值，保留峰值对应的频率及其相邻频率，舍弃剩余的频谱后获得重构频谱，
3)对重构频谱进行反傅里叶变换来计算复杂信号xh rec[n]，然后对xh rec[n]使用解调技术获得相位
4)最后通过公式
得到所需要的频率，
基于频率时间相位回归技术与时间窗口变化技术相结合的算法(FTPR-TWV)的心跳频谱估计方法具体流程是：
1)首先，将区间为T的时间窗内心跳信号分别2fs份，
2)对2fs个时间窗分别进行频谱分析，
3)再把2fs份频谱分析结果组合起来，
4)取最优时间窗，
5)最后用FTPR算法进行心跳信号的估计。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111190185.3A CN113951856A (zh) | 2021-10-13 | 2021-10-13 | 一种基于多普勒雷达心跳检测的频谱估计方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111190185.3A CN113951856A (zh) | 2021-10-13 | 2021-10-13 | 一种基于多普勒雷达心跳检测的频谱估计方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN113951856A true CN113951856A (zh) | 2022-01-21 |
Family
ID=79463674
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111190185.3A Pending CN113951856A (zh) | 2021-10-13 | 2021-10-13 | 一种基于多普勒雷达心跳检测的频谱估计方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113951856A (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115118304A (zh) * | 2022-06-28 | 2022-09-27 | 西安电子科技大学 | 一种用于矢量信号分析仪的跳频信号参数测量方法 |
| CN115844355A (zh) * | 2022-11-22 | 2023-03-28 | 哈尔滨理工大学 | 一种uwb生物雷达回波信号心率提取方法 |
| CN116269260A (zh) * | 2023-03-01 | 2023-06-23 | 亿慧云智能科技（深圳）股份有限公司 | 一种智能手表心率监测方法及系统 |
-
2021
- 2021-10-13 CN CN202111190185.3A patent/CN113951856A/zh active Pending
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115118304A (zh) * | 2022-06-28 | 2022-09-27 | 西安电子科技大学 | 一种用于矢量信号分析仪的跳频信号参数测量方法 |
| CN115118304B (zh) * | 2022-06-28 | 2024-03-22 | 西安电子科技大学 | 一种用于矢量信号分析仪的跳频信号参数测量方法 |
| CN115844355A (zh) * | 2022-11-22 | 2023-03-28 | 哈尔滨理工大学 | 一种uwb生物雷达回波信号心率提取方法 |
| CN116269260A (zh) * | 2023-03-01 | 2023-06-23 | 亿慧云智能科技（深圳）股份有限公司 | 一种智能手表心率监测方法及系统 |
| CN116269260B (zh) * | 2023-03-01 | 2023-07-25 | 亿慧云智能科技（深圳）股份有限公司 | 一种智能手表心率监测方法及系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113951856A (zh) | 一种基于多普勒雷达心跳检测的频谱估计方法 | |
| CN110584631B (zh) | 一种基于fmcw雷达的静态人体心跳和呼吸信号提取方法 | |
| Ahmad et al. | Vital signs monitoring of multiple people using a FMCW millimeter-wave sensor | |
| CN111856452B (zh) | 基于omp的静态人体心跳和呼吸信号分离重构方法 | |
| Mahdiani et al. | Is 50 Hz high enough ECG sampling frequency for accurate HRV analysis? | |
| Jezewski et al. | A novel technique for fetal heart rate estimation from Doppler ultrasound signal | |
| Sahoo et al. | De-noising of ECG signal and QRS detection using Hilbert transform and adaptive thresholding | |
| CN104161509B (zh) | 一种基于幅值谱的心率变异性分析方法及仪器 | |
| Hu et al. | Noncontact accurate measurement of cardiopulmonary activity using a compact quadrature Doppler radar sensor | |
| Sekine et al. | Non-contact heart rate detection using periodic variation in Doppler frequency | |
| Xu et al. | Baseline wander correction in pulse waveforms using wavelet-based cascaded adaptive filter | |
| Ibrahimy et al. | Real-time signal processing for fetal heart rate monitoring | |
| CN113854992A (zh) | 基于77GHz毫米雷达的非接触式精确心率检测方法 | |
| JP2003510722A (ja) | 生理学信号を処理するための方法及び装置 | |
| CN115644840A (zh) | 基于毫米波雷达的生命体征检测方法 | |
| CN109805931A (zh) | 基于太赫兹多普勒雷达的远距离生命微动信号检测方法 | |
| CN108814579B (zh) | 一种基于emd分解的心电、呼吸联合计算心率变异性的方法 | |
| Tan et al. | Non-contact heart rate tracking using Doppler radar | |
| CN111012319A (zh) | 皮肤血流和血管的监测成像方法、系统及存储介质 | |
| Ji et al. | Systematic heartbeat monitoring using a fmcw mm-wave radar | |
| CN117838083A (zh) | 一种基于毫米波雷达的体征快速精确检测方法 | |
| CN105105739B (zh) | 短距离无线心率及心率变异性检测方法 | |
| Pan et al. | A spectrum estimation approach for accurate heartbeat detection using Doppler radar based on combination of FTPR and TWV | |
| CN115813363A (zh) | 一种基于心跳二次谐波重构的心率估计方法 | |
| Feng et al. | Non-invasive monitoring of cardiac function through Ballistocardiogram: An algorithm integrating short-time Fourier transform and ensemble empirical mode decomposition |
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 |