ARMA模型仿真程序.docx

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

2、l = l;c2=l;c3=l;c4=I;%数据输入完成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 二0.7*exp(j*2*pi*0.12),0.7*exp(-j*2*pi*0.12),0.7*ex

3、p(j*2*pi*0.21),0.7*exp(-j*2*pi*0.21)； f2=O.7*exp(j*2*pi*O12),O.7*exp(-j*2*pi*O12)07*exp(j*2*pi*O21),O.7*exp(-j*2*pi*O.21)'f3=0.98*exp(j*2*pi*0ll),0.98*exp(j*2*pi*0ll),0.98*exp(j*2*pi*0.14),0.98*exp(j*2*pi*0.14)'；f4=0.98*exp(j*2*pi*0.1 l),0.98*exp(-j*2*pi*0.1 l),0.98*exp(j*2*pi*0.14),0.98*exp

4、(-j*2*pi*0.14)'；gl=0.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)'50次估计功率谱密度曲线14000.10.20.30.4频率f( Hz)50次估计功率谱密度曲线3o o o o O2 2 4 6 _ _ _ mPMdm豆母50次估计功率谱密度曲线200.10.20.30.4频率f( Hz)50次估计功率谱密度曲线4o o

5、 O2 2 (cop£dss>500.10.20.30.4频率f( Hz)(cnpgd®鄙择00.10.20.30.4频率f( Hz)(CDpedffl阱抑5.最后估计的50次平均功率谱与真实功率谱比较，真实功率谱曲线用黑虚线表示，估计的 平均功率谱曲线用蓝实线表示。比较谱图一比较谱图二00.10.20.3频率f( Hz) 比较谱图三0.20.30.4Q(U6040200-2000.20.30.4频率f( Hz)-406040200-20频率f( Hz)真实功率谱曲线 顼率f( Hz) 平均功率谱曲线比较谱图四0.20.30.406.经测试当改b3=l,0.800,0

6、.900时，第三个ARMA模型估计的效果很好，见下图(COP)(一)d 热耕旨-40-20-3000.050.10.150.20.250.30.350.40.450.5频率f( Hz)50403020100107.实验分析经过对程序的测试，发现4个ARMA模型AR部分的参数估计的很精确，能精确到0.001,效果很好。但是MA部分的参数波动较大，估计的效果不是很好，误差较大。我认为，AR部分的参数估计就是涉及到先根据样本数据估计自相关函数，然后再利用奇异值分 解求解总体最小二乘矩估计问题，涉及的误差积累较少，所以AR部分的参数估计较精确。 而对MA部分的参数估计，先是通过ARMA模型与高阶AR模

7、型的等价性质，我们选的是 50阶的等价AR模型，本身有一定的不合理，再估计高阶AR模型的参数用的是burg法， 再利用高阶AR模型的参数估计逆自相关函数，然后我们可得C-P方程，再对C-P方程求解 总体最小二乘矩估计，也是利用奇异值分解求得MA部分的参数估计，这一求解过程涉及 估计的量有很多，估计中误差积累较多，所以MA部分估计的参数波动大是必然的，我们 需要改进的就是减少估计中间变量,这样能减少误差。对ARMA模型参数的适当改进，使 模型稳定，由模型产生的样本数据估计的AR部分和MA部分的参数更精确，由它们画出的 功率谱与真实功率谱拟合的更好，比如第6幅图拟合的效果就很好。g4=0.2*ex

8、p(j*2*pi*0.17),0.2*exp(-j*2*pi*0.17)' figure(2);subplot(2,2,l);zplane(gl,fl);title('图一') legendC零点？极点)； subplot(2,2,2);zplane(g2,f2);title('图二')subplot(2,2,3)； zplane(g3,f3);title('图三')；subplot(2,2,4);zplane(g4,f4);title。图四);由给定的参数真值绘制完毕figure(3)subplot(2,2,l);glp(al,bl,c

9、l);%定义的函数的调用 title(，功率谱密度曲线)； 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);%定义的函数的调用title(功率谱密度曲线4);%由给定的参数真实功率谱函数曲线绘制完成for i=l:50figure(4);subplot(2,2,l);xl(i,:)=pro(al,bl);cl=z

10、xg(xl(i,:);dl=nxg(al,bl,xl(i,:);cl=l,cl;dl=l,dl;glp(cl,dl,l);title(f5O次估计功率谱密度曲线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=l,cl; dl=l,dl;glp(cl,dl,l);title(*5O次估计功率谱密度曲线2); hold on;endfor i=l:50figure(4);subplot(2,2,3)； xl(i,:)=pro(a3

11、,b3); cl=zxg(xl(i,:); dl=nxg(al,bl,xl(i,:); cl=l,cl;dl=l,dl;glp(cl,dl,l);titleCO次估计功率谱密度曲线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;dl=l,dl；glp(cl,dl,l);title(?5O次估计功率谱密度曲线4) hold on;endfigure(5);subplot(2,2,l);glpp3(al,bl,cl);%定

12、义的函数的调用 hold on;title(此较谱图一);legend。真实功率谱曲线7平均功率谱曲线 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);%定义的函数的调用 title(比较谱图四)附录：自定义的6个脚本函数%第一个自定义函数%产生白噪声激励的模型输出数据function y=pro(a,b) x=rand(l,4)；%产生

13、4个数据u=randn( 1,256);%产生一列含256个数据的白噪声for n=5:256xl=0;for k=2:length(a)xl=-a(k)*x(nk+l)+xl;endfor l=l:length(b)xl=xl +b(l)*u(n-l+1);endx(n)=xl;endy=x;%有ARMA模型控制的输出数据%第二个自定义函数%由ARMA模型真实参数画出真实功率谱function pxp=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*wnl'yl=abs(xl).

14、A2;k=O:length(a)-l;wn2=exp(j.*w'*k);x2二a*wn2；y2=abs(x2).A2;py=yi/y2*c；px=10*logl0(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_max4-20);%第三个自定义函数%有样本数据估计出ARMA模型的AR部分参数function cl=zxg(x)N=length(x);P=N;x=x；for k=0:p-l;xl(k+

15、l)=0;for n=l:N-kx 1 (k+l)=x 1 (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(ij)=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);cl=V(2,5),V(3,5),V(4,5),V(5,5)/V(l,5);% 估计 AR 部分参数%第四个自定义函数%产生50次的平均功率谱function B=glpp

16、(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;wnl=exp(-j.*w'*l);xl=dl*wnl'yl=abs(xl).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,:)=py;endB=0;for i=l:50B=A(i,:)+B;endB=B/50

17、;%第五个自定义函数%估计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-kxl (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(l);for i=l:50h=h+e(i)*conj(r(i+1);endP=4;q=2;for k=0:50-lrix(k+l)=0;for l=l:m-krix(k+1 )=rix

18、(k+1 )+conj(e(l)*e(l+k);endrix(k+1 )=rix(k+1)+1 *e(k+1);rix(k+1 )=rix(k+1 )/h;endfor i=l:10for m=l:2RIX(i,m)=rix(i-m+4+l);endri(i)=rix(i+5);endri=ri'U,S,V=svd(ri,RIX);dl司V(2,3),V(3,3)/V(1,3);%第六个自定义函数%对真实功率谱与平均后的功率谱进行比较function pxp=glpp3(a,b,c)k=0;w=linspace(0,0.5,50)*2*pi;l=O:length(b)-l;wn 1 =

19、exp(j.*w'*l);xl=b*wnl'yl=abs(xl).A2; k=O:length(a)-l;wn2=exp(-j.*w'*k);x2=a*wn2'y2=abs(x2).A2;py=yl./y2;N=length(w);B=glpp(a,b,c);py=10*ogl0(py); plot(w/(2*pi),pyk-);%真实功率谱曲线用黑虚线表示 hold on;B=10*logl0(B);plot(w/(2*pi),Bb-)%估计的平均功率谱曲线用蓝实线表示 h 1 =min(min(py),min(B)-10;h2=max(max(py),max(B)+10; xlabel(濒率 f (Hz)'); ylabel('功率谱 p(f)(dB)');a