遞推的極大似然法辨識程序

1、遞推的極大似然法辨識程序例5.2系統(tǒng)模型如圖5.5所示。試用遞推的極大似然法對系統(tǒng)辨識的參數(shù)集 圖5.5 例5.2系統(tǒng)模型u(k)+ e(k)y(k)+ z(k)v(k) u v(k)隨機(jī)信號，輸入信號為幅值為的M序列或隨機(jī)信號，要求畫出程序流程圖，打印出程序（程序中帶有注釋）和辨識中的參數(shù)、誤差曲線。解: 首先解釋編程所用的部分字母：由于在MATLAB語言中無法用希臘字母描述、無法用上標(biāo)及下標(biāo)，故用 o和o1表示；令；產(chǎn)生M序列時，a(i),b(i),c(i),d(i)表示四級移位寄存器的第1，2，3，4級寄存器的輸出; 編程如下(光盤上該程序：FLch5RMLeg2.m，可在MATLAB6

2、.I下直接運(yùn)行):編程如下：clear %清零a(1)=1;b(1)=0;c(1)=1;d(1)=0;u(1)=d(1);z(1)=0;z(2)=0; %初始化for i=2:1200 %產(chǎn)生m序列u(i) a(i)=xor(c(i-1),d(i-1); b(i)=a(i-1); c(i)=b(i-1); d(i)=c(i-1); u(i)=d(i);endu; %若取去；可以在程序運(yùn)行中觀測到m序列v=randn(1200,1); %產(chǎn)生正態(tài)分布隨機(jī)數(shù)V=0; %計算噪聲方差for i=1:1200 V=V+v(i)*v(i);endV1=V/1200; for k=3:1200 %根據(jù)v和

3、u計算z z(k)=1.2*z(k-1)-0.6*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);endo1=0.001*ones(6,1);p0=eye(6,6); %賦初值 zf(1)=0.1;zf(2)=0.1;vf(2)=0.1;vf(1)=0.1;uf(2)=0.1;uf(1)=0.1; %迭代計算參數(shù)值和誤差值for k=3:1200 h=-z(k-1);-z(k-2);u(k-1);u(k-2);v(k-1);v(k-2); hf=h; K=p0*hf*inv(hf'*p0*hf+1); p=eye(6,6)-K*hf

4、9;*p0; v(k)=z(k)-h'*o1; o=o1+K*v(k) ; p0=p; o1=o; a1(k)=o(1); a2(k)=o(2); b1(k)=o(3); b2(k)=o(4); d1(k)=o(5); d2(k)=o(6); e1(k)=abs(a1(k)+1.2); e2(k)=abs(a2(k)-0.6); e3(k)=abs(b1(k)-1.0); e4(k)=abs(b2(k)-0.5); e5(k)=abs(d1(k)+1.0); e6(k)=abs(d2(k)-0.2); zf(k)=z(k)-d1(k)*zf(k-1)-d2(k)*zf(k-2); uf

5、(k)=u(k)-d1(k)*uf(k-1)-d2(k)*uf(k-2); vf(k)=v(k)-d1(k)*vf(k-1)-d2(k)*vf(k-2); hf=-zf(k-1);-zf(k-2);uf(k-1);uf(k-2);vf(k-1);vf(k-2); end o1 %若取去；可以在程序運(yùn)行中觀測到參數(shù) V1 %繪圖 subplot(4,1,1)k=1:1200;plot(k,a1,'k:',k,a2,'b',k,b1,'r',k,b2,'m:',k,d1,'g',k,d2,'k');xl

6、abel('k')ylabel('parameter')legend('a1=-1.2,','a2=0.6','b1=1.0','b2=0.5','d1=-1.0','d2=0.2'); %圖標(biāo)炷title('The parameter idendification of the RML');endsubplot(4,1,2)k=1:1200;plot(k,e1,'k',k,e2,'b',k,e3,'r',k,e4,'m',k,e5,'g',k,e6,'k');xlabel('k')ylabel('error')%title('誤差曲線')endsubplot(4,1,3)k=1:1200;plot(k,u);xlabel('k')ylabel('input')%title('系統(tǒng)輸入信號')endsubplot(4,1,4)k=1:1200;plot(k,v);