DSP 維納濾波器的計算機實現(xiàn)

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

2、1xxxs）。但要求h(n)滿足因果性要求，維納-霍夫方程便是一個難題，這里利用最佳FIR維納濾波方法求解h(n)的近似，這也便于在計算機上實現(xiàn)，公式為：h =R-1xx rxs。實驗中s(n)由信號生成模型：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)。實驗中s(n)是已知的，但實際中如果s(n)已知，維納濾波也就失去意義了，因此實驗純粹是為了理解維納濾波原理而設(shè)計。2、 維納一步純預(yù)測問題S(n)的生成模型：s(n)+a1(n-1)+aps(n-p)=w(n),

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

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

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

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

7、后所得(n)與實際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)取值的個數(shù)，并通過觀察并結(jié)合N=10的情況可知，N越大(n)與h(n)越接近。從最終均方誤差的比較可知，N越大，濾波效果越好。4、 固定N=10,改變L=1000、5000，根據(jù)實驗結(jié)果，觀察并記錄L的大小對(n)的精度和濾波效果的影響。實驗結(jié)果：圖10為L=1000時估計(n)與理想h(n)的對比圖。圖11為L=5000時估計(n)與理想h

8、(n)的對比圖。圖12為L=1000的FIR濾波后所得(n)與實際S(n)后100位的比較圖。圖13為L=5000的FIR濾波后所得(n)與實際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ù)測 1） 畫出信號生成模型參數(shù)估計的流程圖開始輸入信號生產(chǎn)模型的階數(shù)

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

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

11、(1,2p)即可實現(xiàn)最小均方誤差。五、 思考題答案1、推導(dǎo)公式,驗證式 推導(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=,而由給出的理論，解方程即可得到估計值；用估計值代入方程即可得到估計值。 六、 源程序見附錄1、2附錄1clear allL=i

12、nput('請輸入信號樣本個數(shù)L=');N=input('請輸入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)估計值 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)估計值r_xs = p(L : L+K); r_xs_t=a.0:K;%x(n)與s(n)互相關(guān)理論值%檢測實際值與理論值的近似程度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%同一坐標繪制x(n)自相關(guān)函數(shù)理論值與實際值figure(1),clfstem(r_xx(L-K:L+K),'r') hold onstem(r_xx_t,'k') title('r_xx的實際值（紅色）及理論值（黑色）')%同一坐標繪制x(n)與s(n)互相關(guān)函數(shù)理論值與實際值figure(2),clfstem(r_xs,'r') hold onstem(r_xs_t,'k') title('r_xs的實際值（紅色）及理論值（黑色）')%同一坐標繪制最后10

15、0個s(n)與x(n)figure(3),clfstem(xn(L-99:L),'r')hold onstem(sn(L-99:L),'k')title('最后100個s(n)（黑色）和x(n)（紅色）');%計算h(n)的估計值并與理想值比較%構(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 %同一坐標繪制h(n)的

16、實際值與理想值figure(4),clfstem(hopt,'-','r');hold onstem(hopt_t,'*','k');title('h(n)估計值（紅色），與真值（黑色）的比較');%同一座標繪制理想維納濾波后的s(n)和最后100個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個sn（黑色）和由維納濾波器得到的sn_w（紅色）的比較');%同一坐標繪制實際維納濾波后s(n)與最后100個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個sn（黑色）和由FIR濾波器得到的sn_w（紅色）的比較');%求并比較各個均方差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)計5000次，最終估計值取平均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時，s(n)=w(n)for smp=2:Lsn=sn,wn(smp)-a1*sn(smp-1);%n>=2時，s(n)=w(n)-a1s(n-1)endr_ss=xcorr(sn,'unbiased');%求s(n)的自相關(guān)est_a1(s)=-r_ss(L