APP下载

基于谐波陷波器改进HMLD 的呼吸心跳分离算法

2023-02-08汪新坤张文艳

智能计算机与应用 2023年1期
关键词:陷波基波谐波

汪新坤,曹 乐,张文艳

(上海工程技术大学 电子电气工程学院,上海 201620)

0 引言

基于毫米波雷达可以不受天气,光线等一些环境因素的干扰进行非接触式生命体征检测,对医疗监护和智能家居健康监测具有重要意义和发展前景[1]。但是,基于毫米波雷达的生命体征信号检测仍存在一些困难和挑战,例如测试环境中的静态杂波、人为肢体运动伪影、呼吸谐波等干扰因素都会影响呼吸和心跳信号的提取[2]。一般可以通过滤波处理削弱杂波和噪声干扰,但强呼吸次谐波有可能在心跳有效的检测频率范围内,且心跳信号运动幅度更小,对心跳信号提取造成一定的困难,因此强呼吸次谐波的干扰已成为毫米波雷达体征信号检测中亟待解决的问题。

为了克服呼吸谐波对心跳信号提取的影响,Van Nguyen 团队[3-4]前后提出一种谐波路径算法(Harmonic Path Algorithm,HAPA)和平均谐波路径算法(Spectrum-Averaged Harmonic Path,SHAPA),但这两种算法存在无法找到谐波路径的问题,并且算法容易受阈值设置的影响,算法效率较低;Zhang[5]提 出 一 种 谐 波 倍 数 循 环 检 测 的 算 法(Harmonic Multiple Loop Detection,HMLD)算法,但该算法只计算到呼吸的二次谐波,未考虑到呼吸的三次谐波带来的影响,若身体发生与呼吸同频的噪声,会增强呼吸谐波的强度,继而影响心跳频率的估计;张怡[6]设置的呼吸和心跳信号检测频段被错开,存在频段检测约束问题,所以该算法没有从真正意义上消除呼吸二次谐波及三次谐波带来的影响。为了解决这一问题,本文提出一种基于陷波器谐波消除改进的HMLD 算法,通过陷波滤波器将邻近心跳基波的呼吸谐波滤除,实现心跳基波频率的精确提取;采用连续小波逆变换对已提取到的呼吸和心跳频率重构时域信号;最后,通过静态人体实验进行分析,验证了基于谐波陷波器改进的HMLD 算法能够有效的将呼吸和心跳信号分离和提取。

1 基于毫米波雷达的生命体征检测原理

毫米波雷达发射电磁脉冲,当接触到人体胸廓时,由于身体的高反射率,发射的脉冲信号被反射并被接收。人体的呼吸和心跳引起胸廓周期性的扩张和收缩,雷达发射天线与胸廓的距离d(t)可由公式(1)表示[7]。

其中,d0为接收天线和人体胸廓的初始距离;dr和dh分别为呼吸信号幅度和心跳信号幅度;fr和fh分别为呼吸频率和心跳信号频率。

存在多个传播信道时,接收信号可由公式(2)表示[8]。

其中,t,τ分别表示雷达回波数据中的慢时间和快时间,Akδ(τ -τk(t))和分别表示胸廓运动和静态目标对快时间和振幅峰值的响应。τk(t)由d(t)决定,可由公式(3)表示。

其中,c表示为光速。

雷达将接收到的信号离散化为二维雷达回波矩阵,可由公式(4)表示。

其中,m,n分别以慢时间和快时间显示采样数;Tf是快时间采样间隔,对应系统ADC 采样速率;Ts是脉冲持续时间,对应系统信号采样频率。

雷达回波矩阵包含目标微动位移信息,Δφ目标信号前后相位差由公式(5)表示[9]。

其中,Δx为胸腔运动位移,λ为毫米波对应的波长。

根据公式(5)对雷达回波矩阵进行数据解析,提取胸廓微动信号。

2 基于谐波陷波器改进HMLD 的呼吸心跳分离

2.1 雷达回波预处理

为了准确提取人体胸廓目标微动信号,需要对雷达回波矩阵进行相应的数据预处理。预处理的算法流程图如图1 所示。

图1 预处理算法流程图Fig.1 Flow chart of preprocessing algorithm

因为雷达天线接收到物体反射的电磁波能量最大,反映在雷达回波的数据中即信号幅值最大。所以在雷达矩阵的快时间轴上检测信号幅值最大的所在的距离单元,即目标所在的位置,可由公式(6)表示[10]。

其中,Δd表示距离单元跨度;n表示所在的距离单元;d0表示初始距离。

根据幅值最大检测原理,在慢时间轴上依次寻找每个脉冲的目标距离单元,并提取该距离单元目标的相位信号,最后对提取的相位信号解析提取微动信号,并进行低通滤波和均值滤波,消除高频和静态噪声。对原始的胸廓微动信号进行滤波处理,提取得到的信号如图2 所示。

图2 胸廓微动信号提取、滤波后波形Fig.2 Thoracic micro motion signal extraction and filtered waveform

2.2 基于谐波陷波器改进的HMLD 原理及方法

HMLD 算法对胸廓微动信号在[0.1,0.8]Hz,[0.8,2.0]Hz 呼吸和心跳检测频段分别进行带通滤波,并在对应的频谱中循环查找信号基频和其他高次谐波频率,进而判断呼吸心跳信号基波频率是否存在,信号基频与谐波频率fN的关系可由公式(7)表示[11]。

其中,fbase为信号的基波频率。

由于正常人的呼吸基频在[0.1,0.5]Hz,心跳基频在[0.8,2.0]Hz,若考虑到呼吸高次谐波的影响,HMLD 算法在心跳信号频率中存在检测效率和检测频段约束问题。而基于谐波陷波器改进的HMLD算法可以消除呼吸高次谐波的影响。为进一步提高HMLD 算法有效性,忽略呼吸四次谐波和心跳三次谐波的存在,将呼吸和心跳检测频段分别设置为[0.1,1.5]Hz,[0.8,4.0]Hz。

改进的算法主要分为两部分:呼吸频率提取,心跳频率提取,其具体算法步骤如下。

呼吸频率提取:

步骤1对雷达回波信号预处理,得到胸腔运动信号x(t);

步骤2在呼吸检测频段对x(t)进行带通滤波和快速傅里叶变换(Fast Fourier Transform,FFT),获取呼吸检测频段的频谱信号;

步骤3在频谱图中查找所有的峰值频率;

步骤4在提取的所有峰值频率中以最大峰值频率作为呼吸信号第一条峰值路径的待定基波频率fR1_1;

步骤5根据公式(5),提取呼吸二次谐波频率fR2_1和3 次谐波频率fR3_1,式(8)和式(9):

步骤6在频谱图中剔除第一条峰值频率路径fR1_1、fR2_1和fR3_1,并且寻找最大的峰值频率作为呼吸信号第二条峰值路径的待定基波频率fR1_2;

步骤7重复步骤3~5,为提高检测效率,防止无法找到待定基波的谐波频率,最多查找3 条峰值频率路径Pn,如公式(10)。

步骤8计算每条峰值路径的平均功率,并且以最大的平均功率的峰值路径的待定基波频率估计呼吸频率fR,如公式(11)。

其中,p表示最大平均功率峰值路径的下标。

心跳频率提取:

步骤9在呼吸频率提取的基础上,胸腔运动信号x(t)通过陷波滤波器将呼吸的基波分量和谐波分量滤除得到x′(t);

步骤10在心跳信号检测频段上对x′(t)进行带通滤波和快速傅里叶变换(Fast Fourier Transform,FFT),获取心跳检测频段的频谱信号;

步骤11在心跳频谱图中查找所有的峰值频率;

步骤12在提取的所有的峰值频率中以最大峰值频率作为心跳信号第一条峰值路径的待定基波频率fH1_1;

步骤13根据公式(5),提取心跳二次谐波频率fH2_1为

步骤14在频谱图中剔除第一条峰值频率路径fH1_1和fH2_1,并且寻找最大的峰值频率作为呼吸信号第二条峰值路径的待定基波频率fH1_2;

步骤15重复步骤10~12,与呼吸检测相同,最多查找两条峰值频率路径,为

步骤16计算每条峰值路径的平均功率,并且以最大的平均功率的峰值路径的待定基波频率估计心跳频率fH,即式(14):

基于谐波陷波器改进的HMLD 算法流程图如图3 所示。

图3 基于谐波陷波器改进的HMLD 算法流程图Fig.3 Flow chart of improved HMLD algorithm based on harmonic notch filter

2.3 谐波陷波器设计

本文数据基于MATLAB 软件平台处理的,所以陷波滤波器的传递函数分子和分母系数可以通过MATLAB 的库函数获取。在实际的信号处理过程中,有效的数据频率在[0.1,2]Hz,而且数据在进行FFT 后存在频谱泄露问题,所以对陷波器的参数设置有一定的要求。经过多次的实验验证,本文中将陷波滤波器阶数设置为8 阶,滤波器的带宽设置为0.04 Hz,并且利用MATLAB 软件上绘制出1.0 Hz 和0.3 Hz 的陷波器幅值响应曲线如图4 所示,可以看到该陷波器陷波深度达到了100 dB,能够实现对单一频率信号滤除,且对频率周围其他信号影响较小。

图4 1.0 Hz 及0.3 Hz 陷波器幅值响应曲线Fig.4 Amplitude response curve of 1.0 Hz and 0.3 Hz notch filter

2.4 基于连续小波逆变换的信号重构

连续小波逆变换可以在信号的任意的时频位置对小波系数进行修改[12]。所以可以对提取的人体呼吸心跳频率在原始的胸腔运动信号中重构呼吸和心跳的时域信号。本文采用Morlet 小波基函数对呼吸和心跳信号进行分析,因为Morlet 小波具有良好的时频局部化特性且形状与心跳信号相近[13],小波变换公式(15):

其中,a表示尺度参数;b为平移参数;w0为Morlet 母小波中心频率。

根据Morlet 小波转换公式进行呼吸和心跳信号重构,具体步骤如下:

(1)根据谐波消除改进的HMLD 算法提取的呼吸或心跳频率f,设置小波变换的频率范围[f -BW,f +BW],BW为陷波滤波器的带宽;

(2)根据频率转换尺度的关系,计算尺度参数a,公式(16):

其中,w0为Morlet 母小波中心频率,fs为信号采样频率,f∈[f -BW,f +BW],将频率范围等间隔离散化成[f1,f2,…,fn],即可获得尺度序列S =[a1,a2,…,an]。

(3)根据不同尺度的Morlet 小波函数对胸腔运动信号x(t)进行小波变换,得到小波系数矩阵C;

(4)通过系数矩阵C进行小波逆变换,获得重构信号r(t)。

3 实验及结果分析

3.1 实验设计

本文采用自主设计的单通道的雷达信号采集系统进行实验信号采集。系统由模拟前端、主控模块及信号调理电路组成。其中模拟前端采用脉冲相干式A111 毫米波雷达芯,该芯片工作频率为60 GHz,波长λ=5 mm,检测范围为[0.2,0.8]m,距离单元跨度Δx =0.48 mm,系统的采样率为50 Hz。实验选取3 名平均年龄为24 岁的成年人作为实验对象,在初始位置d0=0.6 m 采集原始雷达回波信号,并用小米手环的实时监测心率值作为实验的心跳参考信号,每次采集时间为2 min。实验场景如图5 所示。

图5 实验场景图Fig.5 Experimental scene

3.2 呼吸心跳分离

选取测试对象1 的雷达回波数据对提取到的胸腔运动信号进行呼吸心跳信号分离。根据谐波陷波器改进的HMLD 算法,先对呼吸的检测频段进行带通滤波,并进行FFT 得到呼吸检测频谱,如图6 所示。

图6 呼吸检测频谱Fig.6 Respiratory detection spectrum

在呼吸检测频谱中查找平均功率最大的呼吸峰值路径为P =[0.300 5,0.608 40,0.908 5],因此此时呼吸频率为峰值路径下的基波频率fR =0.300 5 Hz。对于心跳频率的提取,先将呼吸的基波频率、二次谐波和三次谐波消除掉,然后在心跳检测频段进行带通滤波和FFT,得到心跳信号检测频谱图如图7 所示。

图7 心跳检测频谱Fig.7 Heartbeat detection spectrum

根据提取的呼吸频率fR和心跳频率fH,通过连续小波逆变换进行呼吸和心跳时域信号重构,结果如图8 所示,可以观察到重构后的呼吸和心跳信号整体趋势平稳、光滑。

图8 呼吸心跳时域信号重构Fig.8 Time domain signal reconstruction of respiration and heartbeat

3.3 结果分析

本文采用变分经验模态分解(Variational Mode Decomposition,VMD)算法对每名测试对象进行呼吸心跳分离,并对提取的结果进行信噪比(Single Nosie Ratio,SNR)比较,公式(17):

其中,s2(l)表示提取信号频谱峰值的功率,∑s2(f)为信号频谱的功率和。

基于谐波陷波器改进的HMLD 算法与VMD 算法提取心跳的结果见表1。

通过表1 可知,3 名测试对象采用改进的HMLD 算法提取心跳频率误差率分别为4.61%、11.5%、7.78%,提取心跳信号信噪比相比于VMD 平均提高了2.66 dB。,因此采用谐波陷波器改进的HMLD 算法提取的呼吸和心跳信号更为准确。

表1 基于陷波器改进HMLD 与VMD 算法的实验结果Tab.1 Experimental results of improved HMLD and VMD algorithm based on notch filter

4 结束语

本文针对基于毫米波雷达生命体征信号检测中存在呼吸谐波导致呼吸心跳难以分离的问题,提出了一种基于谐波陷波器改进的HMLD 算法在信号频谱中提取呼吸和心跳频率。首先,对3 名测试对象进行雷达回波预处理提取胸廓微动信号;其次,采用HMLD 算法提取呼吸频率,采用谐波陷波器消除呼吸谐波,再根据HMLD 算法提取心跳频率;最后,对提取的呼吸和心跳频率进行连续小波逆变换,对时域信号进行重构。结果表明,提取的心跳频率误差率在11.5%以内,心跳信噪比相比于VMD 平均提高了2.66 dB,说明基于谐波陷波器改进HMLD 算法能够准确的提取呼吸心跳信号,为基于毫米波雷达的高精度生命体征信号检测提供了一种有效的研究方法。

猜你喜欢

陷波基波谐波
基于SWT的电力系统基波检测
基于数字递归陷波的多通道瞬变电磁法周期噪声去除研究
SFC谐波滤波器的设计及应用
电力系统谐波检测研究现状及发展趋势
自适应的谐波检测算法在PQFS特定次谐波治理中的应用
电力系统谐波状态估计研究综述
最近电平逼近调制的基波谐波特性解析计算
温度对陷波网络阻抗的影响
卫星导航接收机基于IIR陷波器的单频干扰抑制性能分析
全相位FRM陷波原理及其DSP Builder实现