自適應(yīng)濾波器LMS算法

1、用于消除工頻干擾自適應(yīng)濾波器的設(shè)計(jì)與仿真一、背景及意義腦科學(xué)研究不僅是一項(xiàng)重要的前沿性基礎(chǔ)研究，而且是一項(xiàng)對(duì)人類(lèi)健康有重要實(shí)際意義的應(yīng)用研究。隨著社會(huì)的發(fā)展、人類(lèi)壽命的延長(zhǎng)，因腦衰老、紊亂或損傷而引起的腦疾患，對(duì)社會(huì)財(cái) 富消耗和家庭的負(fù)擔(dān)日益增大。 許多國(guó)家紛紛將腦科學(xué)的研究列 入國(guó)家規(guī)劃，并且制訂長(zhǎng)遠(yuǎn)的研究計(jì)劃。人們把21世紀(jì)看成是腦科學(xué)研究高潮的時(shí)代。在腦電信號(hào)的實(shí)際檢測(cè)過(guò)程中，往往含有心電、眼動(dòng)偽跡、 肌電信號(hào)、50Hz工頻干擾以及其它干擾源所產(chǎn)生的干擾信號(hào)， 這給腦電分析以及腦電圖的臨床應(yīng)用帶來(lái)了很大的困難。因此如何從腦電中提取出有用的信息是非常具有挑戰(zhàn)性，且又很有學(xué)術(shù)價(jià)值、實(shí)用價(jià)值的

2、研究課題。本論文從信號(hào)處理的角度出發(fā)， 采集腦電波，使得在強(qiáng)干擾 背景下的腦電信號(hào)得以提取，還原出干凈的腦電波，用于臨床醫(yī) 學(xué)、家庭保健等。醫(yī)生可以利用所采集到的腦電波來(lái)進(jìn)行對(duì)病人 神經(jīng)松弛訓(xùn)練，通過(guò)腦電生物反饋技術(shù)實(shí)現(xiàn)自我調(diào)節(jié)和自我控 制。運(yùn)用生物反饋療法，就是把求治者體內(nèi)生理機(jī)能用現(xiàn)代電子 儀器予以描記，并轉(zhuǎn)換為聲、光等反饋信號(hào)，因而使其根據(jù)反饋 信號(hào)，學(xué)習(xí)調(diào)節(jié)自己體內(nèi)不遂意的內(nèi)臟機(jī)能及其他軀體機(jī)能、達(dá)到防治身心疾病的目的。這種反饋療法是在一定程度上發(fā)掘人體 潛能的一種人一機(jī)反饋方法。 有研究表明腦電生物反饋對(duì)多種神經(jīng)功能失調(diào)疾病有明顯療效。 對(duì)于有腦障礙或腦疾病的人，也可以隨時(shí)監(jiān)測(cè)其腦電

3、信號(hào)，及早地發(fā)現(xiàn)問(wèn)題，避免不必要的損失。二、腦電數(shù)字信號(hào)處理的研究現(xiàn)狀腦電的監(jiān)護(hù)設(shè)備在國(guó)內(nèi)外品種繁多，高新技術(shù)含量高，技術(shù)附加值高，相比而言，我國(guó)的產(chǎn)品較國(guó)際高水平產(chǎn)品落后10-15年。但近年來(lái)，國(guó)內(nèi)產(chǎn)品也逐步利用高新技術(shù)使產(chǎn)品向自動(dòng)化、 智能化、小型化、產(chǎn)品結(jié)構(gòu)模塊化方向發(fā)展。國(guó)內(nèi)產(chǎn)品在抗干擾、 數(shù)字處理、實(shí)時(shí)傳輸數(shù)據(jù)等方面已有很大進(jìn)展，使腦電檢測(cè)不再是只能在屏蔽室進(jìn)行。目前，腦電信號(hào)的數(shù)字濾波從原理上來(lái)看，主要有FIR濾波器和IIR濾波器。FIR濾波器可以提供線性濾波，但存在階數(shù)較 高，運(yùn)算較為復(fù)雜的缺點(diǎn)11;而IIR濾波器是一種非線性濾波 器，它可以用較少的階數(shù)實(shí)現(xiàn)性能良好的濾波，是目前

4、運(yùn)用較 廣泛的一種濾波器10。然而無(wú)論FIR還是IIR濾波器，都是固定 濾波器，其靈活性較差。而本文 提出的自適應(yīng)濾波器能提供易于控制的帶寬、極其深的零點(diǎn)、 以及具有自適應(yīng)地跟蹤干擾頻率和相位的能力1213。三、研究的主要內(nèi)容 本論文研究的主要內(nèi)容有：（1）熟悉腦電信號(hào)基本概念。（2）提出將自適應(yīng)濾波器運(yùn)用于腦電信號(hào)的處理中。（3）對(duì)自適應(yīng)濾波算法進(jìn)行理論分析，研究其規(guī)律，通過(guò)軟件 仿真，驗(yàn)證算法的可靠性。以 50 Hz陷波器為例。(4) 輸入腦電信號(hào)數(shù)據(jù)，觀察濾除 50 Hz工頻干擾之后的腦電 信號(hào)。四、腦電信號(hào)采集原理4.1、腦電圖的一般性質(zhì)人們利用電極把腦細(xì)胞電活動(dòng)的電位與時(shí)間的關(guān)系記錄

5、下來(lái) 的波形稱(chēng)為腦電圖，如圖2.1所示，頭皮上二點(diǎn)之間或頭皮與無(wú) 關(guān)電極(即參考電極)之間的電位差可以用E+A E來(lái)表示，這里E為直流電位，普遍認(rèn)為是由電極極化電壓引起的；E為電位的變化量，人們所觀察到的腦電圖就是厶E的波形，這種電現(xiàn)象是生物電現(xiàn)象的一種表現(xiàn)。 生物一旦死亡，生物電現(xiàn)象隨即很快 消失，不同種類(lèi)的動(dòng)物其腦電活動(dòng)亦各不相同18。圖2.1腦電圖Fig.2.1 Electroe ncephalogram4.2、腦電測(cè)量中的偽差腦電圖的偽差又稱(chēng)偽跡或干擾，是指來(lái)自腦外的電位活動(dòng)在 腦電圖中的反映。偽差的出現(xiàn)常給閱讀、分析、判斷腦電圖造成 困難，尤其是某些偽差與癲癇波很相似， 臨床上很容易

6、造成誤診， 因此正確識(shí)別和排除偽差是很重要的。引起偽差的因素很多，表現(xiàn)也多種多樣，主要有：（1）工頻干擾人體所處環(huán)境有大量的電器設(shè)備，這些設(shè)備周?chē)臻g的電磁 場(chǎng)的頻率主要是50Hz。由于人體自身的分布電容， 所以工頻信號(hào) 極易引入人體，而且強(qiáng)度遠(yuǎn)比腦電信號(hào)大。工頻干擾包括50 Hz正弦信號(hào)及其諧波，在一定的環(huán)境下，頻率相對(duì)穩(wěn)定，但有微小 的波動(dòng)。（2）肌電干擾肌電干擾來(lái)自于人體的肌肉顫抖，人體活動(dòng)頻率范圍大，肌 肉的活動(dòng)可以產(chǎn)生毫伏級(jí)的電勢(shì)，它很可能淹沒(méi)掉有用的腦電信號(hào)。（3）電極接觸噪聲電極接觸噪聲是瞬時(shí)干擾信號(hào)，來(lái)源于電極與肌膚的不良接 觸，如病人運(yùn)動(dòng)和抖動(dòng)造成的檢測(cè)電極接觸不好。電極接觸

7、噪聲是隨機(jī)的階躍信號(hào)，在檢測(cè)的過(guò)程中可能會(huì)多次發(fā)生。這種信號(hào)幅值大，頻率低。（4）人為運(yùn)動(dòng)人為運(yùn)動(dòng)造成電極與皮膚接觸阻抗變大，容易造成基線漂移。五、自適應(yīng)濾波器原理由Glover和Widrow等人2425提出的自適應(yīng)陷波器主要用 于濾除固定頻率的干擾信號(hào)。因此數(shù)字部分的50 Hz陷波濾波器均選用自適應(yīng)干擾對(duì)消器。首先研究一下自適應(yīng)50 Hz陷波器的原理。如圖3.6所示，這是單一頻率的自適應(yīng)陷波濾波器，它等 效于有一個(gè)復(fù)權(quán)的噪聲對(duì)消系統(tǒng)，即用兩個(gè)實(shí)權(quán)達(dá)到同時(shí)調(diào)整單一頻率正弦波的幅度和相位，以達(dá)到消除干擾的目的。 假定原始輸入信號(hào)的類(lèi)型是任意的，即可以是隨機(jī)的、確知的、連續(xù)的或 瞬態(tài)的，也可以是各

8、種類(lèi)型的組合。而參考輸入則是頻率為f0的純正弦波。六、自適應(yīng)濾波器的設(shè)計(jì)七、自適應(yīng)濾波器對(duì)腦電信號(hào)的處理1傳輸函數(shù)零極點(diǎn)位置-1-0.50Real Part0.518642 0 2 4 6 8 - - - - trap vdnraa ml濾波器特性5 o-50 5 0112004XISeoaeHr esanp5050-50100150200250300350400450500Frequency (Hz)o0050100150200250300350400450500Frequency (Hz)0.50-0.5-11002003004005006007008009001000有噪腦電信號(hào)時(shí)域波形

9、204060801001201400%有噪信號(hào)%濾波器輸出值%誤差值%系數(shù)調(diào)整3020101、濾波器軟件設(shè)計(jì)程序：Ts=0.004;%采樣間隔n=4000;%采樣點(diǎn)數(shù)N=10;%濾波器抽頭數(shù)u=0.03125;%調(diào)整步長(zhǎng)C=0.4484;%幅度f(wàn)0=50;%參考輸入頻率 k=1: n;x仁C*cos(2*pi*f0*k*Ts);% 參考輸入x2=C*si n( 2*pi*f0*k*Ts);% 參考輸入d(k)=0.4*x1(k)+(10A(-3)*data(k+7000);y=zeros(1, n);%輸出初始值w仁zeros(1,N);%濾波系數(shù)初始值w2=zeros(1,N);forn=1

10、: ny( n)=w1( n)*x1( n)+w2( n)*x2( n); e( n)=d( n)-y( n);w1( n+1)=w1( n)+2*u.*e( n) .*x1( n);w2( n+1)=w2( n)+2*u.*e (n ).*x2( n);figure(1);subplot(2,1,1);plot(k,d);axis(O 1000 -1 1);xlabel(有噪腦電信號(hào)時(shí)域波形);ylabel(幅度);subplot(2,1,2);plot(k,e);axis(0 1000 -1 1);xlabel(濾去50Hz工頻信號(hào)時(shí)域波形);ylabel(幅度);figure(2);su

11、bplot(2,1,1);D=fft(d,4096);f=250*(0:2047)/4096;plot(f,abs(D(1:2048);axis(0 100 0 200);xlabel(有噪腦電信號(hào)頻域波形);ylabel(幅度);subplot(2,1,2);E=fft(e,4096);f=250*(0:2047)/4096;plot(f,abs(E(1:2048);axis(0 100 0 200);xlabel(濾去50Hz工頻信號(hào)頻域波形);ylabel(幅度);figure(3);D=fft(d,4096);pD=D.*conj(D)/4096;%計(jì)算功率譜密度subplot(2,1

12、,1);f=250*(0:2047)/4096;plot(f,pD(1:2048);xlabel(有噪腦電信號(hào)功率譜波形)；ylabel( (dB/Hz);E=fft(e,4096);pE=E.*conj(E)/4096;%計(jì)算功率譜密度subplot(2,1,2);plot(f,pE(1:2048);xlabel(濾去50Hz信號(hào)功率譜波形);ylabel( (dB/Hz);2、濾波器特性及其零極點(diǎn)位置顯示程序：u=0.03125;%調(diào)整步長(zhǎng)C=0.4484;%幅度wO=O.1*pi;b=1, -2*cos(w0), 1;a=1, -2*cos(w0)+2*u*cos(w0)*(C.A2),1-2*u*(C.A2);figu