ARMA模型仿真程序.docx

《ARMA模型仿真程序.docx》由会员分享，可在线阅读，更多相关《ARMA模型仿真程序.docx（11页珍藏版）》请在第一文库网上搜索。

1、ARMA模型仿真程序实验%利用总体最小二乘矩估计估计ARMA模型的功率谱%这是一个主函数程序clc,clear;close all;%通过模型参数调用函数zplaneO画出零极点图与直接画出真实的零极点图比较al= 11.352,1.338,-0.662,0.240;a2= 11.352,1.338,-0.662,0.240;a3= 1 ,-2.760,3.809,-2.654,0.924;a4=l,-2.760,3.809,-2.654,0.924;bl=l,-0.200,0.040J;b2= 1,0.900,0.8100;b3= 1,0.900,0.81();b4= 1,0.200,0.0

2、40;cl = l;c2=l;c3=l;c4=l;%数据输入完成figure(l);subplot(2,2,l);zplane(bl,al);title(图一)legendC零点？极点力subplot(2,2,2);zplane(b2,a2);title(图二)subplot(2,2,3)；zplane(b3,a3);title(图三)；subplot(2,2,4)；zplane(b4,a4);title，图四上由ARMA模型绘制完毕fl=().7*exp(j*2*pi*().12),0.7*exp(-j*2*pi*0.12),0.7*exp(j*2*pi*().21),().7*exp(-j

3、*2*pi*0.21)T;f2=0.7*exp(j*2*pi*0.12),0.7*exp(-j*2*pi*0.12),0.7*exp(j*2*pi*0.21 ),0.7*exp(-j*2*pi*0.21)1;f3=0.98 * exp(j*2*pi*0.11 ),0.98*exp(-j*2*pi*0.11 ),0.98*exp(j*2*pi*0.14),0.98*exp(-j*2*pi*0.14)；由二0.98*exp(j*2*pi*0.11 ),0.98*exp(j*2*pi*0.11 ),0.98*exp(j*2*pi*0.14),0.98*exp(j*2*pi*0.14)1,；gl =0

4、.2*exp(j*2*pi*0.17),0.2*exp(-j*2*pi*0.17);g2=0.9*exp(j*2*pi*0.17),0.9*exp(-j*2*pi*0.17);g3=0.9*exp(j*2*pi*0.17),0.9*exp(-j*2*pi*0.17);g4=f0.2*exp(j*2*pi*0.17),0.2*exp(-j*2*pi*().17)l，;figure(2);subplot(2,2,l);zplane(gl,fl);title(图一)legendC零点？极点)subplot(2,2,2);zplane(g2,f2);title (图二)subplot(2,2,3)；z

5、plane(g3,f3);title(图三)；subplot(2,2,4)；zplane(g4,f4);title,图四);由给定的参数真值绘制完毕figure(3)subplot(2,2,l);glp(al,bl,cl);%定义的函数的调用titlef功率谱密度曲线1);hold on;subplot(2,2,2);glp(a2,b2,c2);%定义的函数的调用titleC功率谱密度曲线2);hold on;subplot(2,2,3)；glp(a3,b3,c3);%定义的函数的调用titlef功率谱密度曲线3);hold on;subplot(2,2,4)；glp(a4,b4,c4);%定

6、义的函数的调用titlef功率谱密度曲线C);%由给定的参数真实功率谱函数曲线绘制完成for i=l:50figure(4);subplot(2,2,l);xl(i,:)=pro(al,bl);cl=zxg(xl(i,:);dl=nxg(al,bl,xl(i,:);cl=l,cl;dl=l,dl；glp(cl,dl,l);title(5O次估计功率谱密度曲线I1);hold on;endfor i=l:50figure(4);subplot(2,2,2);xl(i,:)=pro(a2,b2);cl=zxg(xl(i,:);dl=nxg(al,bl,xl(i,:);cl=hcl;dl=l,dl;

7、glp(cLdl,l);title(5O次估计功率谱密度曲线h);hold on;endfor i=l:50figure(4);subplot(2,2,3);xl(i,:)=pro(a3,b3);cl=zxg(xl(i,:);d 1 =nxg(a 1 ,b 1 ,x 1 (i,:);cl=l,cl;dl=l,dl;glp(cl,dl,l);title(5O次估计功率谱密度曲线3);hold on;endfor i=l:50figure(4);subplot(2,2,4);xl(i,:)=pro(a4,b4);cl=zxg(xl(i,:);dl=nxg(al,bl,xl(i,:);cl=l,cl

8、;dl=l,dl；glp(cl,dl,l);出1贝50次估计功率谱密度曲线4);hold on;endfigure(5);subplot(2,2,l);glpp3(al,bl,cl);%定义的函数的调用hold on;titleC比较谱图一上legend。真实功率谱曲线？平均功率谱曲线);subplot(2,2,2);glpp3(a2,b2,c2);%定义的函数的调用hold on;titleC比较谱图二。subplot(2,2,3)；glpp3(a3,b3,c3);%定义的函数的调用hold on;title。比较谱图三上subplot(2,2,4)；glpp3(a4,b4,c4);%定义的

9、函数的调用titleC比较谱图四)附录：自定义的6个脚本函数%第一个自定义函数%产生白噪声激励的模型输出数据function y=pro(a,b)x=rand(l,4)；%产生4个数据u=randn( 1,256);%产生一列含256个数据的白噪声for n=5:256xl=0;for k=2:length(a)x 1 =-a(k)*x(n-k+ l)+x 1;endfor l=l:length(b)x 1 =x 1 +b(l)*u(n-l+1);endx(n)=xl;endy=x;%有ARMA模型控制的输出数据%第二个自定义函数%由ARMA模型真实参数画出真实功率谱function pxp=

10、glp(a,b,c)k=0;w=(0:0.01:0.5)*2*pi;l=O:length(b)-l;wnl 二exp(-j.*w*l);xl=b*wn；yl=abs(xl).A2;k=O:length(a)-l;wn2 = exp(-j. * w*k);x2 = a*wn2;y2=abs(x2).A2;py=yl./y2*c;px=10*logl0(py);px_min=min(px);px_max=max(px);plot(w/(2*pi),px);xlabelC频率 f (Hz) 1);ylabelC功率谱 p(f)(dB),);axis(0,0.5,px_min-20,px_max+20

11、);%第三个自定义函数%有样本数据估计出ARMA模型的AR部分参数function cl=zxg(x)N=length(x);P=N;x=x;for k=0:p-l;xl(k+l)=0;for n=l :N-kx l(k+ l)=x l(k+1 )+x(n)*conj(x(n+k);endr(k+l)=xl(k+l)/N;end %由数据输出数组x的自相关函数for i=l:8for j= 1:4if i+2-j+l=0R(i,j)=conj(conj(r(2);elseR(i,j)=r(i+2-j+l);endendrx(i)=r(i+3);endrx=rx;U,S,V=svd(rx,R);

12、c 1 =V(2,5),V(3,5),V(4,5),V(5,5)/V( 1,5);% 估计 AR 部分参数%第四个自定义函数%产生50次的平均功率谱function B=glpp(a,b,c)for i=l:50xl=pro(a,b);cl=zxg(xl);dl=nxg(a,b,xl);cl=l,cl;dl=l,dl;k=0;w=linspace(0,0.5,50)*2*pi;l=0:length(dl)-l;wn 1 =exp(-j.*w*l);xl=dl*wn；yl=abs(xl).A2;k=():length(cl)-l;wn2=exp(-j.*w*k);x2=cl*wn2;y2=abs

13、(x2).A2;py=yl./y2;N=length(w);A(i,:)=py;endB=0;for i=l:50B=A(i,:)+B;endB=B/50;%第五个自定义函数%估计MA部分的参数function dl=nxg(a,b,x)N=length(x);P=N;x=x;for k=0:p-l;xl(k+l)=0;for n=l:N-kx l(k+ l)=x l(k+1 )+x(n)*conj(x(n+k);endr(k+l)=xl(k+l)/N;end %由数据输出数组x的自相关函数m=50;e=arburg(x,50);e=e(2:end);h=r(l);for i=l:50h=h+

14、e(i)*conj(r(i+1);endP=4;q=2;for k=0:50-lrix(k+l)=0;for l=l:m-krix(k+ l)=rix(k+1 )+conj(e(l)*e(l+k);endrix(k+1 )=rix(k+1)+1 *e(k+1);rix(k+l)=rix(k+l)/h;endfor i=l:10for m=l:2RIX(i,m)=rix(i-m+4+1);endri(i)=rix(i+5);endri=ri*;U,S,V=svd(ri,RIX);dl=V(2,3),V(3,3)/V(l,3)；%第六个自定义函数%对真实功率谱与平均后的功率谱进行比较function pxp=glpp3(a,b,c)k=0;w=linspace(0,0.5,50)*2*pi;l=O:length(b)-l;wn 二exp(j.*w*l);xl=b*wnT;yl=abs(xl).A2;k=O:length(a)-l;wn2=exp(