腦電信號的去噪作業(yè)

1、基于算法的腦電信號的去噪基于算法的腦電信號的去噪 在腦電（在腦電（EEGEEG）信號的采集過程中，不可避免地）信號的采集過程中，不可避免地會受到各種各樣的干擾，且噪聲類型具有很大的不可會受到各種各樣的干擾，且噪聲類型具有很大的不可預(yù)知性，因而采集的信號一般是非平穩(wěn)的，這給去噪預(yù)知性，因而采集的信號一般是非平穩(wěn)的，這給去噪帶來了困難，但采用帶來了困難，但采用RLSRLS濾波器對其濾波能有效地解濾波器對其濾波能有效地解決這個問題，利用前一時刻已獲得的濾波器參數(shù)，自決這個問題，利用前一時刻已獲得的濾波器參數(shù)，自適應(yīng)濾波器自動地調(diào)節(jié)現(xiàn)時刻的濾波器參數(shù)，以適應(yīng)適應(yīng)濾波器自動地調(diào)節(jié)現(xiàn)時刻的濾波器參數(shù)，以適

2、應(yīng)信號和噪聲未知的或隨時間變化的統(tǒng)計(jì)特性，實(shí)現(xiàn)最信號和噪聲未知的或隨時間變化的統(tǒng)計(jì)特性，實(shí)現(xiàn)最優(yōu)濾波，得出濾波數(shù)據(jù)；再進(jìn)行匹配濾波，以提高輸優(yōu)濾波，得出濾波數(shù)據(jù)；再進(jìn)行匹配濾波，以提高輸出的信噪比；最后使用出的信噪比；最后使用Matlab7Matlab7軟件繪圖，直觀地反軟件繪圖，直觀地反映濾波結(jié)果。映濾波結(jié)果。自適應(yīng)濾波器及其算法自適應(yīng)濾波器及其算法 RLS算法算法，即遞歸最小二乘方算法，它是，即遞歸最小二乘方算法，它是FIR維納濾波器維納濾波器的一種時間遞歸算法，并且是嚴(yán)格以最小二乘方準(zhǔn)則為依據(jù)的的一種時間遞歸算法，并且是嚴(yán)格以最小二乘方準(zhǔn)則為依據(jù)的算法。算法。優(yōu)點(diǎn)：優(yōu)點(diǎn)：收斂速度快；在快

3、速信道均衡，實(shí)時系統(tǒng)辨識和時間收斂速度快；在快速信道均衡，實(shí)時系統(tǒng)辨識和時間序列分析中得到廣泛的應(yīng)用。序列分析中得到廣泛的應(yīng)用。缺點(diǎn)：缺點(diǎn)：每次迭代計(jì)算量很大。每次迭代計(jì)算量很大。但在生物醫(yī)學(xué)應(yīng)用中數(shù)據(jù)的計(jì)算量不是很大，所以這種算法很但在生物醫(yī)學(xué)應(yīng)用中數(shù)據(jù)的計(jì)算量不是很大，所以這種算法很容易在小型計(jì)算機(jī)上實(shí)現(xiàn)。容易在小型計(jì)算機(jī)上實(shí)現(xiàn)。自適應(yīng)濾波器及其算法自適應(yīng)濾波器及其算法 RLS算法算法的關(guān)鍵是的關(guān)鍵是用二乘方的時間平均的最小化準(zhǔn)則取代用二乘方的時間平均的最小化準(zhǔn)則取代最小均方準(zhǔn)則，并按時間進(jìn)行迭代計(jì)算最小均方準(zhǔn)則，并按時間進(jìn)行迭代計(jì)算。具體說來，就是要對。具體說來，就是要對初始時刻所有誤差

4、的平方進(jìn)行平均并使其最小化，再按這一準(zhǔn)初始時刻所有誤差的平方進(jìn)行平均并使其最小化，再按這一準(zhǔn)則確定則確定FIR濾波器的權(quán)系數(shù)矢量濾波器的權(quán)系數(shù)矢量w，即依據(jù)的準(zhǔn)則是：，即依據(jù)的準(zhǔn)則是：nkken02min)()(其中：其中：)()()(kykdke自適應(yīng)濾波器及其算法自適應(yīng)濾波器及其算法 基本迭代算法如下：基本迭代算法如下：（1 1）在時刻在時刻n n，已獲得，已獲得（n-1n-1）,T(n-1),T(n-1)和和 d(n)d(n)，x(n) x(n) 并存并存儲在濾波器的延時部件中。儲在濾波器的延時部件中。（2 2）利用式（利用式（1 1）- -式（式（4 4）計(jì)算）計(jì)算T T（n n）,

5、(n), k(n), (n), k(n)和和e(n|n-1)e(n|n-1)。) 1.(.)() 1()() 1()()() 1() 1(1)(1)(TnxnTnxnTnxnxnTnTnTnTT自適應(yīng)濾波器及其算法自適應(yīng)濾波器及其算法 )4).()1()()1|()3.(.)()1()()()1()()2.().1|()()1()(nxnndnnenxnTnxnxnTnknnenknnTT（3）得到濾波器輸出響應(yīng)）得到濾波器輸出響應(yīng)y(n)和誤差和誤差e(n)，即：，即：)()()()()(nyndnenxnyT再進(jìn)行第再進(jìn)行第n次迭代。次迭代。匹配濾波器匹配濾波器 匹配濾波器匹配濾波器-輸出

6、信噪比最大的最佳線性濾波器。輸出信噪比最大的最佳線性濾波器。匹配濾波器的匹配濾波器的著眼點(diǎn)著眼點(diǎn)不是盡可能保持信號不失真，而是提高不是盡可能保持信號不失真，而是提高輸出的信噪比。輸出的信噪比。 輸出信噪比輸出信噪比是指輸出在采樣瞬間信噪比。是指輸出在采樣瞬間信噪比。匹配濾波器原理框架圖匹配濾波器原理框架圖x(t)x(t)n(t)n(t)x(t)x(t)H()H()h(t)h(t)y(t)=sy(t)=s0 0(t)+n(t)+n0 0(t)(t)線性濾波器線性濾波器匹配濾波器匹配濾波器 設(shè)給定輸入信號設(shè)給定輸入信號s(t)s(t)，其頻譜是，其頻譜是S(), S(), 輸入噪聲為白噪輸入噪聲為

7、白噪聲聲n(t), pn(t), pn n()=n()=n0 0/2/2，關(guān)鍵是求出與之相匹配的線性濾波器。，關(guān)鍵是求出與之相匹配的線性濾波器。（1 1）匹配濾波器在抽樣時刻）匹配濾波器在抽樣時刻t t0 0時有最大信噪比，即：時有最大信噪比，即：0020max02|)(|nENtsr式中，式中，E為信號為信號s(t)的總能量，即：的總能量，即：dE20| )( |1為使其物理可實(shí)現(xiàn)，為使其物理可實(shí)現(xiàn)，t0是輸出信號是輸出信號s(t)消失的時刻。消失的時刻。匹配濾波器匹配濾波器 （2 2）確定）確定H H（） 出現(xiàn)出現(xiàn)r r0max0max的條件是線性濾波器的傳輸特性的條件是線性濾波器的傳輸特

8、性H（）滿足如）滿足如下關(guān)系，即：下關(guān)系，即：tjeSKH)(*)( H（）就是最佳線性濾波器傳輸特性，由于線性該傳輸）就是最佳線性濾波器傳輸特性，由于線性該傳輸特性與信號頻譜的復(fù)共軛相一致，故又稱匹配濾波器。特性與信號頻譜的復(fù)共軛相一致，故又稱匹配濾波器。（3）匹配濾波器的單位沖激響應(yīng)）匹配濾波器的單位沖激響應(yīng)h(t)為：為：)()(0ttKsthh(t)是信號是信號s(t)的鏡像信號的鏡像信號s(-t) 在時間上再平移在時間上再平移t0得到的。得到的。匹配濾波器匹配濾波器 （4）匹配濾波器的輸出信號）匹配濾波器的輸出信號s s0 0(t)(t)為：為：)()(00ttKRts即匹配濾波器的

9、輸出信號即匹配濾波器的輸出信號s s0 0(t)(t)是輸入信號是輸入信號s(t)s(t)的自相關(guān)函的自相關(guān)函數(shù)的數(shù)的K K倍，所以匹配濾波器可以看成是一個相關(guān)器。倍，所以匹配濾波器可以看成是一個相關(guān)器。實(shí)驗(yàn)結(jié)果實(shí)驗(yàn)結(jié)果-編程實(shí)現(xiàn)編程實(shí)現(xiàn)RLSRLS濾波濾波 用用C語言變成實(shí)現(xiàn)對源含噪信號的自適應(yīng)濾波處理，并用語言變成實(shí)現(xiàn)對源含噪信號的自適應(yīng)濾波處理，并用Matlab 7繪圖輸出結(jié)果。繪圖輸出結(jié)果。下面為實(shí)現(xiàn)下面為實(shí)現(xiàn)RLS自適應(yīng)濾波的程序：自適應(yīng)濾波的程序：# include# include# define A n # define A n /這里這里n n為輸入數(shù)據(jù)數(shù)為輸入數(shù)據(jù)數(shù)/mai

10、n()main() int i,l=0.98; int i,l=0.98; float e, rA=0,0,pA=0,0, wA=0,0, float e, rA=0,0,pA=0,0, wA=0,0, xA=0,0, yA=0,0, dA=0,0,;xA=0,0, yA=0,0, dA=0,0,;/ xA=0,0,/ xA=0,0,零后面輸入含噪數(shù)據(jù)，零后面輸入含噪數(shù)據(jù)， dA=0,0,dA=0,0,零后面輸零后面輸入期望數(shù)據(jù)，入期望數(shù)據(jù)， yA=0,0,yA=0,0,零后面輸入濾波數(shù)據(jù)零后面輸入濾波數(shù)據(jù) /實(shí)驗(yàn)結(jié)果實(shí)驗(yàn)結(jié)果-編程實(shí)現(xiàn)編程實(shí)現(xiàn)RLSRLS濾波濾波 for(i=2;iA;i+)

11、for(i=2;iA;i+) ri=(xi ri=(xi* *xi+lxi+l* *(xi-1(xi-1* *xi-1)+lxi-1)+l* *l l* *(xi-2(xi-2* *xi-2);xi-2); e=di-wi-1 e=di-wi-1* *xi;xi; wi=wi-1+xi wi=wi-1+xi* *e/(ri);e/(ri); yi=wi yi=wi* *xi;xi; ; ;for(i=2;iA;i+)for(i=2;iA;i+)printf(“%f”,yi);printf(“%f”,yi);getch( );getch( ); 實(shí)驗(yàn)結(jié)果實(shí)驗(yàn)結(jié)果-編程實(shí)現(xiàn)編程實(shí)現(xiàn)RLSRLS濾波

12、濾波下面為下面為Matlab繪圖程序：繪圖程序：x= , x= , /為輸入數(shù)據(jù)為輸入數(shù)據(jù)/t=0:0.01:5.11, t=0:0.01:5.11, /為為512512個數(shù)據(jù)的個數(shù)據(jù)的間隔時間間隔時間/plot(t,x), plot(t,x), /顯示顯示/實(shí)驗(yàn)結(jié)果實(shí)驗(yàn)結(jié)果-實(shí)現(xiàn)功率譜估計(jì)實(shí)現(xiàn)功率譜估計(jì) 下面為實(shí)現(xiàn)功率譜估計(jì)的下面為實(shí)現(xiàn)功率譜估計(jì)的Matlab程序：程序： x= , x= , /為輸入數(shù)據(jù)為輸入數(shù)據(jù)/ p=floor(length(x)/3)+1, p=floor(length(x)/3)+1, nfft=1024, nfft=1024, Fs=1000, Fs=1000,

13、xpsd,f=pburg(x,p,nfft,Fs,half) xpsd,f=pburg(x,p,nfft,Fs,half) pmax=max(xpsd), pmax=max(xpsd), xpsd=xpsd/pmax, xpsd=xpsd/pmax, xpsd=10 xpsd=10* *log10(xpsd+0.000001).log10(xpsd+0.000001). subplot(212), subplot(212), plot(f,xpsd) plot(f,xpsd) tilte( tilte(基于基于BurgBurg算法的功率譜估計(jì)算法的功率譜估計(jì)) ) ylabel( ylabel

14、(功率譜密度功率譜密度(dB)(dB) xlabel( xlabel(頻率頻率(Hz)(Hz) grid on grid on實(shí)驗(yàn)結(jié)果實(shí)驗(yàn)結(jié)果-匹配濾波輸出匹配濾波輸出 以下為對輸入數(shù)據(jù)序列取倒置的以下為對輸入數(shù)據(jù)序列取倒置的C程序：程序：IncludeIncludemain()main()int i=0,j=0;int i=0,j=0; float aA= ,bA=0,; float aA= ,bA=0,; for(i=0;iA;i+) bj=aA-j;j+; for(i=0;iA;i+) bj=aA-j;j+;for(i=0;iA;i+)for(i=0;iA;i+) printf(“%f”,bi); printf(“%f”,bi);getch();getch(); /aA/aA為輸入數(shù)據(jù)；為輸入數(shù)據(jù)；bAbA輸輸 出數(shù)據(jù)出數(shù)據(jù)/下面為實(shí)現(xiàn)匹配濾波的下面為實(shí)現(xiàn)匹配濾波的MatlabMatlab程序：程序：x= , x= , /輸入數(shù)據(jù)輸入數(shù)據(jù)/y= , y= , /輸入原數(shù)據(jù)的倒置輸入原數(shù)據(jù)的倒置數(shù)據(jù)數(shù)據(jù)/c=conv(x,y), c=conv(x,y), /求卷積求卷積/t=0:0.01:10.23,t=0:0.01:10.23,plot(t,c), plot(t,c),