递推最小二乘法RLS及模型阶次辨识F-Test

系统辨识课上的实验作业，该压缩文件里包含我的实验报告以及matlab源代码，包括由M序列产生白噪声，递推最小二乘法以及使用F-Test法进行模型阶次辨识，比较适合学生。

close all;
clear all;

%%%%%%%%%%%人机对话%%%%%%%%%%%%%%%%%
lambda=input（‘噪声标准差为（0.0或0.1或0.5或1.0）? ‘）;
L=input（‘数据长度为（100或300或500）? ‘）;
miu=input（‘遗忘因子为（01之间）? ‘）;
%%%%%%%%设定模型阶次从2到6，即Nbeg=2Nend=6（假定na=nb）%%%%%%%%%

%%%%%%%%%%%%%%%赋初值并生成M序列和噪声序列%%%%%%%%%%
%%%%%%%%%%%%%%生成a=1P=4的M序列%%%%%%%%%%%%%%%%%
M=[0 1 0 1 0 ];
Np=2^4-1;
for k=1:534
     if M（5）==0 
        u（k）=1;end;
     if M（5）==1  
        u（k）=-1;end;
     for i=5:-1:2
        M（i）=M（i-1）;
     end
     M（1）=M（2）+M（5）;
     if M（1）==2 
        M（1）=0;end;
end
%%%%%%%%%%%%%%%生成白噪声v（k）%%%%%%%%%%%%%%%%%%%%
A=179;  xi=11;  M=32768;
for k=1:534       
    ksai=0;
    for i=1:12
        xi=A*xi;
        xi=mod（xiM）;
        ksai=ksai+（xi/M）;
    end
v（k）=ksai-6;
end

%%%%%%%%%%%%%%%%%过程仿真%%%%%%%%%%%%%%%%
z（1）=0;
z（2）=0;
for k=3:length（u）
    z（k）=1.5*z（k-1）-0.7*z（k-2）+u（k-1）+0.5*u（k-2）+lambda*v（k）;%
end

%%%%%%%%%%%%%%%%%%模型参数估计，RLS算法%%%%%%%%%%%%
for n=2:6
    for i=1:2*n
        theta（in）=eps;
    end
    P=eye（2*n）;
    for i=1:L
        for j=1:n
            if （i-j）<=0
                H（ji）=0;
            else
            H（ji）=-z（i-j）;
            end
        end 
        for m=（n+1）:2*n
            if （i-m+n）<=0
                H（mi）=0;
            else
            H（mi）=u（i-m+n）;
            end
        end
    end
    J（n）=0;
    
    for k=1:L
        K=P*H（:k）*inv（H（:k）‘*P*H（:k）+miu）;
        J（n）=miu*（J（n）+（（z（k）-H（:k）‘*theta（:n））^2/（H（:k）‘*P*H（:k）+miu）））;
        theta（:n）=theta（:n）+K*（z（k）-H（:k）‘*theta（:n））;
        P=（（eye（2*n）-K*H（:k）‘）*P）/miu;
    end
end

%%%%%%%%%%%%%%%模型阶次识别%%%%%%%%%%%%%%%%  
if L==100
    for n=1:5
        if （J（n）-J（n+1））*（L-2*n-2）/（2*J（n+1））<=3.09
           N=n+1
           break
        end
    end
end
if L==300
     for n=1:5
        if （J（n）-J（n+1））*（L-2*n-2）/（2*J（n+1））<=3.03
           N=n+1
           break
        end
     end
end
if L==500
     for n=1:5
        if （J（n）-J（n+1））*（L-2*n-2）/（2*J（n+1））<=3.01
           N=n+1
           break
        end
     end
end

%%%%%%%%%%%%%%%输出结果%%%%%%%%%%%%%%%%
if N==2
    A=theta（1:2N）
    B=theta（3:4N）
end
if N==3
    A=theta（1:3N）
    B=theta（4:6N）
end   
if N==4
    A=theta（1:4N）
    B=theta（5:8N）
end
if N==5
    A=theta（1:5N）
    B=theta（6:10N）
end
if N==6
    A=theta（1:6N）
    B=theta（7:12N）
end

%%%%%%%%%%%%%%%%%%%计算性能指标%%%%%%%%%%%%%%%%%%
lam1=sqrt（J（N）/（L-2*N））%%%%%%%%%%%%%%噪声标准差的估计
lam2=sum（B）/（1+sum（A））%%%%%%%%%%%%%模型静态增益估计
a0=1;a1=-1.5;a2=0.7;b0=0;b1=1;b2=0.5;
thetazhen=[a1a2b1b2]‘;
chazhi=thetazhen-theta（1:42）;
lam3=sqrt（sum（（chazhi./thetazhen）.^2））%%%%%%参数估计平方相对偏差
lam4=sqrt（sum（chazhi.^2）/sum（thetazhen.^2））%%%%%%%%%参数估计平方根偏差
Kzhen=（1+0.5）/（1-1.5+0.7）;
lam5=sqrt（abs（lam2-Kzhen）/Kzhen）%%%%%%%%%静态增益估计相对偏差
a02=a0;a12=a1;a22=a2;
a01=（a02*a02-a22*a22）/a02;
a11=（a02*a12-a22*a12）/a02;
a00=（a01*a01-a11*a11）/a01;
b02=b0;b12=b1;b22=b2;
b11=

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件     418613  2011-10-18 15:57  实验二报告.docx

     文件       3310  2011-04-26 21:04  main2.m

----------- ---------  ---------- -----  ----

               421923                    2