www.gusucode.com > MA模型算法仿真程序源码 > MA模型算法仿真程序源码/x/MA.m
clc; clear; close all; N=456; B1=[1 0.3544 0.3508 0.1736 0.2401]; A1=[1 -1.3817 1.5632 -0.8843 0.4096]; w=linspace(0,pi,512); H1=freqz(B1,A1,w);%产生信号的频域响应 Ps1=abs(H1).^2; SPy11=0;%20次AR(4) SPy12=0;%20次AR(8) SPy13=0;%20次AR周期图 SPy14=0;%20次ARMA(4,4) SPy15=0;%20次ARMA(8,8) VSPy11=0;%20次AR(4) VSPy12=0;%20次AR(8) VSPy13=0;%20次AR周期图 VSPy14=0;%20次ARMA(4,4) VSPy15=0;%20次ARMA(8,8) % for k=1:20 %采用自协方差法对AR模型参数进行估计% %gA1:AR模型的参数;gE1:激励白噪声的方差% y1=filter(B1,A1,randn(1,N)).*[zeros(1,200),ones(1,256)]; [Py11,F]=pcov(y1,4,512,1);%AR(4)的估计% [Py12,F]=pcov(y1,8,512,1);%AR(8)的估计% [Py13,F]=periodogram(y1,[],512,1); SPy11=SPy11+Py11; SPy12=SPy12+Py12; SPy13=SPy13+Py13; VSPy11=VSPy11+abs(Py11).^2; VSPy12=VSPy12+abs(Py12).^2; VSPy13=VSPy13+abs(Py13).^2; figure(1) plot(w./(2*pi),Ps1,F,Py11); legend('真实功率谱','20次AR(4)估计图'); hold on; figure(2) plot(w./(2*pi),Ps1,F,Py12); legend('真实功率谱','20次AR(8)估计图'); hold on; figure(3) plot(w./(2*pi),Ps1,F,Py13); legend('真实功率谱','20次周期图法估计图'); hold on; %------------ARMA模型---------------% y=zeros(1,256); for i=1:256 y(i)=y1(200+i); end ny=[0:255]; z=fliplr(y);nz=-fliplr(ny); nb=ny(1)+nz(1);ne=ny(length(y))+nz(length(z)); n=[nb:ne]; Ry=conv(y,z); R4=zeros(8,4);%ARMA(4,4)的R r4=zeros(8,1);%ARMA(4,4)的r for i=1:8 r4(i,1)=-Ry(260+i); for j=1:4 R4(i,j)=Ry(260+i-j); end end R4%R矩阵 r4%r矩阵 a4=inv(R4'*R4)*R4'*r4%利用最小二乘法得到的a的估计参数 %----------------ARMA(4,4)对MA的参数b(1)-b(4)进行估计----------------------% A1 A14=[1,a4']%AR的参数a(1)-a(4)的估计值 B14=fliplr(conv(fliplr(B1),fliplr(A14)));%MA模型的分子 y24=filter(B14,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)];%由估计出的MA模型产生数据 %---因为(q=4)<<L<<N=256,所以选取L=32---% [Ama4,Ema4]=arburg(y24,32),%利用数据y2估计AR(32)的参数 B1 b4=arburg(Ama4,4)%求出MA模型的参数 %---求功率谱---% w=linspace(0,pi,512); %H1=freqz(B1,A1,w) H14=freqz(b4,A14,w);%产生信号的频域响应 %Ps1=abs(H1).^2;%真实谱 Py14=abs(H14).^2;%估计谱 %if Py14>200 % PPy14=200; %elseif Py14<200 % PPy14=Py14; %end SPy14=SPy14+Py14; VSPy14=VSPy14+abs(Py14).^2; figure(4) plot(w./(2*pi),Ps1,w./(2*pi),Py14); legend('真实功率谱','20次ARMA(4,4)的估计图'); hold on; R8=zeros(16,8);%ARMA(8,8)的R r8=zeros(16,1);%ARMA(8,8)的r for i=1:16 r8(i,1)=-Ry(264+i); for j=1:8 R8(i,j)=Ry(264+i-j); end end R8%R矩阵 r8%r矩阵 a8=inv(R8'*R8)*R8'*r8%利用最小二乘法得到的a的估计参数 %----------------ARMA(8,8)对MA的参数b(1)-b(8)进行估计----------------------% A1 A18=[1,a8']%AR的参数a(1)-a(8)的估计值 B18=fliplr(conv(fliplr(B1),fliplr(A18)))%MA模型的分子 y28=filter(B18,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)];%由估计出的MA模型产生数据 %---因为(q=8)<<L<<N=256,所以选取L=48---% [Ama8,Ema8]=arburg(y28,48),%利用数据y2估计AR(32)的参数 B1 b8=arburg(Ama8,8)%求出MA模型的参数 %---求功率谱---% %w=linspace(0,pi,512); %H1=freqz(B1,A1,w) H18=freqz(b8,A14,w);%产生信号的频域响应 %Ps1=abs(H1).^2;%真实谱 Py15=abs(H18).^2;%估计谱 %if Py15>200 % PPy15=200; %elseif Py15<200 % PPy15=Py15; %end SPy15=SPy15+Py15; VSPy15=VSPy15+abs(Py15).^2; figure(5) plot(w./(2*pi),Ps1,w./(2*pi),Py15); legend('真实功率谱','20次ARMA(8,8)的估计图'); hold on; %-----------------------------------% % end V1=VSPy11/20-abs(SPy11/20).^2; V2=VSPy12/20-abs(SPy12/20).^2; V3=VSPy13/20-abs(SPy13/20).^2; V4=VSPy14/20-abs(SPy14/20).^2; V5=VSPy15/20-abs(SPy15/20).^2; figure(6) plot(w./(2*pi),Ps1,F,SPy11/20); legend('真实功率谱','20次AR(4)估计的均值'); figure(7) plot(w./(2*pi),Ps1,F,SPy12/20); legend('真实功率谱','20次AR(8)估计的均值'); figure(8) plot(w./(2*pi),Ps1,F,SPy13/20); legend('真实功率谱','20次周期图估计的均值'); figure(9) plot(w./(2*pi),Ps1,w./(2*pi),SPy14/20); legend('真实功率谱','20次ARMA(4,4)估计的平均值'); figure(10) plot(w./(2*pi),Ps1,w./(2*pi),SPy15/20); legend('真实功率谱','20次ARMA(8,8)估计的平均值'); figure(12) plot(F,V1,F,V2,F,V3,w./(2*pi),V4,w./(2*pi),V5); legend('AR(4)方差','AR(8)方差','周期图方差','ARMA(4,4)方差','ARMA(8,8)方差'); axis([0 0.5 0 2000]);