實(shí)驗(yàn)一利用相關(guān)函數(shù)辨識(shí)脈沖響應(yīng)

1、北京工商大學(xué)系統(tǒng)辨識(shí)課程實(shí)驗(yàn)報(bào)告（2014-2015 1 學(xué)期）課程名稱 ：系統(tǒng)辨識(shí)題 目：利用相關(guān)分析法辨識(shí)脈沖響應(yīng)專業(yè)班級(jí) ：控制工程學(xué)生姓名： 指導(dǎo)教師：文U文U成 績(jī)：2015年1月18日、實(shí)驗(yàn)?zāi)康耐ㄟ^(guò)仿真實(shí)驗(yàn)掌握利用相關(guān)分析法辨識(shí)脈沖響應(yīng)的原理和方法。、實(shí)驗(yàn)內(nèi)容圖1為本實(shí)驗(yàn)的原理框圖。過(guò)程傳遞函數(shù)為G(s),其中K 120, 8.3Sec,T2 6.2Sec； u(k)和z(k)分別為過(guò)程的輸入和輸出變量；v(k)為過(guò)程測(cè)量白噪聲，服從正態(tài)分布，均值為零，方差為v2，記作v(k)N(o, v2)；?g0(k)為過(guò)程的脈沖響應(yīng)理論值，g(k)為過(guò)程脈沖響應(yīng)估計(jì)值，(k)為過(guò)程脈 沖響應(yīng)

2、估計(jì)誤差。過(guò)程的輸入驅(qū)動(dòng)采用M序列，輸出受到白噪聲v(k)的污染。根據(jù)過(guò)程的輸入Q和輸出數(shù)據(jù)u(k),z(k)，利用相關(guān)分析算法根據(jù)輸出過(guò)程的脈沖響應(yīng)值 g(k)，并 與過(guò)程脈沖響應(yīng)理論值 go(k)比較，得到過(guò)程脈沖響應(yīng)估計(jì)誤差值 g(k)，當(dāng)k 時(shí),圖1相關(guān)分析法辨識(shí)脈沖響應(yīng)原理框圖三、實(shí)驗(yàn)要求進(jìn)行方案設(shè)計(jì)，模擬過(guò)程傳遞函數(shù)，獲得輸出數(shù)據(jù)，用 時(shí)列作為辨識(shí)的輸 入信號(hào)，噪聲采用標(biāo)準(zhǔn)正態(tài)分布的白噪聲， 計(jì)算互相關(guān)函數(shù)，不同 值的脈沖響 應(yīng)估計(jì)值、脈沖響應(yīng)理論值和脈沖響應(yīng)估計(jì)誤差， 計(jì)算信噪比，畫出實(shí)驗(yàn)流程圖， 用MATLA編程實(shí)現(xiàn)。四、實(shí)驗(yàn)原理1、采用串聯(lián)傳遞函數(shù)G(s)仿真G(s)T2 s

3、 1 /T1 s 1/T22、來(lái)描述g( c)則 y(t) 0 g( )x(t )d上式兩端同乘x(t)，進(jìn)而取時(shí)間均值，有1 T 1 TTim T 0 y(t)x(t )dt 0 g( )Tm 。x(t )x(t )dtd 則Rxy( )0 g( )Rx()d這就是著名的維納霍夫積分方程。如果輸入是白噪聲，這時(shí)x(t)的自相關(guān)函數(shù)為Rx( ) k ( ),Rx() k ()Rxy()k這樣，只要記錄x(t)y(t)的值，并計(jì)算它們的互相關(guān)函數(shù)，即可求得脈則根據(jù)維納霍夫積分方程可得Rxy( )0 g()Rx()dkg()或者g()沖響應(yīng)函數(shù)g( T )。而在系統(tǒng)有正常輸入的情形下，辨識(shí)脈沖響應(yīng)

4、的原理圖如下圖所示止席輸人止席輸人五、實(shí)驗(yàn)框圖止席輸人幵殆L產(chǎn)生21丿產(chǎn)少擴(kuò)磁21序句作為 牛俞A»u(k>aR ;Wx( k)求徉v5)產(chǎn)土廣"I聚聲I 口憾彈f%v(町吾 力11 彳Ef 3lJz(k)信1邑0占古i十-L吳delta_gTV 釜1( 結(jié)瑕 )六、實(shí)驗(yàn)代碼function ex2clc;clear all ;close all ;%創(chuàng)建M序列Np=63; %循環(huán)周期delta_T = 1;%時(shí)鐘節(jié)拍a=1; %幅度M(1)=1;M(2)=0;M(3)=0;M(4)=1;M(5)=1;M(6)=0; %初始化M序列M_XuLie(Np) = 0;fo

5、r n = 1 : Nptemp = xor(M(6), M(5);if (temp = 0)M_XuLie(n) = a;elseM_XuLie(n) = -a;endM(6) = M(5);M(5) = M(4);M(4) = M(3);M(3) = M(2);M(2) = M(1);M(1) = temp;end%生成M序列完畢r=3; %周期數(shù)u=repmat(M_XuLie,1,r+1);%將M序列賦給輸入，作為輸入信號(hào)%第一步，從 u(k) 得到 x(k),y(k)K = 120;T0 = 1; % 采樣時(shí)間T1 = 8.3;T2 = 6.2;K1=K/(T1*T2);%初始化 X

6、(k),Y(k) 為0K2=1x(63)=0;y(63)=0for k = 2 : 63*4 %取得 x(k) 序列 x(k)=exp(-T0/T1)*x(k-1)+T1*K1*(1-exp(-T0/T1)*u(k-1)+T1*K1 *(T1*(exp(-T0/T1)-1)+T0)*(u(k)-u(k-1)/T0%取得 y(k) 序列 y(k)=exp(-T0/T2)*y(k-1)+T2*K2*(1-exp(-T0/T2)*x(k-1)+T2*K2 *(T2*(exp(-T0/T2)-1)+T0)*(x(k)-x(k-1)/T0 end%獲取沒(méi)有白噪聲時(shí)候輸出完畢 %作圖figure(1);

7、plot(u, 'r' ); hold on;plot(x, 'k' ); plot(y, 'b' );legend( 'u(k)' , 'x(k)' , 'y(k)' ); %第二步，將白噪聲添加入輸出信號(hào) %產(chǎn)生白噪聲信號(hào) v fangcha = 0.5; %隨意指定的方差 v = fangcha * randn(1,63*4);%信號(hào)疊加 , 輸出實(shí)際信號(hào) z(k) z = y + v;figure(2); %打印無(wú)白噪聲污染信號(hào) plot(y, 'b' );hold on;

8、%打印白噪聲信號(hào) plot(v, 'm' );%打印白噪聲污染后的信號(hào) plot(z, 'k' );legend( 'y(k)' , 'v(k)' , 'z(k)' );%計(jì)算 Rmz(k)for k = 1 : NpRmz(k)=0; %初始化為 0 for i = (Np + 1) : (r+1)*Np) Rmz(k)=Rmz(k) + u(i-k)*z(i); endRmz(k)=Rmz(k)/(r*Np);end%十算c c=-Rmz(Np - 1);%十算脈沖響應(yīng)估十值 g1g仁Np*(Rmz+c)/(N

9、p+1)*aA2*delta_T);%十算理論脈沖 g0for k = 1: NpgO(k)=K/(T1-T2)*(exp(-k*delta_T/T1)-exp(-k*delta_T/T2);end%十算脈沖響應(yīng)估計(jì)誤差delta_gdelta_g=sqrt(sum(g0-g1).A2)/sum(g0.A2);figure(3);plot(gO, 'k');hold on;plot(g1, 'r');%axis(0,100,0,10);legend('脈沖響應(yīng)理論值g0(k)','脈沖響應(yīng)估計(jì)值g1');七、實(shí)驗(yàn)結(jié)果1、輸入u(k),中間輸入x(k)，無(wú)干擾輸入(k)2、白噪聲標(biāo)準(zhǔn)差為1.5時(shí)，理想輸出y(k)，帶干擾的輸出z(k)，干擾v(k)3、輸入白噪聲標(biāo)準(zhǔn)差為1.5,周期數(shù)r為3時(shí)，脈沖響應(yīng)理論值與估計(jì)值脈沖晌電理論倩010203040606070脈沖響應(yīng)估計(jì)誤差：0.0