DSP 維納濾波器的計(jì)算機(jī)實(shí)現(xiàn)

1、維納濾波器的計(jì)算機(jī)實(shí)現(xiàn)一、 實(shí)驗(yàn)?zāi)康?、 利用計(jì)算機(jī)編程實(shí)現(xiàn)加性噪聲信號(hào)的維納濾波。2、 將計(jì)算機(jī)模擬實(shí)驗(yàn)結(jié)果與理論分析結(jié)果相比較，分析影響維納濾波效果的各種因素，從而加深對(duì)維納濾波的理解。3、 利用維納濾波一步純預(yù)測(cè)方法實(shí)現(xiàn)對(duì)信號(hào)生成模型的參數(shù)估計(jì)。二、 實(shí)驗(yàn)原理1、 維納濾波器是一種從噪聲中提取信號(hào)的最佳線性估計(jì)方法，假定一個(gè)隨機(jī)信號(hào)形式為：x(n)=s(n)+v(n),其中s(n)為有用信號(hào)，v(n)為噪聲信號(hào)。而維納濾波的作用就是讓x(n)通過(guò)一個(gè)系統(tǒng)h(n)盡可能濾掉噪聲，提取近似s(n),h(n)的選擇以最小均方誤差為準(zhǔn)則。由維納-霍夫方程知，只要求出xx 及xs就可求出h（h=-

2、1xxxs）。但要求h(n)滿足因果性要求，維納-霍夫方程便是一個(gè)難題，這里利用最佳FIR維納濾波方法求解h(n)的近似，這也便于在計(jì)算機(jī)上實(shí)現(xiàn)，公式為：h =R-1xx rxs。實(shí)驗(yàn)中s(n)由信號(hào)生成模型：s(n)=as(n-1)+w(n)確定，其中a=0.95,w(n)是均值為0，方差為w2=1的高斯白噪聲，v(n)為均值為0，方差為1的高斯白噪聲，且s(n)與v(n)不相關(guān)。實(shí)驗(yàn)中s(n)是已知的，但實(shí)際中如果s(n)已知，維納濾波也就失去意義了，因此實(shí)驗(yàn)純粹是為了理解維納濾波原理而設(shè)計(jì)。2、 維納一步純預(yù)測(cè)問(wèn)題S(n)的生成模型：s(n)+a1(n-1)+aps(n-p)=w(n),

3、已知xx（n）,利用Yule-walker方程即可得到信號(hào)生成模型參數(shù)ai(i=1,2p)和2w 。三、 實(shí)驗(yàn)步驟及結(jié)果分析1、 根據(jù)維納濾波原理繪制程序流程圖開始 輸入樣本個(gè)數(shù)L,FIR濾波器階數(shù)N產(chǎn)生L個(gè)v(n),w(n),s(n)和x(n)，利用L個(gè)s(n)和x(n)，估計(jì)RSS和rxs檢驗(yàn)產(chǎn)生序列x(n的自相關(guān)和互相關(guān)函數(shù)是否與理論值相符 N在同一坐標(biāo)內(nèi)繪出x(n)自相關(guān)函數(shù)的理論值和實(shí)際值 Y在同一坐標(biāo)內(nèi)繪出最后100個(gè)s(n)和x(n)。調(diào)矩陣求逆子程序計(jì)算,將N個(gè)理想的h(n)和估計(jì)的h(n) 繪于同一坐標(biāo)內(nèi)進(jìn)行理想的維納濾波得L個(gè)SI (n),和最后100個(gè)s(n)繪制于同一坐

4、標(biāo) 對(duì)x(n)進(jìn)行過(guò)濾得L個(gè)SR(n)，和最后100個(gè)s(n)和繪于同一坐標(biāo)內(nèi)結(jié)束L個(gè)x(n),s(n), SI (n), SR(n),統(tǒng)計(jì)ex2,eI2,eR22、 根據(jù)流程圖編寫程序（見附錄1）并分析運(yùn)行結(jié)果：選擇L=5000,N=10觀察并記錄、分析實(shí)驗(yàn)結(jié)果。1） 與s(n)相比，信號(hào)x(n)在維納濾波前后效果比較： 圖1 圖2圖1為維納濾波后的s(n)與最后100個(gè)s(n)比較圖圖2為未經(jīng)維納濾波的x(n)與最后100個(gè)s(n)比較圖。分析：顯然與s(n)相比，x(n)在維納濾波前與s(n)相差很大，維納濾波后較接近s(n),可見濾波效果比較好。2） 估計(jì)(n)與理想h(n)的比較：

5、圖3圖3為估計(jì)(n)與理想h(n)的對(duì)比圖。分析：由圖可見，二者近似程度除最后幾個(gè)點(diǎn)外，其他近似度還是滿高的，總體而言，近似效果不錯(cuò)。3） 理想的維納濾波與FIR維納濾波效果對(duì)比： 圖 4 圖5圖4為理想維納濾波效果，圖5為FIR維納濾波效果分析：直接從圖形觀察，差異太小，無(wú)法觀察其精度。只能通過(guò)最小均方差來(lái)比較其差異，結(jié)果為：理想維納濾波ei= 0.2287,FIR維納濾波ef=0.2254?？梢?，理想維納濾波效果要好過(guò)FIR維納濾波。4） 自相關(guān)與互相關(guān)數(shù)據(jù)判斷對(duì)效果的影響分析：若去掉流程圖中自相關(guān)與互相關(guān)數(shù)據(jù)判斷步驟，可能會(huì)得到理想維納濾波不如FIR濾波的效果，如其中一個(gè)結(jié)果：ei= 0

6、.2503，ef= 0.2495。這里的判斷步驟就是為了檢測(cè)實(shí)際產(chǎn)生序列的自相關(guān)或互相關(guān)特性與理論值的近似程度，若誤差很小且通過(guò)我們?cè)O(shè)定的某一下限則認(rèn)為二者近似，所以最終的濾波效果才很近似。如果沒有這里的判斷，實(shí)際自相關(guān)或互相關(guān)則是任意的，完全有可能出現(xiàn)比理想維納濾波更好的效果。3、 固定L=5000,分別取N=3、20，根據(jù)實(shí)驗(yàn)結(jié)果，觀察N的大小對(duì)(n)的估計(jì)和濾波效果的影響并記錄實(shí)驗(yàn)結(jié)果。實(shí)驗(yàn)結(jié)果：圖6為N=3時(shí)估計(jì)(n)與理想h(n)的對(duì)比圖。圖7為N=20時(shí)估計(jì)(n)與理想h(n)的對(duì)比圖。圖8為N=3的FIR濾波后所得(n)與實(shí)際S(n)后100位的比較圖。圖9為N=20的FIR濾波

7、后所得(n)與實(shí)際S(n)后100位的比較圖。其均方誤差分別為：ei= 0.3175(N=3), 0.2500（N=20） ef= 0.2762(N=3), 0.2488 (N=20)圖6 圖7 圖8 圖9分析：由圖6、7可知，N的大小決定(n)與h(n)取值的個(gè)數(shù)，并通過(guò)觀察并結(jié)合N=10的情況可知，N越大(n)與h(n)越接近。從最終均方誤差的比較可知，N越大，濾波效果越好。4、 固定N=10,改變L=1000、5000，根據(jù)實(shí)驗(yàn)結(jié)果，觀察并記錄L的大小對(duì)(n)的精度和濾波效果的影響。實(shí)驗(yàn)結(jié)果：圖10為L(zhǎng)=1000時(shí)估計(jì)(n)與理想h(n)的對(duì)比圖。圖11為L(zhǎng)=5000時(shí)估計(jì)(n)與理想h

8、(n)的對(duì)比圖。圖12為L(zhǎng)=1000的FIR濾波后所得(n)與實(shí)際S(n)后100位的比較圖。圖13為L(zhǎng)=5000的FIR濾波后所得(n)與實(shí)際S(n)后100位的比較圖。其均方誤差分別為：ei=0.2400 (L=1000), 0.2381（L=5000） ef= 0.2390 (L=1000), 0.2375 (L=5000) 圖10 圖11 圖12 圖13分析：由圖10、11可知，L越大(n)與h(n)越接近，(n)的精度越高。由均方誤差可知，L越大，濾波效果越高。這也容易理解，樣本越大，精度自然越高。5、 維納一步純預(yù)測(cè) 1） 畫出信號(hào)生成模型參數(shù)估計(jì)的流程圖開始輸入信號(hào)生產(chǎn)模型的階數(shù)

9、p, AR模型的參數(shù)ai(i=1,2p),w2,信號(hào)s(n)的樣本數(shù)L利用randn函數(shù)產(chǎn)生L個(gè)w(n),并產(chǎn)生L個(gè)s(n)利用Yule-Walker 方程，求出1.p結(jié)束2）根據(jù)流程圖編寫程序（見附錄2）3）運(yùn)行信號(hào)生成模型程序，選擇p=1,a1=-0.6,L=100. 理論值： w2=1-a12=0.6400 a1=-0.6 估計(jì)值： w2= 0.9860 1= -0.5876 相對(duì)誤差：error-w2= -0.0140 error-a1= -0.02064)固定p=1,a1=-0.6, w2=1,改變L=50、500，觀察L的大小對(duì)信號(hào)生成模型參數(shù)估計(jì)精度的影響。實(shí)驗(yàn)結(jié)果：理論值： w

10、2=1-a12=0.6400 a1=-0.6估計(jì)值：L=50 w2= 0.9730 1= -0.5760 估計(jì)值：L=500 w2= 0.9967 1= -0.5965相對(duì)誤差：L=50 error-w2= -0.0270 error-a1= -0.0401相對(duì)誤差：L=500 error-w2= -0.0033 error-a1= -0.0058分析：顯然樣本個(gè)數(shù)L的增大，使得信號(hào)模型參數(shù)精度明顯提高。四、 實(shí)驗(yàn)總結(jié)通過(guò)實(shí)驗(yàn)結(jié)果及分析可得出以下結(jié)論：1、 樣本個(gè)數(shù)越大，參數(shù)精度越高。2、 影響維納濾波效果的因素包括樣本個(gè)數(shù)L、FIR濾波階數(shù)，且均成正比關(guān)系。3、 維納一步純預(yù)測(cè)，只要調(diào)整ai

11、(1,2p)即可實(shí)現(xiàn)最小均方誤差。五、 思考題答案1、推導(dǎo)公式,驗(yàn)證式 推導(dǎo)：已知a= 0.95,w(n)為零均值方差為的高斯白噪聲，v(n)是與s(n)互不相關(guān)的高斯白噪聲，其均值為零，方差。A(z)=,所有S(z)=W(z)A(z)=2、由公式s(n)+(n-1)+s(n-p)=w(n),怎樣得到和？ 分析： 理論w(n)已知，即均值及已知，那么根據(jù)Yule-Walker方程有,其中為(p+1)*(p+1)的s(n)自相關(guān)矩陣,A為（p+1）*1的系數(shù)列向量及A=,而由給出的理論，解方程即可得到估計(jì)值；用估計(jì)值代入方程即可得到估計(jì)值。 六、 源程序見附錄1、2附錄1clear allL=i

12、nput('請(qǐng)輸入信號(hào)樣本個(gè)數(shù)L=');N=input('請(qǐng)輸入FIR濾波器的階數(shù)N=');a = 0.95;K = 50;sigma_a2 = 1-a2;a_ = 1, -a;while(1)% 利用randn()函數(shù)產(chǎn)生白噪聲wn = sqrt(sigma_a2)*( randn(L,1);sn = filter(1, a_, wn);%H(Z)=1/(1-az_-1)vn = randn(L,1);xn = sn + vn;r_xx = xcorr(xn,'unbiased'); %x(n)z自相關(guān)估計(jì)值 r_xx_t = a.abs(-K

13、:K); % x(n)自相關(guān)理論值 r_xx_t(K+1)=r_xx_t(K+1)+1;p = xcorr(sn,xn,'unbiased');%x(n)與s(n)互相關(guān)估計(jì)值r_xs = p(L : L+K); r_xs_t=a.0:K;%x(n)與s(n)互相關(guān)理論值%檢測(cè)實(shí)際值與理論值的近似程度rou_xx = sum(r_xx(L-K:L+K)-r_xx_t').2)/sum(r_xx_t.2);rou_xs = (sum(r_xs-r_xs_t').2)/sum(r_xs_t.2);if rou_xx < 0.03 & rou_xs &l

14、t; 0.01break;endend%同一坐標(biāo)繪制x(n)自相關(guān)函數(shù)理論值與實(shí)際值figure(1),clfstem(r_xx(L-K:L+K),'r') hold onstem(r_xx_t,'k') title('r_xx的實(shí)際值（紅色）及理論值（黑色）')%同一坐標(biāo)繪制x(n)與s(n)互相關(guān)函數(shù)理論值與實(shí)際值figure(2),clfstem(r_xs,'r') hold onstem(r_xs_t,'k') title('r_xs的實(shí)際值（紅色）及理論值（黑色）')%同一坐標(biāo)繪制最后10

15、0個(gè)s(n)與x(n)figure(3),clfstem(xn(L-99:L),'r')hold onstem(sn(L-99:L),'k')title('最后100個(gè)s(n)（黑色）和x(n)（紅色）');%計(jì)算h(n)的估計(jì)值并與理想值比較%構(gòu)造x(n)的N階自相關(guān)矩陣R_xxn=0:N-1;for i=1:N for j=1:N R_xx(i,j)=r_xx(i-j+L); endendhopt=inv(R_xx)*r_xs(1:N);%利用維納霍夫方程求h hopt_t=0.238*(0.724).n; %理想h %同一坐標(biāo)繪制h(n)的

16、實(shí)際值與理想值figure(4),clfstem(hopt,'-','r');hold onstem(hopt_t,'*','k');title('h(n)估計(jì)值（紅色），與真值（黑色）的比較');%同一座標(biāo)繪制理想維納濾波后的s(n)和最后100個(gè)s(n)sn_w=filter(hopt_t,1,xn);%理想維納濾波figure(5),clfstem(sn_w(L-99:L),'*','r');hold onstem(sn(L-99:L),'.','k&#

17、39;)title('最后100個(gè)sn（黑色）和由維納濾波器得到的sn_w（紅色）的比較');%同一坐標(biāo)繪制實(shí)際維納濾波后s(n)與最后100個(gè)s(n)sn_f=filter(hopt,1,xn);figure(6),clfstem(sn_f(L-99:L),'*','r');hold onstem(sn(L-99:L),'.','k')title('最后100個(gè)sn（黑色）和由FIR濾波器得到的sn_w（紅色）的比較');%求并比較各個(gè)均方差e_x=sum(xn(1:L)-sn(1:L).2)/L

18、 e_i=sum(sn_w(1:L)-sn(1:L).2)/Le_f=sum(sn_f(1:L)-sn(1:L).2)/L附錄二clear all;%輸入AR模型的參數(shù)p=1;a1=-0.6;sigma_w2=1;L=500;%s(m)的樣本數(shù)s=1;%統(tǒng)計(jì)5000次，最終估計(jì)值取平均while s<=5000wn=sqrt(sigma_w2)*randn(1,L);%w(n)為均值為0，方差為1的白噪聲sn=wn(1);%s(n)+a1s(n-1)=w(n),n=1時(shí)，s(n)=w(n)for smp=2:Lsn=sn,wn(smp)-a1*sn(smp-1);%n>=2時(shí)，s(n)=w(n)-a1s(n-1)endr_ss=xcorr(sn,'unbiased');%求s(n)的自相關(guān)est_a1(s)=-r_ss(L