wiener濾波程序

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

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

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

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

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

6、4096); %畫(huà)出加入噪聲后信號(hào)的頻譜圖title('加入噪聲后語(yǔ)音信號(hào)的頻域圖形');xlabel('頻率 f');ylabel('FFT');yyhxcorr=xcorr(x(1:4096); %求取信號(hào)的信號(hào)的自相關(guān)函數(shù)size(yyhxcorr); A=yyhxcorr(4096:4595);yyhdcorr=xcorr(d(1:4096),x(1:4096); %求取信號(hào)和理想信號(hào)的互相關(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); %對(duì)維納濾波的結(jié)果進(jìn)行傅立葉變換，看其頻譜變化subplot(3,1,3); f=Fs*(0:4095)/8192;plot(f,abs(fq(1:4096); %畫(huà)出維納濾波后信號(hào)的頻譜圖title('經(jīng)過(guò)維納濾波后語(yǔ)音信號(hào)的頻域圖形');xlabel('頻率 f');ylabel('FFT');求出的頻譜圖如下所示：由上述結(jié)果可見(jiàn)，經(jīng)過(guò)維納濾波后信號(hào)的噪聲減弱，信噪比提高。yyhwiener.mload handel

8、 %sound(5*y,Fs) %原始語(yǔ)音信號(hào),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('原始語(yǔ)音信號(hào)的頻域圖形'); xlabel('頻率 f'); ylabel('FFT'); m,n=size(d); x_noise=randn(1,n); x=d+x_noise; %加入噪聲后的語(yǔ)音信號(hào)，噪聲為(0,1)分布的高斯白噪聲 fq=fft(x,8192); subp

9、lot(3,1,2); plot(f,abs(fq(1:4096); title('加入噪聲后語(yǔ)音信號(hào)的頻域圖形'); 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)過(guò)維納濾波后語(yǔ)音信號(hào)的頻域圖形'); xlabel('頻率 f'); ylabel('FFT');wienerfilter.mfunction y=wienerfilter(x,Rxx,Rxd,N) %進(jìn)行維納濾波 %x是輸入信號(hào)，Rxx是輸入信號(hào)的自相關(guān)向量 %Rxx是輸入信號(hào)和理想信號(hào)的的互相

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

12、xx(1),.,Rxx(M-1) %B是接收信號(hào)和沒(méi)有噪聲干擾信號(hào)的互相關(guān)向量為 Rxd(0),Rxd(1),.,Rxd(M-1) %M是濾波器的長(zhǎng)度 %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(