DSP维纳滤波器的计算机实现

《DSP维纳滤波器的计算机实现》由会员分享，可在线阅读，更多相关《DSP维纳滤波器的计算机实现（20页珍藏版）》请在装配图网上搜索。

1、1、2、3、1 、维纳滤波器的计算机实现实验目的利用计算机编程实现加性噪声信号的维纳滤波。将计算机模拟实验结果与理论分析结果相比较，分析影响维纳 滤波效果的各种因素，从而加深对维纳滤波的理解。利用维纳滤波一步纯预测方法实现对信号生成模型的参数估 计。实验原理维纳滤波器是一种从噪声中提取信号的最佳线性估计方法，假 定一个随机信号形式为：x(n)=s(n)+v(n),其中s(n)为有用信号， v(n)为噪声信号。而维纳滤波的作用就是让x(n)通过一个系统h(n)尽可能滤掉噪声，提取近似s(n),h(n)的选择以最小均方误差 为准则。由维纳-霍夫方程知，只要求出XX及 xs就可求出 h (h= -1

2、xx xs)。但要求h(n)满足因果性要求，维纳-霍夫方 程便是一个难题，这里利用最佳FIR维纳滤波方法求解h(n)的 近似，这也便于在计算机上实现，公式为：h =R-1xx rxs。实验中 s(n)由信 号生成 模型：s(n)=as(n-1)+w(n)确定，其中 a=0.95,w(n)是均值为0,方差为6 w2=1的高斯白噪声，v(n)为 均值为0，方差为1的高斯白噪声，且s(n)与v(n)不相关。实 验中s(n)是已知的，但实际中如果s(n)已知，维纳滤波也就失 去意义了，因此实验纯粹是为了理解维纳滤波原理而设计。2、维纳一步纯预测问题S(n)的生成模型：s(n)+a 1 (n-1)+ +

3、apS(n-p)=w(n), 已知 xx(n)利用Yule-walker方程即可得到信号生成模型参数实验步骤及结果分析1、根据维纳滤波原理绘制程序流程图开始输入样本个数L,FIR滤波器阶数N产生 L 个 v(n),W(n),S(n)和 x(n),在同一坐标内绘出x(n)自相关函数的理论值和实际值在同一坐标内绘出最后100个s(n)和 调矩阵求逆子程序计算，将N个理想的h(n)和 估计的h(n)绘于同一坐标内进行理想的维纳滤波得 L个Si (n),和最后100个s(n)绘制于同一坐标对x(n)进行过滤得L个S (n)，和最后100个根据流程图编写程序(见附录)并分析运行结果：L 个 x(n),S

4、 (n), SI (n), SR(n),统 2 2. 2计 ee-eR选择L=5OOO,N=1。观察并记录、分析实验结果。结束1)与s(n)相比，信号x(n)在维纳滤波前后效果比较:最后100个s(n)(黑色)和x(n)(红色). 411IElD20 3D 4。506070 8D W 100图1为维纳滤波后的s(n)与最后100个s(n)比较图图2为 未经维纳滤波的x(n)与最后100个s(n)比较图。分析：显然与s(n)相比，x(n)在维纳滤波前与s(n)相差很大，维纳滤波后较接近s(n),可见滤波效果比较好。图3为估计？（n）与理想h（n）的对比图。分析：由图可见，二者近似程度除最后几个点

5、外，其他近似度还是满高的，总体而言，近似效果不错。3） 理想的维纳滤波与 FIR维纳滤波效果对比：最后10 （个sn （黑色）和由维纳滤波器得到的sn（红色）的比较最后10C个sn （黑色）和由FIF滤波器得到的sn （红色）的比较 w图4图5图4为理想维纳滤波效果，图5为FIR维纳滤波效果 分 析：直接从图形观察，差异太小，无法观察其精度。只能通过最小均方差来比较其差异，结果为：理想维纳 滤 波ej=0.2287,FIR维纳滤波8=0.2254。可见，理想维纳滤波效果要好过FIR维纳滤波。自相关与互相关数据判断对效果的影响分析：若去掉流程 图中自相关与互相关数据判断步骤，可能会得到理想维纳

6、滤波不如FIR滤波的效果，如其中一个结果：ei= 0.2503, ef= 0.2495。这里的判断步骤就是为了检测实际产生序列的自相关或互相关特性与理论值的近 似程度，若误差很小且通过我们设定的某一下限则 认为二者 近似，所以最终的滤波效果才很近似。如果没 有这里的判 断，实际自相关或互相关则是任意的，完全 有可能出现比理 想维纳滤波更好的效果。3、固定L=5000,分别取N=3、20,根据实验结果，观察N的大小对？(n)的估计和滤波效果的影响并记录实验结果实验结果: 图6为N=3时估计？ (n)与理想h(n)的对比图。图7为N=20时估计？ (n)与理想h(n)的对比图。图8为N=3的FIR

7、滤波后所得S (n)与实际S(n)后100位的比较 图。图9为N=20的FIR滤波后所得S(n)与实际S(n)后100位的比较 图。其均方误差分别为：e 0.3175(N=3), 0.2500(N=20)ef= 0.2762(N=3), 0.2488 (N=20)图7分析：由图6、7可知，N的大小决定？(n)与h(n)取值的个数，并 通过观察并结合N=10的情况可知，N越大？(n)与h(n)越接近。从最终均方误差的比较可知，N越大，滤波效果越好。4、固定N=10，改变L=1000、5000,根据实验结果，观察并记录L的大小对R(n)的精度和滤波效果的影响。实验结果:图10为L=1000时估计？

8、(n)与理想h(n)的对比图。图11为L=5000时估计？(n)与理想h(n)的对比图。图12为L=1000的FIR滤波后所得S(n)与实际S(n)后100位 的比 较图。图13为L=5000的FIR滤波后所得S(n)与实际S(n)后100位 的比 较图。其均方误差分别为：0=0.2400 (L=1000), 0.2381 (L=5000)ef= 0.2390 (L=1000), 0.2375 (L=5000)图10图11开始输入信号生产模型的阶数P, AR模型的参数a(i=1,2P), 6评信号1 s(n)的样本数L利用randn函数产 生s(n)匚L个w(n),并产生L个r利用 Yule-

9、Walker方程，求出自1.？ p最后100Nn (黑色)和由FIF滤波器得到的Sn (红色)的比较w2.5|最后10个sn (黑色)和由F1滤波器得到的srw (红色)的比较2.5-2-2.3jiIIIIII0102030405060 703090 1OD图121.50.5teiB %-0.5-1_1 5J111nr11N0102030405060708090100图13分析：由图10、11可知，L越大？(n)与h(n)越接近，？(n)的精度 越高。由均方误差可知，L越大，滤波效果越高。这也容易理解， 样本越大，精度自然越高。维纳一步纯预测1)画出信号生成模型参数估计的流程图结束2)根据流程

10、图编写程序(见附录2)3）运行信号生成模型程序，选择p=1,a1=-0.6,L=100.理论值：6 w2=1-a12=0.6400 a1=-0.6估计值：%w2= 0.9860 a?1= -0.5876相对误差：error-6 w2= -0.0140error-a1= -0.02064）固定p =1,a1=-0.6, 6 w2=1,改变L=50、500,观察L的大小对信号生成模型参数估计精度的影响。实验结果：理论值： 6 2=1-a12=0.6400a1=-0.6w估计值：L=506?w2= 0.9730a?1= -0.5760估计值：L=500 6?w2= 0.9967 a?1= -0.59

11、65相对误差：L=50 error- 6 w2= -0.0270error-a1= -0.0401相对误差：L=500 error- 6 w2= -0.0033error-a1= -0.0058分析：显然样本个数L的增大，使得信号模型参数精度明显提 高。四、实验总结通过实验结果及分析可得出以下结论：1、样本个数越大，参数精度越高。2、影响维纳滤波效果的因素包括样本个数L、FIR滤波阶数，且均成正比关系。3、维纳一步纯预测，只要调整4(1,2.P)即可实现最小均方误差。五、思考题答案“4 推且分术日“己当士 聆术 口眉卜亡、) s(n) 23 9、于推寸公工式口 心、口由，验垃证工式E 目忙E

12、?( 27推导：已知a= 0.95,w(n)为零均值方差为璧=1-a2的高斯白噪声，v(n)是与s(n)互不相关的高斯白噪声，其均值为零，方差 碱=1。A(z)二厂 m,所有 S(z)=W(z)A(z)= -waz)iI azI azS(z) -az 弋(Z)=W(z)= s(n) =a(n 1) +w(n)打 a=0.95 二 =1 a2 =0.0975x(n) =s(n) +v(n)而且v(n)与s(n)不相关，所有有(m) =Wss(m) +% (m)= xx(z) =ess(z) + vv (z)xx10.952 其中ss如 w (z)A(z)A(z) = (1 _095zW(= 09

13、5Z)xx(Z2331oz+59. -0=:七皿95Z(1-0.95z)(1-0.95Z)曰 2P. 381315化解得 0.95b -2b +0.95 =0= b12, a =1.3122490.723947代入可得xx(Z) =1.312249 Q-V93?470.9。.723947 *1B(Z)B(Z).41.312249,B(Z)二“0B(Z)=，。sx(Z)以(Z)10.952g(z厂(Z)x顷GX|1=GZ2 , G(z )喘0.95Z1(0 723947Z)11 0.952AZ + BZ(1 -0.95Z )(1 -0.723947Z) Z 0.95 Z -0.723947 A

14、=0.312248, B =-0.312248、0.3122481 0.723947Z 1 cccrc,c1 0.723947Z 1g(z) =il zr 中I =0.237948j+I1.312248 b-0.95Z1 -0.72394Z 1-0.95Z10.72394Z 10.723947Zn += 0.95 u(n)+0.723947 u(n1) 1 -0.95Z1 -0.72394Z1 0.723947Z 11 厂 0.237948二 i +=_二 g(z) =Z出-0.95Z1 -0.72394Zjf 1 0.95z1-0.95Z所有 H t(z) =g=O.237948optB(Z

15、) 11 f(不723947 Z SHo(nn= Ee (n) min 旨呵 cds(z)Ho四舍五入取四位小数点就是所求公式不dzPt(z)0xs(z );0.237948 dzZ一 z22 丙c (1 0.90)951 -0.95z) . 1 10.72394z(z0.95 )(z0.723947) dz=-0.0782106 X=0.23795,四舍五入即为所求值0.723947 0.952/. Ee (n)min=-0.0782106X =0.2379500.723947 0.952、由公式 s(n)+ a1(n-1)+ apS(n-p)=w(n),怎样得到印ap和 6w ?分析：理论

16、w(n)已知，即均值及 抚已知，那么根据Yule-Walker方 程有R AY%，其中R为(P +1)*(P+1)的s(n)自相关矩阵危为卜+1) *1的 系数列向量及 A= 1,a1.apT,而& =1,0.0t由给出的理论6：，解方程即可得到估计值a1ap ；用估计值 a1-a代入方程即可得到估计值dW。八源程序见附录1、2、附录1 clear allL=input(请输入信号样本个数L=);N=input(请输入FIR滤波器的阶数N=);a = 0.95;K = 50;sigma_a2 = 1-a 八2;a_ = 1, -a;while(1)%利用randn()函数产生白噪声wn = s

17、qrt(sigma_a2)*( randn(L,1);sn = filter(1, a_, wn);%H(Z)=1/(1-az_T) vn = randn(L,1);xn = sn + vn;r_xx = xcorr(xn,unbiased); %x(n)z 自相关估计值r_xx_t = 2八 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. A O:K;%x(n)与 s(n)互相关理论值%检测实际值

18、与理论值的近似程度 rou_xx = sum(r_xx(L- K:L+K)-r_xx_t).A2)/sum(r_xx_t.八 2);rouxs =(sum(r_xs-r_xs_t)八 2)/sum(r_ xs_t.八 2);if rou_xx 0.03 & rou_xs 0.01 break;end end %同一坐标绘制x(n)自相关函数理论值与实际值 figure(1),clf stem(r_xx(L-K:L+K),r) hold on stem(r_xx_t,k)title(r_xx的实际值(红色)及理论值(黑色)%同一坐标绘制x(n)与s(n)互相关函数理论值与实际值 figure(2

19、),clf stem(r_xs,r) hold on stem(r_xs_t,k)title(r_xs的实际值(红色)及理论值(黑色)%同一坐标绘制最后100个s(n)与x(n) figure(3),clfstem(xn(L-99:L),r) hold on stem(sn(L-99:L),k) titleC 最后 100 个 s(n)(黑色)和x(n)(红色)；%计算h(n)的估计值并与理想值比较构造x(n)的N阶自相关矩阵R_xx n=0:N-1;for i=1:Nfor j=1:NR_xx(i,j)=r_xx(i-j+L);end end hopt=inv(R_xx)*r_xs(1:N)

20、;% 利用维纳一霍夫方程求 h hopt_t=0.238*(0.724ln; %理想h %同一坐标绘制h(n)的实际值与理想 值 figure(4),clf stem(hopt,-,r);hold on stem(hopt_t,*,k);title(h(n)估计值(红色)，与真值(黑色)的比较)；%同一座标绘制理想维纳滤波后的s(n)和最后100个s(n)sn_w=filter(ho pt_t,1,x n);% 理想维纳滤波figure(5),clf stem(sn_w(L-99:L),*,r);hold on stem(sn(L-99:L),.,k) title (最后 100 个 sn (

21、黑色)和 由维纳滤波器得到的sn_w (红色)的 比较);%同一坐标绘制实际维纳滤波后s(n)与最后100个s(n) sn_f=filter(hopt,1,xn);figure(6),clf stem(sn_f(L-99:L),*,r);hold on stem(sn(L-99:L),.,k) titleC 最后 100 个 sn (黑色)和由 FIR滤波器得到的sn_w (红色)的 比较);%求并比较各个均方差e_x=sum(x n(1:L)-s n(1:L)八2)几 e_i=sum(sn_w(1:L)-s n(1:1)八 2)/L e_f=sum(sn_f(1:L)-s n(1:L)A 2

22、)/L附录二 clear all;%输入AR模型的参数p=1；a1=-0.6;sigma_w2=1;L=500;%s(m)的样本数s=1;%统计5000次，最终估计值取平均while s=2 时，s(n)-二 w(n)-a1s(n-1)end r_ss=xcorr(sn,unbiased);%求 s(n) 的 自相关est_a1(s)=-r_ss(L+1)/r_ss(L);%由 r_ss(l)+a1r_ss(l-1)=0 可得 a1 的 估计%由 r_ss(l)+a1r_ss(l-1)=sigma_w2 可得 sigma_w2 的估计值est_sigma(s)=r_ss(L)+est_a1(s)*r_ss(L-1);s=s+1;end %估计值取统计平均值 estimate_a1=mean(est_a1)estimate_sigma=mean(est_sigma) %求估计值与理论的相对误差 error_a1=(estimate_a1-a1)/a1 error_sigma=(estimate_sigma-sigma_w2)/sigma_w2/、21 打(10.95 )x 0.762052 51