维纳滤波原理及其matlab实现.doc

维纳滤波原理及其matlab实现.doc 摘要 本文介绍了维纳滤波的原理及其matlab实现，以案例的形式展示FIR维纳滤波的特性。 关键字:FIR维纳滤波 Matlab 1.引言 滤波技术是信号分析、处理技术的重要分支，无论是信号的获取、传输，还是信号的处理和交换都离不开滤波技术，它对信号安全可靠和有效灵活地传递是至关重要的。信号分析检测与处理的一个十分重要的就是从噪声中提取信号，实现这种功能的有效手段之一是一种具有最佳线性过滤特性的滤波器，当伴有噪声的信号通过这种滤波器的时候，它可以将信号尽可能精确地重现或对信号做出尽可能精确的估计，而对所伴随噪声进行最大限度地抑制。维纳滤波器就是这种滤波器的典型代表之一。 2.维纳滤波概述 维纳(Wiener)是用来解决从噪声中提取信号的一种过滤(或滤波)。这种线性滤波问，可以看做是一种估计问题或一种线性估计问题。 一个线性系统，如果它的单位样本响应为，当输入一个随机信号，且 h(n)x(n) (1) x(n),s(n)，v(n) 其中表示信号，)表示噪声，则输出为 x(n)v(n)y(n) y(n),h(m)x(n,m) (2) ,m 我们希望通过线性系统后得到的尽量接近于，因此称为的x(n)h(n)y(n)s(n)y(n)s(n) ^ 估计值，用表示，即 s(n) ^ (3) y(n),s(n) 则维纳滤波器的输入—输出关系可用下面图1表示。 图1 实际上，式(2)所示的卷积形式可以理解为从当前和过去的观察值x(n)，x(n,1)， ^ „，„来估计信号的当前值。因此，用进行过滤问题实际上是x(n,2)x(n,m)h(n)s(n) 一种统计估计问题。 一般地，从当前的和过去的观察值x(n)，x(n,1)，x(n,2)„估计当前的信号值 ^ 成为过滤或滤波;从过去的观察值，估计当前的或者将来的信号值y(n),s(n) ^ 称为外推或预测;从过去的观察值，估计过去的信号值y(n),s(n，N)(N,0) ^ 称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预y(n),s(n,N)(N,1) 测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准则的。 ^ 如果我们分别以与表示信号的真实值与估计值，而用表示他们之间的误s(n)e(n)s(n) 差，即 ^ (4) e(n),s(n),s(n) 显然可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来e(n) 表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小: 2 (5) ,(n),E[(n)],mine 采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概 率的描述。 3.维纳-霍夫方程的求解 为了按(5)式所示的最小均方误差准则来确定维纳滤波器的冲激响应，令对h(n),(n) 的导数等于零，即可得 h(j) R(m),h(i)R(m,i),,m (6) ,xsxxi 式中，是与的互相关函数，是的自相关函数，分别定义为 R(m)s(n)x(n)R(m)x(n)xsxx R,E[x(n)s(n，m)]xs R,E[x(n)x(n，m)]xx 式(6)称为维纳滤波器的方程或维纳-霍夫(Wiener-Hopf)方程。如果已知和R(m)xs ，那么解此方程即可求的维纳滤波器的冲激响应。 R(m)xx 式(6)所示标准方程右端的求和范围即的取值范围没有具体标明，实际上有三种情i 况: (1) 有限冲激响应(FIR)维纳滤波器，从到取得有限个整数值; i0N,1 ,,，,(2) 非因果无限冲激响应(非因果IIR)维纳滤波器，从到取所有整数值; i ，,(3) 因果无限冲激响应(因果IIR)维纳滤波器，从到取正整数值。 i0 上述三种情况下标准方程的解法不同，本文只描述FIR维纳滤波器的求解。 设滤波器冲激响应序列的长度为，冲激响应矢量为 N T (7) h,[h(0)h(1)....h(N,1)] 滤波器输入数据矢量为 T (8) x(n),[x(n)x(n,1)...x(n,N，1)] 则滤波器的输出为 ^TT (9) y(n),s(n),x(n)h,hx(n) 这样，式(6)所示的维纳-霍夫方程可写成 TT 或 (10) P,hRP,Rh 其中 (11) P,E[x(n)s(n)] 与的互相关函数，它是一个维列矢量;是的自相关函数，是阶方是Rs(n)x(n)x(n)NN阵 T (12) R,E[x(n)x(n)] 利用求逆矩阵的方法直接求解式(10)，得 ,1 (13) h,RPopt 这里表示“最佳”，这就是FIR维纳滤波器的冲激响应。 opt 4.FIR维纳滤波器的matlab实现 4.1问题描述 假设一个点目标在x，y平面上绕单位圆做圆周运动，由于外界干扰，其运动轨迹发生了偏移。其中，x方向的干扰为均值为0，方差为0.05的高斯噪声;y方向干扰为均值为0，方差为0.06的高斯噪声。 1) 产生满足要求的x方向和y方向随机噪声500个样本; 2) 明确期望信号和观测信号; ,13) 试设计一FIR维纳滤波器，确定最佳传递函数:，并用该滤波器处理观hRR,optxxxs 测信号，得到其最佳估计。(注:自行设定误差判定阈值，根据阈值确定滤波器的阶 数或传递函数的长度)。 4) 分别绘制出x方向和y方向的期望信号、噪声信号、观测信号、滤波后信号、最小均 方误差信号的曲线图; 5) 在同一幅图中绘制出期望信号、观测信号和滤波后点目标的运动轨迹。 4.2 Matlab仿真及运行结果 用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); for i=1:100 for j=1:100 mrxx(i,j)=rxx(500-i+j); end end xd=cos(sita); %产生维纳滤波中x方向上观测信号与期望信号的互相关矩阵 rxd=xcorr(x,xd); for i=1:100 mrxd(i)=rxd(499+i); end hoptx=inv(mrxx)*mrxd';%由维纳-霍夫方程得到的x方向上的滤波器最优解 fx=conv(x,hoptx);%滤波后x方向上的输出 nx=sum(abs(xd).^2); eminx=nx-mrxd*hoptx;%x方向上最小均方误差 %产生维纳滤波中y方向上观测信号的自相关矩阵 ryy=xcorr(y); for i=1:100 for j=1:100 mryy(i,j)=ryy(500-i+j); end end yd=sin(sita); %产生维纳滤波中y方向上观测信号与期望信号的互相关矩阵 ryd=xcorr(y,yd); for i=1:100 mryd(i)=ryd(499+i); end hopty=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'); hold on; plot(x,y,'b:'); hold on; plot(fx,fy,'g-'); title('最终结果'); 运行结果如下: 图2 x方向及y方向的期望信号、噪声信号、观测信号以及滤波后的最小均方误差如上图2 所示。 图3 滤波后的到的信号与原始信号和噪声信号的对比如上图3所示，滤波后的结果与期望信号还是很接近的，整体上达到了最优滤波的效果。 参考文献: 1. 姚天仁 孙洪。现代数字信号处理。华中科技大学出版社，2005 2. 崔晓杰，维纳滤波的应用研究，2006