相關(guān)分析法辯識(shí)離散線性系統(tǒng)的脈沖響應(yīng)特性

1、相關(guān)分析法辯識(shí)離散線性系統(tǒng)的脈沖響應(yīng)特性朱熙S201002153指導(dǎo)老師：余春暄2010年11月13日北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153 原理介紹辯識(shí)系統(tǒng)如圖1所示。圖1 辨識(shí)系統(tǒng)設(shè)其在靜態(tài)工作點(diǎn)(U0,Y0)附近作局部線性化所得動(dòng)態(tài)模型為：y(k)+a1(k-1)+a2u(k-2)=b1u(k-1-d)+b2u(k-2-d)，式中參數(shù)建議值為a1=-1，a2=0.5，a3=-1，a4=0.5。為了進(jìn)行系統(tǒng)沖激響應(yīng)特性的辯識(shí)，在對(duì)象輸入端施加電平幅值V為1的M序列信號(hào)（其波形見圖2）作為動(dòng)態(tài)測試信號(hào)u(k)。本次實(shí)驗(yàn)的基本任務(wù)就是在輸入、輸出的實(shí)際觀測值z(mì)(k)=y0(k)+y(

2、k)+(k)和u(k)的基礎(chǔ)上，用批量或遞推算法估計(jì)出脈沖響應(yīng)序列h(k)。t圖2 M序列在相關(guān)法中，采用偽隨機(jī)信號(hào)M序列對(duì)系統(tǒng)模型進(jìn)行在線辯識(shí)。由于M序列的自相關(guān)函數(shù)很接近于周期性的函數(shù)，可近似地當(dāng)作白噪聲信號(hào)疊加在系統(tǒng)輸入信號(hào)上，按其與輸出信號(hào)間的互相關(guān)函數(shù)，可得到系統(tǒng)的脈沖響應(yīng)曲線。設(shè)對(duì)象輸入u(t)為具有平穩(wěn)遍歷性的隨機(jī)信號(hào)，其自相關(guān)函數(shù)為- 1 -北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153Ru()=Eu(t)u(t+)=u(t)u(t+)=ru()。當(dāng)系統(tǒng)初始過渡過程結(jié)束后，輸入和輸出之間的互相關(guān)函數(shù)為Ruz()=Eu(t)z(t+)=u(t)z(t+)=ruz()。其中，E.

3、表示統(tǒng)計(jì)平均，.表示時(shí)間平均。由維納-霍甫方程有ruy()=h()ru(-v)d0TpTph()r(-v)d。 u當(dāng)u(t)為偽隨機(jī)信號(hào)M序列時(shí)，自相關(guān)函數(shù)的離散值為V2,l=0,±Np,±2Np, 。 ru(l)=ru()|=l=V2-,otherwiseNp式中，V、和Np分別為M序列的電平值、基本電平時(shí)間和離散周期長度，為時(shí)差的離散值。相應(yīng)的功率譜密度函數(shù)為sin()Su(j)=S(-i0)。 0i=-2式中， 2=1+S0若將M序列的自相關(guān)函數(shù)近似為 (1Np)V2。ruy()S0并取其周期大于對(duì)象的抽樣時(shí)間Ts，即 i=-h(-iTp)Tp=Np>Ts則有h

4、(),0<T。 ruy()=S0p這樣就可以得到脈沖響應(yīng)曲線,0<T。 h()=ruy()S0p- 2 -北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153 批量算法（Batch algorithm）Np-1ruy'(l)=ruy'()|=l=Npk=0u(k-l)y'(k)uy'h()=遞推算法（Recursive algorithm） Np(Np+1)V2r()-ruy'(0) hu(m)(0)m 1 ±2Npu(m-1)h(1)mhm=hm-1+y'(m)-hm-1 m+1 (Np+1)V uh(m-Np+1)(N-1)

5、mpM序列又稱最大長度偽隨機(jī)二進(jìn)制序列，它實(shí)際上是一種特殊的周期序列，其電平值m(k)只取+V和-V（或者說其邏輯電平值u(k)只取1或0）。m(k)=m(k-k0)m(k-n)。M序列的當(dāng)前值可以由過去的值進(jìn)行異或運(yùn)算得到。選取恰當(dāng)?shù)膋0，便可以產(chǎn)生最長周期Np為(2n-1)的M序列。很容易由計(jì)算機(jī)程序來產(chǎn)生M序列。如取n=4，k0=3，m(0)=m(1)=m(2)=m(3)=1，便可按m(k)=m(k-3)m(k-4)依次推出m(4)，m(5)最后得出圖2所示的M序列。L序列又稱逆重復(fù)M序列，可由M序列和t(k)=(-1)序列進(jìn)行異或運(yùn)算得到： kl(k)=m(k)t(k)。- 3 -北京

6、工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153圖3 遞推算法流程圖- 4 -北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153圖4 批量算法流程圖程序代碼及實(shí)驗(yàn)結(jié)果clear all;close all;clc;a1=-1;a2=0.5;b1=1;b2=0.5;U0=0.1;Y0=0.2;%系統(tǒng)參數(shù)d=0;%系統(tǒng)延時(shí)Ts=0.02;%系統(tǒng)穩(wěn)定時(shí)間% M 序列的產(chǎn)生Np=15;fm=1;r=4;m_length=(r+1)*Np;%total length of M;X1=1;X2=1;X3=1;X4=1;a=1;%M序列幅值for i=1:1:m_lengthY4=X4;Y3=X3;Y2=X2;Y

7、1=X1;X4=Y3;X3=Y2;X2=Y1;X1=xor(Y3,Y4);U(i)=Y4;- 5 -北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153if Y4=0M(i)=-a;elseM(i)=a;endend%draw figurei=1:1:m_length;plot(i,M)所產(chǎn)生M序列如下圖%L序列的產(chǎn)生for i=1:1:m_length;tt=(-1)i;if tt=-1t_l(i)=0;elset_l(i)=1;endl_L(i)=xor(U(i),t_l(i);if l_L(i)=0;L(i)=-a;else- 6 -北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153L(i)=

8、a;endendi=1:1:m_length;plot(i,L)所產(chǎn)生L序列如下圖%脈沖序列impulse=1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;%白噪聲noise=zeros(1,m_length);for i=1:1:m_lengthnoise=noise +0.5*rands(1,m_length);endnoise=noise/12;%figurei=1:1:m_length;plot(i,noise)所產(chǎn)生白噪聲如下圖- 7 -北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153%formulationy(1)=0;y(2)=0;yp(1)=0;yp(2)=0;

9、Y(1)=Y0+y(1)+noise(1);Y(2)=Y0+y(2)+noise(2);for k=3:1:m_lengthy(k)=b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2);Y(k)=Y0+y(k)+noise(k);endfor k =3:1:Npyp(k)=b1*impulse(k-1-d)+b2*impulse(k-2-d)-a1*yp(k-1)-a2*yp(k-2);end%draw figurei=1:1:m_length;plot(i,Y)%figurei=1:1:Np;%m_lengthplot(i,yp(1:Np)所產(chǎn)生圖像為-

10、8 -北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153%批量法U=M;%或U=L;mid1=ones(Np);mid2=eye(Np);mid=mid1 + mid2;- 9 -北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153 A(1,:)=U(1:(r*(Np+1);for k=2:1:Np+1%移位tt(1)=U(m_length)for i=1:1:(m_length-1)tt(i+1)=U(i);endU=tt;A(k,:)=U(1:(r*(Np+1);endA1=A(1:Np,1:(r*Np-1);y1=Y(1:(r*Np-1);g=mid*A1*y1'/(a2*r*(Np+1);%draw figurei=1:1:Np;plot(i,g)所產(chǎn)生批量法結(jié)果如下圖%遞推算法U=M;%或U=L;m=10;% m次h=zeros(Np,1);for i=1:1:m- 10 -北京工業(yè)大學(xué)電控學(xué)院 朱熙 S201002153for j=1:1:NpUm(j)=U(i+1-j+Np+1)