CN114869294B - 基于vmd分解和let模型的粒子滤波运动伪迹抑制方法 - Google Patents
基于vmd分解和let模型的粒子滤波运动伪迹抑制方法 Download PDFInfo
- Publication number
- CN114869294B CN114869294B CN202210479424.5A CN202210479424A CN114869294B CN 114869294 B CN114869294 B CN 114869294B CN 202210479424 A CN202210479424 A CN 202210479424A CN 114869294 B CN114869294 B CN 114869294B
- Authority
- CN
- China
- Prior art keywords
- electrocardiosignal
- particle
- amplitude
- value
- sampling point
- 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
- 239000002245 particle Substances 0.000 title claims abstract description 98
- 238000000034 method Methods 0.000 title claims abstract description 76
- 230000001629 suppression Effects 0.000 title claims abstract description 23
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 22
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 60
- 238000001914 filtration Methods 0.000 claims abstract description 55
- 230000003068 static effect Effects 0.000 claims abstract description 29
- 238000005070 sampling Methods 0.000 claims description 56
- 239000011159 matrix material Substances 0.000 claims description 28
- 238000005259 measurement Methods 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000012546 transfer Methods 0.000 claims description 9
- 230000002829 reductive effect Effects 0.000 claims description 8
- 230000005764 inhibitory process Effects 0.000 claims description 6
- 230000001133 acceleration Effects 0.000 claims description 4
- 235000003642 hunger Nutrition 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 4
- 238000007781 pre-processing Methods 0.000 claims description 4
- 230000037351 starvation Effects 0.000 claims description 4
- 238000012549 training Methods 0.000 claims description 3
- 238000000819 phase cycle Methods 0.000 claims description 2
- 238000012216 screening Methods 0.000 abstract description 3
- 230000000694 effects Effects 0.000 description 19
- 230000006870 function Effects 0.000 description 9
- 238000001228 spectrum Methods 0.000 description 9
- 230000009467 reduction Effects 0.000 description 6
- 238000005452 bending Methods 0.000 description 5
- 238000012545 processing Methods 0.000 description 4
- 238000011084 recovery Methods 0.000 description 4
- 238000012935 Averaging Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000012806 monitoring device Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000000284 resting effect Effects 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 238000012952 Resampling Methods 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000006793 arrhythmia Effects 0.000 description 1
- 206010003119 arrhythmia Diseases 0.000 description 1
- 230000000747 cardiac effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 208000029078 coronary artery disease Diseases 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 208000019622 heart disease Diseases 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 210000003205 muscle Anatomy 0.000 description 1
- 208000010125 myocardial infarction Diseases 0.000 description 1
- 230000003183 myoelectrical effect Effects 0.000 description 1
- 238000005312 nonlinear dynamic Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000000926 separation method Methods 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/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/318—Heart-related electrical modalities, e.g. electrocardiography [ECG]
- A61B5/346—Analysis of electrocardiograms
- A61B5/349—Detecting specific parameters of the electrocardiograph cycle
- A61B5/352—Detecting R peaks, e.g. for synchronising diagnostic apparatus; Estimating R-R interval
-
- 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
- A61B5/7207—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal of noise induced by motion artifacts
-
- 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/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Molecular Biology (AREA)
- Animal Behavior & Ethology (AREA)
- Pathology (AREA)
- Physics & Mathematics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Cardiology (AREA)
- Surgery (AREA)
- Biophysics (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Signal Processing (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physiology (AREA)
- Psychiatry (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明公开了一种基于VMD分解和LET模型的粒子滤波运动伪迹抑制方法,首先提取采集心电信号的周期,然后将静态心电信号和采集心电信号拼接得到组合心电信号,采用VMD算法对组合心电信号进行分解,筛选出参考分量和含噪分量,提取参考分量和含噪分量的包络信号,采用LET模型对含噪分量进行重构,再将重构的含噪分量与其余模态分量进行叠加,得到重构的采集心电信号,构建心电动态模型并估计其参数,最后根据心电动态模型和重构采集心电信号对采集心电信号进行粒子滤波,得到运动伪迹抑制后的采集心电信号。本发明通过采用VMD分解算法和LET模型对采集心电信号进行预处理,提高运动伪迹抑制效果。
Description
技术领域
本发明属于心电信号处理技术领域,更为具体地讲,涉及一种基于VMD分解和LET模型的粒子滤波运动伪迹抑制方法。
背景技术
近年来,可穿戴监测设备得益于传感器技术的革新,得到了极大的发展。在可穿戴监测设备所监测的数据中,心电信号是一项十分重要的生理指标。心电信号在临床上常常通过心电图(ECG)的多项特征对冠心病、心率失常、心肌梗死等心脏疾病进行诊断和预测。
不同于传统医院场景中较为静止的心电检测环境,在可穿戴的场景中,尤其是人体处于运动状态时,采集到的心电图常常因身体与电极之间的相对运动引入噪声,导致采集的心电信号波形失真,无法用于心电特征获取。因此,为了提升采集到的动态心电质量,有必要在可穿戴监测仪的软件总体系统中,对动态心电进行有效降噪处理。心电采集设备中经常引入的噪声包括基线漂移、工频干扰、肌电噪声和运动伪迹等,其中运动伪迹噪声因其产生原因与身体运动及皮肤形变相关,在可穿戴生理监测场景中尤为明显。同时,运动伪迹的频率分布与心电信号的主要频率存在重叠,无法直接使用如数字滤波等常见降噪方法达到较好的抑制效果。因此,对心电信号降噪算法进行深入研究,尤其是实现运动伪迹的降噪算法,对于提升可穿戴监测仪的功能性,恢复心电特征提高采集的心电信号可用性有重要意义。
在运动伪迹噪声抑制方法中,基于心电动态模型的粒子滤波算法是一种表现较为优良的算法。图1是基于心电动态模型的粒子滤波算法的流程图。如图1所示,在基于心电动态模型的粒子滤波算法中,首先对采集心电信号进行R波定位,以R波峰点为核心提取心电周期,对齐各心电周期数据后,获取心电平均信号ECGmean,然后从心电平均信号ECGmean中提取出心电动态模型参数,建立状态方程,通过平均心电周期组成与采集心电信号等长的心电序列以得到合成心电信号ECGsynth,最后根据状态方程和合成心电信号对采集心电信号进行粒子滤波实现对采集心电信号的运动伪迹噪声抑制。该算法的详细说明可以参见文献“Hesar H D,Mohebbi M.Muscle artifact cancellation in ECG signal using adynamical model and particle filter[C]//2015 22nd Iranian Conference onBiomedical Engineering(ICBME).IEEE,2015.”
虽然基于心电动态模型的粒子滤波算法可以对运动伪迹噪声实现较好的抑制,但是处理后的信号在部分原心电信号信噪比较低的位置,滤波效果较差;同时,平均心电信号方法最早用于受高斯白噪声影响的心电信号处理,在心率变化以及信号中存在非高斯噪声时平均心电信号与心电信号相关性将下降,导致算法整体效果降低。
此外,粒子滤波是在贝叶斯估计方法的基础上引入了序贯蒙特卡洛方法,使用粒子群来表示概率,并以权值来衡量各粒子的质量。在现有基于心电动态模型的粒子滤波算法中,粒子滤波算法在计算权值时利用了各粒子点与合成心电信号以及采集心电信号两者的马氏距离,使得粒子滤波在信噪比较低时能够根据估计的心电信号获得较好的滤波效果。但是在采集心电信号质量不同时,权值计算由于缺乏指导无法根据心电质量自适应变化,导致滤波性能下降,对运动伪迹噪声抑制造成不良影响。
发明内容
本发明的目的在于克服现有技术的不足,提供一种基于VMD分解和LET模型的粒子滤波运动伪迹抑制方法,通过采用VMD分解算法和LET模型对采集心电信号进行预处理,提高运动伪迹抑制效果。
为了实现上述发明目的,本发明基于VMD分解和LET模型的粒子滤波运动伪迹抑制方法包括以下步骤:
S1:对含有运动伪迹噪声的采集心电信号Xd(n)进行R波定位,提取心电信号周期,n=1,2,…表示采样点序号;
S2:将Q个完整心电信号周期的静态心电信号Xs(n)与采集心电信号Xd(n)进行拼接,得到组合心电信号X(n);
S3:采用VMD算法对组合心电信号X(n)进行分解得到K个模态分量IMFk(n),k=1,2,…,K,K表示模态总数;对于第1层至第K层模态分量,将中心频率和预设频率阈值最为接近的模态分量作为参考分量,记其序号为k*,将第1层至第k*-1层作为含噪分量;
S6:采用LET模型对含噪分量进行重构,具体方法如下:
S7:将步骤S6中拟合还原的采集心电分量IMF′d,k″(n)和其余模态分量叠加,得到重构的采集心电信号X′d(n):
S8:构建如下心电动态模型:
对采集心电信号X′d(t)使用最小二乘法估计得到心电动态模型的关键参数,包括采集心电信号中P/Q/R/S/T五部分波形对应高斯项的振幅ai、宽度bi、中心参数θi,i∈{P,Q,R,S,T};
S9:根据步骤S8确定的心电动态模型和重构采集心电信号X′d(n),采用粒子滤波算法对采集心电信号Xd(n)进行粒子滤波,得到运动伪迹抑制后的采集心电信号X* d(n)。
本发明基于VMD分解和LET模型的粒子滤波运动伪迹抑制方法,首先提取采集心电信号的周期,然后将静态心电信号和采集心电信号拼接得到组合心电信号,采用VMD算法对组合心电信号进行分解,筛选出参考分量和含噪分量,提取参考分量和含噪分量的包络信号,采用LET模型对含噪分量进行重构,再将重构的含噪分量与其余模态分量进行叠加,得到重构的采集心电信号,构建心电动态模型并估计其参数,最后根据心电动态模型和重构采集心电信号对采集心电信号进行粒子滤波,得到运动伪迹抑制后的采集心电信号。
本发明先根据运动伪迹的频域特性,使用受影响小的模态拟合并替换被噪声严重污染的模态,在尽可能保存心电周期特征的情况下实现采集心电信号的重构,使重构得到的采集心电信号更加准确,以提高粒子滤波的性能,提高对于运动伪迹的抑制效果。
附图说明
图1是基于心电动态模型的粒子滤波算法的流程图;
图2是本发明基于VMD分解和LET模型的粒子滤波运动伪迹抑制方法的具体实施方式流程图;
图3是LET模型的模块化结构图;
图4是本发明中基于LET模型重构含噪模态分量的流程图;
图5是本实施例中粒子滤波算法的流程图;
图6是步行状态下原始心电信号以及本发明和五种对比方法实现运动伪迹抑制后的波形对比图;
图7是步行状态下静态心电与本发明和五种对比方法滤波前后的心电频谱对比图;
图8是扩胸状态下原始心电信号以及本发明和五种对比方法实现运动伪迹抑制后的波形对比图;
图9是扩胸状态下静态心电与本发明和五种对比方法滤波前后的心电频谱对比图;
图10是弯腰状态下原始心电信号以及本发明和五种对比方法实现运动伪迹抑制后的波形对比图;
图11是弯腰状态下静态心电与本发明和五种对比方法滤波前后的心电频谱对比图。
具体实施方式
下面结合附图对本发明的具体实施方式进行描述,以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是,在以下的描述中,当已知功能和设计的详细描述也许会淡化本发明的主要内容时,这些描述在这里将被忽略。
实施例
图2是本发明基于VMD分解和LET模型的粒子滤波运动伪迹抑制方法的具体实施方式流程图。如图2所示,本发明VMD分解和LET模型的粒子滤波运动伪迹抑制方法的具体步骤包括:
S201:提取采集心电信号周期:
对含有运动伪迹噪声的采集心电信号Xd(n)进行R波定位,提取心电信号周期,n=1,2,…表示采样点序号。
S202:组合心电信号:
将Q个完整心电信号周期的静态心电信号Xs(n)与采集心电信号Xd(n)进行拼接,得到组合心电信号X(n)。为了降低对运动伪迹频域特征的影响,Q的优选取值范围为[3,5]。
S203:VMD分解:
采用VMD(Variational Mode Decomposition,变分模态分解)算法对组合心电信号X(n)进行分解得到K个模态分量IMFk(n),k=1,2,…,K,K表示模态总数,对于第1层至第K层模态分量,将中心频率和预设频率阈值最为接近的模态分量作为参考分量,记其序号为k*,将第1层至第k*-1层作为含噪分量。根据研究,运动伪迹噪声频率主要在1~10Hz范围内,因此本实施例中设置频率阈值为20Hz,从而保证选取的参考分量受运动伪迹影响较小。
VMD算法是一种常用的信号分解方法,对非平稳、非线性信号具有良好处理效果。该算法假设各本征模态都存在一个中心频率,当前模态在该频率附近会最紧凑。VMD算法的具体原理和步骤可以参考文献“K.Dragomiretskiy and D.Zosso,"Variational ModeDecomposition,"in IEEE Transactions on Signal Processing,vol.62,no.3,pp.531-544,Feb.1,2014.”。
在VMD算法中,需要设置的参数包括模态总数K,惩罚系数α和收敛容差ε。模态总数K决定了VMD算法分解出的本征模态层数,若K值偏小,会出现分解不充分的欠分解情况,导致模态混叠;若K值偏大,会出现分解成分超过信号中有用成分的过分解情况,导致出现虚假分量。
惩罚系数α决定IMF各分量的带宽大小同时影响分解精度,若α值偏大,会导致IMF带宽偏大,使得同一个IMF分量中包含不同成分信号即模态混叠;若α值偏小,会导致IMF带宽偏小,使得被分解信号中部分信号的丢失。
收敛容差ε是决定优化停止的准则之一,通常取值范围为10-6~5×10-6。
对于上述三个参数,由于本发明针对的是心电信号,因此本实施例中根据经验设置收敛容差ε=10-6,对于模态总数K和惩罚系数α采用以下方法确定:
基于模态总数K过大会出现虚假模态以及心电信号主要能量集中在0.5-40Hz内这两条性质,在固定惩罚系数α的情况下,在预设的模态总数取值范围内对各个取值进行搜索,计算相邻模态总数的能量差值,选取其中能量差值最大的模态总数取值作为最终使用的模态总数K。能量差值的计算公式如下:
其中,EK、EK-1分别表示模态总数取值为K和K-1时的所有模态分量的能量总和。
在确定模态总数K后,根据运动伪迹噪声频谱与心电局部重合这一性质,在固定模态总数K的情况下,在预设的惩罚系数取值范围内对各个取值进行搜索,计算每个惩罚系数取值下的模态带宽和其中BWk表示第k个模态分量带宽,在模态带宽和小于预设阈值的惩罚系数取值中选取一个作为惩罚系数α。本实施例中,在满足阈值要求的惩罚系数取值中选择计算时间复杂度最小的惩罚系数取值。
S204:提取模态分量的包络:
常用的包络提取方法包括希尔伯特变换法、平方检波滤波法、香农能量法等,根据VMD分解中模态分量是调频-调幅信号这一性质,本实施例中选择希尔伯特变换方法对各模态进行包络提取。
S205:包络信号划分:
S206:对含噪分量进行重构:
Volterra(沃尔特拉)核估计方法是目前常用的一种非参数建模方法之一,核心是Volterra模型,该模型的离散时间公式定义如下:
其中,x(n)、y(n)分别是输入信号和输出信号,Ts代表采样时间,kp(m1,...,mp)代表p阶Volterra核函数,M代表记忆长度。
Volterra核估计方法的最佳实现是使用Ogura等人提出的拉盖尔展开技术(Laguerre Expansion Technique,LET),利用离散Laguerre函数(Discrete Laguerrefunctions,DLFs)的标准正交集在离散时间下压缩Volterra核数量,从而实现更紧凑的Volterra核扩展形式,构成基于拉盖尔展开技术的沃尔特拉模型,也称为LET模型。图3是LET模型的模块化结构图。如图3所示,在LET模型中,将输入信号采用线性滤波器组进行滤波,再由多输入稳态非线性函数进行处理,得到输出信号。其中线性滤波器组由L个DLFs组成,DLFs与其对应输出如下式所示:
其中,bj(m)代表第j个DLFs,j=1,2,…,L;m为正整数,其取值范围在0至M-1,与/>分别表示从m和j中取k个值的组合计算得到系数/>与/>α取值范围在0和1之间,决定滤波器的衰减速率;vj(n)代表第j个滤波器的输出信号。
在确定模型阶数Q和滤波器总数L后,可以对滤波器组的输出信号vj(n)进行求解。为了降低计算复杂度,可以采用文献“Marmarelis V Z.Nonlinear dynamic modeling ofphysiological systems[M].Wiley-Interscience,2004.”中的递归求解方法对滤波器输出进行求解,求解公式为:
Volterra核函数扩展形式如下式所示,其中kp(m1,...,mp)代表p阶Volterra核函数:
将公式(3)和(5)代入(2)中得到修正的Volterra模型形式如下:
其中,{cp}是拉盖尔扩展后的核系数,p=1,2,…,Q,Q表示阶数,ε(n)是误差项,包含了可能的模型截断误差和数据中的噪声。LET模型的矩阵形式如下式所示:
y(n)=Vc+ε(n) (7)
C表示拉盖尔扩展后的核系数矩阵,表达式如下式所示:
核系数矩阵C的估计方法根据滤波器输出矩阵V是否满秩可以分为两种,当滤波器组输出矩阵V满秩时,使用最小二乘法进行估计:
使用递归公式(4)确定各级滤波器输出vj(n)后,Volterra核估计的任务就简化为对核扩展系数{cp}的估计,将滤波器输出矩阵V代入公式(10)或(11)求出扩展系数估计矩阵后,通过公式(7)计算获得输出y(n)。
基于以上分析可知,可以采用LET模型来对各个模态分量中的含噪分量进行重构。图4是本发明中基于LET模型重构含噪模态分量的流程图。如图4所示,本发明中基于LET模型重构含噪模态分量的具体步骤包括:
S401:计算滤波器输出矩阵:
S402:计算扩展系数矩阵:
由于各个含噪分量的特性有所不同,因此在实际应用中,可以针对各个含噪分量分别设置LET模型的参数,从而提高扩展系数矩阵的准确性。
S403:重构待拟合分量的采集心电包络信号:
S404:还原采集心电信号分量:
S207:重构采集心电信号:
将步骤S206中拟合还原的采集心电分量IMF′d,k″(n)和其余模态分量叠加,得到重构的采集心电信号X′d(n):
S208:估计心电动态模型参数:
构建心电动态模型,表达式如下:
其中,与zn分别是心电信号在采样点n的相位值与幅值,/>与zn-1分别是心电信号在采样点n-1的相位值与幅值,ω表示心电信号在x-y平台内的投影环的角速度,由心率大小决定,δ是采样间隔,η表示模拟噪声,用于模拟模型的不确定性,可以根据需要由不同的噪声模拟函数构成,mod表示求取余数。
对采集心电信号Xd′(t)使用最小二乘法估计得到心电动态模型的关键参数,包括采集心电信号中P/Q/R/S/T五部分波形对应高斯项的振幅ai、宽度bi、中心参数θi,i∈{P,Q,R,S,T}。
S209:粒子滤波:
根据步骤S208确定的心电动态模型和重构采集心电信号Xd′(n),采用粒子滤波算法对采集心电信号Xd(n)进行粒子滤波,得到运动伪迹抑制后的采集心电信号X* d(n)。
在粒子滤波算法中,粒子权值是一个非常重要的参数。现有的粒子滤波算法中,粒子权值计算由于缺乏指导无法根据心电质量自适应变化,导致滤波性能下降。因此,本实施例中提出了一种使用同步采集的运动数据指导权值计算的方法,对粒子滤波算法进行了改进。图5是本实施例中粒子滤波算法的流程图。如图5所示,本实施例中粒子滤波算法的具体步骤包括:
S501:构建状态空间模型:
在粒子滤波中,首先需要构建状态空间模型,包括状态转移方程和量测方程。由于本发明所针对的是心电信号,因此其状态转移方程采用心电动态模型,量测方程表示式如下:
S502:同步采集运动数据:
在对采集心电信号Xd(n)进行采集的同时,同步采集包括三轴加速度、加速度混合数据与三轴角速度,对于每个运动数据采用预设方法进行预处理后,再与采集心电信号Xd(t)进行采样速率同步,记得到的第p路运动数据为Sp(n),p=1,2,…,7。
本实施例中,运动数据的预处理方法为移动平均滤波,以去除噪声,采样速度同步的方法采用三次样条插值方法。
S503:筛选指导运动数据:
经过研究发现,当用户处于不同运动情况(如步行、坐下与起身、弯腰、扩胸)时,七路运动数据与采集心电信号的相关性有所不同。因此为了提高后续粒子权值计算的准确度,首先分别计算每路运动数据Sp(n)和采集心电信号Xd(n)之间的相关性,选择相关性最大的运动数据进行归一化后作为指导运动数据S*(n)。本实施例中相关性采用皮尔逊相关系数。
S504:计算指导参数数据:
采用如下公式计算得到指导参数数据γ(n):
γ(n)=τ+sigmoid(abs(ρS*(n))) (15)
其中,sigmoid()函数作为约束算子,τ表示预设的正常数,配合约束算子实现权值范围调整,取值范围为[0.8,1.2],ρ表示预设的运动数据调整系数,abs()表示求取绝对值。本实施例中所采用的约束算子sigmoid()的表达式如下:
其中,a、b是调整约束算子的两个参数,根据实验分别取0.1与0.2时可以获得较好的效果。
S505:粒子滤波初始化:
本实施例中,心电幅值粒子的初值由重构采集心电信号X′d(n)在采样点0的幅值与[0,0.5]间的随机数相加得到,相位粒子的初值/>由采集心电信号Xd(n)在采样点0的相位与[0,0.1]间的随机数相加得到。
S506:更新粒子权值:
其中,上标T表示转置,上标-1表示求逆,表示根据第g个粒子在上一采样点n-1的相位/>由状态转移方程计算得到的当采样点n该粒子的相位估计值,/>表示根据第g个粒子在上一采样点n-1的幅值/>由状态转移方程计算得到的当前采样点n该粒子的幅值估计值,当n=0,/> 表示采集心电信号Xd(n)在当前采样点n的相位,R表示重构采集心电信号X′d(n)协方差矩阵,其表达式为:
S507:计算幅值的量测估计值:
S508:评估是否出现粒子匮乏现象,即判断是否满足以下公式:
其中,λ表示预设阈值,本实施例中λ=2G/3。如果满足,则判断出现粒子匮乏现象,进入步骤S509,否则不作任何操作,直接进入步骤S510。
S509:粒子重采样:
S510:计算状态转移方程估计值:
S511:令n=n+1,返回步骤S506。
为了更好地说明本发明的技术效果,采用一个具体实例对本发明进行实验验证。本次实验验证中采集了4位25-26岁健康男性在日常生活中模仿三种不同运动状态下的心电信号与运动数据。测试者进行的3种运动模拟包括步行、弯腰与扩胸,每一组实验流程与之前相同,包括静止20s,持续动作20s以及恢复静止20s共60s。本次实验验证中选取了五种对比方法,包括最小均方算法(简称LMS算法),基于单通道的盲源分离算法(简称CCPC算法)、扩展卡尔曼算法(简称EKF算法),基于PCA与蒙特卡洛滤波算法(简称PMP算法)、基于心电模型的粒子滤波算法(简称DPF算法)。
图6是步行状态下原始心电信号以及本发明和五种对比方法实现运动伪迹抑制后的波形对比图。如图6所示,LMS算法降噪效果一般,这可能与参考信号选取有关;CCPC算法对于QRS波段恢复效果较好,但P、T波段恢复不佳;EKF算法、PMF算法、DPF算法与本发明这四种基于模型的降噪算法对波形恢复较好,尤其是对于P、T波段恢复效果明显,但可以看到EKF、DPF与PMF在采样点500、700以及550左右仍存在未滤除噪声,并且PMF算法R波峰幅值普遍较低,这与其预处理时提取的心电主成分与原始心电R波存在差异有关。
图7是步行状态下静态心电与本发明和五种对比方法滤波前后的心电频谱对比图。如图7所示,从频域看,步行MA噪声主要集中在0-10Hz内,本发明相比于其它五种对比方法在滤波后噪声更少,频谱更趋近于静态心电信号。
图8是扩胸状态下原始心电信号以及本发明和三种对比方法实现运动伪迹抑制后的波形对比图。如图8所示,在扩胸状态下MA噪声幅度较小,LMS与CCPC算法效果较步行数据有所提升,但波形仍存在较多毛刺;基于模型的四种算法中,EKF算法可以看出仍存在MA引起的漂移噪声,DPF算法在550采样点处仍存在未被完全去除的MA噪声,本发明与PMF算法总体效果均较为良好。。
图9是扩胸状态下静态心电与本发明和五种对比方法滤波前后的心电频谱对比图。如图9所示,在扩胸状态下从频域上看MA噪声主要集中在0~8Hz内,基于模型的四种方法对信号改善效果均较明显。
图10是弯腰状态下原始心电信号以及本发明和五种对比方法实现运动伪迹抑制后的波形对比图。如图10所示,LMS与CCPC算法处理后心电信号含有噪声最多,而EKF与DPF算法在300与850采样点附近存在少量噪声,PMF与本发明效果相似均较优,从R波幅值看本发明恢复效果更优
图11是弯腰状态下静态心电与本发明和五种对比方法滤波前后的心电频谱对比图。如图11所示,从弯腰状态下各方法滤波后频谱可以看出,看PMF与本发明改善效果较明显,更接近静态心电信号。
表1是三种运动状态下本发明和五种对比方法滤波前后的信号指标对比表。
表1
如表1所示,步行与弯腰状态下由于身体整体幅度变化较大,MA噪声较强,各算法滤波后信号质量提升较大,本发明相较于对比方法信噪比SNR提升>12%,均方根误差RMSE提升>9%;扩胸状态下,由于运动数据与MA相关性相比其它两种动作更好,因此本算法效果提升最佳,本发明相较于对比方法信噪比SNR提升>20%,均方根误差RMSE提升>18%,而本发明相关系数与PMF算法近似,较其它四种对比方法平均提升>2%。
综上所述,本发明在抑制心电信号的运动伪迹噪声方面是有效的,在部分运动状态下具有较好技术优势。
尽管上面对本发明说明性的具体实施方式进行了描述,以便于本技术领域的技术人员理解本发明,但应该清楚,本发明不限于具体实施方式的范围,对本技术领域的普通技术人员来讲,只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内,这些变化是显而易见的,一切利用本发明构思的发明创造均在保护之列。
Claims (8)
1.一种基于VMD分解和LET模型的粒子滤波运动伪迹抑制方法,其特征在于,包括以下步骤:
S1:对含有运动伪迹噪声的采集心电信号Xd(n)进行R波定位,提取心电信号周期,n=1,2,…表示采样点序号;
S2:将Q个完整心电信号周期的静态心电信号Xs(n)与采集心电信号Xd(n)进行拼接,得到组合心电信号X(n);
S3:采用VMD算法对组合心电信号X(n)进行分解得到K个模态分量IMFk(n),k=1,2,…,K,K表示模态总数;对于第1层至第K层模态分量,将中心频率和预设频率阈值最为接近的模态分量作为参考分量,记其序号为k*,将第1层至第k*-1层作为含噪分量;
S6:采用LET模型对含噪分量进行重构,具体方法如下:
S7:将步骤S6中拟合还原的采集心电分量IMF′d,k″(n)和其余模态分量叠加,得到重构的采集心电信号X′d(n):
S8:构建如下心电动态模型:
对采集心电信号X′d(t)使用最小二乘法估计得到心电动态模型的关键参数,包括采集心电信号中P/Q/R/S/T五部分波形对应高斯项的振幅ai、宽度bi、中心参数θi,i∈{P,Q,R,S,T};
S9:根据步骤S8确定的心电动态模型和重构采集心电信号X′d(n),采用粒子滤波算法对采集心电信号Xd(n)进行粒子滤波,得到运动伪迹抑制后的采集心电信号X* d(n)。
2.根据权利要求1所述的粒子滤波运动伪迹抑制方法,其特征在于,所述步骤S2中Q的取值范围为[3,5]。
3.根据权利要求1所述的粒子滤波运动伪迹抑制方法,其特征在于,所述步骤S3中频率阈值为20Hz。
4.根据权利要求1所述的粒子滤波运动伪迹抑制方法,其特征在于,所述步骤S3中VMD算法的模态总数K,惩罚系数α和收敛容差ε采用如下方法确定:
收敛容差ε根据经验设置,然后固定惩罚系数α,在预设的模态总数取值范围内对各个取值进行搜索,计算相邻模态总数的能量差值,选取其中能量差值最大的模态总数取值作为最终使用的模态总数K,能量差值的计算公式如下:
其中,EK、EK-1分别表示模态总数取值为K和K-1时的所有模态分量的能量总和;
5.根据权利要求1所述的粒子滤波运动伪迹抑制方法,其特征在于,所述步骤S9中粒子滤波的具体方法如下:
S9.1:构建状态空间模型,包括状态转移方程和量测方程,其中状态转移方程采用心电动态模型,量测方程表示式如下:
S9.2:在对采集心电信号Xd(n)进行采集的同时,同步采集包括三轴加速度、加速度混合数据与三轴角速度,对于每个运动数据采用预设方法进行预处理后,再与采集心电信号Xd(t)进行采样速率同步,记得到的第p路运动数据为Sp(n),p=1,2,…,7;
S9.3:分别计算每路运动数据Sp(n)和采集心电信号Xd(n)之间的相关性,选择相关性最大的运动数据进行归一化后作为指导运动数据S*(n);
S9.4:采用如下公式计算得到指导参数数据γ(n):
γ(n)=τ+sigmoid(abs(ρS*(n)))
其中,sigmoid()函数作为约束算子,τ表示预设的正常数,ρ表示预设的运动数据调整系数;
其中,上标T表示转置,上标-1表示求逆,表示根据第g个粒子在上一采样点n-1的相位/>由状态转移方程计算得到的当采样点n该粒子的相位估计值,/>表示根据第g个粒子在上一采样点n-1的幅值/>由状态转移方程计算得到的当前采样点n该粒子的幅值估计值,当n=0,/> 表示采集心电信号Xd(n)在当前采样点n的相位,R表示重构采集心电信号Xd′(n)协方差矩阵,其表达式为:
S9.8:评估是否出现粒子匮乏现象,即判断是否满足以下公式:
其中,λ表示预设的阈值,如果满足,则判断出现粒子匮乏现象,进入步骤S9.9,否则不作任何操作,直接进入步骤S9.10;
S9.11:令n=n+1,返回步骤S9.6。
8.根据权利要求5所述的粒子滤波运动伪迹抑制方法,其特征在于,所述步骤S9.8中阈值λ=2G/3。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210479424.5A CN114869294B (zh) | 2022-05-05 | 2022-05-05 | 基于vmd分解和let模型的粒子滤波运动伪迹抑制方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210479424.5A CN114869294B (zh) | 2022-05-05 | 2022-05-05 | 基于vmd分解和let模型的粒子滤波运动伪迹抑制方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114869294A CN114869294A (zh) | 2022-08-09 |
CN114869294B true CN114869294B (zh) | 2023-05-30 |
Family
ID=82673555
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210479424.5A Active CN114869294B (zh) | 2022-05-05 | 2022-05-05 | 基于vmd分解和let模型的粒子滤波运动伪迹抑制方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114869294B (zh) |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110392550A (zh) * | 2016-07-25 | 2019-10-29 | 撒哈拉·凯瑟林·帕奇 | 用声波测量进行辐射束射程验证的系统和方法 |
CN111904406A (zh) * | 2020-08-25 | 2020-11-10 | 上海交通大学 | 一种生理信号运动伪迹抑制装置与方法 |
CN113066502A (zh) * | 2021-03-11 | 2021-07-02 | 电子科技大学 | 基于vmd和多小波的心音分割定位方法 |
Family Cites Families (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103690160B (zh) * | 2013-11-18 | 2015-04-29 | 浙江大学 | 一种基于非高斯时序模型的脑电特征提取和状态识别方法 |
CN104605841A (zh) * | 2014-12-09 | 2015-05-13 | 电子科技大学 | 可穿戴心电信号监测装置及方法 |
CN204909435U (zh) * | 2015-08-04 | 2015-12-30 | 浙江好络维医疗技术有限公司 | 植入式心脏监护系统 |
EP3219356A1 (en) * | 2016-03-14 | 2017-09-20 | Max-Planck-Gesellschaft zur Förderung der Wissenschaften e.V. | Apparatus for applying electric pulses to living myocardial tissue |
CN107951484A (zh) * | 2017-12-01 | 2018-04-24 | 电子科技大学 | 一种可拆卸的抑制运动伪迹织物心电有源干电极 |
CN107981859A (zh) * | 2017-12-05 | 2018-05-04 | 电子科技大学 | 具有运动伪迹抑制功能的心电监测胸带 |
CN111329467A (zh) * | 2018-12-18 | 2020-06-26 | 上海图灵医疗科技有限公司 | 一种基于人工智能的心脏疾病辅助检测方法 |
CN109998524A (zh) * | 2019-03-29 | 2019-07-12 | 山东理工大学 | 一种基于变分模态分解理论和k最近邻算法的心电信号分类方法 |
CN110292374B (zh) * | 2019-05-31 | 2022-05-17 | 辽宁师范大学 | 基于奇异谱分析和变分模态分解的心电信号去基线漂移方法 |
CN110575164B (zh) * | 2019-09-20 | 2022-04-12 | 桂林电子科技大学 | 脑电信号伪迹去除方法及计算机可读存储介质 |
CN114403820A (zh) * | 2022-01-13 | 2022-04-29 | 武汉烽火凯卓科技有限公司 | 一种移动目标的生命体征检测方法及系统 |
-
2022
- 2022-05-05 CN CN202210479424.5A patent/CN114869294B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110392550A (zh) * | 2016-07-25 | 2019-10-29 | 撒哈拉·凯瑟林·帕奇 | 用声波测量进行辐射束射程验证的系统和方法 |
CN111904406A (zh) * | 2020-08-25 | 2020-11-10 | 上海交通大学 | 一种生理信号运动伪迹抑制装置与方法 |
CN113066502A (zh) * | 2021-03-11 | 2021-07-02 | 电子科技大学 | 基于vmd和多小波的心音分割定位方法 |
Also Published As
Publication number | Publication date |
---|---|
CN114869294A (zh) | 2022-08-09 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Rahman et al. | Efficient sign based normalized adaptive filtering techniques for cancelation of artifacts in ECG signals: Application to wireless biotelemetry | |
Sameni et al. | A nonlinear Bayesian filtering framework for ECG denoising | |
Zhang et al. | QRS detection based on multiscale mathematical morphology for wearable ECG devices in body area networks | |
Rahman et al. | An efficient cardiac signal enhancement using time–frequency realization of leaky adaptive noise cancelers for remote health monitoring systems | |
Pant et al. | Compressive sensing of electrocardiogram signals by promoting sparsity on the second-order difference and by using dictionary learning | |
Rahman et al. | Noise cancellation in ECG signals using computationally simplified adaptive filtering techniques: Application to biotelemetry | |
JPH0798346A (ja) | 高分解能スペクトル分析方法 | |
Sharma et al. | ECG denoising using weiner filter and adaptive least mean square algorithm | |
Prashar et al. | Semiautomatic detection of cardiac diseases employing dual tree complex wavelet transform | |
CN114648048B (zh) | 基于变分自编码和PixelCNN模型的心电信号降噪方法 | |
Yadav et al. | Denoising and SNR improvement of ECG signals using wavelet based techniques | |
Popescu et al. | High resolution ECG filtering using adaptive Bayesian wavelet shrinkage | |
Yang et al. | Removal of pulse waveform baseline drift using cubic spline interpolation | |
Elbuni et al. | ECG parameter extraction algorithm using (DWTAE) algorithm | |
Rahman et al. | Adaptive noise removal in the ECG using the block LMS algorithm | |
Qureshi et al. | Performance analysis of adaptive algorithms for removal of low frequency noise from ECG signal | |
CN114869294B (zh) | 基于vmd分解和let模型的粒子滤波运动伪迹抑制方法 | |
Gowri et al. | An efficient variable step size least mean square adaptive algorithm used to enhance the quality of electrocardiogram signal | |
Priyadharsini et al. | An Efficient method for the removal of ECG artifact from measured EEG Signal using PSO algorithm | |
Charleston et al. | Interference cancellation in respiratory sounds via a multiresolution joint time-delay and signal-estimation scheme | |
Olmos et al. | ECG signal compression plus noise filtering with truncated orthogonal expansions | |
Rahul et al. | Baseline correction of ECG using regression estimation method | |
Mou et al. | Noise removal and QRS detection of ECG signal | |
Castaño et al. | Autoregressive models of electrocardiographic signal contaminated with motion artifacts: benchmark for biomedical signal processing studies | |
Chauhan et al. | Optimal choice of thresholding rule for denoising ECG using DWT |
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 |