ARMA(2，2)模型的RLS-RELS参数估计

最新推荐文章于 2022-07-17 20:58:17 发布

ZC_lianshuang

最新推荐文章于 2022-07-17 20:58:17 发布

阅读量3k

点赞数 2

分类专栏：系统辨识的MATLAB实现文章标签： RLS-RELS ARMA 改进RELS

本文链接：https://blog.csdn.net/ZC_lianshuang/article/details/90719330

版权

系统辨识的MATLAB实现专栏收录该内容

4 篇文章 4 订阅

订阅专栏

%*****ARMA(2，2)模型的RLS-RELS参数估计*******
clear;
clc;
close all;
N=2000;
n=2;
lamda=1;
rand('seed',10);
e=0.6*randn(1,N+1);%a(t)
a=-[-1.7000    0.7200];%poly([0.8 0.9])
d=[-0.3 0.02];
y(1)=e(1);
y(2)=a(1)*y(1)+e(2)+d(1)*e(1);

for t=3:N+1
    y(t)=a(1)*y(t-1)+a(2)*y(t-2)+e(t)+d(1)*e(t-1)+d(2)*e(t-2);
end

%-------RLS算法-----------------
n0=25;
fai2(:,1)=zeros(1,n0)';
for i=2:n0
    for j=1:i-1
        fai2(j,i)=-y(i-j);      
        %fai2(j+n0,i)=e(i-j);     
    end
end
sita2(:,1)=zeros(n0,1);
p2(:,:,1)=10^5*eye(n0);
for i=2:N
    if i>n0
       for j=1:n0
        fai2(j,i)=-y(i-j);
        %fai2(j+n0,i)=e(i-j);
       end
    end 
   sita2(:,i)=sita2(:,i-1)+p2(:,:,i-1)*fai2(:,i)/(1+fai2(:,i)'*p2(:,:,i-1)*fai2(:,i))...
        *(y(i)-fai2(:,i)'*sita2(:,i-1));
    p2(:,:,i)=p2(:,:,i-1)-p2(:,:,i-1)*fai2(:,i)*fai2(:,i)'*p2(:,:,i-1)/...
         (1+fai2(:,i)'*p2(:,:,i-1)*fai2(:,i));
     ejian(i)=y(i)-fai2(:,i)'*sita2(:,i-1);
end
%---------RELS算法-----------------
sita(:,2)=zeros(2*n,1);
P(:,:,2)=10^5*eye(2*n);
ee(1)=y(1);
fai(:,2)=[y(1) 0 ee(1) 0];
for t=2:N
    ee(t)=y(t)-fai(:,t)'*sita(:,t-1);
    fai(:,t+1)=[y(t) y(t-1) ee(t) ee(t-1)]';%add '
    sita(:,t+1)=sita(:,t)+P(:,:,t)*fai(:,t+1)/(lamda+fai(:,t+1)'*P(:,:,t)*fai(:,t+1))*...
        [y(t+1)-fai(:,t+1)'*sita(:,t)];
    P(:,:,t+1)=1/lamda*(P(:,:,t)-P(:,:,t)*fai(:,t+1)/(lamda+fai(:,t+1)'*P(:,:,t)*fai(:,t+1))*...
        fai(:,t+1)'*P(:,:,t));
end   

ee(1)=y(1);
taoe(1)=ee(1)^2;%y(0)的初值为0
for t=2:N
    %对噪声方差的估计
    ee(t)=y(t)-fai(:,t)'*sita(:,t-1);
    taoe(t)=taoe(t-1)+1/t*[ee(t)^2-taoe(t-1)];
end
%----------作图-----------
t=1:N;
subplot(3,2,1);
plot(t,sita(1,t),'r');
line([0,N],[a(1),a(1)])
axis([0 N 0 3]);
title('a1');
subplot(3,2,2);
plot(t,sita(2,t),'r');
line([0,N],[a(2),a(2)])
axis([0 N -2 0]);
title('a2');
subplot(3,2,3);
plot(t,sita(3,t),'r');
line([0,N],[d(1),d(1)])
axis([0 N -0.5 0.5]);
title('d1');
subplot(3,2,4);
plot(t,sita(4,t),'r');
line([0,N],[d(2),d(2)])
axis([0 N -0.5 0.5]);
title('d2');
subplot(3,2,5);
plot(t,taoe,'r');
line([0,N],[0.36,0.36])
axis([0 N 0 1]);
title('taoe');