%随机性问题梯度校正参数辨识
%z(k)=1.75*z(k-1)-0.8*z(k-2)+1.4*u(k-1)+0.2*u(k-2)+s(k);
%========================================
clc
clear
% 产生4位m序列作为输入数据
L=800; % L输入数据长度
y1=1;y2=0;y3=1;y4=0; %4个移位寄存器的输出初始值
for i=1:L;
x1=xor(y3,y4); %第一个移位寄存器的输入信号
x2=y1; %第二个移位寄存器的输入信号
x3=y2; %第三个移位寄存器的输入信号
x4=y3; %第四个移位寄存器的输入信号
y(i)=y4; %第四个移位寄存器的输出信号，M序列, 幅值"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
s=normrnd(0,0.3,1,800); %生成输入噪声，均值为0，方差为0.1，500个
% 产生输出数据
z(2)=0; z(1)=0; %系统最初几个输出的初值
for k=3:L; %开始求K
z(k)=1.75*z(k-1)-0.8*z(k-2)+1.4*u(k-1)+0.2*u(k-2)+s(k);
%系统在M序列输入下的输出信号
end
%赋初值
Theita{3}=[1 1 1 1]'; %参数矩阵始值（充分小）
Theita1=[-1.75 0.8 1.4 0.2];
R=eye(4);
for k=3:L; %辅助变量法，递推循环
h{k}=[-z(k-1),-z(k-2),u(k-1),u(k-2)]'; %计算h(k)矩阵
a=h{k};
b=a(1,1)^2+a(2,1)^2+a(3,1)^2+a(4,1)^2;
c=1/b;
Theita{k+1}=Theita{k}+c*inv(R)*h{k}*(z(k)-h{k}'*Theita{k}); %计算参数矩阵
R=R+c*(h{k}*h{k}'-R);
re(:,k)=Theita{k}; %记录参数变化过程
end
%结果输出
disp('辨识参数')
Theita{k}
disp('参数误差')
E=abs((Theita{k}-Theita1')./Theita1')
figure(1); %绘制参数变化轨迹
a1=re(1,:); a2=re(2,:); b1=re(3,:); b2=re(4,:); %分离出a1、 a2、 b1、 b2
i=1:L;plot(i,a1,'r',i,a2,'b',i,b1,'black',i,b2,'m:')
title('随机性问题辨识过程曲线')
