f-检验法部分程序%els.m

主要内容： 【1】随机序列产生程序 【2】白噪声产生程序 【3】M序列产生程序 【4】二阶系统一次性完成最小二乘辨识程序 【5】实际压力系统的最小二乘辨识程序 【6】递推的最小

% 【14】F-检验法部分程序%ELS
clear
clc
L=1024; 
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)=0; %M序列的值为"0"时,辨识的输入信号取“0”
end
y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号作准备
end

A=6; %初始化；
x0=1; M=255;
for k=1:L %乘同余法递推100次；
x2=A*x0; %x2和x0分别表示xi和xi-1；
x1=mod (x2,M); %将x2存储器的数除以M，取余数放x1（xi）中；
v1=x1/256; %将x1存储器的数除以256得到小于1的随机数放v1中；
v2(:,k)=v1; % 将v1中的数（ ）存放在矩阵存储器v的第k列中，v(:,k)
%表示行不变、列随递推循环次数变化； 
x0=x1; %xi-1= xi；
v0=v1;
end %递推100次结束；
v=v2; 

z=zeros(4,1024);
z(2)=0;z(1)=0;z(3)=0;
%增广递推最小二乘辨识
c0=[0.001 0.001 0.001 0.001 0.001 0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(6,6);%直接给出初始状态P0，即一个充分大的实数单位矩阵
E=5.0e-15;%取相对误差E
c=[c0,zeros(6,1023)];%被辨识参数矩阵的初始值及大小
e=zeros(6,1024); %相对误差的初始值及大小

for k=4:1024
z(k)=1.8*z(k-1)-1.3*z(k-2)+0.4*z(k-3)+1.1*u(k-1)+0.288*u(k-2)+v(k);
end

for k=301:1024; %开始求K 
z(k)=1.8*z(k-1)-1.3*z(k-2)+0.4*z(k-3)+1.1*u(k-1)+0.288*u(k-2)+v(k);%系统在M序列输入下的输出采样信号 
h1=[-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),v(k)]';%为求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 
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%循环结束
for i=1:724
for j=1:5
a=301+i-j;
f(i,j)=z(a);

end
for j=5:10
b=301+i-j-4;
f(i,j)=u(b);
end
end
t=[-1.8 1.3 -0.4 1.1 0.288 0]';
fy_n0=f*t;
s=[1,2,3,4,5];
fy_n1=f(:,s)*[-1.8 1.3 -0.4 1.1 0.288]';
xl1=fy_n0-fy_n1;
J1=xl1'*xl1/724+1