%LSeg3
clear
L=60;%四位移位積存器產生的M序列的周期
y1=1;y2=1;y3=1;y4=0;%四個移位積存器的輸出初始值
for i=1:L;
    x1=xor(y3,y4);%第一個移位積存器的輸入信號
    x2=y1;%第二個移位積存器的輸入信號
    x3=y2;%第三個移位積存器的輸入信號
    x4=y3;%第四個移位積存器的輸入信號
    y(i)=y4;%第四個移位積存器的輸出信號，幅值"0"和"1"
    if y(i)>0.5,u(i)=-1;%M序列的值為"1"時,辨識的輸入信號取“-1”
    else u(i)=1;%M序列的值為"0"時,辨識的輸入信號取“1”
    end
    y1=x1;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號作準備
end
v=randn(1,60); %產生一組60個正態分布的隨機噪聲
R=corrcoef(u,v);%計算輸入信號與隨機噪聲信號的相關系數
r=R(1,2);%取出互相關系數
u%顯示輸入型號
v%顯示噪聲型號
z=zeros(7,60);zs=zeros(7,60);zm=zeros(7,60);zmd=zeros(7,60);%輸出采樣、不考慮噪聲時系統輸出、不考慮噪聲時模型輸出、模型輸出矩陣的大小
z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;zm(2)=0;zm(1)=0;zmd(2)=0;zmd(1)=0;%給輸出采樣、不考慮噪聲時系統輸出、不考慮噪聲時模型輸出、模型輸出賦初值
%增廣遞推最小二乘辨識
c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001]';%直接給出被辨識參數的初始值,即一個充分小的實向量
p0=10^6*eye(7,7);%直接給出初始狀態P0，即一個充分大的實數單位矩陣
E=5.0e-15;%取相對誤差E
c=[c0,zeros(7,59)];%被辨識參數矩陣的初始值及大小
e=zeros(7,60);%相對誤差的初始值及大小

for k=3:60; %開始求K 
    z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);%系統在M序列輸入下的輸出采樣信號    
    h1=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]';%為求K(k)作準備
    x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %K
    d1=z(k)-h1'*c0; c1=c0+k1*d1;%辨識參數c 
    zs(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%系統在M序列的輸入下不考慮擾動時的輸出響應
    zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[c1(1);c1(2);c1(3);c1(4)];%模型在M序列的輸入下不考慮擾動時的的輸出響應
    zmd(k)=h1'*c1;%模型在M序列的輸入下的的輸出響應
    e1=c1-c0; 
    e2=e1./c0; %求參數的相對變化
    e(:,k)=e2;
    c0=c1;%給下一次用
    c(:,k)=c1;%把辨識參數c 列向量加入辨識參數矩陣 
    p1=p0-k1*k1'*[h1'*p0*h1+1];%find p(k)
    p0=p1;%給下次用
    if e2<=E break;%若收斂情況滿足要求，終止計算
    end%判斷結束
end%循環結束
c, e, %顯示被辨識參數及參數收斂情況
z,zs ,zmd %顯示輸出采樣值、系統實際輸出值、模型輸出值