(完整word版)多种最小二乘算法分析+算法特点总结

第一部分：程序设计思路、辨识结果分析和算法特点总结 (2)一：RLS遗忘因子法 (2)RLS遗忘因子法仿真思路和辨识结果 (2)遗忘因子法的特点： (3)二：RFF遗忘因子递推算法 (4)仿真思路和辨识结果 (4)遗忘因子递推算法的特点： (5)三：RFM限定记忆法 (5)仿真思路和辨识结果 (5)RFM限定记忆法的特点： (7)四：RCLS偏差补偿最小二乘法 (7)仿真思路和辨识结果 (7)RCLS偏差补偿最小二乘递推算法的特点： (9)五：增广最小二乘法 (9)仿真思路和辨识结果 (9)RELS增广最小二乘递推算法的特点： (11)六：RGLS广义最小二乘法 (12)仿真思路和辨识结果 (12)RGLS广义最小二乘法的特点： (14)七：RIV辅助变量法 (14)仿真思路和辨识结果 (14)RIV辅助变量法的特点： (16)八：Cor-ls相关最小二乘法（二步法） (17)仿真思路和辨识结果 (17)Cor-ls相关最小二乘法（二步法）特点： (18)九：MLS多级最小二乘法 (19)仿真思路和辨识结果 (19)MLS多级最小二乘法的特点： (22)十：yule_walker辨识算法 (23)仿真思路和辨识结果 (23)yule_walker辨识算法的特点： (24)第二部分：matlab程序 (24)一：RLS遗忘因子算法程序 (24)二：RFF遗忘因子递推算法 (26)三：RFM限定记忆法 (28)四：RCLS偏差补偿最小二乘递推算法 (31)五：RELS增广最小二乘的递推算法 (33)六;RGLS 广义最小二乘的递推算法 (36)七：Tally辅助变量最小二乘的递推算法 (39)八：Cor-ls相关最小二乘法（二步法） (42)九：MLS多级最小二乘法 (45)十yule_walker辨识算法 (49)第一部分：程序设计思路、辨识结果分析和算法特点总结一：RLS遗忘因子法RLS遗忘因子法仿真思路和辨识结果仿真对象如下：其中， v(k )为服从N(0,1)分布的白噪声。

输入信号u(k)采用M 序列，幅度为 1。

M 序列由 9 级移位寄存器产生，x(i)=x(i-4)⊕x(i-9)。

选择如下辨识模型：加权阵取Λ=I。

衰减因子β = 0.98，数据长度 L = 402。

辨识结果与理论值比较，基本相同。

辨识结果可信：Estimate =-1.46660.65030.97360.3035遗忘因子法的特点：对老数据加上遗忘因子，以降低老数据对辨识的影响，相对增加新数据对辨识的影响，不会出现“数据饱和”现象。

如模型噪声是有色噪声，则Ø是有偏估计量。

常用作其他辨识方式的起步，以获得其他方式的初始值。

二：RFF遗忘因子递推算法仿真思路和辨识结果辨识模型与遗忘因子法所用模型相同。

其中， 0 ≤µ≤1为遗忘因子，此处取0.98。

数据长度L=402，初始条件：参数a1 a2 b1 b2的估计值:ans =-1.49770.68631.19030.4769待估参数变化过程如图所示：遗忘因子递推算法的特点：从上面两个例子可以看出对于相同的仿真对象，一次算法和递推算法结果基本一致，但递推算法可以实现在线实时辨识，而且可以减少计算量和存储量。

三：RFM限定记忆法仿真思路和辨识结果辨识模型与遗忘因子法所用模型相同。

辨识结果与理论值比较，基本相同。

辨识结果可信：参数 a1 a2 b1 b2 的估计值为：Theta_a =-1.51280.70990.83930.4416待估参数的过渡过程如下：RFM限定记忆法的特点：辨识所使用的数据长度保持不变，每增加一个新数据就抛掉一个老数据，使参数估计值始终只依赖于有限个新数据所提供的新消息，克服了遗忘因子法不管多老的数据都在起作用的缺点，因此该算法更能有效的克服数据饱和现象。

四：RCLS偏差补偿最小二乘法仿真思路和辨识结果辨识模型与遗忘因子法所用模型相同。

辨识结果与理论值比较，基本相同。

辨识结果可信：参数a1 a2 b1 b2的估计值为：ans =-1.49160.70051.03650.4271RCLS偏差补偿最小二乘递推算法的特点：算法思想：：在最小二乘参数估计值的基础上，引进补偿项σW2C-1D Ø0，则获得了参数的无偏估计。

针对模型噪声来说，RCLS算法的适应能力比RLS更好。

五：增广最小二乘法仿真思路和辨识结果考虑如下仿真对象：其中，为服从N(0,1)分布的白噪声。

输入信号采用 M 序列，幅度为 1。

M 序列由 9 级移位寄存器产生，x(i)=x(i-4)⊕x(i-9)。

选择如下的辨识模型：观测数据长度取L =402 。

加权阵取Λ=I。

辨识结果与理论值比较，基本相同，同时又能获得噪声模型的参数估计。

辨识结果可信：参数a1、a2、b1、b2、d1、d2估计结果：ans =-1.50000.70001.00010.5002-0.99990.2000RELS增广最小二乘递推算法的特点：增广最小二乘的递推算法对应的噪声模型为滑动平均噪声，扩充了参数向量和数据向量H（k）的维数，把噪声模型的辨识同时考虑进去。

最小二乘法只能获得过程模型的参数估计，而增广最小二乘法同时又能获得噪声模型的参数估计，若噪声模型为平均滑动模型，，则只能用RELS算法才能获得无偏估计。

当数据长度较大时，辨识精度低于极大似然法。

六：RGLS广义最小二乘法仿真思路和辨识结果模型结构选择：模型结构选用：其中，各个参数的真值为：广义最小二乘算法为：辨识结果与理论值比较，基本相同，同时又能获得噪声传递系数的参数估计。

辨识结果可信：参数a1 a2 b1 b2的估计结果：ans =-1.50580.69720.93160.4833噪声传递系数c1 c2的估计结果：ans =0.62030.2210RGLS广义最小二乘法的特点：该算法用于自回归输入模型，是一种迭代的算法。

其基本思想是基于对数据先进行一次滤波处理，后利用普通最小二乘法对滤波后的数据进行辨识，进而获得无偏一致估计。

但是当过程的输出信噪比比较大或模型参数较多时，这种数据白色化处理的可靠性就会下降，辨识结果往往会是有偏估计。

数据要充分多，否则辨识精度下降。

模型阶次不宜过高。

初始值对辨识结果有较大影响。

七：RIV辅助变量法仿真思路和辨识结果辨识模型与遗忘因子法所用模型相同,只不过此处噪声为有色噪声，产生过程为：e(k)=v(k)+0.5v(k-1)+0.2v(k-2),v(k)为0均值的不相关随机噪声。

按照Tally法选取辅助变量x(k)=z(k-n d), n d为误差传递函数的阶数，此处为2.则有辅助变量法的递推公式可写成：辨识结果与理论值比较，基本相同。

辨识结果可信：参数a1 a2 b1 b2的估计结果：ans =-1.53140.74610.99990.4597RIV辅助变量法的特点：适当选择辅助变量，使之满足相应条件，参数估计值就可以是无偏一致。

估计辅助变量法的计算量与最小二乘法相当，但辨识效果却比最小二乘法好的多。

尤其当噪声是有色的，而噪声的模型结构又不好确定时，增广最小二乘法和广义最小二乘法一般都不好直接应用，因为他们需要选用特定的模型结构，而辅助变量法不需要确定噪声的模型结构，因此辅助变量法就显得更为灵活，但辅助变量法不能同时获得噪声模型的参数估计。

八：Cor-ls相关最小二乘法（二步法）仿真思路和辨识结果辨识模型与遗忘因子法所用模型相同：，e(k)=v(k)+0.5v(k-1)+0.2v(k-2),v(k)为0均值的不相关随机噪声。

Cor-ls的递推公式可写成:其中：，M（k）为输入M序列。

初始条件：辨识结果与理论值比较，基本相同，辨识结果可信：参数a1 a2 b1 b2的估计结果：ans =-1.48960.68581.01680.4362Cor-ls相关最小二乘法（二步法）特点：把辨识分成两步进行：第一步：利用相关分析法获得对象的非参数模型（脉冲响应或相关函数）；第二步：利用最小二乘法、辅助变量法或增广最小二乘法等，进一步求的对象的参数模型。

如果模型噪声与输入无关，则Cor-ls相关最小二乘法（二步法）可以得到较好的辨识结果。

Cor-ls相关最小二乘法（二步法）实质上是先对数据进行一次相关分析，滤除了有色噪声的影响，再利用最小二乘法必然就会改善辨识结果。

能适应较宽广的噪声范围，计算量不大，初始值对辨识结果影响较小。

但要求输入信号与噪声不相关。

九：MLS多级最小二乘法仿真思路和辨识结果仿真对象如下:其中，u (k)是输入变量，此处为 M 序列；v (k ) 是零均值、方差为 1 的不相关随机噪声，通过控制λ的大小来控制信噪比。

辨识模型结构选用:其中，辨识过程如下：第一级，辅助模型参数辨识原模型可写为:利用最小二乘法可获得辅助模型的参数无偏一致估计值:数据长度 L=400，第二级，过程模型参数辨识:根据最小二乘算法可以获得过程模型的参数估计值为:第三级，噪声模型参数辨识:根据最小二乘算法可以获得过程模型的参数估计值为辨识结果与理论值比较，基本相同。

辨识结果可信：第一级辅助模型参数 e1 e2 e3 e3 e4 f1 f2 f3 f4 辨识结果：E =1.90621.44540.52790.0613-0.00260.7988-0.8694-1.3037-0.6318第二级过程模型参数 a1 a2 a3 b1 b2 辨识结果：E2 =0.93040.15960.01130.7998-1.6502第三级噪声模型参数 c1 c2 辨识结果：E3 =0.97500.3824MLS多级最小二乘法的特点：当信噪比较大时，采用广义最小二乘法可能会出现多个局部收敛点，解决这个问题的方法可用多级最小二乘法，一般来说多级最小二乘法包含三级辨识过程。

利用输入输出数据，通过多级最小二乘法，可分别求的辅助模型，过程模型和噪声模型的参数估计值。

在高噪声的情况下，多级最小二乘法明显优于广义最小二乘法，其收敛点唯一。

十：yule_walker辨识算法仿真思路和辨识结果仿真对象如下：，z (k)是可观测变量；v (k )是均值为零，方差为 1 的不相关随机噪声；数据长度取 L=1024。

相关函数按下式计算：参数的估计算法按下式计算：辨识结果与理论值比较，基本相同，同时又能获得噪声模型的参数估计。

辨识结果可信：辨识结果为：Theta =0.85970.2955-0.0034d =1.0025yule_walker辨识算法的特点：yule_walker辨识算法可以方便的辨识形如的参数估计值。

第二部分：matlab程序一：RLS遗忘因子算法程序clearclc%==========================================%最小二乘法辨识对象% Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k)%==========产生M序列作为输入===============x=[0 1 0 1 1 0 1 1 1]; %初始值n=403; %n为脉冲数目M=[]; %存放M序列for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end;%产生高斯白噪声v=randn(1,400);z=[];z(1)=-1;z(2)=0;u=0.98;% 遗忘因子L=400;for i=3:402z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); zstar(i)=z(i)*u^(L-i+2);endH=zeros(400,4);for i=1:400H(i,1)=-z(i+1)*u^(L-i);H(i,2)=-z(i)*u^(L-i);H(i,3)=M(i+1)*u^(L-i);H(i,4)=M(i)*u^(L-i);endEstimate=inv(H'*H)*H'*(zstar(3:402))'二：RFF遗忘因子递推算法%最小二乘遗忘因子的递推算法仿真对象%Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k)%========================================clearclc%==========400 个产生M序列作为输入===============x=[0 1 0 1 1 0 1 1 1]; %initial valuen=403; %n为脉冲数目M=[]; %存放M 序列for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end%===========产生均值为0，方差为1 的高斯白噪声=========v=randn(1,400);%==============产生观测序列z=================z=zeros(402,1);z(1)=-1;z(2)=0;for i=3:402z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); end%==============递推求解=================P=10*eye(4); %估计方差Theta=zeros(4,401); %参数的估计值，存放中间过程估值Theta(:,1)=[0.001;0.001;0.001;0.001];K=zeros(4,400); %增益矩阵K=[10;10;10;10];u=0.98; %遗忘因子for i=3:402h=[-z(i-1);-z(i-2);M(i-1);M(i-2)];K=P*h*inv(h'*P*h+u);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2)); P=(eye(4)-K*h')*P/u;end%==========================输出结果及作图=============================disp('参数a1 a2 b1 b2的估计值:')Theta(:,401)i=1:401;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) title('待估参数过渡过程')三：RFM限定记忆法%限定记忆最小二乘的递推算法辨识对象%Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k)%========================================clearclc%==========产生M序列作为输入===============x=[0 1 0 1 1 0 1 1 1]; %initial valuen=403; %n为脉冲数目M=[]; %存放M 序列for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end%===========产生均值为0，方差为1 的高斯白噪声=========v=randn(1,402);%==============产生观测序列z=================z=zeros(402,1);z(1)=-1;z(2)=0;for i=3:402z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i);end%递推求解P_a=100*eye(4); %估计方差Theta_a=[3;3;3;3];L=20; %记忆长度for i=3:L-1 %利用最小二乘递推算法获得初步参数估计值和P 阵h=[-z(i-1);-z(i-2);M(i-1);M(i-2)];K=P_a*h*inv(h'*P_a*h+1);Theta_a=Theta_a+K*(z(i)-h'*Theta_a);P_a=(eye(4)-K*h')*P_a;endfor k=0:380hL=[-z(k+L-1);-z(k+L-2);M(k+L-1);M(k+L-2)];%增加新数据的信息K_b=P_a*hL*inv(1+hL'*P_a*hL);Theta_b=Theta_a+K_b*(z(k+L)-hL'*Theta_a);P_b=(eye(4)-K_b*hL')*P_a;hk=[-z(k+L);-z(k+L-1);M(k+L);M(k+L-1);];%去掉老数据的信息K_a=P_b*hk*inv(1+hk'*P_b*hk);Theta_a=Theta_b-K_a*(z(k+L+1)-hk'*Theta_b);P_a=(eye(4)+K_a*hk')*P_b;Theta_Store(:,k+1)=Theta_a;end%========================输出结果及作图===========================disp('参数 a1 a2 b1 b2 的估计值为：')Theta_ai=1:381;figure(1)plot(i,Theta_Store(1,:),i,Theta_Store(2,:),i,Theta_Store(3, :),i,Theta_Store(4,:))title('待估参数过渡过程')四：RCLS偏差补偿最小二乘递推算法%偏差补偿最小二乘的递推算法辨识对象%Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k)%========================================clearclc%==========产生M序列作为输入===============x=[0 1 0 1 1 0 1 1 1]; %initial valuen=403; %n为脉冲数目M=[]; %存放M 序列for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end%===========产生均值为0，方差为1 的正态分布噪声========= v=random('Normal',0,1,1,400);%==============产生观测序列z=================z=zeros(402,1);z(1)=-1;z(2)=0;for i=3:402z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2);end%===================递推求解==================%赋初值P=100*eye(4); %估计方差Theta=zeros(4,401); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];K=[10;10;10;10]; %增益J=0;ThetaC=zeros(4,401); %偏差补偿后的估计值ThetaC(:,1)=[2;3;1;3.5];D=[1 0 0 0;0 1 0 0;0 0 0 0;0 0 0 0];for i=3:402h=[-z(i-1);-z(i-2);M(i-1);M(i-2)];J=J+(z(i-1)-h'*Theta(:,i-1))^2/(1+h'*P*h);K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(4)-K*h')*P;endes=J/((i-1)*(1+(ThetaC(:,i-2))'*D*Theta(:,i-1)));ThetaC(:,i-1)=Theta(:,i-1)+(i-1)*es*P*D*ThetaC(:,i-2);%==============输出参数估计结果及作图================disp('参数a1 a2 b1 b2的估计值为：')Theta(:,401)i=1:401;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) title('待估参数过渡过程')五：RELS增广最小二乘的递推算法%增广最小二乘的递推算法辨识对象%Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)-v(k+1)+0.2*v(k) %========================================clearclc%==========产生M序列作为输入===============x=[0 1 0 1 1 0 1 1 1]; %initial valuen=403; %n为脉冲数目M=[]; %存放M 序列for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end%===========产生均值为0，方差为1 的高斯白噪声=========v=randn(1,402);%==============产生观测序列z=================z=zeros(402,1);z(1)=-1;z(2)=0;for i=3:402z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)-v(i-1)+0.2*v(i -2);end%递推求解P=100*eye(6); %估计方差Theta=zeros(6,401); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3;3;3];% K=zeros(4,400); %增益矩阵K=[10;10;10;10;10;10];for i=3:402h=[-z(i-1);-z(i-2);M(i-1);M(i-2);v(i-1);v(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(6)-K*h')*P;end%========================================================== =============disp('参数a1、a2、b1、b2、d1、d2估计结果：')Theta(:,401)i=1:401;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i, Theta(5,:),i,Theta(6,:))title('待估参数过渡过程')六;RGLS 广义最小二乘的递推算法%广义最小二乘的递推算法仿真模型%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)%========================================clearclc%==========400 个产生M序列作为输入===============x=[0 1 0 1 1 0 1 1 1]; %initial valuen=403; %n为脉冲数目M=[]; %存放M 序列for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end%===========产生均值为0，方差为1 的高斯白噪声========= v=randn(1,400);e=[]; e(1)=v(1); e(2)=v(2);for i=3:400e(i)=0*e(i-1)+0*e(i-2)+v(i);end%==============产生观测序列z=================z=zeros(400,1);z(1)=-1;z(2)=0;for i=3:400z(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:400zf(i)=z(i)-0*z(i-1)-0*z(i-2);end%变换后的输入序列uf=[]; uf(1)=M(1); uf(2)=M(2);for i=3:400uf(i)=M(i)-0*M(i-1)-0*M(i-2);end%赋初值P=100*eye(4); %估计方差Theta=zeros(4,400); %参数的估计值，存放中间过程估值Theta(:,2)=[3;3;3;3];K=[10;10;10;10]; %增益PE=10*eye(2);ThetaE=zeros(2,400);ThetaE(:,2)=[0.5;0.3];KE=[10;10];%递推Thetafor i=3:400h=[-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;endhe=[-e(i-1);-e(i-2)];%递推ThetaEKE=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(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) title('待估参数过渡过程')七：Tally辅助变量最小二乘的递推算法%Tally辅助变量最小二乘的递推算法%Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+e(k)，e(k)为有色噪声%e(k)=v(k)+0.5*v(k-1)+0.2*v(k-2),v(k)为零均值的不相关随机噪声%========================================clearclc%==========产生M序列作为输入===============x=[0 1 0 1 1 0 1 1 1]; %initial valuen=403; %n为脉冲数目M=[]; %存放M 序列for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end%===========产生均值为0，方差为1 的高斯白噪声========= v=randn(1,400);e=[];e(1)=0.3;e(2)=0.5;for i=3:400e(i)=v(i)+0.5*v(i-1)+0.2*v(i-2);end%==============产生观测序列z=================z=zeros(402,1);z(1)=-1;z(2)=0;for i=3:400z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+e(i); end%递推求解P=100*eye(4); %估计方差Theta=zeros(4,400); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];Theta(:,2)=[3;3;3;3];Theta(:,3)=[3;3;3;3];Theta(:,4)=[3;3;3;3];% K=zeros(4,400); %增益矩阵K=[10;10;10;10];for i=5:400h=[-z(i-1);-z(i-2);M(i-1);M(i-2)];hstar=[-z(i-2-1);-z(i-2-2);M(i-1);M(i-2)]; %辅助变量 %递推算法K=P*hstar*inv(h'*P*hstar+1);Theta(:,i)=Theta(:,i-1)+K*(z(i)-h'*Theta(:,i-1));P=(eye(4)-K*h')*P;end%==================结果输出及作图=================== disp('参数a1 a2 b1 b2的估计结果：')Theta(:,400)i=1:400;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) title('待估参数过渡过程')八：Cor-ls相关最小二乘法（二步法）%两步法的递推算法%Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+e(k)，e(k)为零均值的不相关随机噪声%e(k)=v(k)+0.5*v(k-1)+0.2*v(k-2)%========================================clearclc%==========产生M序列作为输入===============x=[0 1 0 1 1 0 1 1 1]; %initial valuen=403; %n为脉冲数目M=[]; %存放M 序列for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end%===========产生均值为0，方差为1 的高斯白噪声========= v=randn(1,400);e=[];e(1)=0.3;e(2)=0.5;for i=3:400e(i)=v(i)+0.5*v(i-1)+0.2*v(i-2);end%==============产生观测序列z===========z=zeros(402,1);z(1)=-1;z(2)=0;for i=3:400z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+e(i);end%递推求解P=100*eye(4); %估计方差Theta=zeros(4,400); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];Theta(:,2)=[3;3;3;3];Theta(:,3)=[3;3;3;3];Theta(:,4)=[3;3;3;3];K=zeros(4,400); %增益矩阵K=[10;10;10;10];for i=5:400h=[-z(i-1);-z(i-2);M(i-1);M(i-2)];hstar=[M(i-1);M(i-2);M(i-3);M(i-4)]; %辅助变量%递推K=P*hstar*inv(h'*P*hstar+1);Theta(:,i)=Theta(:,i-1)+K*(z(i)-h'*Theta(:,i-1));P=(eye(4)-K*h')*P;end%==================结果输出及作图===================disp('参数a1 a2 b1 b2的估计结果：')Theta(:,400)i=1:400;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) title('待估参数过渡过程')九：MLS多级最小二乘法clearclc%==========================================%Z(k+3)=-0.9*Z(k+2)-0.15*Z(k+1)-0.02*z(k)+0.7*u(k+2)-1.5*u(k +1)+e(k)%e(k+2)+1.0*e(k+1)+0.41*e(k)=r*v(k+2)%==========产生M 序列作为输入===============x=[0 1 0 1 1 0 1 1 1]; %initial valuen=405; %n为脉冲数目M=[]; %存放M 序列for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end%===========产生均值为0，方差为1 的高斯白噪声============= v=randn(1,405);e=[];e(1)=0.3;e(2)=0.7;r=0.9; %控制信噪比for i=3:405e(i)=-1.0*e(i-1)-0.41*e(i-2)+r*v(i);end%=================产生观测序列===================z=[];z(1)=-1;z(2)=0;z(3)=1.5;for i=4:405z(i)=-0.9*z(i-1)-0.15*z(i-2)-0.02*z(i-3)+0.7*M(i-1)-1.5*M(i -2)+e(i);end%================第一级辨识辅助模型参数辨识==================H=zeros(400,9);for i=1:400H(i,1)=-z(i+4);H(i,2)=-z(i+3);H(i,3)=-z(i+2);H(i,4)=-z(i+1);H(i,5)=-z(i);H(i,6)=M(i+4);H(i,7)=M(i+3);H(i,8)=M(i+2);H(i,9)=M(i+1);enddisp('第一级辅助模型参数 e1 e2 e3 e3 e4 f1 f2 f3 f4 辨识结果：')E=inv(H'*H)*H'*(z(6:405))'e1=E(1);e2=E(2);e3=E(3);e4=E(4);e5=E(5);f1=E(6);f2=E(7);f3=E(8);f4=E(9);%=================第二级辨识过程模型参数辨识====================z2=[f1;f2;f3;f4;0;0;0];H2=[ 0 0 0 1 0;-f1 0 0 e1 1;-f2 -f1 0 e2 e1;-f3 -f2 -f1 e3 e2;-f4 -f3 -f2 e4 e3;0 -f4 -f3 e5 e4;0 0 -f4 0 e5;];disp('第二级过程模型参数 a1 a2 a3 b1 b2 辨识结果：') E2=inv(H2'*H2)*H2'*z2a1=E2(1);a2=E2(2);a3=E2(3);b1=E2(4);b2=E2(5);%================第三级辨识噪声模型参数辨识=======================z3=[e1-a1;e2-a2;e3-a3;e4;e5;f2-b2;f3;f4];H3=[1 0;a1 1;a2 a1;a3 a2;0 a3;b1 0;b2 b1;0 b2;];disp('第三级噪声模型参数 c1 c2 辨识结果：')E3=inv(H3'*H3)*H3'*z3十yule_walker辨识算法%Yule-Walker 辨识算法%辨识模型：z(k)=-0.9*z(k-1)-0.36*z(k-2)-0.054*z(k-3)+v(k) %========================================================== ==%产生随机噪声v=random('Normal',0,1,1,1024); %均值为零，方差为 1%产生观测序列z=[];z(1)=0;z(2)=1;z(3)=1.5;for i=4:1024z(i)=-0.9*z(i-1)-0.36*z(i-2)-0.054*z(i-3)+v(i); end%计算 z(k)的自相关函数Rz0=0;Rz1=0;Rz2=0;Rz3=0;for i=1:1024Rz0=Rz0+z(i)^2;endRz0=Rz0/1024;for i=1:1023Rz1=Rz1+z(i+1)*z(i);endRz1=Rz1/1024;for i=1:1022。