%廣義最小二乘的遞推算法仿真模型
%Z（k+2）=1.5*Z（k+1）-0.7*Z（k）+u（k+1）+0.5*u（k）+e（k）
%e（k+2）+2.1*e（k+1）-2.5*e（k）=v（k+2）
%========================================
clear?
clc
%==========400?個產生M序列作為輸入===============
x=[0?1?0?1?1?0?1?1?1];??%initial?value?
n=403;?%n為脈沖數目
M=[];	%存放M?序列
for?i=1:n
????temp=xor（x（4）x（9））;?
????M（i）=x（9）;
for?j=9:-1:2
x（j）=x（j-1）;
?end
?x（1）=temp;
end?
%===========產生均值為0，方差為1?的高斯白噪聲=========
v=randn（1400）;?
e=[];?e（1）=v（1）;?e（2）=v（2）;
for?i=3:400
e（i）=0*e（i-1）+0*e（i-2）+v（i）;
end
%==============產生觀測序列z=================
z=zeros（4001）;
z（1）=-1;
z（2）=0;
for?i=3:400
z（i）=1.5*z（i-1）-0.7*z（i-2）+M（i-1）+0.5*M（i-2）+e（i）;
end
%變換后的觀測序列
zf=[];
zf（1）=-1;
zf（2）=0;
for?i=3:400
????zf（i）=z（i）-0*z（i-1）-0*z（i-2）;
end
%變換后的輸入序列
uf=[];?uf（1）=M（1）;?uf（2）=M（2）;?
for?i=3:400
uf（i）=M（i）-0*M（i-1）-0*M（i-2）;
end
%賦初值
P=100*eye（4）;	%估計方差
Theta=zeros（4400）;??%參數的估計值，存放中間過程估值?
Theta（:2）=[3;3;3;3];
K=[10;10;10;10];??%增益?
PE=10*eye（2）;?
ThetaE=zeros（2400）;?
ThetaE（:2）=[0.5;0.3];
KE=[10;10];
%遞推Theta
for?i=3:400
h=[-zf（i-1）;-zf（i-2）;uf（i-1）;uf（i-2）];?
K=P*h*inv（h‘*P*h+1）;
Theta（:i）=Theta（:i-1）+K*（z（i）-h‘*Theta（:i-1））;?
P=（eye（4）-K*h‘）*P;
end
he=[-e（i-1）;-e（i-2）];
%遞推ThetaE
KE=PE*he*inv（1+he‘*PE*he）;
ThetaE（:i）=ThetaE（:i-1）+KE*（e（i）-he‘*ThetaE（:i-1））;?
PE=（eye（2）-KE*he‘）*PE;
%=====================輸出結果及作圖=========================
disp（‘參數a1?a2?b1?b2的估計結果：‘）
Theta（:400）
disp（‘噪聲傳遞系數c1?c2的估計結果：‘）?
ThetaE（:400）
i=1:400;
figure（1）?
plot（iTheta（1:）iTheta（2:）iTheta（3:）iTheta（4:））?
title（‘待估參數過渡過程‘）