作業一遞推最小二乘法引數辨識

設被辨識系統的數學模型由下式描述：

式中x(k)為方差為0.1的白噪聲，要求：

當輸入信號u(k)是方差為1的白噪聲序列時，利用系統的輸入輸出資料在線辨識上述模型的引數；
當輸入信號u(k)是幅值為1的逆M序列時，利用系統的輸入輸出資料在線辨識上述模型的引數；
當輸入信號u(k)是單位階躍信號時，利用系統的輸入輸出資料在線辨識上述模型的引數；

分析比較在不同輸入信號作用下，對系統模型引數辨識精度的影響，

解： $\tiny y(k)=-[-1.5y(k-1)+0.7y(k-2)+0.1y(k-3)]+u(k-2)+2u(k-3)+1.5u(k-4)+\xi (k)$

代碼待完善，還沒修改完

【參考代碼】

①方差為1的白噪聲序列

%遞推最小二乘引數估計（RLS）
clear all; close all;
 
a=[1 -1.5 0.7 0.1]'; 
b=[1 2 1.5]'; d=2; %物件引數
na=length(a)-1; nb=length(b)-1; %na=3、nb=2為A、B階次
 
L=400; %仿真長度
uk=zeros(d+nb,1); %輸入初值：uk(i)表示u(k-i)
yk=zeros(na,1); %輸出初值
u=randn(L,1); %輸入采用白噪聲序列
xi=sqrt(1)*randn(L,1); %白噪聲序列  方差=1
 
theta=[a(2:na+1);b]; %物件引數真值
 
thetae_1=zeros(na+nb+1,1); %thetae初值
P=10^6*eye(na+nb+1); 
for k=1:L
    phi=[-yk;uk(d:d+nb)]; %此處phi為列向量
    y(k)=phi'*theta+xi(k); %采集輸出資料
   
    %遞推最小二乘法
    K=P*phi/(1+phi'*P*phi);
    thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);
    P=(eye(na+nb+1)-K*phi')*P;
    
    %更新資料
    thetae_1=thetae(:,k);
    
    for i=d+nb:-1:2
        uk(i)=uk(i-1);
    end
    uk(1)=u(k);
    
    for i=na:-1:2
        yk(i)=yk(i-1);
    end
    yk(1)=y(k);
end
plot([1:L],thetae); %line([1,L],[theta,theta]);
xlabel('k'); ylabel('引數估計a、b');
legend('a_1','a_2','a_3','b_0','b_1','b_2'); axis([0 L -2 1.5]);

②幅值為1的逆M序列

%遞推最小二乘引數估計（RLS）
clear all; close all;
 
a=[1 -1.5 0.7 0.1]'; 
b=[1 2 1.5]'; d=2; %物件引數
na=length(a)-1; nb=length(b)-1; %na=3、nb=2為A、B階次
 
L=400; %仿真長度
uk=zeros(d+nb,1); %輸入初值：uk(i)表示u(k-i)
yk=zeros(na,1); %輸出初值
%u=randn(L,1); %輸入采用白噪聲序列
xi=sqrt(1)*randn(L,1); %白噪聲序列  方差=1

theta=[a(2:na+1);b]; %物件引數真值
 
thetae_1=zeros(na+nb+1,1); %thetae初值
P=10^6*eye(na+nb+1); 

%M序列及逆M序列的產生
%M序列長度 L=400; %仿真長度
x1=1; x2=1; x3=1; x4=0; %移位暫存器初值xi-1、xi-2、xi-3、xi-4
S=1; %方波初值

for k=1:L
    M(k)=xor(x3,x4); %進行異或運算,產生M序列
    IM=xor(M(k),S); %進行異或運算,產生逆M序列
    if IM==0
        u(k)=-1;
    else
        u(k)=1;
    end
    
    S=not(S); %產生方波
    
    x4=x3; x3=x2; x2=x1; x1=M(k); %暫存器移位
	
	
	phi=[-yk;uk(d:d+nb)]; %此處phi為列向量
    y(k)=phi'*theta+xi(k); %采集輸出資料
   
    %遞推最小二乘法
    K=P*phi/(1+phi'*P*phi);
    thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);
    P=(eye(na+nb+1)-K*phi')*P;
    
    %更新資料
    thetae_1=thetae(:,k);
    
    for i=d+nb:-1:2
        uk(i)=uk(i-1);
    end
    uk(1)=u(k);
    
    for i=na:-1:2
        yk(i)=yk(i-1);
    end
    yk(1)=y(k);
	
end

plot([1:L],thetae); %line([1,L],[theta,theta]);
xlabel('k'); ylabel('引數估計a、b');
legend('a_1','a_2','a_3','b_0','b_1','b_2'); axis([0 L -2 4]);

③單位階躍信號

%遞推最小二乘引數估計（RLS）
clear all; close all;
 
a=[1 -1.5 0.7 0.1]'; 
b=[1 2 1.5]'; d=2; %物件引數
na=length(a)-1; nb=length(b)-1; %na=3、nb=2為A、B階次
 
L=400; %仿真長度
uk=zeros(d+nb,1); %輸入初值：uk(i)表示u(k-i)
yk=zeros(na,1); %輸出初值

u=ones(L,1);%輸入采用單位階躍信號

xi=sqrt(1)*randn(L,1); %白噪聲序列  方差=1
 
theta=[a(2:na+1);b]; %物件引數真值
 
thetae_1=zeros(na+nb+1,1); %thetae初值
P=10^6*eye(na+nb+1); 
for k=1:L
    phi=[-yk;uk(d:d+nb)]; %此處phi為列向量
    y(k)=phi'*theta+xi(k); %采集輸出資料
   
    %遞推最小二乘法
    K=P*phi/(1+phi'*P*phi);
    thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);
    P=(eye(na+nb+1)-K*phi')*P;
    
    %更新資料
    thetae_1=thetae(:,k);
    
    for i=d+nb:-1:2
        uk(i)=uk(i-1);
    end
    uk(1)=u(k);
    
    for i=na:-1:2
        yk(i)=yk(i-1);
    end
    yk(1)=y(k);
end
plot([1:L],thetae); %line([1,L],[theta,theta]);
xlabel('k'); ylabel('引數估計a、b');
legend('a_1','a_2','a_3','b_0','b_1','b_2'); axis([0 L -2 1.5]);