wiener滤波程序.docx

《wiener滤波程序.docx》由会员分享，可在线阅读，更多相关《wiener滤波程序.docx（9页珍藏版）》请在第一文库网上搜索。

1、基于最小均方误差(MMSE)估计的因果维纳滤波的实现一.功能简介基于最小均方误差(MMSE)估计的因果维纳滤波的Mat1ab实现，用莱文森-德宾(1eVinSon-DUrbin)算法求解维纳-霍夫方程(YU1e-Wa1ker)方程，得到滤波器系数，进行维纳滤波。二.维纳滤波简介信号处理的实际问题，常常是要解决在噪声中提取信号的问题，因此，我们需要寻找一种所谓有最佳线性过滤特性的滤波器，这种滤波器当信号与噪声同时输入时，在输出端能将信号尽可能精确地重现出来，而噪声却受到最大抑制。维纳(Wiener)滤波就是用来解决这样一类从噪声中提取信号问题的一种过滤(或滤波)方法。一个线性系统，如果它的单位样

2、本响应为人()，当输入一个随机信号x()，且x(n)=s(n)+(h)其中$()表示信号，表示噪声，则输出y()为y(n)=Zh(m)x(n-m)m我们希望通过线性系统励)后得到的M)尽量接近于s()，因此称刈)为s()的估计值,用假)表示，即y(n)=s(n)Jh(n)-%x(n)=s(n)+(n)维纳滤波器的输入一输出关系如上图所示。这个线性系统z()称为对于s(n)的一种估计器。如果我们以S与分别表示信号的真值与估计值，而用e()表示它们之间的误差，即e(n)=s5)s(n)显然，式)可能是正的，也可能是负的，并且它是一个随机变量。因此，用它的均方值来表达误差是合理的，所谓均方误差最小即

3、它的平方的统计平均值最小：Eys)=eg一丹最小已知希望输出为：NTy(n)=s(n)=Zh(m)x(n-ni)/W=O误差为：NTe(n)=s(n)-s(n)=s(n)-Zh(m)x(n-tn)w=0均方误差为：N-七/()=七($()一(m)x(一n)21W=O.上式对力(z)m=0,1,N-1求导得到：,-12E(5(n)-Xhtyp1(ni)x(n-ni)xn-j)=Oj=0,1,2N-I_n=0_进一步得：N-Is()x(-J)=Zh0pt(m)Ex(n-m)x(n-j)j=0,1,N-I/W=O从而有：Ar-IR式j)=*(m)Rx(j-m)7=0,1,2,.,7V-1m=0于是就

4、得到N个线性方程：7=0RKO)=MO)RKo)+%(1)RK1)+(NT)RKN-I)7=1RK1)=(0)段41)+/川/式0)+以-1)/?式-2)J=TV-IRxS(N-1)=MO)RXX(N-1)+MDRxX(N-2)+VN-I)RKo)写成矩阵形式为:RU(O)此RKN-“h(0)RKo)RK1)RKo)Ra(N-2)=RK1)AMDRKN-2)鼠(0)M(N-1)A(N力简化形式：RRH=RXS其中：H=h(O)h(1)h(N-1)是滤波器的系数RXS=凡(0),4(1),&(NT)T是互相关序列-鼠(0)&RnNFRRK1)RKO)RKN-2)是白柏至矩阵Rxx=.是自相关矩阵

5、A(NT)RxX(N-2).&(0)由上可见，设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳-霍夫(Wiener-Hopf)方程。另外，设计维纳滤波器要求已知信号与噪声的相关函数。三.程序求解过程由上述可见，本程序实现的关键是在已知输入信号的自相关函数和输入信号和理想输出信号的互相关函数的情况下，求解维纳-霍夫(Wiener-Hopf)方程，从而得到滤波器系数,再进行维纳滤波。求解步骤：1 .初始化值4(0)=力(0)/噎(0)仇0)=力/仁(0)2 .对于j=1,2,M-1,进行如下计算：rxj)SJ-i)a(i)tempi=Wi=0a(

6、)=a(i)-tempib(i)i=O,1,j-1a(j)=tempitempi=rxx)-rxx(j-iW)i=Qb(i)=b(i-1)-tempib(j-i)i=1,j伙0)=tempi3 .滤波器系数为：h(i)a(i)i=0,1,.M-14 .利用上面的得到的滤波器对输入信号进行维纳滤波，得到输出信号。四.函数说明函数使用方法：y=wienerfiIter(x,Rxx,Rxd,M)参数说明：X是输入信号，RXX是输入信号的自相关向量，RXd是输入信号和理想信号的的互相关向量，M是维纳滤波器的长度，输出y是输入信号通过维纳滤波器进行维纳滤波后的输出。具体程序见Ma11ab的.m文件。五.

7、程序示例加载MatIab中的语音数据hande1,人为地加入高斯白噪声，分别计算加入噪声后信号的自相关RXX和加入噪声后信号和理想信号的互相关Rxd，取滤波器的长度为M=500,将以上参数代入函数中进行维纳滤波，得到输出。程序如下：1oadhande1%加载语音信号d=y；d=d*8;%增强语音信号强度d=d;fq=fft(d,8192);%进行傅立叶变换得到语音信号频频subp1ot(3,1,1);f=Fs*(0:4095)/8192;p1ot(f,abs(fq(1:4096);%画出频谱图IiUeC原始语音信号的频域图形);XIabeIC频率f)；y1abe1(FFT,);m,n=size

8、(d);x_noise=randn(1,n);%(01)分布的高斯白噪声x=d+x_noise;%加入噪声后的语音信号fq=fft(x,8192);%对加入噪声后的信号进行傅立叶变换，看其频谱变化subp1ot(3,1,2);p1ot(f,abs(fq(1:4096);%画出加入噪声后信号的频谱图tit1e，加入噪声后语音信号的频域图形)；X1abe1C频率)；y1abe1(FFT,);yyhxcorr=xcorr(x(1:4096);%求取信号的信号的自相关函数size(yyhxcorr);A=yyhxcorr(4096:4595);yyhdcorr=xcorr(d(1:4096),x(1:

9、4096);%求取信号和理想信号的互相关函数size(yyhdcorr);B=yyhdcorr(4096:4595);M=500;yyhresu1t=wienerfi11er(x,A,B,M);%进行维纳滤波yyhresu1t=yyhresu11(300:8192+299);fq=fft(yyhresu1t);%对维纳滤波的结果进行傅立叶变换，看其频谱变化SUbP1ot(3,1,3);f=Fs*(0:4095)/8192;p1ot(f,abs(fq(1:4096);%画出维纳滤波后信号的频谱图IiUec经过维纳滤波后语音信号的频域图形力X1abeIe频率0；y1abe1(FFT,);求出的频谱

10、图如下所示：原始语音信号的频域图形频率f加入噪声后语音信号的频域图形频率f由上述结果可见，经过维纳滤波后信号的噪声减弱，信噪比提高。yyhwiener.m1oadhande1%sound(5*y,FS)%原始语音信号，Fs=9192d=y;CI=C1*8;d=d,;fq=fft(d,8192);subp1ot(3,1,1);f=Fs*(0:4095)/8192;p1ot(f,abs(fq(1:4096);tit1e(,原始语音信号的频域图形);x1abe1(,频率f,);y1abe1(,FFT);m,n=size(d);Xnoise=randn(1,n);x=c1+xnoise;%加入噪声后的

11、语音信号，噪声为(0,1)分布的高斯白噪声fq=fft(x,8192);subp1ot(3,1,2);p1ot(f,abs(fq(1:4096);tit1e(,加入噪声后语音信号的频域图形)；x1abe1(,频率f,);y1abe1(,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;yyhresu1t=wienerfiIter(x,A,B,M);%用维纳滤

12、波进行去噪yyhresu1t=yyhresu1t(300:8192+299);fq=fft(yyhresu1t);subp1ot(3,1,3);f=Fs*(0:4095)/8192;p1ot(f,abs(fq(1:4096);tit1e(,经过维纳滤波后语音信号的频域图形);x1abe1(,频率f);y1abe1(,FFT);wienerfi1ter.mfunctionenerfi1ter(x,Rxx,Rxd,N)%进行维纳滤波%是输入信号，Rxx是输入信号的自相关向量%R是输入信号和理想信号的的互相关向量，N是维纳滤波器的长度%输出y是输入信号通过维纳滤波器进行维纳滤波后的输出h=yu1ew

13、a1ker(Rxx,Rxd,N)求解维纳滤波器系数t=conv(x,h);%进行滤波1h=Iength(h);%得到滤波器的长度1x=Iength(x);%得到输入信号的长度y=t(doub1e(uint16(1h2):1+(10此1,6(111111:16(5/2)-1);%输出序列丫的长度和输入序列X的长度相同%以下是维纳滤波器系数的求解functionh=yu1ewa1ker(,B,M)%求解Yu1e-Wa1ker方程%是接收信号的自相关向量为RXX(O),Rxx(1),Rxx(M-I)%B是接收信号和没有噪声干扰信号的互相关向量为Rxd(0),Rxd(I),Rxd(M-I)%M是滤波器

14、的长度%h保存滤波器的系数%例如A=6,5,4,3,2,1;B=100,90,120,50,80,200;M=6;%求解出h=26.4286-20.000050.0000-50.0000-45.000081.4286T1=Zeros(IJI);%T1存放中间方程的解向量T2=zeros(1,M);%T2存放中间方程的解向量T1(I)=B(I)ZA(I);T2(1)=A(2)A(1);X=zeros(1,M);fori=2:M-1temp1=0;temp2=0;forj=1:i-1temp1=temp1+A(i-j1)*T1(j);temp2=temp2+A(i-j+1)*T2(j);endX(i)=(B(i)-temp1)/(A(1)-temp2);forj=1:i-1X(j)=T1(j)-X(i)*T2(j);endforj=1:iT1(j)=X(j);endtemp1=0;temp2=0;forj=1:i-1temp1=temp1+A(j+1)*T2(j);temp2=temp2+A