实验报告册数字图形图像处理维纳滤波器matlab实现学院：人民武装学院学院专业：计算机科学与技术班级： 11级计科班学号： 1120070544 学生姓名：苏靖指导教师：维纳滤波的原理及其matlab 实现，以案例的形式展示FIR 维纳滤波的特性。

2.维纳滤波概述维纳（Wiener ）是用来解决从噪声中提取信号的一种过滤（或滤波）方法。

这种线性滤波问题，可以看做是一种估计问题或一种线性估计问题。

一个线性系统，如果它的单位样本响应为)(n h ，当输入一个随机信号)(n x ，且)()()(n v n s n x += （1） 其中)(n x 表示信号，)(n v )表示噪声，则输出)(n y 为∑-=mm n x m h n y )()()( （2）我们希望)(n x 通过线性系统)(n h 后得到的)(n y 尽量接近于)(n s ，因此称)(n y 为)(n s 的估计值，用^)(n s 表示，即^)()(n s n y = （3） 则维纳滤波器的输入—输出关系可用下面图1表示。

图1实际上，式（2）所示的卷积形式可以理解为从当前和过去的观察值)(n x ，)1(-n x ，)2(-n x …)(m n x -，…来估计信号的当前值^)(n s 。

因此，用)(n h 进行过滤问题实际上是一种统计估计问题。

一般地，从当前的和过去的观察值)(n x ，)1(-n x ，)2(-n x …估计当前的信号值^)()(n s n y =成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值)0)(()(^≥+=N N n s n y 称为外推或预测；从过去的观察值，估计过去的信号值)1)(()(^>-=N N n s n y 称为平滑或内插。

因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。

这里所谓的最佳与最优是以最小均方误差为准则的。

如果我们分别以)(n s 与^)(n s 表示信号的真实值与估计值，而用)(n e 表示他们之间的误差，即)()()(^n s n s n e -= （4） 显然)(n e 可能是正值，也可能是负值，并且它是一个随机变量。

因此，用它的均方误差来表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小：min )]([)(2==n E n e ξ （5） 采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概率的描述。

4.FIR 维纳滤波器的matlab 实现4.1问题描述假设一个点目标在x ，y 平面上绕单位圆做圆周运动，由于外界干扰，其运动轨迹发生了偏移。

其中，x 方向的干扰为均值为0，方差为0.05的高斯噪声；y 方向干扰为均值为0，方差为0.06的高斯噪声。

1) 产生满足要求的x 方向和y 方向随机噪声500个样本；2) 明确期望信号和观测信号；3) 试设计一FIR 维纳滤波器，确定最佳传递函数：1opt xx xs h R R -=，并用该滤波器处理观测信号，得到其最佳估计。

（注：自行设定误差判定阈值，根据阈值确定滤波器的阶数或传递函数的长度）。

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:100for j=1:100mrxx(i,j)=rxx(500-i+j);endend%产生维纳滤波中x方向上观测信号与期望信号的互相关矩阵rxd=xcorr(x,xd);for i=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);for i=1:100for j=1:100mryy(i,j)=ryy(500-i+j);endendyd=sin(sita);%产生维纳滤波中y方向上观测信号与期望信号的互相关矩阵ryd=xcorr(y,yd);for i=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方向期望信号');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('最终结果');运行结果如下：图2x方向及y方向的期望信号、噪声信号、观测信号以及滤波后的最小均方误差如上图2所示。

图3滤波后的到的信号与原始信号和噪声信号的对比如上图3所示，滤波后的结果与期望信号还是很接近的，整体上达到了最优滤波的效果。

维纳滤波器的局限维纳滤波复原法存在着几个实质性的局限。

第一，最有标准是基于最小均方误差的且对所有误差等权处理，这个标准在数学上可以接受，但却是个不适合人眼的方式，原因在于人类对复原错误的感知在具有一致灰度和亮度的区域中更为严重，而对于出现在暗的和高梯度区域的误差敏感性差得多。

第二，空间可变的退化不能用维纳滤波复原法复原，而这样的退化是常见的。

第三，维纳滤波不能处理非平稳信号和噪声。

四、模拟仿真运行结果运行程序代码clear;I=imread('img_orignal.tif');figure;subplot(2,2,1);imshow(I);title('原图像');[m,n]=size(I);F=fftshift(fft2(I));k=0.005;for u=1:mfor v=1:nH(u,v)=exp((-k)*(((u-m/2)^2+(v-n/2)^2)^(5/6)));endendG=F.*H;I0=real(ifft2(fftshift(G)));I1=imnoise(uint8(I0),'gaussian',0,0.001)subplot(2,2,2);imshow(uint8(I1));title('模糊退化且添加高斯噪声的图像');F0=fftshift(fft2(I1));K=0.1;for u=1:mfor v=1:nH(u,v)=exp(-k*(((u-m/2)^2+(v-n/2)^2)^(5/6)));H0(u,v)=(abs(H(u,v)))^2;H1(u,v)=H0(u,v)/(H(u,v)*(H0(u,v)+K));endendF2=H1.*F0;I2=ifft2(fftshift(F2));subplot(2,2,3);imshow(uint8(I2));title('维纳滤波复原图');五、结论与心得体会通过这个实验，使我们更加深刻和具体的了解到了维纳滤波的原理，功能以及在图像处理方面的应用。

维纳滤波器是对噪声背景下的信号进行估计，它是最小均方误差准则下的最佳线性滤波器。

在实验的过程中，我发现采用维纳滤波复原可以得到比较好的效果，这个算法可以使估计的点扩散函数值更加接近它的真实值。

但实现维纳滤波的要求是①输入过程是广义平稳的；②输入过程的统计特性是已知的。

根据其他最佳准则的滤波器也有同样的要求。

然而，由于输入过程取决与外界信号，干扰环境，这种环境的统计特性常常是未知的，变化的，因而这两个要求很难满足，这就促使人们研究自适应滤波器。

附：维纳滤波器的设计方法维纳-霍夫方程维纳滤波器的设计，实际上就是在最小均方误差条件下探索和确定滤波器的冲激函数h(n)或系统函数H （z ），也就是求解维纳-霍夫方程的问题。

对于物理可实现系统，由（1）式得∑∞=-==0)()()(ˆ)(y m m n x m h n sn 式（6） 它实现的是将当前的及过去的诸输入值作相应的加权后的求和运算。

故维纳滤波的设计则是确定均方误差}])()()({[)]([E 202∑∞=--=m m n x m h n s E n e 式（7） 最小意义下的冲激响应h(n)。

为便于得出矩阵表达式，我们将（6）式改写成i x ∑∞==1i i h (n)sˆ 式（8） 式中1()()1()(1m ,1i +-=-=-==-=+=i n x m n x x i h m h h i m i i 或 式（9）因此])[()]([E 212∑∞=-=i i i x h s E n e 式（10） 为求得])[(21∑∞=-i ii x h s E 最小时的{h i }，我们将（10）式对h i 求偏导，得1},x )]([2{)]([E i 22112≥-⋯++-=∂∂i x h x h s E h n e i）（ 式（11） 再令其为零，即1,0])[(E 1≥=-∑∞=i x x h s i i i i或，1,0][E ≥=j ex j 式（12）从而可以确定我们所需要的{h i }。