实验三-AR模型的参数估计

生物医学信号处理实验指导书 2011-3实验三 AR模型的参数估计一、 设计目的1. 利用维纳预测方法实现对AR模型的参数估计。2. 实现AR模型参数的自适应估计二、 设计原理与方法1. 利用维纳预测方法来估计AR模型的参数实验1中如果已知s(n)，维纳滤波也就没有多少意义了。因此，实验一纯粹是为了理解维纳滤波原理而设计的。下面我们考虑利用维纳预测方法来估计AR模型的参数。 假定s(n)是一个p阶AR模型，即 (3-1) 其中w(n)是均值为零，方差等于的高斯白噪声。在已知准确自相关函数 的情况下，由下面Yule-Walker方程可以得到AR模型参数和 (3-2) 其中为的自相关矩阵，其意义类似于(1-9)式，只是将N换成, 换成而已，A为的系数列向量，定义为 (3-3) 为的单位列向量，除第一个元素等于1外，其余元素均为零，即 (3-4)2. 利用LMS算法实现AR模型参数的估计自适应信号处理方法的应用十分广泛，其中一个非常重要的方面是用来进行参数估计。我们已经知道，如果信号为一个M阶的AR模型，即 (3-5) 通过解Yule-Walker方程可以得到AR模型的参数估计，同样，利用LMS算法，我们也可以对AR模型的参数估计进行自适应估计，其算法如下： (3-6) (3-7) （3-8） 这种算法的实现框图如图3.1所示。 图3.1同样可以证明，只要步长值选择合适，当时，上述自适应算法得到的也收敛于AR模型的参数。三、 设计步骤1. 仔细阅读有关维纳预测原理，弄清AR模型参数估计与维纳预测的关系，根据框图3.2编制AR模型参数估计程序。2. 运行AR模型参数估计程序，选择p=1, =-0.6, L=100, =1, 观察并记录与的最佳估值，与理论值进行比较。3. 固定p=1, =-0.6, =1, 改变L=50,500, 观L的大小对AR模型参数估计精度的影响。(3-1)图3.2 4. 仔细阅读有关自适应系统仿真的内容，按照图3.3给出的框图,编制自适应AR模型参数估计程序(将实验二中自适应滤波程序稍作修改即可)。5. 运行自适应AR模型参数估计程序,选择 观察并记录的收敛情况及和6. 利用100个，通过实验一解Yule-Walker方程的方法，得和估计和，与步骤7中的和比较，有什么差别？为什么？7. 改变噪声的方差，其它条件同步骤7，观察的方差对自适应算法的收敛性，收敛速度以及失调量的影响。(3-5)图3.3 四、 设计报告要求1 简述设计目的和原理。 2 按设计步骤附主要结果。 3 根据结果总结主要结论。 4 如果使用自编程序，附上源程序。 5. 实验感想。 五、 附录 参考程序（1） 利用维纳预测方法实现对AR模型的参数估计。clear all;%输入：AR模型阶数p，AR模型参数ai，i1，.p，dw2，信号s(n)样本数L%L=input(L=);R=input(重复次数=);ac=0;dwc=0;a1=-0.6;dw2=1;for(z=1:R);w=randn(1,L);%利用MATLAB的rand函数产生Lw(n),根据(120)产生L个s(n)%s(1)=w(1);for(i=2:L); s(i)=w(i)-a1*s(i-1);end;%根据公式(1-18)(将x换成s)，由L个s(n)估计p1个自相关函数，解方程(121)得a11，.，app%fss=zeros(1,2); %计算fssfor(i=1:2); for(k=0:(L-i-1); fss(i)=(1/(L-i)*s(k+1)*s(k+i)+fss(i); end;end;for(i=1:2); %生成Rss矩阵 for(k=1:2); Rss(k,i)=fss(abs(k-i)+1); end;end;a11=-(Rss(2,1)/Rss(2,2);dw22=Rss(1,1)+a11*Rss(1,2);temp1(z)=a11;temp2(z)=dw22;ac=a11+ac;dwc=dw22+dwc;end;a1ac=ac/Rdw2dwc=dwc/R（2） 实现AR模型参数的自适应估计%自适应AR滤波clear;m=2;u=0.01;v=0.01; %方差l=100;am=zeros(m,l);cstep=l; %收敛步数p=2;a=zeros(m,1);a(1)=-1.3;a(2)=0.8;w=randn(1,l);meanw=sum(w)/l;covw=sum(w-meanw)*(w-meanw)/l;w=sqrt(v)*(w-meanw)/sqrt(covw);y=zeros(1,l);y(1)=w(1);y(2)=-a(1)*y(1)+w(2);for i=3:l y(i)=-a(1)*y(i-1)-a(2)*y(i-2)+w(i); %y(n)endym=zeros(1,l);e=zeros(1,l);ym(2)=-am(1,2)*y(1);e(2)=y(2)-ym(2);am(1,3)=am(1,2)-2*u*e(2)*y(1);em=0;for i=3:l-1 for j=1:m ym(i)=ym(i)-am(j,i)*y(i-j); end e(i)=y(i)-ym(i); for j=1:m am(j,i+1)=am(j,i)-2*u*e(i)*y(i-j); end if (am(:,i+1)-a)*(am(:,i+1)-a)/(a*a)=0.01 cstep=i;break,end %收敛速度endfor i=cstep+1:l-1 for j=1:m ym(i)=ym(i)-am(j,i)*y(i-j); end e(i)=y(i)-ym(i); em=em+e(i)2; for j=1:m am(j,i+1)=am(j,i)-2*u*e(i)*y(i-j); end endem=em/(l-2-cstep);m=(em-v)/vt=1:l;x=zeros(1,l)+a(1);plot(t,am(1,:),r,t,x,-b);title(a1(