%-----CAR(2)--------
clear;
clc;
close all;
N=4000; %步长
rand('seed',10); %改变种子改变数据
e=0.6*randn(1,N+1); %白噪声序列
a=-[-0.5000 0.0400]; %poly([0.1 0.4])
b=[3 1];
rand('seed',4); %改变种子改变数据
u=2*rand(1,N+1)-1; %u(t)在[-1,1]之间服从均匀分布
y(1)=e(1); %初值y(1),y(2),y(3)
y(2)=a(1)*y(1)+b(1)*u(1)+e(2);
for t=3:N+1
y(t)=a(1)*y(t-1)+a(2)*y(t-2)+b(1)*u(t-1)+b(2)*u(t-2)+e(t);%AR(3)的RLS算法实现
end
n=2;
lamda=0.998;%遗忘因子
sita(:,2)=zeros(2*n,1);%初值,均为0
P(:,:,2)=10^5*eye(2*n);%P的初值,选取足够大的单位阵
for t=2:N
fai(:,t+1)=[y(t) y(t-1) u(t) u(t-1)]';%最小二乘格式的fai(t)的表达式
%参数估值
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 1]);
title('a(1)');
subplot(3,2,2);
plot(t,sita(2,t),'r');
line([0,N],[a(2),a(2)])
axis([0 N -1 1]);
title('a(2)');
subplot(3,2,3);
plot(t,sita(3,t),'r');
line([0,N],[b(1),b(1)])
axis([0 N 0 5]);
title('b(1)');
subplot(3,2,4);
plot(t,sita(4,t),'r');
line([0,N],[b(2),b(2)])
axis([0 N 0 3]);
title('a(2)');
subplot(3,2,5);
plot(t,taoe,'r');
axis([0 N 0 1]);
line([0,N],[0.36,0.36])
title('taoe');
subplot(3,2,6);
plot(t,u(t),'r');
title('u(t)');
其中,e代表,为白噪声序列
lamda代表,为遗忘因子
sita代表,估值参数之一
fai代表,为估值参数之一
ee代表
taoe代表
其中一次的测试结果如下图所示: