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

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=-1x

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

3、x（n）,利用Yule-walker方程即可得到信号生成模型参数ai(i=1,2p)和2w 。三 实验步骤及结果分析1. 仔细阅读维纳滤波原理，根据图 1.1 给出的框图编制维纳滤波程序。(程序见附录)开始 输入样本个数L,FIR滤波器阶数N产生L个v(n),w(n),s(n)和x(n)，利用L个s(n)和x(n)，估计RSS和rxs检验产生序列x(n的自相关和互相关函数是否与理论值相符NY在同一坐标内绘出x(n)自相关函数的理论值和实际值 在同一坐标内绘出最后100个s(n)和x(n)。调矩阵求逆子程序计算,将N个理想的h(n)和估计的h(n) 绘于同一坐标内进行理想的维纳滤波得L个SI (

4、n),和最后100个s(n)绘制于同一坐标 对x(n)进行过滤得L个SR(n)，和最后100个s(n)和绘于同一坐标内结束L个x(n),s(n), SI (n), SR(n),统计ex2,eI2,eR22. 运行维纳滤波程序，选择L=5000,N=10，观察并记录实验结果，分析比较下列三个问题： 与s(n)比较，信号x(n)在维纳滤波前后有何差别？滤波效果如何？（注意：比较噪声方差时应取多次实现的平均值,在本实验中我们统一取100 次实现的平均）可知滤波前后x(n)围绕s(n)的波动比较大，这种变化是由滤波前有很大噪声造成的。滤波后x(n)变得比较光滑，与s(n)更为接近。这说明维纳滤波器对广

5、义平稳输入的滤波效果是相当明显。从波形上可以看出经过维纳滤波获得的信号比原信号有一定的滞后，这是维纳滤波器的因果性造成的。 估计出的h(n)和理想的 h(n)比较，近似程度如何？由图可见，二者近似程度除最后几个点外，其他近似度还是满高的，总体而言，近似效果不错。 理想的维纳滤波和FIR 维纳滤波效果有何差异？从仿真100次后的均方误差值上看可知FIR维纳滤波的均方误差比理想维纳滤波的均方误差非常接近。由于每次产生的信号都不一样，仿真结果会有一些差异，但FIR维纳滤波效果总是逼近与维纳滤波。 若去掉流程图1.1 中判断数据自相关和互相关特性的步骤，可能得出理想维纳滤波效果不如FIR 滤波的结论，

6、思考原因。实际上在估计有限长因果序列h(n)时，由于是利用有限个x(n)和s(n)来估计的，每次仿真x(n)和s(n)也都不同，因此h(n)会有差异。这样我们利用维纳滤波器的得到了信号会于原来的输入信号x(n)误差较大，只是滤波效果很差。3.固定 L=5000, 分别取 N=3、20,根据实验结果，观察 N 的大小对h(n)的估计和滤波效果的影响，记录实验结果。N=3时：N=20时：N的大小决定(n)与h(n)取值的个数，并通过观察并结合N=10的情况可知，N越大(n)与h(n)越接近。从最终均方误差的比较可知，N越大，滤波效果越好。4.固定 N=10, 改变 L=10000,50000,根据

7、实验结果，观察并记录 L 的大小对h(n)的精度和滤波效果的影响。L=10000时：L=50000时：L越大(n)与h(n)越接近，(n)的精度越高。由均方误差可知，L越大，滤波效果越高。这也容易理解，样本越大，精度自然越高。5.仔细阅读有关维纳一步纯预测原理，弄清信号生成模型参数估计与维纳预测的关系，根据框图1.2 编制信号生成模型参数估计程序。（程序见附录）开始输入信号生产模型的阶数p, AR模型的参数ai(i=1,2p),w2,信号s(n)的样本数L利用randn函数产生L个w(n),并产生L个s(n)利用Yule-Walker 方程，求出1.p结束6.运行信号生成模型参数估计程序，选择

8、观察并记录的最佳估值，与理论值进行比较。蓝色表示信号sn，红色表示噪声wn：估计值与理论值之比有一定的误差。估计值与理论值之比有一定的误差。7固定, 改变L=50,500, 观L 的大小对信号生成模型参数估计精度的影响。L=50:L=500:显然样本个数L的增大，使得信号模型参数精度明显提高。四 实验总结通过实验结果及分析可得出以下结论：1.增加输入信号样本个数L和滤波器的阶数N，可以显著提高维纳滤波器的性能。2.维纳一步纯预测，只要调整ai(1,2p)即可实现最小均方误差。五 思考题答案1. 推导公式,验证： 推导：已知a= 0.95,w(n)为零均值方差为的高斯白噪声，v(n)是与s(n)

9、互不相关的高斯白噪声，其均值为零，方差。A(z)=,所有S(z)=W(z)A(z)=2.由公式s(n)+(n-1)+s(n-p)=w(n),怎样得到和？ 分析： 理论w(n)已知，即均值及已知，那么根据Yule-Walker方程有,其中为(p+1)*(p+1)的s(n)自相关矩阵,A为（p+1）*1的系数列向量及A=,而由给出的理论，解方程即可得到估计值；用估计值代入方程即可得到估计值。 附录clc;clear all;close all;L = 5000; % 维纳滤波数据长度 La = 0.95; % w(n)的方差(sigma_w)2K = 50; % 滤波器滤波长度 K N=10;si

10、gma_a2 = 1-a2; % wn的方差a_ = 1, -a; %x(n)d的生成序列while(1) wn = sqrt(sigma_a2)*( randn(L,1); % w(n)是零均值方差为(sigma_w)2的高斯白噪声 sn = filter(1, a_, wn); % s(n)=as(n-1)+w(n) vn = randn(L,1); % v(n)是零均值方差为1的高斯白噪声 xn = sn + vn; r_xx = xcorr(xn,'unbiased'); % 实际的x的自相关函数 r_xx_t = a.abs(-K:K); % 理想的s的自相关函数r_

11、xx_t(K+1)=r_xx_t(K+1)+1; % 理想的x的自相关函数 p = xcorr(sn,xn,'unbiased'); % 实际的x与s的互相关函数 r_xs = p(L : L+K); % 实际的x与s的互相关函数的后k个值rou_xx = sum(r_xx(L-K:L+K)-r_xx_t').2)/sum(r_xx_t.2); rou_xs = (sum(r_xs-a.0:K').2)/sum(a.0:K.2); if rou_xx < 0.03 & rou_xs < 0.01 break; endend figure(9)

12、,clfstem(r_xx(L-K:L+K),'*','r')hold onstem(r_xx_t,'b')figure(1),clfstem(sn(L-99:L),'*','r')hold onstem(xn(L-99:L),'b')legend('最后100个s(n)','最后100个x(n)') % 将h（n）的理想值和实际值进行比较for i=1:N r_xx1(:,i)=r_xx(L-i+1:L-i+N);end hn=inv(r_xx1)*p(L:L+N-

13、1); % h(n)的估计值hn1=0.2379*0.7239.0:N-1; % h(n)的理想值figure(2),clfstem(hn,'*','r')hold onstem(hn1,'b')legend('h(n)的实际值','h(n)的理论值') %将s（n）的理想值和实际值进行比较s=filter(0.2379, 1,-0.7239, xn);figure(3),clfstem(sn(L-99:L),'*','r')hold onstem(s(L-99:L),'b&

14、#39;)legend('最后100个s(n)','理想维纳滤波后的最后100个s(n)') s_1=conv(hn,xn);figure(4),clfstem(sn(L-99:L),'*','r')hold onstem(s(L-99:L),'v','g')hold onstem(s_1(L-99:L),'b')legend('最后100个s(n)','理想维纳滤波后的最后100个s(n)','FIR维纳滤波后的最后100个s(n)'

15、)% legend('最后100个s(n)','FIR维纳滤波后的最后100个s(n)') ex=(1/5000)*sum(xn-sn).2)ei=(1/5000)*sum(s-sn).2)er=(1/5000)*sum(s_1(1:L)-sn).2) %维纳预测sigma_3=1;L1=50;p=1;a=1 -0.6;wn1 = sqrt(sigma_3)*( randn(L1,1);sn1 = filter(1, a, wn1);r_sn1sn1 = xcorr(sn1,'unbiased'); %实际的s1的自相关函数figure(6),clfstem(wn1,'r')hold onstem(sn1,'b') %生成Rfor i=1:p+1 r_ss2(:,i) = r_sn1sn1(L1-i+1:L1-i+p+1);endyipurou=zeros(p+1,1);yipurou(1)=1;A=inv(r_ss2)*yipurou*sigma_3;