example9_6.m

系统辨识的M文件

clear all;
x=[ 9050, 9970,11760, 8010, 8910, 7100, 7450, 9060 , 14480,7950, 7790, 7464,  10800, 7870, 9410, 6230, 9300, 7980, 7580, 9720, 7810, 9000, 9920, 9450, 9230, 7850, 9080];
x=[x, 8270,  10560, 6950, 9680, 8060,  9990, 6920, 8600, 6820, 8460, 6940, 9700, 9590, 3870,  10880, 6070,  10720, 7200, 12080, 7310, 7940, 7410,  10000, 9440, 9170,  12880, 9640];
x=[x, 7900, 9880];     N=56;
AIC=zeros(4,5);   FPE=zeros(4,5);
xavr=sum(x)/N;     w=x-xavr;                     %零均值化
z=[w',w'];   r=covf(z, 17);      rho=r/r(1);     %求自相关函数
Q=w;    Q=[Q', [] ];                             %观测数据
for i=1:4                                         % ARMA(na,nb)模型中的na
    for j=0:i                                     % ARMA(na,nb)模型中的nb
        order=[i  j];
    Model=armax(Q,order);                         %辨识模型参数
    AIC(i,j+1)=aic(Model);                        %辨识阶次的AIC准则
    FPE(i,j+1)=fpe(Model);                        %辨识阶次的FPE准则
    end
end
AIC                                               % 输出AIC矩阵
for k=1:15
    order=[k  0];
    Model=armax(Q,order)                          %求偏自相关函数
end
et=zeros(N,1);   xt=x;
for t=2:N
    xt(t)=16187-0.8381*x(t-1)+0.6162*et(t-1);     %一步预报 
    et(t)=x(t)-xt(t);                             %预报误差
end
plot(x,'r*');   hold on;                          %作图显示
plot(x,'r');    plot(ones(N,1)*xavr,'g');
plot(xt, 'b.');    plot(xt,'b--');   hold off;