维纳滤波原理及其matlab实现

摘要本文介绍了维纳滤波的原理及其matlab实现，以案例的形式展示FIR维纳滤波的特性。关键字：FIR维纳滤波Matlab1.引言滤波技术是信号分析、处理技术的重要分支，无论是信号的获取、传输，还是信号的处理和交换都离不开滤波技术，它对信号安全可靠和有效灵活地传递是至关重要的。信号分析检测与处理的一个十分重要的内容就是从噪声中提取信号，实现这种功能的有效手段之一是设计一种具有最佳线性过滤特性的滤波器，当伴有噪声的信号通过这种滤波器的时候，它可以将信号尽可能精确地重现或对信号做出尽可能精确的估计，而对所伴随噪声进行最大限度地抑制。维纳滤波器就是这种滤波器的典型代表之一。2.维纳滤波概述维纳（Wiener）是用来解决从噪声中提取信号的一种过滤（或滤波）方法。这种线性滤波问题，可以看做是一种估计问题或一种线性估计问题。一个线性系统，如果它的单位样本响应为)(nh，当输入一个随机信号)(nx，且)()()(nvnsnx（1）其中)(nx表示信号，)(nv)表示噪声，则输出)(ny为mmnxmhny)()()(（2）我们希望)(nx通过线性系统)(nh后得到的)(ny尽量接近于)(ns，因此称)(ny为)(ns的估计值，用^)(ns表示，即^)()(nsny（3）则维纳滤波器的输入—输出关系可用下面图1表示。图1实际上，式（2）所示的卷积形式可以理解为从当前和过去的观察值)(nx，)1(nx，)2(nx…)(mnx，…来估计信号的当前值^)(ns。因此，用)(nh进行过滤问题实际上是一种统计估计问题。一般地，从当前的和过去的观察值)(nx，)1(nx，)2(nx…估计当前的信号值^)()(nsny成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值)0)(()(^NNnsny称为外推或预测；从过去的观察值，估计过去的信号值)1)(()(^NNnsny称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准则的。如果我们分别以)(ns与^)(ns表示信号的真实值与估计值，而用)(ne表示他们之间的误差，即)()()(^nsnsne（4）显然)(ne可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小：min)]([)(2nEne（5）采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概率的描述。3.维纳-霍夫方程的求解为了按（5）式所示的最小均方误差准则来确定维纳滤波器的冲激响应)(nh，令)(n对)(jh的导数等于零，即可得mimRihmRixxxs,)()()(（6）式中，)(mRxs是)(ns与)(nx的互相关函数，)(mRxx是)(nx的自相关函数，分别定义为)]()([mnsnxERxs)]()([mnxnxERxx式（6）称为维纳滤波器的标准方程或维纳-霍夫（Wiener-Hopf）方程。如果已知)(mRxs和)(mRxx，那么解此方程即可求的维纳滤波器的冲激响应。式（6）所示标准方程右端的求和范围即i的取值范围没有具体标明，实际上有三种情况：（1）有限冲激响应（FIR）维纳滤波器，i从0到1N取得有限个整数值；（2）非因果无限冲激响应（非因果IIR）维纳滤波器，i从到取所有整数值；（3）因果无限冲激响应（因果IIR）维纳滤波器，i从0到取正整数值。上述三种情况下标准方程的解法不同，本文只描述FIR维纳滤波器的求解。设滤波器冲激响应序列的长度为N，冲激响应矢量为TNhhhh)]1()....1()0([(7)滤波器输入数据矢量为TNnxnxnxnx)]1()...1()([)(（8）则滤波器的输出为)()()()(^nxhhnxnsnyTT（9）这样，式（6）所示的维纳-霍夫方程可写成RhPTT或RhP（10）其中)]()([nsnxEP（11）是)(ns与)(nx的互相关函数，它是一个N维列矢量；R是)(nx的自相关函数，是N阶方阵)]()([nxnxERT（12）利用求逆矩阵的方法直接求解式（10），得PRhopt1（13）这里opt表示“最佳”，这就是FIR维纳滤波器的冲激响应。4.FIR维纳滤波器的matlab实现4.1问题描述假设一个点目标在x，y平面上绕单位圆做圆周运动，由于外界干扰，其运动轨迹发生了偏移。其中，x方向的干扰为均值为0，方差为0.05的高斯噪声；y方向干扰为均值为0，方差为0.06的高斯噪声。1)产生满足要求的x方向和y方向随机噪声500个样本；2)明确期望信号和观测信号；3)试设计一FIR维纳滤波器，确定最佳传递函数：1optxxxshRR，并用该滤波器处理观测信号，得到其最佳估计。（注：自行设定误差判定阈值，根据阈值确定滤波器的阶数或传递函数的长度）。4)分别绘制出x方向和y方向的期望信号、噪声信号、观测信号、滤波后信号、最小均方误差信号的曲线图；5)在同一幅图中绘制出期望信号、观测信号和滤波后点目标的运动轨迹。4.2Matlab仿真及运行结果用Matlab实现FIR滤波器，并将先前随机产生的500个样本输入，得到最佳估计。具体程序如下：clear;clf;sita=0:pi/249.5:2*pi;xnoise=sqrt(0.05)*randn(1,500);%产生x轴方向噪声ynoise=sqrt(0.06)*randn(1,500);%产生y轴方向噪声x=cos(sita)+xnoise;%产生x轴方向观测信号y=sin(sita)+ynoise;%产生y轴方向观测信号%产生维纳滤波中x方向上观测信号的自相关矩阵rxx=xcorr(x);fori=1:100forj=1:100mrxx(i,j)=rxx(500-i+j);endendxd=cos(sita);%产生维纳滤波中x方向上观测信号与期望信号的互相关矩阵rxd=xcorr(x,xd);fori=1:100mrxd(i)=rxd(499+i);endhoptx=inv(mrxx)*mrxd';%由维纳-霍夫方程得到的x方向上的滤波器最优解fx=conv(x,hoptx);%滤波后x方向上的输出nx=sum(abs(xd).^2);eminx=nx-mrxd*hoptx;%x方向上最小均方误差%产生维纳滤波中y方向上观测信号的自相关矩阵ryy=xcorr(y);fori=1:100forj=1:100mryy(i,j)=ryy(500-i+j);endendyd=sin(sita);%产生维纳滤波中y方向上观测信号与期望信号的互相关矩阵ryd=xcorr(y,yd);fori=1:100mryd(i)=ryd(499+i);endhopty=inv(mryy)*mryd';%由维纳-霍夫方程得到的y方向上的滤波器最优解fy=conv(y,hopty);%滤波后y方向上的输出ny=sum(abs(yd).^2);eminy=ny-mryd*hopty;%y方向上最小均方误差subplot(2,4,1)plot(xd);title('x方向期望信号');subplot(2,4,2)plot(xnoise);title('x方向噪声信号');subplot(2,4,3)plot(x);title('x方向观测信号');subplot(2,4,4)n=0:500;plot(n,eminx);title('x方向最小均方误差');subplot(2,4,5)plot(yd);title('y方向期望信号');subplot(2,4,6)plot(ynoise);title('y方向噪声信号');subplot(2,4,7)plot(y);title('y方向观测信号');subplot(2,4,8)plot(n,eminy);title('y方向最小均方误差');figure;plot(xd,yd,'k');holdon;plot(x,y,'b:');holdon;plot(fx,fy,'g-');title('最终结果');运行结果如下：图2x方向及y方向的期望信号、噪声信号、观测信号以及滤波后的最小均方误差如上图2所示。图3滤波后的到的信号与原始信号和噪声信号的对比如上图3所示，滤波后的结果与期望信号还是很接近的，整体上达到了最优滤波的效果。参考文献：1.姚天仁孙洪。现代数字信号处理。华中科技大学出版社，20052.崔晓杰，维纳滤波的应用研究，2006