滤除ECG中肌电和宽频率范围工频干扰的小波算法

1、第3期赵捷,等.滤除ECG中肌电和宽频率范围工频干扰的小波算法225通道两种。在单通道中的自适应泸波器可使用频率估计,来跟踪在一定范围内跟踪工频干扰的频率变化5。在多通道中常用的方法是Widrow提出的自适应噪音滤波方法,该方法的前提是假设心电信号和干扰信号不相关,在此基础上利用自适应滤波方法,自动调整滤波器系数,以跟踪输入过程,实现工频干扰的抵消6。但自适应泸波器的缺点是频率跟踪范围较窄;3):2dentcomponent)波器虽然能通过平滑抑制噪音,但同时也会使信号的边沿模糊。高通滤波器可以使边沿更加陡峭,但背景噪音同时被加强。小波变换(wavelettransform,WT)是一种时频局

2、部化、多分辨分析方法,具有自动“变焦距”的功能。以WT为理论基,能,因此在临床。多resolutionanalysis)又称多尺度分析。它的基本思想是利用正交小波基函数的多尺度特性将信号在不同尺度下展开并加以比较,以得到有用的信息。多分辨率分析不但为正交小波基函数的构造提供了一种简单的方法,而且还为正交小波变换的快速算法(Mallat算法)提供了理论依据。离散二进小波变换的Mallat算法为:mm-1(n-2m-1k)x(n)=h0(k)xk相互统计独立的信源经线性组合而产生的混合信号,最终从混合信号中提出各自的独立的信号分量,使用ICA可从多通道心电信号中分出工频干扰信号7。ICA的缺点是使

3、用多通道,不适合仅有一两个通道的各类心电监护仪使用。另外工频干扰其中还包含丰富的谐波分量,以及因电网不稳定造成的其它噪声干扰。因此对工频干扰需要设计多吸收点,以消除工频干扰的基频和谐波成分。上面3种算法都不能很好地满足对谐波成分的吸收。心电信号中肌电干扰主要由骨骼肌产生。它的频率范围很宽,一般在2500Hz之间,表现为不规则的快速变化波形。由于肌电的频谱较宽,而且又与心电的频谱混在一起,故很难用一般常规的方法将其与心电分开。常规的心电图及一般使用转折频率为37Hz的低通模拟滤波器,这样在滤除部分肌电干扰的同时,也使QRS波群的高频成分受到了衰减,表现为R波幅度的降低。mm-1(n-2m-1k)

4、d(n)=h1(k)xk其中m,nZ,m:为尺度,xm(n)是分辨率为2m的信号,称为近似(approximations)。x0(n)为被处理的原始数字信号x(n),dm(n)为信号m-1(n)的二进小波变换,称为细节(details)。x01h(n)和h(n)为有限冲击响应滤波器,称为分解滤波器。信号xm-1(n)可由如下Mallat重构算法来恢复:xm-1(n)=g0(k)xm(n-2k)+g1(k)kkd(n-2k)m这里g0(n)和g1(n)为有限冲击响应滤波器,称为重构滤波器。这里选取双正交小波(biorthogonal)中的小波函数(bior2.2)。双正交小波的主要特征是具有线性

5、相位、对称性和正规性,这点在心电处理中尤为重要,因为相位失真将导致心电波形的失真。小波函数(bior2.2)分解滤波器系数如下:h0(-2,-1,0,1,2)小波除噪算法从广义上来说,含干扰的心电信号可表示为x(n)=s(n)+v(n),其中:x(n)是实测的心电信号,s(n)是不含噪声的心电信号,v(n)是各种干扰信号总和,如50/60Hz工频干扰、肌电干扰等。x(n)具有较强的随机性和背景噪声,而且属于非线性、非平稳的微弱信号。此类信号如果采用经典频谱分析方法,则即难以有效去噪,又难以揭示出信号的动态变化过程。传统的建立在付氏变换基础上的滤波方法在提高信噪比和提高分辨率之间存在矛盾。低通滤

6、=0,0,0.3536,0.7071,0.3536;h1(-2,-1,0,1,2)=0.1768,0.3536,-1.0607,0.3536,0.1768。小波函数(bior2.2)重构滤波器系数如下:g0(-2,-1,0,1,2)=-0.1768,0.3536,1.0607,0.3536,-0.1768;第2/5页 226航天医学与医学工程第17卷g1(-2,-1,0,1,2)=0,0.3536,-0.7071,0.3536,0。根据上述分析可知,基于小波变换分析的多分辨分析即相当于对信号进行低通和高通滤波,可将信号分解为位于不同频带和时段内的各个成分。因此,通过Mallat算法将信号分解后

7、,就可根据先验知识,引入门限来作为甄别受到噪声污染的小波系数。由噪声产生,小波系数,再由Mallat来重建信号,从而即获得滤除噪声后的信号,又不致于引起重建结果的明显失真。这就是非线性小波方法用于从噪声中恢复信号的实质。要用小波方法很好地实现信噪分离,关键的问题是如何设计出好的门限8。设Wf(a,b)为信号f(t)的小波变换。取任意点(a0,b0)使得当b属于b0的左邻域或右邻域时,有成立,则称Wf(a0,b0)为小波变换的极大模。极大模是由信号中奇异点及噪声产生的,并且信号和噪声的二进小波变换的极大值的形态将随着分辨级数的不断增加表现出不同的状态。由信号的奇异点产生的的模极大值具有沿尺度传递

8、的性质。即在各个尺度上,信号的奇异点附近都具有模极大值点,而且这些模极大值点的符号保持不变,而且随着分辨级数的增加而缓慢增加。由噪声产生的极大模,却随着分辨级数的增加而迅速减小9因此,换,。如图1(bior2.2),采样频率为200Hz。其中:d3,d2,d1是尺度为3,2,1的小波变换(细节),a3是尺度为3的近似;s是含噪音的心电信号。由噪声产生的极大模,却随着分辨级数的增加而迅速减小,在d3上基本上只有心电信号中奇异点产生的极大模,噪声产生的极大模很小。而在d3上极大模产生处,对应与心电信号中QRS波群。我们在前面指出:QRS波群含有的最高频率成分较高(>50Hz),而T波和P波含

9、有的最高频率成分只有十几赫兹,根据这一特点,设计了小波去除噪声方法。方法如下:先取23s的尺度为3,2,1小波变换(至少包含一个心跳周期),令:d3=d3(n),d2=d2(n),d1=d1(n),a3=a3(n),d3max=max|d3(n)|。在一个心跳周期以内可以近似图1含噪音的心电信号小波分解Fig.1ThewaveletsofECGcontaminatedbynoise第3/5页 第3期赵捷,等.滤除ECG中肌电和宽频率范围工频干扰的小波算法22认为d1(n)的方差为:d1=Ed1-Ed1,由于小波变换相当于通过一带通滤波器,故:Ed1=0,22所以:d1=EEd1同理d2和d3的

10、方差为:2222d2=EEd2,d3=EEd3227由于d3(n)上的干扰较小(图1),小波变换的极大值在QRS波群中,从有利于防止干扰及防止漏掉小幅度的QRS波群的角度,取门限thqrs=d3max/3。该门限值适应于噪声的小波分解在3上的最大值小于dd3(nmd1(n)=先定位一段在80ms内|d3(n)|<thqrs的电作为起始点,由此点向后搜索d3(n),如果当|d3(n)|>thqrs,即找到QRS波群的起点,记为QRSs,找到起点后。根据不应期理论,再后跳200ms,此点向前搜索d3(n),如果当|d3(n)|>thqrs,即找到QRS,记为,thd13,d2=2

11、d2,thd3=(n)-thd1当d1(n)thd1与(QRSsnQRSe)d1(n)-thd1当d1(n)-thd1与(QRSsnQRSe)0,(其它n)d2(n)-thd2当d2(n)thd2与(QRSsnQRSe)md2(n)=d2(n)-thd2当d2(n)-thd2与(QRSsnQRSe)md3(n)=0,(其它n)d3(n),当(|d3(n)|thd3)和61Hz工频干扰已消除。ST段、T波、P波段的0当(|d3(n)|<thd3)肌电信号也已基本消除。重构系数的md1,md2,md3代替原有的d1,d2,d3与a3重新构成消噪的心电。由于d1(n)和d2(n)的信噪(SNR

12、)较低,故采用SoftThresh2olding法进行去噪,以较好地平滑信号。而d3(n)的SNR较高,故采用HardThresholding法进行去噪,以较好地保持信号的细节。d3max应逐跳更新。上面小波消噪方法,可以看成一个时变低通滤波器。该滤波器对于QRS波群的截止频率较高,而对T和P波的截止频率较低。这样在保留QRS波群较高的频率同时,消去了噪音。小波消噪算法验证图2为受肌电干扰的心电信号,选自美国MIT2BIH心律失常数据库106信号,原抽样频率为360Hz,再抽样频率为200Hz,A/D为12bit,分辨率为4.88V。从图可见ST段,T波,P波段有明显的肌电信号。下图为经过上节

13、所述小波消噪算法后输出的心电。可见经小波消噪算法后ST段,T波,P波段的肌电信号已基本消除。图3和4所示为受肌电干扰的心电信号又分别叠加幅度为20%的49Hz和61Hz工频干扰的心电信号。下图为经过上节所述小波消噪算法后输出的心电。可见经小波消噪算法后49Hz图2受肌电干扰的心电信号的消噪Fig.2De2noisedECGsignalcontaminatedbymuscleelec2tricity图3受肌电干扰并叠加幅度为20%的49Hz工频干扰心电信号的消噪Fig.3De2noisedECGsignalcontaminatedbymuscleelec2tricityand20%49Hzpow

14、erline第4/5页 228航天医学与医学工程第17卷化并不敏感,对于50/60Hz的工频干扰可用同样的算法。参考文献1AlsteJAV.ECGbaselinewanderreductionusinglinerphasefilterJ.ComputerandBiomedicalsearch,1986,19:41712BadiliniF,JM,Edwardsplinebaselineesti2inthemeasurementofSTC.In:Annualin2inMedicineand,(2):584258513PS.Acomparisonofadaptiveandnonadaptivefil

15、tersforreductionofpowerlineinterferenceintheECGJ.IEEETransactiononBiomedicalEngineering,1996,43(1):105210914YinDY,YangYX.High2speedreal2timeprocessinginmulti2channelsynchronousECGJ.SpaceMedicine&MedicalEn2gineering,2002,15(1):6426715ZiaraniAK,KonradA.ANonlinearadaptivemethodofelimi2nationofpower

16、lineinterferenceinECGsignalsJ.IEEETransactionsonBiomedicalEngineering,2002,49(6):540254716YangFS,GaoSK.BiomedicalSignalProcessingM.Beijing:HigherEducationPress,1989.7WuXP,FengHQ,ZhouHQ,etal.Removalofpowereinter2ferencefromdigitalsignalsbyusingindependentcomponentanalysisJ.JournalofChinaUniversityofS

17、cienceandTech2nology,2000,30(6):671267618XuL,ZhengXX.AnovelmethodformedicalsignaldenoisingbasedonwavelettransformJ.ComputerApplication,1998,6:3323419LiYH,YangXJ,MiYG.FilteringofsampledatabasedonwavelettransformationJ.JournalofNanjingUniversityofChemicalTechnology,1999,21(2):562581作者简介:赵捷,男,博士,教授,研究方向为信号处理、通信工程图4%的Fig.4De2noisedcontaminatedbymuscleelec2tricityand2