基于最小均方誤差(MMSE)估計(jì)的因果維納濾波的實(shí)現(xiàn)

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

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

3、陣形式為：簡化形式：其中：是濾波器的系數(shù)是互相關(guān)序列是自相關(guān)矩陣由上可見，設(shè)計(jì)維納濾波器的過程就是尋求在最小均方誤差下濾波器的單位脈沖響應(yīng)或傳遞函數(shù)的表達(dá)式，其實(shí)質(zhì)就是解維納霍夫（WienerHopf）方程。另外，設(shè)計(jì)維納濾波器要求已知信號(hào)與噪聲的相關(guān)函數(shù)。三程序求解過程由上述可見，本程序?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ù)說明函數(shù)使用方

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

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

6、4096); %畫出加入噪聲后信號(hào)的頻譜圖title('加入噪聲后語音信號(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,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); %畫出維納濾波后信號(hào)的頻譜圖title('經(jīng)過維納濾