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

维纳滤波器的计算机实现一、 实验目的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=-1xxxs）。但要求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),已知xx（n）,利用Yule-walker方程即可得到信号生成模型参数ai(i=1,2p)和2w 。三、 实验步骤及结果分析1、 根据维纳滤波原理绘制程序流程图开始 输入样本个数L,FIR滤波器阶数N产生L个v(n),w(n),s(n)和x(n)，利用L个s(n)和x(n)，估计RSS和rxs检验产生序列x(n的自相关和互相关函数是否与理论值相符 N在同一坐标内绘出x(n)自相关函数的理论值和实际值 Y在同一坐标内绘出最后100个s(n)和x(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), SI (n), SR(n),统计ex2,eI2,eR22、 根据流程图编写程序（见附录1）并分析运行结果：选择L=5000,N=10观察并记录、分析实验结果。1） 与s(n)相比，信号x(n)在维纳滤波前后效果比较： 图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)的比较： 图3图3为估计(n)与理想h(n)的对比图。分析：由图可见，二者近似程度除最后几个点外，其他近似度还是满高的，总体而言，近似效果不错。3） 理想的维纳滤波与FIR维纳滤波效果对比： 图 4 图5图4为理想维纳滤波效果，图5为FIR维纳滤波效果分析：直接从图形观察，差异太小，无法观察其精度。只能通过最小均方差来比较其差异，结果为：理想维纳滤波ei= 0.2287,FIR维纳滤波ef=0.2254。可见，理想维纳滤波效果要好过FIR维纳滤波。4） 自相关与互相关数据判断对效果的影响分析：若去掉流程图中自相关与互相关数据判断步骤，可能会得到理想维纳滤波不如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滤波后所得(n)与实际S(n)后100位的比较图。图9为N=20的FIR滤波后所得(n)与实际S(n)后100位的比较图。其均方误差分别为：ei= 0.3175(N=3), 0.2500（N=20） ef= 0.2762(N=3), 0.2488 (N=20)图6 图7 图8 图9分析：由图6、7可知，N的大小决定(n)与h(n)取值的个数，并通过观察并结合N=10的情况可知，N越大(n)与h(n)越接近。从最终均方误差的比较可知，N越大，滤波效果越好。4、 固定N=10,改变L=1000、5000，根据实验结果，观察并记录L的大小对(n)的精度和滤波效果的影响。实验结果：图10为L=1000时估计(n)与理想h(n)的对比图。图11为L=5000时估计(n)与理想h(n)的对比图。图12为L=1000的FIR滤波后所得(n)与实际S(n)后100位的比较图。图13为L=5000的FIR滤波后所得(n)与实际S(n)后100位的比较图。其均方误差分别为：ei=0.2400 (L=1000), 0.2381（L=5000） ef= 0.2390 (L=1000), 0.2375 (L=5000) 图10 图11 图12 图13分析：由图10、11可知，L越大(n)与h(n)越接近，(n)的精度越高。由均方误差可知，L越大，滤波效果越高。这也容易理解，样本越大，精度自然越高。5、 维纳一步纯预测 1） 画出信号生成模型参数估计的流程图开始输入信号生产模型的阶数p, AR模型的参数ai(i=1,2p),w2,信号s(n)的样本数L利用randn函数产生L个w(n),并产生L个s(n)利用Yule-Walker 方程，求出1.p结束2）根据流程图编写程序（见附录2）3）运行信号生成模型程序，选择p=1,a1=-0.6,L=100. 理论值： w2=1-a12=0.6400 a1=-0.6 估计值： w2= 0.9860 1= -0.5876 相对误差：error-w2= -0.0140 error-a1= -0.02064)固定p=1,a1=-0.6, w2=1,改变L=50、500，观察L的大小对信号生成模型参数估计精度的影响。实验结果：理论值： w2=1-a12=0.6400 a1=-0.6估计值：L=50 w2= 0.9730 1= -0.5760 估计值：L=500 w2= 0.9967 1= -0.5965相对误差：L=50 error-w2= -0.0270 error-a1= -0.0401相对误差：L=500 error-w2= -0.0033 error-a1= -0.0058分析：显然样本个数L的增大，使得信号模型参数精度明显提高。四、 实验总结通过实验结果及分析可得出以下结论：1、 样本个数越大，参数精度越高。2、 影响维纳滤波效果的因素包括样本个数L、FIR滤波阶数，且均成正比关系。3、 维纳一步纯预测，只要调整ai(1,2p)即可实现最小均方误差。五、 思考题答案1、推导公式,验证式 推导：已知a= 0.95,w(n)为零均值方差为的高斯白噪声，v(n)是与s(n)互不相关的高斯白噪声，其均值为零，方差。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=,而由给出的理论，解方程即可得到估计值；用估计值代入方程即可得到估计值。 六、 源程序见附录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 = 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)与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.01break;endend%同一坐标绘制x(n)自相关函数理论值与实际值figure(1),clfstem(r_xx(L-K:L+K),r) hold onstem(r_xx_t,k) title(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)（红色）);%计算h(n)的估计值并与理想值比较%构造x(n)的N阶自相关矩阵R_xxn=0:N-1;for i=1:N for j=1:N R_xx(i,j)=r_xx(i-j+L); endendhopt=inv(R_xx)*r_xs(1:N);%利用维纳霍夫方程求h hopt_t=0.238*(0.724).n; %理想h %同一坐标绘制h(n)的实际值与理想值figure(4),clfstem(hopt,-,r);hold onstem(hopt_t,*,k);title(h(n)估计值（红色），与真值（黑色）的比较);%同一座标绘制理想维纳滤波后的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),*,r);hold onstem(sn(L-99:L),.,k)title(最后100个sn（黑色）和由FIR滤波器得到的sn_w（红色）的比较);%求并比较各个均方差e_x=sum(xn(1:L)-sn(1:L).2)/L e_i=sum(sn_w(1:L)-sn(1:L).2)/Le_f=sum(sn_f(1:L)-sn(1:L).2)/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)endr_ss=xcorr(sn,unbiased);%求s(n)的自相关est_a1(s)=-r_ss(L+1)/r_ss(L);%由r_ss(l)+a