维纳滤波器设计

1.设计要求Sequences(n)ofN=2000pointsisgeneratedbyAR(1)model:s(n)=as(n-1)+w(n),inwhicha=0.8,w(n)iswhitenoisesequence,themeanandvarianceofw(n)is0wm，20.36w.Themeasurementmodelisx(n)=s(n)+v(n),inwhichwhitenoisesequencev(n)andw(n)isnotrelated,themeanandvarianceofv(n)is0vm,21m.Requirements:(1)DesignIIRcausalWienerfilter,calculatethefilteredsequenceandmeansquareerror;(2)DesignFIRWienerfilter,calculatethefilteredsequenceandmeansquareerror;(3)Displayrawdata,noisedataandfiltereddataonthesamegraph,comparethemeansquareerrorbetweenthetwocasesanddrawaconclusion.2.设计原理2.1维纳滤波原理概述维纳（Wiener）是用来解决从噪声中提取信号的一种过滤（或滤波）方法。这种线性滤波问题，可以看做是一种估计问题或一种线性估计问题。一个线性系统，如果它的单位样本响应为)(nh，当输入一个随机信号)(nx，且)()()(nvnsnx（1）其中)(nx表示信号，)(nv)表示噪声，则输出)(ny为mmnxmhny)()()(（2）我们希望)(nx通过线性系统)(nh后得到的)(ny尽量接近于)(ns，因此称)(ny为)(ns的估计值，用^)(ns表示，即^)()(nsny（3）则维纳滤波器的输入—输出关系可用下面图1表示。图1实际上，式（2）所示的卷积形式可以理解为从当前和过去的观察值)(nx，)1(nx，)2(nx…)(mnx，…来估计信号的当前值^)(ns。因此，用)(nh进行过滤问题实际上是一种统计估计问题。一般地，从当前的和过去的观察值)(nx，)1(nx，)2(nx…估计当前的信号值^)()(nsny成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值)0)(()(^NNnsny称为外推或预测；从过去的观察值，估计过去的信号值)1)(()(^NNnsny称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准则的。如果我们分别以)(ns与^)(ns表示信号的真实值与估计值，而用)(ne表示他们之间的误差，即)()()(^nsnsne（4）显然)(ne可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小：min)]([)(2nEne（5）采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概率的描述。2.2维纳-霍夫方程的求解为了按（5）式所示的最小均方误差准则来确定维纳滤波器的冲激响应)(nh，令)(n对)(jh的导数等于零，即可得mimRihmRixxxs,)()()(（6）式中，)(mRxs是)(ns与)(nx的互相关函数，)(mRxx是)(nx的自相关函数，分别定义为)]()([mnsnxERxs)]()([mnxnxERxx式（6）称为维纳滤波器的标准方程或维纳-霍夫（Wiener-Hopf）方程。如果已知)(mRxs和)(mRxx，那么解此方程即可求的维纳滤波器的冲激响应。式（6）所示标准方程右端的求和范围即i的取值范围没有具体标明，实际上有三种情况：（1）有限冲激响应（FIR）维纳滤波器，i从0到1N取得有限个整数值；（2）非因果无限冲激响应（非因果IIR）维纳滤波器，i从到取所有整数值；（3）因果无限冲激响应（因果IIR）维纳滤波器，i从0到取正整数值。上述三种情况下标准方程的解法不同，本文将描述因果IIR维纳滤波器和FIR维纳滤波器的求解。3维纳滤波器的设计与实现3.1因果IIR维纳滤波器的求解解维纳-霍夫方程得，因果IIR维纳滤波器的传输函数为：211sxcSzHzBzBz（7）现要设计一个因果IIR维纳滤波器，对)(nx进行处理，以得最佳估计。为分析简单起见，作如下假设：（1）wn是方差等于Q的白噪声，其自相关函数和功率谱分别为1,0,nininiEwnwiQni（2）vn是方差等于R的白噪声，其自相关函数和功率谱分别为niEvnviR，且wwvvSzQSzR（3）vn与()sn不相关，与wn也不相关，即00EvnsiEvnwi我们需要求出如下几个功率谱：111ssQSzazaz111sxcQSzazaz22111xxssvvcQSzcSzSzRazaz对出入信号功率谱进行因式分解，得12111,,11xxfzSzBzBzBzfaz（8）继续推导得到2222221111fcQaRfaRfcQaR（9）联立两个方程，得22,0aRPQPPRcP（10）求得因果IIR滤波器系统函数211121111111111ccQHzGzGfaBzfzazfzaz（11）其中，维纳增益2cPGRcP（12）21aRfacGRcP（13）3.2FIR维纳滤波器的求解设滤波器冲激响应序列的长度为N，冲激响应矢量为TNhhhh)]1()....1()0([（14）滤波器输入数据矢量为TNnxnxnxnx)]1()...1()([)(（15）则滤波器的输出为)()()()(^nxhhnxnsnyTT（16）这样，式（6）所示的维纳-霍夫方程可写成RhPTT或RhP（17）其中)]()([nsnxEP（18）是)(ns与)(nx的互相关函数，它是一个N维列矢量；R是)(nx的自相关函数，是N阶方阵)]()([nxnxERT（19）利用求逆矩阵的方法直接求解式（10），得PRhopt1（20）这里opt表示“最佳”，这就是FIR维纳滤波器的冲激响应。3.3维纳滤波器的matlab实现由设计要求，可得维纳滤波器的数据模型和测量模型参数：a=0.8,c=1,Q=0.36,R=1;由3.1节，我们可得IIR维纳滤波器的设计步骤如下：（1）根据一直参数a,c,R,Q求解式（10），得到正解P=0.6;（2）由式（12）计算维纳增益G=0.375；（3）由式（13）计算滤波器系数f=0.5；（4）将G和f值带入式（11），得到IIR维纳滤波器的传输函数10.37510.5cHzz，于是单位脉冲响应为()0.375*(0.5)()nhnun。而由3.2节的计算推导，我们可以设计FIR维纳滤波器的程序流程图为：开始产生L个v(n),w(n),s(n),x(n)计算Rxx和Rxs计算h(n)将最后100个x(n),s(n),y(n)绘于同一坐标利用卷积求y(n)结束计算均方误差MSE图1FIR维纳滤波器算法流程图由此，具体Matlab代码如下：clcclearallL=2000;%signallengthN=20;%lengthoftheFIRfiltera=0.8;%whitenoisewithmeanof0andvarof0.36w=sqrt(12*0.36)*(rand(1,L)-0.5);%truesignal:s(n)=a*s(n-1)+w(n)s=zeros(1,L);s(1)=w(1);forii=2:Ls(ii)=a*s(ii-1)+w(ii);end%whitenoisewithmeanof0andvarof1v=sqrt(12)*(rand(1,L)-0.5);%receivedsignal:x=s+vx=s+v;%r_xxistheautocorrelationofxr_xx=xcorr(x);%R_xxistheN-dimentionalautocorrelationmatrixofxR_xx=zeros(N);forii=1:NR_xx(:,ii)=r_xx(L+1-ii:L+1-ii+N-1)';end%r_xsisthecross-correlationofxandsr_xs=xcorr(x,s);r_xs=r_xs(L:L+N-1)';%accordingtoR_xx*h_FIR=R_xs,wecancalculatetheh_FIRh_FIR=inv(R_xx)*r_xs;h_FIR=h_FIR';%y_FIRisthereslutsignalfilteredbyh_FIRy_FIR=cconv(h_FIR,x,L);%y_IIRistheresultsignalfilteredbyIIRcausalWienerfilterh_IIRn=0:L-1;h_IIR=0.375*(0.5.^n);y_IIR=cconv(h_IIR,x,L);%%plotthetruesignals(n),thereceivedsignalx(n),theresultsignaly_FIR(n)%%filteredbyh_FIR,theresultsignaly_IIR(n)filteredbyh_IIRofthelast100pointst=L-99:L;figure(1);subplot(2,2,1);plot(t,x(t),'k');title('Therawdatas(n)');xlabel('n');ylabel('SignalAmplitude');subplot(2,2,2);plot(t,s(t),'r');title('thenoisedatax(n)');xlabel('n');ylabel('SignalAmplitude');subplot(2,2,3);plot(t,y_FIR(t),'b');title('ThedatafilteredbyFIRfilter');xlabel('n');ylabel('SignalAmplitude');subplot(2,2,4);plot(t,y_IIR(t),'g');title('ThedatafilteredbyIIRfilter');xlabel('n');ylabel('SignalAmplitude');figure(2);plot(t,x(t),'--k',t,s(t),'r',t,y_FIR(t),'g',t,y_IIR(t),'b');legend('rawdatas(n)','noisedatax(n)','datafilteredbyFIRfilter','datafilteredbyIIRfilter',0);xlabel('n');ylabel('SignalAmplitude');%计算均方误差e_f=sum((y_FIR-s(1:L)).^2)/Le_i=sum((y_IIR-s(1:L)).^2)/L4.实验结果与结论运行代码，所得原始数据，噪声数据，FIR维纳滤波数据和因果IIR维纳滤波数据如图2所示：图2(a)图2(b)图2实验数据和结果运行结果中，不论是FIR维纳滤波还是因果IIR维纳滤波滤波后的到的信号，与原始信号和噪声信号的对比可以看出，滤波后的结果与期望信号还是很接近的，整体上达到了最优滤波的效果。FIR滤波后得到的信号与因果II