FIR維納濾波的設(shè)計(jì)

FIR維納濾波器的設(shè)計(jì)在信號(hào)處理的許多實(shí)際應(yīng)用中，人們往往不能直接獲得所需要的有用信號(hào)，需要從噪聲中提取信號(hào)。比如，在信號(hào)傳輸過(guò)程中，由于存在信道噪聲等干擾，在接收端觀測(cè)到的信號(hào)必然與原始信號(hào)不同。為了從觀測(cè)數(shù)據(jù)中盡可能精確地重現(xiàn)原始信號(hào)，而最大成都地抑制噪聲，需要設(shè)計(jì)一種濾波器，其輸出盡可能逼近原始信號(hào)，成為原始信號(hào)的最佳估計(jì)。這種濾波器成為最佳濾波器。維納(Wiener)濾波器就是用來(lái)解決這樣一類(lèi)問(wèn)題的一種濾波器。本文將應(yīng)用MATLAB并結(jié)合實(shí)例介紹FIR維納濾波器的設(shè)計(jì)方法。一、維納濾波的原理維納濾波的本質(zhì)是一種最佳估計(jì)問(wèn)題，采用的是最小均方誤差準(zhǔn)則。一個(gè)線(xiàn)性系統(tǒng)，其單位樣本響應(yīng)為h(n)，當(dāng)輸入一個(gè)隨機(jī)信號(hào)其中s(n)表示信號(hào)，表示噪聲，則輸出y(n)為 (1)系統(tǒng)是通過(guò)y(n)來(lái)估計(jì)s(n)，因此將其稱(chēng)為s(n)的估計(jì)值，用表示，即 (2)圖1維納濾波器基本框圖圖1所示為維納濾波器的基本框圖。式(1)為一卷積，可以理解為從當(dāng)前和過(guò)去的觀察值x(n)，x(n-1)，x(n-2)…x(n-m)，…來(lái)估計(jì)信號(hào)的當(dāng)前值。維納濾波器一般有三種用途。用當(dāng)前的和過(guò)去的觀察值x(n)，x(n-1)，x(n-2)，…來(lái)估計(jì)當(dāng)前的信號(hào)值稱(chēng)為濾波；用過(guò)去的觀察值來(lái)估計(jì)當(dāng)前的或?qū)?lái)的信號(hào)值稱(chēng)為預(yù)測(cè)；用全部數(shù)據(jù)來(lái)估計(jì)過(guò)去的信號(hào)值稱(chēng)為平滑。維拉濾波采用的是最小均方誤差準(zhǔn)則，以與分別表示信號(hào)的真值與估計(jì)值，而用e(n)表示它們之間的誤差，即 (3)那么最小均方誤差表示為 = (4)為了得到使式(4)最小時(shí)的維納濾波器沖擊響應(yīng)，令對(duì)的導(dǎo)數(shù)等于零，即(5)由式(5)可得(6)式(6)稱(chēng)為正交方程，表明任何時(shí)刻的估計(jì)誤差都與用于估計(jì)的所有數(shù)據(jù)正交。將式(1)，式(3)代入式(6)可得(7)式(7)稱(chēng)為維納-霍夫（Wiener-Hopf）方程，其中是與的互相過(guò)函數(shù)，是的自相關(guān)函數(shù)，其表達(dá)式如下(8)(9)維納濾波器有三種情況，一是FIR維納濾波器；二是非因果IIR維納濾波器；三是因果的IIR維納濾波器。本文僅討論FIR維納濾波器的求解。二、FIR維納濾波器求解設(shè)維納濾波器的沖擊響應(yīng)序列長(zhǎng)度為N，則沖擊響應(yīng)矢量為h=(10)濾波器輸入數(shù)據(jù)矢量為=(11)著濾波器輸出為(12)這樣，維納-霍夫方程可表示成P=Rh(13)其中P=，R=。求解式(13)可得h=三、實(shí)例設(shè)計(jì) 結(jié)合維拉濾波器的原理，利用MATLAB設(shè)計(jì)濾波器。 將隨機(jī)信號(hào)X(n)看成是由典型白噪聲序列源W(n)激勵(lì)一個(gè)線(xiàn)性系統(tǒng)產(chǎn)生，用一個(gè)差分方程來(lái)描述。進(jìn)行Z變換得到，那么均值為１的高斯白噪聲序列W（n）可以用randn函數(shù)產(chǎn)生，再利用函數(shù)X=filter(B,A,W)產(chǎn)生隨機(jī)信號(hào)x(n)。 這里將濾波器的階數(shù)設(shè)為101，根據(jù)維納－霍夫方程：，其中是觀測(cè)信號(hào)的自相關(guān)函數(shù)，是觀測(cè)信號(hào)和期望信號(hào)的互相關(guān)函數(shù)。定義維納濾波的模型，最后帶入filter。具體MATLAB程序代碼如下：maxlag=100;N=100; %采樣次數(shù)為100x=zeros(N,1);y=zeros(N,1);var=1; %%%%%%%%%%%%%%%%%%列出狀態(tài)方程x(1)=randn(1,1);%令x(-1)=x(-2)=x(-3)=x(-4)=0x(2)=randn(1,1)+1.352*x(1);x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);forn=5:Nx(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);end;v=randn(N,1);y=x+v; %z_x為觀測(cè)樣本值=真值+噪聲%%%%%%%%%%%%%%%%%%%%%%%濾波x=x';y=y';xk_s(1)=y(1); %賦初值xk_s(2)=y(2);xk_s(3)=y(3);xk_s(4)=y(4);xk=[y(1);y(2);y(3);y(4)];%%%%%%%%%%%%%%%%%%計(jì)算觀測(cè)函數(shù)%維納濾波器的生成[rx,lags]=xcorr(y,maxlag,'biased');%觀測(cè)信號(hào)的自相關(guān)函數(shù)rx1=toeplitz(rx(101:end));%對(duì)稱(chēng)化自相關(guān)函數(shù)矩陣使之成為方陣，濾波器的階數(shù)為101階rx2=xcorr(x,y,maxlag,'biased');%觀測(cè)信號(hào)與期望信號(hào)的互相關(guān)函數(shù)rx2=rx2(101:end);h=inv(rx1)*rx2'; %維納-霍夫方程xk_s=filter(h,1,y);%加噪信號(hào)通過(guò)濾波器后的輸出%%%%%%%%%%%%%%%%%%%%%計(jì)算誤差e_x=0;eq_x=0;e_x1=N:1;%計(jì)算濾波的均值,計(jì)算濾波誤差的均值fori=1:Ne_x(i)=x(i)-xk_s(i);%誤差=真實(shí)值-濾波估計(jì)值end%%%%%%%%%%%%%%%%%%%%%%%作圖t=1:N;figure(1);plot(t,x,'r',t,y,'g',t,xk_s,'b');legend('