濾除ECG中肌電和寬頻率范圍工頻干擾的小波算法.doc

第3期趙捷,等.濾除ECG中肌電和寬頻率范圍工頻干擾的小波算法225通道兩種。在單通道中的自適應(yīng)瀘波器可使用頻率估計(jì),來(lái)跟蹤在一定范圍內(nèi)跟蹤工頻干擾的頻率變化5。在多通道中常用的方法是Widrow提出的自適應(yīng)噪音濾波方法,該方法的前提是假設(shè)心電信號(hào)和干擾信號(hào)不相關(guān),在此基礎(chǔ)上利用自適應(yīng)濾波方法,自動(dòng)調(diào)整濾波器系數(shù),以跟蹤輸入過(guò)程,實(shí)現(xiàn)工頻干擾的抵消6。但自適應(yīng)瀘波器的缺點(diǎn)是頻率跟蹤范圍較窄;3):2dentcomponent)波器雖然能通過(guò)平滑抑制噪音,但同時(shí)也會(huì)使信號(hào)的邊沿模糊。高通濾波器可以使邊沿更加陡峭,但背景噪音同時(shí)被加強(qiáng)。小波變換(wavelettransform,WT)是一種時(shí)頻局部化、多分辨分析方法,具有自動(dòng)“變焦距”的功能。以WT為理論基,能,因此在臨床。多resolutionanalysis)又稱多尺度分析。它的基本思想是利用正交小波基函數(shù)的多尺度特性將信號(hào)在不同尺度下展開并加以比較,以得到有用的信息。多分辨率分析不但為正交小波基函數(shù)的構(gòu)造提供了一種簡(jiǎn)單的方法,而且還為正交小波變換的快速算法(Mallat算法)提供了理論依據(jù)。離散二進(jìn)小波變換的Mallat算法為:mm-1(n-2m-1k)x(n)=h0(k)xk相互統(tǒng)計(jì)獨(dú)立的信源經(jīng)線性組合而產(chǎn)生的混合信號(hào),最終從混合信號(hào)中提出各自的獨(dú)立的信號(hào)分量,使用ICA可從多通道心電信號(hào)中分出工頻干擾信號(hào)7。ICA的缺點(diǎn)是使用多通道,不適合僅有一兩個(gè)通道的各類心電監(jiān)護(hù)儀使用。另外工頻干擾其中還包含豐富的諧波分量,以及因電網(wǎng)不穩(wěn)定造成的其它噪聲干擾。因此對(duì)工頻干擾需要設(shè)計(jì)多吸收點(diǎn),以消除工頻干擾的基頻和諧波成分。上面3種算法都不能很好地滿足對(duì)諧波成分的吸收。心電信號(hào)中肌電干擾主要由骨骼肌產(chǎn)生。它的頻率范圍很寬,一般在2500Hz之間,表現(xiàn)為不規(guī)則的快速變化波形。由于肌電的頻譜較寬,而且又與心電的頻譜混在一起,故很難用一般常規(guī)的方法將其與心電分開。常規(guī)的心電圖及一般使用轉(zhuǎn)折頻率為37Hz的低通模擬濾波器,這樣在濾除部分肌電干擾的同時(shí),也使QRS波群的高頻成分受到了衰減,表現(xiàn)為R波幅度的降低。mm-1(n-2m-1k)d(n)=h1(k)xk其中m,nZ,m:為尺度,xm(n)是分辨率為2m的信號(hào),稱為近似(approximations)。x0(n)為被處理的原始數(shù)字信號(hào)x(n),dm(n)為信號(hào)m-1(n)的二進(jìn)小波變換,稱為細(xì)節(jié)(details)。x01h(n)和h(n)為有限沖擊響應(yīng)濾波器,稱為分解濾波器。信號(hào)xm-1(n)可由如下Mallat重構(gòu)算法來(lái)恢復(fù):xm-1(n)=g0(k)xm(n-2k)+g1(k)kkd(n-2k)m這里g0(n)和g1(n)為有限沖擊響應(yīng)濾波器,稱為重構(gòu)濾波器。這里選取雙正交小波(biorthogonal)中的小波函數(shù)(bior2.2)。雙正交小波的主要特征是具有線性相位、對(duì)稱性和正規(guī)性,這點(diǎn)在心電處理中尤為重要,因?yàn)橄辔皇д鎸?dǎo)致心電波形的失真。小波函數(shù)(bior2.2)分解濾波器系數(shù)如下:h0(-2,-1,0,1,2)小波除噪算法從廣義上來(lái)說(shuō),含干擾的心電信號(hào)可表示為x(n)=s(n)+v(n),其中:x(n)是實(shí)測(cè)的心電信號(hào),s(n)是不含噪聲的心電信號(hào),v(n)是各種干擾信號(hào)總和,如50/60Hz工頻干擾、肌電干擾等。x(n)具有較強(qiáng)的隨機(jī)性和背景噪聲,而且屬于非線性、非平穩(wěn)的微弱信號(hào)。此類信號(hào)如果采用經(jīng)典頻譜分析方法,則即難以有效去噪,又難以揭示出信號(hào)的動(dòng)態(tài)變化過(guò)程。傳統(tǒng)的建立在付氏變換基礎(chǔ)上的濾波方法在提高信噪比和提高分辨率之間存在矛盾。低通濾=0,0,0.3536,0.7071,0.3536;h1(-2,-1,0,1,2)=0.1768,0.3536,-1.0607,0.3536,0.1768。小波函數(shù)(bior2.2)重構(gòu)濾波器系數(shù)如下:g0(-2,-1,0,1,2)=-0.1768,0.3536,1.0607,0.3536,-0.1768;第2/5頁(yè) 226航天醫(yī)學(xué)與醫(yī)學(xué)工程第17卷g1(-2,-1,0,1,2)=0,0.3536,-0.7071,0.3536,0。根據(jù)上述分析可知,基于小波變換分析的多分辨分析即相當(dāng)于對(duì)信號(hào)進(jìn)行低通和高通濾波,可將信號(hào)分解為位于不同頻帶和時(shí)段內(nèi)的各個(gè)成分。因此,通過(guò)Mallat算法將信號(hào)分解后,就可根據(jù)先驗(yàn)知識(shí),引入門限來(lái)作為甄別受到噪聲污染的小波系數(shù)。由噪聲產(chǎn)生,小波系數(shù),再由Mallat來(lái)重建信號(hào),從而即獲得濾除噪聲后的信號(hào),又不致于引起重建結(jié)果的明顯失真。這就是非線性小波方法用于從噪聲中恢復(fù)信號(hào)的實(shí)質(zhì)。要用小波方法很好地實(shí)現(xiàn)信噪分離,關(guān)鍵的問(wèn)題是如何設(shè)計(jì)出好的門限8。設(shè)Wf(a,b)為信號(hào)f(t)的小波變換。取任意點(diǎn)(a0,b0)使得當(dāng)b屬于b0的左鄰域或右鄰域時(shí),有成立,則稱Wf(a0,b0)為小波變換的極大模。極大模是由信號(hào)中奇異點(diǎn)及噪聲產(chǎn)生的,并且信號(hào)和噪聲的二進(jìn)小波變換的極大值的形態(tài)將隨著分辨級(jí)數(shù)的不斷增加表現(xiàn)出不同的狀態(tài)。由信號(hào)的奇異點(diǎn)產(chǎn)生的的模極大值具有沿尺度傳遞的性質(zhì)。即在各個(gè)尺度上,信號(hào)的奇異點(diǎn)附近都具有模極大值點(diǎn),而且這些模極大值點(diǎn)的符號(hào)保持不變,而且隨著分辨級(jí)數(shù)的增加而緩慢增加。由噪聲產(chǎn)生的極大模,卻隨著分辨級(jí)數(shù)的增加而迅速減小9因此,換,。如圖1(bior2.2),采樣頻率為200Hz。其中:d3,d2,d1是尺度為3,2,1的小波變換(細(xì)節(jié)),a3是尺度為3的近似;s是含噪音的心電信號(hào)。由噪聲產(chǎn)生的極大模,卻隨著分辨級(jí)數(shù)的增加而迅速減小,在d3上基本上只有心電信號(hào)中奇異點(diǎn)產(chǎn)生的極大模,噪聲產(chǎn)生的極大模很小。而在d3上極大模產(chǎn)生處,對(duì)應(yīng)與心電信號(hào)中QRS波群。我們?cè)谇懊嬷赋?QRS波群含有的最高頻率成分較高(50Hz),而T波和P波含有的最高頻率成分只有十幾赫茲,根據(jù)這一特點(diǎn),設(shè)計(jì)了小波去除噪聲方法。方法如下:先取23s的尺度為3,2,1小波變換(至少包含一個(gè)心跳周期),令:d3=d3(n),d2=d2(n),d1=d1(n),a3=a3(n),d3max=max|d3(n)|。在一個(gè)心跳周期以內(nèi)可以近似圖1含噪音的心電信號(hào)小波分解Fig.1ThewaveletsofECGcontaminatedbynoise第3/5頁(yè) 第3期趙捷,等.濾除ECG中肌電和寬頻率范圍工頻干擾的小波算法22認(rèn)為d1(n)的方差為:d1=Ed1-Ed1,由于小波變換相當(dāng)于通過(guò)一帶通濾波器,故:Ed1=0,22所以:d1=EEd1同理d2和d3的方差為:2222d2=EEd2,d3=EEd3227由于d3(n)上的干擾較小(圖1),小波變換的極大值在QRS波群中,從有利于防止干擾及防止漏掉小幅度的QRS波群的角度,取門限thqrs=d3max/3。該門限值適應(yīng)于噪聲的小波分解在3上的最大值小于dd3(nmd1(n)=先定位一段在80ms內(nèi)|d3(n)|thqrs,即找到QRS波群的起點(diǎn),記為QRSs,找到起點(diǎn)后。根據(jù)不應(yīng)期理論,再后跳200ms,此點(diǎn)向前搜索d3(n),如果當(dāng)|d3(n)|thqrs,即找到QRS,記為,thd13,d2=2d2,thd3=(n)-thd1當(dāng)d1(n)thd1與(QRSsnQRSe)d1(n)-thd1當(dāng)d1(n)-thd1與(QRSsnQRSe)0,(其它n)d2(n)-thd2當(dāng)d2(n)thd2與(QRSsnQRSe)md2(n)=d2(n)-thd2當(dāng)d2(n)-thd2與(QRSsnQRSe)md3(n)=0,(其它n)d3(n),當(dāng)(|d3(n)|thd3)和61Hz工頻干擾已消除。ST段、T波、P