最小二乘法系统辨识

clear all;
close all;

r=35;
Np=15;
N=（r+1）*Np; %数据长度
%% ------生成M序列和输入数据
P=4;%寄存器长度
M=[0 1 0 0 1];
u=zeros（1N）;
a=1; %M序列幅度
s=1; %方波
for k=1:N
    M（1）=M（2）+M（5）;
    if M（1）==2
        M（1）=0;
    end 
    IM=M（1）+s; %逆M序列
    if IM==2 || IM==0
        u（k）=a;
    else
         u（k）=-a;       
    end    
    for i=P+1:-1:2 %移位
        M（i）=M（i-1）;
    end 
    s=~s;
end
%% ------生成白噪声
A=179;
m=32768;
sigma=1; %白噪声标准差
v=zeros（1N）; %白噪声
xi=11;
for k=1:N
    ksai=0;
    for i=1:12   
          xi=A*xi;
          xi=mod（xim）;
          ksai=ksai+xi/m;
    end
    v（k）=sigma*（ksai-6.0）;
end
%% ------生成噪声
LAMBDA=0.1; %噪声标准差
e=zeros（1N）;
for k=3:N
    e（k）=1.5*e（k-1）-0.7*e（k-2）+LAMBDA*v（k）;
end
%% ------生成输出数据
y=zeros（1N）; %不含噪声
for k=3:N
    y（k）=1.5*y（k-1）-0.7*y（k-2）+u（k-1）+0.5*u（k-2）;
end
z=y+e; %含噪声
%% ------遗忘因子法
Mu=0.98; %遗忘因子
L=300; %数据长度
EPSILON=1e-3;
a2=1e6;
n=1; % 初始阶次
nend=10; %最大阶次
J=zeros（1nend）; %损失函数
if L==100
    talpha=3.09; %阈值
else
    if L==300
        talpha=3.03;
    else
        if L==500
            talpha=3.01;
        end
    end
end
flag=1;
t=talpha+1;
while （flag && n    THETA=EPSILON*ones（n+n1）; %参数向量初始值
    P0=a2*eye（n+n）; %协方差矩阵初始值
    h=zeros（n+n1）; %数据向量   
    THETA0=THETA（:1）;
    for k=n+1:L
       h=[-z（k-1:-1:k-n）u（k-1:-1:k-n）]‘; %输入输出数据
       K=P0*h/（h‘*P0*h+Mu）;