实验题目维纳滤波器的计算机实现

1、实验题目：维纳滤波器的计算机实现专业：姓名：学号：指导老师：2012-11-26维纳滤波器的计算机实现一、实验目的1、利用计算机编程实现加性噪声信号的维纳滤波。2、将计算机模拟实验结果与理论分析结果相比较，分析影响维纳滤波效果的各种因素，从而加深对维纳滤波的理解。3、利用维纳滤波一步纯预测方法实现对信号生成模型的参数估计。二、实验原理1、维纳滤波器是一种从噪声中提取信号的最佳线性估计方法，假定一个随机信号形式为：x(n)=s(n)+v(n), 其中 s(n)为有用信号，v(n)为噪声信号。而维纳滤波的作用就是让x(n)通过一个系统h(n)尽可能滤掉噪声，提取近似 s(n),h(n)的选择以最小

2、均方误差为准则。由维纳 -霍夫方程知，只要求出 xx 及xs 就可求出h（h=-1 xxxs ）。但要求 h(n) 满足因果性要求，维纳 - 霍夫方程便是一个难题， 这里利用最佳FIR 维纳滤波方法求解h(n) 的近似，这也便于在计算机上实现，公式为：h =R-1 xx r xs 。实验中s(n) 由 信 号 生 成 模 型 ： s(n)=as(n-1)+w(n)确 定 ， 其 中2a=0.95,w(n) 是均值为 0，方差为 w =1 的高斯白噪声， v(n) 为均值为 0，方差为 1 的高斯白噪声，且 s(n) 与 v(n) 不相关。实验中 s(n) 是已知的，但实际中如果 s(n) 已知

3、，维纳滤波也就失去意义了，因此实验纯粹是为了理解维纳滤波原理而设计。2、维纳一步纯预测问题S(n) 的生成模型：s(n)+a 1(n-1)+ +aps(n-p)=w(n), 已知 xx（n）,利用 Yule-walker 方程即可得到信号生成模型参数 ai (i=1,2 p) 和2 w 。三、实验步骤及结果分析1、根据维纳滤波原理绘制程序流程图开始输入样本个数L,FIR 滤波器阶数N产生 L个 v(n),w(n),s(n) 和 x(n) ，利用 L 个 s(n)和 x(n) ，估计 RSS和 r xsN检验产生序列 x(n 的自相关和互相关函数Y在同一坐标内绘出x(n) 自相关函数的理论值和实

4、际值在同一坐标内绘出最后100个 s(n) 和调矩阵求逆子程序计算, 将N个理想的 h(n) 和估计的 h(n)绘于同一坐标内进行理想的维纳滤波得L 个 SI (n),和最后 100 个 s(n)绘制于同一坐标对 x(n) 进行过滤得 L个 SR(n) ，和最后 100 个s(n)和绘于同一坐标内L 个 x(n),s(n), S I (n), SR(n),2,e I2,2统计 ex,e R结束2、根据流程图编写程序（见附录1）并分析运行结果：选择 L=5000,N=10 观察并记录、分析实验结果。1）与 s(n)相比，信号 x(n)在维纳滤波前后效果比较：最 后 100个 sn（黑色 ） 和由

5、维纳 滤波 器 得到 的 snw （ 红色 ）的 比 较最 后100个 s(n)（ 黑色 ）和 x(n)（红色 ）251.54130.5201-0.50-1-1-1.5-2-2-3-2.5-401020304050607080901000102030405060708090100图1图2图 1 为维纳滤波后的 s(n)与最后 100 个 s(n)比较图图 2 为未经维纳滤波的 x(n)与最后 100 个 s(n)比较图。分析：显然与 s(n)相比，x(n)在维纳滤波前与 s(n)相差很大，维纳滤波后较接近 s(n),可见滤波效果比较好。2）估计 ?(n)与理想 h(n)的比较：hh(n) 估

6、计 值 （ 红 色 ） ， 与 真 值 （ 黑 色 ） 的 比 较0.250.20.150.10.05012345678910图 3图 3 为估计 ?(n)与理想 h(n)的对比图。h分析：由图可见，二者近似程度除最后几个点外，其他近似度还是满高的，总体而言，近似效果不错。3）理想的维纳滤波与FIR 维纳滤波效果对比：最后100个sn（黑色）和由维纳滤波器得到的sn （红色）的比较最后100个sn（黑色）和由 FIR滤波器得到的sn（红色）的比较2w2w1.51.5110.50.500-0.5-0.5-1-1-1.5-1.5-2-2-2.5-2.501020304050607080901000

7、102030405060708090100图4图5图 4 为理想维纳滤波效果，图 5 为 FIR 维纳滤波效果分析：直接从图形观察，差异太小，无法观察其精度。只能通过最小均方差来比较其差异，结果为：理想维纳滤波 ei= 0.2287,FIR 维纳滤波 ef=0.2254。可见，理想维纳滤波效果要好过 FIR 维纳滤波。4）自相关与互相关数据判断对效果的影响分析：若去掉流程图中自相关与互相关数据判断步骤，可能会得到理想维纳滤波不如FIR 滤波的效果，如其中一个结果： ei= 0.2503，ef= 0.2495。这里的判断步骤就是为了检测实际产生序列的自相关或互相关特性与理论值的近似程度，若误差很

8、小且通过我们设定的某一下限则认为二者近似，所以最终的滤波效果才很近似。如果没有这里的判断，实际自相关或互相关则是任意的，完全有可能出现比理想维纳滤波更好的效果。3、固定 L=5000,分别取 N=3、 20，根据实验结果，观察N 的大小对 ? (n)的估计和滤波效果的影响并记录实验结果。h实验结果：图 6 为 N=3 时估计 ?(n)与理想 h(n)的对比图。h图 7 为 N=20 时估计 ? (n)与理想 h(n)的对比图。h图 8 为 N=3 的 FIR 滤波后所得 s(n)与实际 S(n)后 100 位的比较图。图 9 为 N=20 的 FIR 滤波后所得 s(n)与实际 S(n)后 1

9、00 位的比较图。其均方误差分别为： ei= 0.3175(N=3), 0.2500（N=20）ef = 0.2762(N=3), 0.2488(N=20)0.35h(n)估计 值（ 红 色） ，与 真 值（ 黑色 ）的 比 较0.3h(n)估计值（红 色） ，与 真值 （黑 色） 的比 较0.30.250.250.20.20.150.150.10.10.050.0500-0.0511.21.41.61.822.22.42.62.8302468101214161820图6图7最后100个sn（黑色）和由 FIR滤波器得到的snw（红色）的比较最后 100个sn（黑 色） 和由 FIR滤波 器得

10、 到的 snw （红色）的比较1.5211.50.5100.50-0.5-0.5-1-1-1.5-1.5-2102030405060708090100-210203040506070809010000图 8图 9分析：由图 6、7 可知，N 的大小决定?(n)与取值的个数，并hh(n)通过观察并结合 N=10 的情况可知， N 越大?(n)与 h(n)越接近。h从最终均方误差的比较可知，N 越大，滤波效果越好。4、固定 N=10,改变 L=1000 、5000，根据实验结果，观察并记录 L 的大小对 ? (n)的精度和滤波效果的影响。h实验结果：图 10 为 L=1000 时估计 ? (n)与

11、理想 h(n)的对比图。h图 11 为 L=5000 时估计 ? (n)与理想 h(n)的对比图。h图 12 为 L=1000 的 FIR 滤波后所得 s(n)与实际 S(n)后 100 位的比较图。图 13 为 L=5000 的 FIR 滤波后所得 s(n)与实际 S(n)后 100 位的比较图。其均方误差分别为： ei=0.2400 (L=1000), 0.2381（L=5000）ef = 0.2390 (L=1000), 0.2375(L=5000)h(n)估计值（红色），与真值（黑色）的比较h(n)估计值（红色），与真值（黑色）的比较0.250.250.20.20.150.150.10

12、.10.050.05023456789100234567891011图 10图 11最后 100个sn（黑色）和由FIR滤波器得到的snw（红色）的比较最后100个sn（黑色）和由 FIR滤波器得到的 snw （红色）的比较2.52.5221.51.510.5100.5-0.50-1-0.5-1.5-2-1-2.5102030405060708090100-1.510203040506070809010000图 12图 13?与越接近，?的精度分析：由图 10、11 可知， L 越大 h(n)h (n)h(n)越高。由均方误差可知， L 越大，滤波效果越高。 这也容易理解，样本越大，精度自然越

13、高。5、维纳一步纯预测1）画出信号生成模型参数估计的流程图开始输入信号生产模型的阶数p, AR 模型的参数2, 信号 s(n) 的样本数 Lai (i=1,2 p), w利用 randn 函数产生L 个 w(n), 并产生 L 个s(n)利用 Yule-Walker方程，求出a? 1 . a?p结束2）根据流程图编写程序（见附录2）3）运行信号生成模型程序，选择p=1,a1=-0.6,L=100.理论值：w2=1-a12=0.6400a1=-0.6估计值：? 2= 0.9860a1= -0.5876w?相对误差： error-w 2= -0.0140error-a1= -0.02064)固定

14、p=1,a1=-0.6, w 2=1,改变 L=50、500，观察 L 的大小对信号生成模型参数估计精度的影响。实验结果：理论值： w2=1-a12=0.6400a1=-0.6估计值： L=50? 2a1= -0.5760w = 0.9730?估计值： L=500? 2a1= -0.5965w = 0.9967?相对误差： L=50error-w2= -0.0270error-a1= -0.0401相对误差： L=500 error-w2= -0.0033error-a1= -0.0058分析：显然样本个数L 的增大，使得信号模型参数精度明显提高。四、实验总结通过实验结果及分析可得出以下结论：

15、1、样本个数越大，参数精度越高。2、影响维纳滤波效果的因素包括样本个数L 、 FIR滤波阶数，且均成正比关系。3、维纳一步纯预测， 只要调整ai(1,2 p)即可实现最小均方误差。五、思考题答案1、推导公式0.23792?2H ( Z )1 0.7239 Z 1 ,验证式 Ee(n)E s(n)s(n)0.2379推导：已知 a= 0.95,w(n)为零均值方差为22的高斯白噪声，w1 av(n)是与 s(n)互不相关的高斯白噪声，其均值为零，方差21。v1W (z)A(z)= 1 az1 ,所有 S(z)=W(z)A(z)=1 az 1S(z)1S( z)W (z)s( n) a(n 1)

16、w(n)aza0.95220.0975w1 ax(n)xx (m)其中s(n) v( n)而且 v( n)与 s( n)不相关，所有有ss ( m)vv ( m)xx( z)ss( z)vv( z)12ss(z)10.95, vv(Z ) 1ww( z)A(Z ) A( Z )0.95Z1)(1 0.95Z1(1)xx( Z )10.9521a(1bz 1 )(1bz)(10.95Z1 )(10.95Z1 )(10.95z 1 )(10.95z)化解得 0.95b 22b0.950b1, 21.3813151.312249, a0.723947代入可得xx ( Z)1.312249(10.72

17、3947 Z10.723947 Z )w1 B( Z ) B(Z).)(121(10.95Z1)(10.95Z )11 0.72394Z21.312249, B( Z )1 0.723947 Z, B(Z1，sx (Z )ss (Z )w110.95Z 1)10.95Z21Y(Z )G(Z)11Y(Z) G(Z)X(Z)，G(Z)0.95B(Z )H (Z)B(Z)2(10.95Z1)(1 0.723947Z )X(Z )w110.95 2AZBZ(110.723947 Z)Z0.95Z0.72394710.95Z )(1A0.312248 , B0.312248G(Z)0.31224810.

18、723947 Z0.23794810.723947 Z1.31224810.95Z 110.72394Z0.95Z 110.72394 Z1110.723947 Z0.95n u(n)0.723947 n u (n1)1 0.95Z10.72394 Z10.723947 Z1G(Z )0.23794810.95Z110.72394Z1110.95Z10.95z所有 H opt ( z)G(Z)0.2379481 , 四舍五入取四位小数点就是所求公式B(Z)10.723947 z由H opt ( z)2E e (n )minE e2 (n) min1jc ss (z)H opt ( z)xs (

19、 z 1) dz2z110.952.10.237948dz2jc(10.95 z1)(11 z0.95z)10.72394 z1c(10.952） 0.7620521dz2j0.95.10.723947)(z0.95 )( z0.0782106110.23795，四舍五入即为所求值0.7239470.95E e2 (n )min0.078210610.95 10.2379500.7239472、由公式 s(n)+a1 (n-1)+ a p s(n-p)=w(n),怎样得到 a1 ap 和 w2 ？分析： 理论 w(n) 已知，即均值及w2 已知，那么根据 Yule-Walker方程有RSS A

20、W2, 其中 RSS 为(p+1)*(p+1) 的s(n)自相关矩阵 ,A 为（ p+1）*1 的系数列向量及A= 1,a1 .ap T,而1,0.0T由给出的理论 w2 ，解方程即可得到估计值 a1 a p ；用估计值 a1 a p 代入方程即可得到估计值 w2 。六、源程序见附录1、 2附录 1clear allL=input(' 请输入信号样本个数L=');N=input(' 请输入 FIR 滤波器的阶数 N=');a = 0.95;K=50;sigma_a2 = 1-a2;a_ = 1, -a;while(1)% 利用 randn()函数产生白噪声wn

21、= sqrt(sigma_a2)*( randn(L,1);sn = filter(1, a_, wn);%H(Z)=1/(1-az_-1)vn = randn(L,1);xn = sn + vn;r_xx = xcorr(xn,'unbiased'); %x(n)z 自相关估计值r_xx_t = a.abs(-K:K); % x(n)自相关理论值 r_xx_t(K+1)=r_xx_t(K+1)+1;p = xcorr(sn,xn,'unbiased');%x(n)与 s(n)互相关估计值 r_xs = p(L : L+K);r_xs_t=a.0:K;%x(n)

22、 与 s(n)互相关理论值%检测实际值与理论值的近似程度rou_xx = sum(r_xx(L-K:L+K)-r_xx_t').2)/sum(r_xx_t.2);rou_xs = (sum(r_xs-r_xs_t').2)/sum(r_xs_t.2); if rou_xx < 0.03 & rou_xs < 0.01 break;endend%同一坐标绘制 x(n)自相关函数理论值与实际值figure(1),clfstem(r_xx(L-K:L+K),'r')hold onstem(r_xx_t,'k')title('

23、r_xx 的实际值（红色）及理论值（黑色）')%同一坐标绘制 x(n)与 s(n)互相关函数理论值与实际值figure(2),clfstem(r_xs,'r')hold onstem(r_xs_t,'k')title('r_xs 的实际值（红色）及理论值（黑色）')%同一坐标绘制最后100 个 s(n)与 x(n)figure(3),clfstem(xn(L-99:L),'r')hold onstem(sn(L-99:L),'k')title(' 最后 100 个 s(n)（黑色）和 x(n)（红色

24、） ');%计算 h(n)的估计值并与理想值比较%构造 x(n)的 N 阶自相关矩阵 R_xxn=0:N-1;for i=1:Nfor j=1:NR_xx(i,j)=r_xx(i-j+L);endendhopt=inv(R_xx)*r_xs(1:N);% 利用维纳霍夫方程求hhopt_t=0.238*(0.724).n; %理想 h%同一坐标绘制 h(n)的实际值与理想值figure(4),clfstem(hopt,'-','r');hold onstem(hopt_t,'*','k');title('h(n) 估

25、计值（红色），与真值（黑色）的比较');%同一座标绘制理想维纳滤波后的s(n)和最后 100 个 s(n)sn_w=filter(hopt_t,1,xn);% 理想维纳滤波figure(5),clfstem(sn_w(L-99:L),'*','r');hold onstem(sn(L-99:L),'.','k')title(' 最后 100 个 sn（黑色）和由维纳滤波器得到的sn_w（红色）的比较 ');%同一坐标绘制实际维纳滤波后s(n)与最后 100 个 s(n)sn_f=filter(hopt,1,xn);figure(6),clfstem(sn_f(L-99:L),'*','