% 【10】Bayes辨識程序

%FLch6BAeg1
clear
L=60;%四位移位積存器產(chǎn)生的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;%為下一次的輸入信號作準(zhǔn)備
end
figure(1);%畫第一個圖形
subplot(2,1,1); %畫第一個圖形的第一個子圖
stem(u),grid on%畫出M序列輸入信號
v=randn(1,60); %產(chǎn)生一組60個正態(tài)分布的隨機(jī)噪聲
subplot(2,1,2); %畫第一個圖形的第二個子圖
plot(v),grid on;%畫出隨機(jī)噪聲信號
R=corrcoef(u,v);%計算輸入信號與隨機(jī)噪聲信號的相關(guān)系數(shù)
r=R(1,2)%取出互相關(guān)系數(shù)
rv=std(v)*std(v)%計算隨機(jī)噪聲的方差
u%顯示輸入型號
v%顯示噪聲型號
z=zeros(1,60);zmd=zeros(1,60);%定義輸出采樣矩陣與模型輸出矩陣的大小
z(2)=0;z(1)=0;zmd(2)=0;zmd(1)=0;%輸出采樣、系統(tǒng)實際輸出、模型輸出賦初值
%增廣遞推最小二乘辨識
c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001]';%直接給出被辨識參數(shù)的初始值,即一個充分小的實向量
p0=10^6*eye(7,7);%直接給出初始狀態(tài)P0，即一個充分大的實數(shù)單位矩陣
E=5.e-9;%相對誤差E=0.000000005
c=[c0,zeros(7,59)];%被辨識參數(shù)矩陣的初始值及大小
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);%系統(tǒng)在M序列輸入下的輸出采樣信號 
h1=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]';%為求K(k)作準(zhǔn)備
x=h1'*p0*h1+rv;
x1=inv(x);
k1=p0*h1*x1;%K
d1=z(k)-h1'*c0;%開始求被辨識參數(shù)c
c1=c0+k1*d1;%辨識參數(shù)c
zmd(k)=h1'*c1%模型在M序列的輸入下的輸出響應(yīng)
e1=c1-c0; 
e2=e1./c0; %求參數(shù)的相對變化
e(:,k)=e2;
c0=c1;%給下一次用
c(:,k)=c1;%把辨識參數(shù)c 列向量加入辨識參數(shù)矩陣 
p1=p0-k1*k1'*[h1'*p0*h1+1];%find p(k)
p0=p1;%給下次用
if e2<=E break;%若收斂情況滿足要求，終止計算
end%判斷結(jié)束
end%循環(huán)結(jié)束
c,e, z, zmd %顯示被辨識參數(shù)、誤差情況、輸出采樣值、模型輸出值
%分離賦值
a1=c(1,; a2=c(2,; b1=c(3,; b2=c(4,; %分離出a1、a2、b1、b2
d1=c(5,; d2=c(6,; d3=c(7,; %分離出d1、 d2、 d3
ea1=e(1,; ea2=e(2,; eb1=e(3,; eb2=e(4,; %分離出a1、 a2、b1、b2的收斂情況
ed1=e(5,; ed2=e(6,; ed3=e(7,; %分離出d1、 d2、 d3的收斂情況
figure(2);%畫第二個圖形
i=1:60;
plot(i,a1,'r',i,a2,'r:',i,b1,'b',i,b2,'b:',i,d1,'g',i,d2,'g:',i,d3,'g+')%畫出各個被辨識參數(shù)
title('Parameter Identification with Recursive Least Squares Method')%標(biāo)題
figure(3);%畫出第三個圖形
i=1:60;
plot(i,ea1,'r',i,ea2,'r:',i,eb1,'b',i,eb2,'b:',i,ed1,'g',i,ed2,'g:',i,ed2,'r+')%畫出各個參數(shù)收斂情況
title('Identification Precision')%標(biāo)題
figure(4);
subplot(2,1,1);i=1:60;plot(i,z(i),'r')%畫出系統(tǒng)的采樣輸出
subplot(2,1,2);i=1:60;plot(i,zmd(i),'b')%畫出模型的輸出