ARMA模型仿真程序

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

2、c1=1;c2=1;c3=1;c4=1;%数据输入完成figure(1);subplot(2,2,1);zplane(b1,a1); title(图一) legend(零点,极点)； 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模型绘制完毕f1=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(-

3、j*2*pi*0.21); 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); 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 ),;f4=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 ),;g1=0.2

4、*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=0.2*exp(j*2*pi*0.17),0.2*exp(-j*2*pi*0.17);figure(2);subplot(2,2,1);zplane(g1,f1);title(图一)legend(零点,极点)；subplot(2,2,2);zplane(g2,f2);title(图二)subplot(2,2,3);zplan

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

6、的调用title(功率谱密度曲线4);%由给定的参数真实功率谱函数曲线绘制完成for i=1:50figure(4);subplot(2,2,1);x1(i,:)=pro(a1,b1);c1=zxg(x1(i,:);d1=nxg(a1,b1,x1(i,:);c1=1,c1;d1=1,d1;glp(c1,d1,1);title(50次估计功率谱密度曲线1);hold on;endfor i=1:50figure(4);subplot(2,2,2);x1(i,:)=pro(a2,b2);c1=zxg(x1(i,:);d1=nxg(a1,b1,x1(i,:);c1=1,c1;d1=1,d1;glp(

7、c1,d1,1);title(,50次估计功率谱密度曲线2);hold on;endfor i=1:50figure(4);subplot(2,2,3);x1(i,:)=pro(a3,b3);c1=zxg(x1(i,:);d1=nxg(a1,b1,x1(i,:);c1=1,c1;d1=1,d1;glp(c1,d1,1);title(50次估计功率谱密度曲线3);hold on;endfor i=1:50figure(4);subplot(2,2,4);x1(i,:)=pro(a4,b4);c1=zxg(x1(i,:);d1=nxg(a1,b1,x1(i,:);c1=1,c1;d1=1,d1;g

8、lp(c1,d1,1);title(50次估计功率谱密度曲线4);hold on;endfigure(5);subplot(2,2,1);glpp3(a1,b1,c1);%定义的函数的调用hold on;title(比较谱图一)；legend(真实功率谱曲线,平均功率谱曲线);subplot(2,2,2);glpp3(a2,b2,c2);%定义的函数的调用hold on;title(比较谱图二)subplot(2,2,3);glpp3(a3,b3,c3);%定义的函数的调用hold on;title(比较谱图三)；subplot(2,2,4);glpp3(a4,b4,c4);%定义的函数的调用

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

10、w=(0:0.01:0.5)*2*pi;l=0:length(b)-1;wn1=exp(-j.*w*l);x1=b*wn1;y1=abs(x1).A2;k=0:length(a)-1;wn2=exp(-j.*w*k);x2=a*wn2;y2=abs(x2).A2;py=y1./y2*c;px=10*log10(py);px_min=min(px);px_max=max(px);plot(w/(2*pi),px);xlabel(频率 f (Hz);ylabel(功率谱 p(f)(dB),);axis(0,0.5,px_min-20,px_max+20);%第三个自定义函数%有样本数据估计出ARM

11、A模型的AR部分参数function c1=zxg(x)N=length(x);p=N;x=x;for k=0:p-1;x1(k+1)=0;for n=1:N-kx1(k+1)=x1(k+1)+x(n)*conj(x(n+k);endr(k+1)=x1(k+1)/N;end %由数据输出数组x的自相关函数for i=1:8for j=1:4if i+2-j+1=0R(i,j)=conj(conj(r(2);elseR(i,j)=r(i+2-j+1);endendrx(i)=r(i+3);endrx=rx;U,S,V=svd(rx,R);c1=V(2,5),V(3,5),V(4,5),V(5,5

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

13、,:)=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+l0;for n=l:N-kx l(k+1 )=x 1 (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 ；for i=l:50h=h+e(i) *conj (r(i+1);endp=4;q=2;for

14、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+ l)=rix(k+1)+1 *e(k+1);rix(k+ l)=rix(k+ l)/h;endfor i=l:10for m=1:2RIX(i,m)=rix(i-m+4+1);endri(i)=rix(i+5);endri=ri;U,S,V=svd(ri,RIX);d1=V(2,3),V(3,3)/V(1,3);%第六个自定义函数%对真实功率谱与平均后的功率谱进行比较function pxp=glpp3(a,b,c)k=0;w=linsp

15、ace(0,0.5,50)*2*pi;l=0:length(b)-1;wn1=exp(-j.*w*l);x1=b*wn1;y1=abs(x1).A2;k=0:length(a)-1;wn2=exp(-j.*w*k);x2=a*wn2;y2=abs(x2).A2;py=y1./y2;N=length(w);B=glpp(a,b,c);py=10*log10(py);plot(w/(2*pi),py,k-);%真实功率谱曲线用黑虚线表示hold on;B=10*log10(B);plot(w/(2*pi),B,b-)%估计的平均功率谱曲线用蓝实线表示 h1=min(min(py),min(B)-1

16、0;h2=max(max(py),max(B)+10;xlabel(频率 f(Hz);ylabel(功率谱 p(f)(dB),);axis(0,0.5,h1,h2);图二1 5 0 5 1 - - 00- +tar VKanuuam1.直接利用ARMA模型的参数画出零级点图图5 0 5 1 00- trnp VKanuuam零点极点50t r aoIP三 ea图,5R0 -11-0.500.5Real Part-1 -0.500.51Real Part图四1 5 0 5 1 - 00- trnp v/an.oam1 5 0 5 1 - - 00- +tar vtan-oam1 .5 0 0.5

17、 0. -.利用零极点的位置画出4个ARMA模型的零极点图1图二1 5 0 5 1 - - 00-图一trar vtan-oam1 5 0 5 100 -trar vtan-oam零点极点10.50-0.5-1-0.500.5Real Part 图四-1 -0.500.51Real Part-1 -0.500.51Real Part1 5 0 5 1 - - 00- trnp VKanuuam.由给定的ARMA模型参数绘制的真实功率谱曲线功率谱密度曲线100-20c(dB5 率功0c(dB5 率功功率谱密度曲线200002 1-10.10.20.30.4频率f ( Hz)功率谱密度曲线30.1

18、0.20.30.4频率f( Hz)功率谱密度曲线40 0 0 P4 2-2 0c(dB5 率功0c(dB5 率功00006 4 20.10.20.30.40.10.20.30.4频率f（ Hz）00.10.20.30.4频率f（ Hz）.利用估计的ARMA模型参数分别画出50次功率谱曲线4050次估计功率谱密度曲线1060000 2 -2 -4 0c(dB5 率功50次估计功率谱密度曲线2000 22- 0c(dB5 率功00.10.20.30.4频率f( Hz)50次估计功率谱密度曲线300.10.20.30.4频率f ( Hz)50次估计功率谱密度曲线400-50c(dB5 率功00 50c(dB5 率功估计的比较谱图二o o o O 2 1 -1 DO-)p5 率功比较谱图真实功率谱曲线平均功率谱曲线00.10.20.30.4估计的比较谱图二o o o O 2 1 -1 DO-)p5 率功比较谱图真实功率谱曲线平均功率谱曲线00.10.20.30.4频率f( Hz)-4000.10.20.30.4频率f( Hz)1000.10.20.30.400.10.20.30.4频率f (