维纳滤波(带程序)

维纳滤波器的计算机实现专业：信息与通信工程实验一维纳滤波器的计算机实现一、实验目的1．MATLAB编程实现加性干扰信号的维纳滤波。2．仿真比较实验结果与理论分析结果，分析影响维纳滤波效果的各种因素，从而加深对维纳滤波的理解。3．利用维纳预测方法实现对AR模型的参数估计。二、实验原理及方法维纳滤波实际上就是在最小均方误差条件下探索和确定滤波器的冲激响应h(n)或系统函数H(z)，也就是求解维纳－霍夫方程。假定一个随机信号x(n)具有以下形式：x(n)=s(n)+v(n)其中，s(n)为有用信号，v(n)为噪声干扰，将其输入一个单位脉冲响应为h(n)的线性系统，则其输出为：()()()()mynsnhmxnm我们希望x(n)在经过系统h(n)后得到y(n)，即s(n)的估计值能尽可能接近s(n)，按照最小均方误差准则，h(n)应满足下面的正则方程：()()()xsxxmkhmkm其中，()xsk是x(n)与s(n)的互相关函数，()xxm是x(n)的自相关函数。在h(n)满足因果性的条件下，，求解维纳-霍夫方程是一个典型的难题。虽然目前有几种求解h(n)的解析方法，但它们在计算机上实现起来非常困难。因此，本实验中，利用近似方法，即最佳FIR维纳滤波方法，在计算机上实现随机信号的维纳滤波。设h(n)为一因果序列，其长度为N，则10()()()Nmynhmxnm同样利用最小均方误差准则，h(n)满足下面正则方程：xxxsRhrxxR为信号x(n)的N阶自相关矩阵，xsr为x(n)与s(n)的互相关函数向量。当xxR为满秩矩阵时，可得1xxxshRr由此可见，利用有限长的h(n)实现维纳滤波器，只要已知xxR和xsr就可以按上式解得满足因果性的h。虽然它不同于真正的维纳滤波器，但是只要N选择的足够大，它就可以很好地逼近真正的维纳滤波器。利用维纳预测方法来估计AR模型的参数。即假定s(n)是一个p阶AR模型。1()(1)()()psnasnasnpwn其中w(n)是均值为零，方差等于2w的高斯白噪声。在已知准确自相关函数()ssn的情况下，由下面Yule-Walker方程可以得到AR模型参数ia(i=1，…p)和2w2sswRAA为(p+1)╳1的系数列向量，定义为1[1,,]TpAaaε为(p+1)×1的单位列向量，除第一个元素等于1外，其余元素均为零，即[1,0,0]T三、仿真结果及分析（1）、L=500，N=10；运行维纳滤波程序得到如下结果：图1100个是S(n)与X(n)比较图2图3其误差为：不经滤波：TheE_xis0.98152理伦滤波：TheE_Iis0.77672维纳滤波：TheE_Ris0.809841、通过比较图1与图3可以看出在滤波后，估计值更近似等于S(n)，其数据点围绕s(n)的曲线做很小的扰动，其值在s(n)的上下波动。2、因为只用了有限个数据S(n)和X(n)来求自相关函数和互相关函数，由此得到的估计值不精确，从而带来了实验误差。（2）、固定L=500，N=3与N=20的比较结果如下：N=3N=20N=3N=20N=3N=20N=3N=20分析比较可发现，L=500固定不变，N值越大，误差越小，滤波后的信号更加接近原来的信号，即滤波效果得到改善。（3）、固定N=10，L=200与L=1000比较：L=200L=1000L=200L=1000L=200L=1000L=200L=1000当N保持不变，而将数据长度L从200增大到1000，可以发现，均方误差都有增大，而且h(n)其理论值与实验值的差别也变大，这说明单方面增加数据长度将会降低滤波效果。（4）、AR模型参数估计：1、P=2，a=[0.5,0.6],sigma2=1,L=5002、P=2，a=[0.5,0.6],sigma2=1,L=50从上述两幅图中可以看出样本个数L取的越多对参数a及噪声功率估计的越准确。（5）、程序代码：维纳滤波：clear;L=input('PleaseinputthevalueofL:L=');N=input('PleaseinputthevalueofN:N=');a=0.95;w(1)=0;w=rand(1,L);w=w-mean(w)*ones(1,L);w=w*sqrt(12*(1-a^2));s(1)=1;forn1=2:Ls(n1)=a*s(n1-1)+w(n1-1);endv=rand(1,L);v=v-mean(v)*ones(1,L);v=v*12^(1/2);forn2=1:Lx(n2)=s(n2)+v(n2);endm=100;j=(L-m+1):1:L;sj=s((L-m+1):1:L);xj=x((L-m+1):1:L);figure(1)exp=zeros(1,m);plot(j,sj,'or',j,xj,'*b',j,exp,'-.k');title('100samplesofs(n)andx(n)');xlabel('variablei');ylabel('s(i)andx(i)');legend('s(i)','x(i)');figure(2)PHIxx=xcorr(x,x,'unbiased');PHIxx=PHIxx(L:(L+m-1));Rxx=toeplitz(PHIxx);PHIxs=xcorr(x,s,'unbiased');rxs=PHIxs(L:(L+m-1));inv_Rxx=inv(Rxx);rxs=rxs';h=inv_Rxx*rxs;h=h';formm=0:(m-1)hd(mm+1)=0.238*(0.724)^mm;endi=0:(m-1);plot(i,h,'or',i,hd,'xb');title('估计h(n)与理想hd(n)比较');xlabel('数值i');ylabel('h(i)andhd(i)');legend('h(i)','hd(i)');gridon;EST_sI=circonv(hd,x,L);EST_sIj=EST_sI((L-m+1):1:L);figure(3)subplot(211)plot(j,sj,'.:r',j,EST_sIj,'x-b');title('真实值s(n)与通过理想滤波器估计值s(n)的比较');xlabel('数值i');ylabel('s(i)与EST_sI(i)');legend('s(i)','EST_sI(i)');gridon;subplot(212)EST_s=circonv(h,x,L);EST_sj=EST_s((L-m+1):1:L);plot(j,sj,'.:r',j,EST_sj,'+-g');title('真实值s(n)and通过维纳滤波器估计值s(n)的比较');xlabel('数值i');ylabel('s(i)与EST_s(i)');legend('s(i)','EST_s(i)');gridon;E_x=sum((x-s).^2)/L;E_I=sum((x-EST_sI).^2)/L;E_R=sum((x-EST_s).^2)/L;fprintf('TheE_xis%8.5f\n',E_x);fprintf('TheE_Iis%8.5f\n',E_I);fprintf('TheE_Ris%8.5f\n',E_R);functiony=circonv(x,h,N)x=[x,zeros(1,N-length(x))];h=[h,zeros(1,N-length(h))];m=[0:N-1];hm=h(mod(-m,N)+1);H=toeplitz(hm,[0,h(2:N)]);y=x*H;AR模型：clear;order=input('请输入模型阶数：P=');a=input('请输入系数变量：a=');if~(order==length(a))dips('thelengthofcoefficientsvectorshouldbeequaltoorder!')endvar=input('请输入噪声功率：var=');L=input('请输入样本个数：L=');s=est_AR(order,a,var,L);PHIss=xcorr(s,s,'unbiased');PHIssp=PHIss(L:(L+order));Rss=toeplitz(PHIssp);Rss_l=Rss((2:order+1),(2:(order+1)));Rss_r=Rss((2:order+1),1);Rss_r=(-1)*Rss_r;est_a=inv(Rss_l)*Rss_r;est_var=Rss(1,:)*[1;est_a];est_a=est_a';est_var=est_var';fprintf('\n');fprintf('结果如下：\n');fprintf('a的估计值为：%8.5f\n',est_a);E_a=((a-est_a).^2)/order;fprintf('误差为：%8.5f\n',E_a);fprintf('\n');fprintf('白噪声功率估计为：%8.5f\n',est_var);E_var=sum(var-est_var).^2;fprintf('误差为：%8.5f\n',E_var);functions=est_AR(order,a,var,L)order=2;var=1;L=500;a=[0.5,0.6];a=a';w=randn(1,L);w=w*(var^(1/2));forn1=1:orders(n1)=w(n1);endforn2=(n1+1):Lsj=[s(n2-1),s(n2-order)];s(n2)=(-1)*(sj*a)+w(n2);endend