wiener濾波程序

1、基于最小均方誤差(MMSE)估計(jì)的因果維納濾波的實(shí)現(xiàn)一功能簡介 基于最小均方誤差(MMSE)估計(jì)的因果維納濾波的Matlab實(shí)現(xiàn)，用萊文森-德賓(Levinson-Durbin)算法求解維納-霍夫方程(Yule-wa1ker)方程，得到濾波器系數(shù)，進(jìn)行維納濾波。二維納濾波簡介信號處理的實(shí)際問題，常常是要解決在噪聲中提取信號的問題，因此，我們需要尋找一種所謂有最佳線性過濾特性的濾波器，這種濾波器當(dāng)信號與噪聲同時輸入時，在輸出端能將信號盡可能精確地重現(xiàn)出來，而噪聲卻受到最大抑制。維納(Wiener)濾波就是用來解決這樣一類從噪聲中提取信號問題的一種過濾(或?yàn)V波)方法。一個線性系統(tǒng)，如果它的單位樣本

2、響應(yīng)為h(n)，當(dāng)輸入一個隨機(jī)信號x(n)，且其中s(n)表示信號，表示噪聲，則輸出y(n)為我們希望x(n)通過線性系統(tǒng)h(n)后得到的y(n)盡量接近于s(n)，因此稱y(n)為s(n)的估計(jì)值，用表示，即維納濾波器的輸入輸出關(guān)系如上圖所示。這個線性系統(tǒng)稱為對于的一種估計(jì)器。如果我們以分別表示信號的真值與估計(jì)值，而用e(n)表示它們之間的誤差，即顯然，e(n)可能是正的，也可能是負(fù)的，并且它是一個隨機(jī)變量。因此，用它的均方值來表達(dá)誤差是合理的，所謂均方誤差最小即它的平方的統(tǒng)計(jì)平均值最?。鹤钚∫阎Ｍ敵鰹椋赫`差為：均方誤差為：上式對求導(dǎo)得到：進(jìn)一步得：從而有：于是就得到N個線性方程：寫成矩

3、陣形式為：簡化形式：其中：是濾波器的系數(shù)是互相關(guān)序列是自相關(guān)矩陣由上可見，設(shè)計(jì)維納濾波器的過程就是尋求在最小均方誤差下濾波器的單位脈沖響應(yīng)或傳遞函數(shù)的表達(dá)式，其實(shí)質(zhì)就是解維納霍夫（WienerHopf）方程。另外，設(shè)計(jì)維納濾波器要求已知信號與噪聲的相關(guān)函數(shù)。三程序求解過程由上述可見，本程序?qū)崿F(xiàn)的關(guān)鍵是在已知輸入信號的自相關(guān)函數(shù)和輸入信號和理想輸出信號的互相關(guān)函數(shù)的情況下，求解維納霍夫（WienerHopf）方程，從而得到濾波器系數(shù)，再進(jìn)行維納濾波。求解步驟：1. 初始化值2. 對于，進(jìn)行如下計(jì)算：3濾波器系數(shù)為：4利用上面的得到的濾波器對輸入信號進(jìn)行維納濾波，得到輸出信號。四函數(shù)說明函數(shù)使用方

4、法：y=wienerfilter(x,Rxx,Rxd,M)參數(shù)說明：x是輸入信號，Rxx是輸入信號的自相關(guān)向量，Rxd是輸入信號和理想信號的的互相關(guān)向量，M是維納濾波器的長度，輸出y是輸入信號通過維納濾波器進(jìn)行維納濾波后的輸出。具體程序見Matlab的.m文件。五程序示例加載Matlab中的語音數(shù)據(jù)handel，人為地加入高斯白噪聲，分別計(jì)算加入噪聲后信號的自相關(guān)和加入噪聲后信號和理想信號的互相關(guān)，取濾波器的長度為M=500，將以上參數(shù)代入函數(shù)中進(jìn)行維納濾波，得到輸出。程序如下：load handel %加載語音信號d=y; d=d*8; %增強(qiáng)語音信號強(qiáng)度d=d'fq=fft(d,8

5、192); %進(jìn)行傅立葉變換得到語音信號頻頻subplot(3,1,1);f=Fs*(0:4095)/8192;plot(f,abs(fq(1:4096); %畫出頻譜圖title('原始語音信號的頻域圖形');xlabel('頻率 f');ylabel('FFT');m,n=size(d);x_noise=randn(1,n); %（0，1）分布的高斯白噪聲x=d+x_noise; %加入噪聲后的語音信號fq=fft(x,8192); %對加入噪聲后的信號進(jìn)行傅立葉變換，看其頻譜變化subplot(3,1,2);plot(f,abs(fq(1:

6、4096); %畫出加入噪聲后信號的頻譜圖title('加入噪聲后語音信號的頻域圖形');xlabel('頻率 f');ylabel('FFT');yyhxcorr=xcorr(x(1:4096); %求取信號的信號的自相關(guān)函數(shù)size(yyhxcorr); A=yyhxcorr(4096:4595);yyhdcorr=xcorr(d(1:4096),x(1:4096); %求取信號和理想信號的互相關(guān)函數(shù)size(yyhdcorr);B=yyhdcorr(4096:4595);M=500;yyhresult=wienerfilter(x,A,B,

7、M); %進(jìn)行維納濾波yyhresult=yyhresult(300:8192+299);fq=fft(yyhresult); %對維納濾波的結(jié)果進(jìn)行傅立葉變換，看其頻譜變化subplot(3,1,3); f=Fs*(0:4095)/8192;plot(f,abs(fq(1:4096); %畫出維納濾波后信號的頻譜圖title('經(jīng)過維納濾波后語音信號的頻域圖形');xlabel('頻率 f');ylabel('FFT');求出的頻譜圖如下所示：由上述結(jié)果可見，經(jīng)過維納濾波后信號的噪聲減弱，信噪比提高。yyhwiener.mload handel

8、 %sound(5*y,Fs) %原始語音信號,Fs=9192 d=y; d=d*8; d=d' fq=fft(d,8192); subplot(3,1,1); f=Fs*(0:4095)/8192; plot(f,abs(fq(1:4096); title('原始語音信號的頻域圖形'); xlabel('頻率 f'); ylabel('FFT'); m,n=size(d); x_noise=randn(1,n); x=d+x_noise; %加入噪聲后的語音信號，噪聲為(0,1)分布的高斯白噪聲 fq=fft(x,8192); subp

9、lot(3,1,2); plot(f,abs(fq(1:4096); title('加入噪聲后語音信號的頻域圖形'); xlabel('頻率 f'); ylabel('FFT'); yyhxcorr=xcorr(x(1:4096); size(yyhxcorr); A=yyhxcorr(4096:4595); yyhdcorr=xcorr(d(1:4096),x(1:4096); size(yyhdcorr); B=yyhdcorr(4096:4595); M=500; yyhresult=wienerfilter(x,A,B,M); %用維納濾

10、波進(jìn)行去噪 yyhresult=yyhresult(300:8192+299); fq=fft(yyhresult); subplot(3,1,3); f=Fs*(0:4095)/8192; plot(f,abs(fq(1:4096); title('經(jīng)過維納濾波后語音信號的頻域圖形'); xlabel('頻率 f'); ylabel('FFT');wienerfilter.mfunction y=wienerfilter(x,Rxx,Rxd,N) %進(jìn)行維納濾波 %x是輸入信號，Rxx是輸入信號的自相關(guān)向量 %Rxx是輸入信號和理想信號的的互相

11、關(guān)向量，N是維納濾波器的長度 %輸出y是輸入信號通過維納濾波器進(jìn)行維納濾波后的輸出 h=yulewalker(Rxx,Rxd,N);%求解維納濾波器系數(shù) t=conv(x,h); %進(jìn)行濾波 Lh=length(h); %得到濾波器的長度 Lx=length(x); %得到輸入信號的長度 y=t(double(uint16(Lh/2):Lx+double(uint16(Lh/2)-1);%輸出序列y的長度和輸入序列x的長度相同 %以下是維納濾波器系數(shù)的求解 function h=yulewalker(A,B,M) %求解Yule-Walker方程 %A是接收信號的自相關(guān)向量為 Rxx(0),R

12、xx(1),.,Rxx(M-1) %B是接收信號和沒有噪聲干擾信號的互相關(guān)向量為 Rxd(0),Rxd(1),.,Rxd(M-1) %M是濾波器的長度 %h保存濾波器的系數(shù) %例如A=6,5,4,3,2,1;B=100,90,120,50,80,200;M=6; %求解出h=26.4286 -20.0000 50.0000 -50.0000 -45.0000 81.4286' T1=zeros(1,M);%T1存放中間方程的解向量 T2=zeros(1,M);%T2存放中間方程的解向量 T1(1)=B(1)/A(1); T2(1)=A(2)/A(1); X=zeros(1,M); fo

13、r i=2:M-1 temp1=0; temp2=0; for j=1:i-1 temp1=temp1+A(i-j+1)*T1(j); temp2=temp2+A(i-j+1)*T2(j); end X(i)=(B(i)-temp1)/(A(1)-temp2); for j=1:i-1 X(j)=T1(j)-X(i)*T2(j); end for j=1:i T1(j)=X(j); end temp1=0; temp2=0; for j=1:i-1 temp1=temp1+A(j+1)*T2(j); temp2=temp2+A(j+1)*T2(i-j); end X(1)=(A(i+1)-temp1)/(A(1)-temp2); for j=2:i X(j)=T2(j-1)-X(1)*T2(