Burg算法估计模型参数

已知AR(4)过程：x(n)=2.760x(n-1)-3.809x(n-2)+2.654x(n-3)-0.924x(n-4)+w(n)试用Burg算法估计模型参数、阶数和功率谱。实验仿真50次，给出估计结果的方差和偏差（用估计结果均值与真值的差来估计偏差）。用Burg算法对该AR过程估计,50次实验得平均阶数p=4.8，取整数即为P=5，阶数方差var(p)=12.53，阶数偏差为0.8。AR参数（平均）:[-2.06882.2674-1.16920.3671-0.00640.0082-0.01020.00290.0067-0.00550.00090.0007]AR参数方差:AR参数偏差：Burg最大熵谱估计：matlab代码：u=1;whileu=50%50次仿真实验n=10000;w=randn(1,n);%生成零均值，方差为0.2的白噪声w=w-mean(w);w=w/std(w);w=sqrt(0.2)*w;x(1)=1+w(1);x(2)=2+w(2);x(3)=3+w(3);x(4)=4+w(4);%由ar模型生成含噪信号；fori=5:nx(i)=2.760*x(i-1)-3.809*x(i-2)+2.654*x(i-3)-0.924*x(i-4)+w(i);endP(1+1)=1/n*sum(x.^2);%计算预测误差功率的初始值g(1,:)=x;f(1,:)=g(1,:);m=1;k=0;l=0;a(1,1)=1;P(1)=1;whileabs(P(m+1)-P(m))/P(m)=0.01forj=m+1:nk=k+f(m,j)*g(m,j-1);l=l+f(m,j)^2+g(m,j-1)^2;endK(m)=-k/(1/2*l);%求反射系数a(1,1)=1;ifm=2forj=1:m-1%计算前向预测滤波器系数a(m,j)=a(m-1,j)+K(m)*a(m-1,m-j);a(m,m)=K(m);endendP(m+1+1)=(1+abs(K(m))^2)*P(m+1);%计算预测误差功率forj=1:n%计算滤波器输出f(m+1,j)=f(m,j)+K(m)*g(m,n-1);g(m+1,j)=K(m)*f(m,j)+g(m,n-1);endm=m+1;endM(u)=m-1;%记录每次实验的阶数u=u+1;endfori=1:50ar=arburg(x,M(i));%估计AR系数ar_coef(i,1:length(ar))=ar;[Px,fx]=pburg(x,M(i));%估计功率谱Pxx(i,1:length(Px))=10*log10(Px');fxx(i,1:length(fx))=fx';end%%求阶数和AR系数的方差和偏差var_p=var(M);fori=1:max(M)var_ar(i)=var(ar_coef(:,i+1));endbias_p=mean(M)-4;ar_r(max(M))=0;ar_r(1)=-2.76;ar_r(2)=3.809;ar_r(3)=-2.654;ar_r(4)=0.924;fori=1:max(M)mean_ar(i)=mean(ar_coef(:,i+1));bias_ar(i)=mean(ar_coef(:,i+1))-ar_r(i);end%%50次实验功率谱平均fori=1:length(Pxx(1,:))mean_Pxx(i)=mean(Pxx(:,i));end%%作图plot(fxx(1,:),mean_Pxx);%作功率谱图title('burg最大熵估计功率谱');xlabel('频率');ylabel('功率谱密度');figure(2);plot(var_ar);%AR系数估计方差title('AR系数估计方差');xlabel('AR(n)');ylabel('方差');figure(3);plot(bias_ar);%AR系数估计偏差title('AR系数估计偏差');xlabel('AR(n)');ylabel('偏差');disp('阶数方差：');disp(var_p);disp('阶数偏差：');disp(bias_p);disp('AR系数:');disp(mean_ar);