實驗一利用相關函數辨識脈沖響應分析解析

1、北京工商大學系統(tǒng)辨識課程實驗報告（2014-20151學期）課程名稱:系統(tǒng)辨識題目:利用相關分析法辨識脈沖響應專業(yè)班級:控制工程學生姓名:指導教師:劉劉成績:2015年1月18日一、實驗目的通過仿真實驗掌握利用相關分析法辨識脈沖響應的原理和方法。二、實驗內容圖1為本實驗的原理框圖。過程傳遞函數為G(s),其中K=120,T=83Se(T,=62Secu(k)ffz(k)分別為過程的輸入和輸出變量；v(k)12為過程測量白噪聲，服從正態(tài)分布，均值為零，方差為b2,記作v(k)N(0,b2)；vv(k)為過程的脈沖響應理論值，g(k)為過程脈沖響應估計值，(k)為過程脈沖響應估計誤差。過程的輸入驅

2、動采用M序列，輸出受到白噪聲v(k)的污染。根據過程的輸入和輸出數據b(k),z(k),利用相關分析算法根據輸出過程的脈沖響應值g(k)，并與過程脈沖響應理論值g°(k)比較，得到過程脈沖響應估計誤差值(k)，當圖1相關分析法辨識脈沖響應原理框圖三、實驗要求進行方案設計，模擬過程傳遞函數，獲得輸出數據，用M序列作為辨識的輸入信號，噪聲采用標準正態(tài)分布的白噪聲，計算互相關函數，不同九值的脈沖響應估計值、脈沖響應理論值和脈沖響應估計誤差，計算信噪比，畫出實驗流程圖，用MATLAB編程實現。四、實驗原理1、采用串聯傳遞函數G(s)仿真K11TTs+1/Ts+1/T1212令尺令，則G(s)

3、的表達框圖為:1u(k)K1x(k)fa-1y(k)F-s+1/TlilFs+1/T2-2、一個單輸入單輸出線性定常系統(tǒng)的動態(tài)特性可用它的脈沖響應函數g(a)來描述。x(t)線性系統(tǒng)皿g(6則y(t)=jsgQ)x(t-)dc0上式兩端同乘(t-T)，進而取時間均值，有l(wèi)imTTsjT0y(t)x(t-T)dtsg(c)limljT0TTsT0x(t-c)x(t-T)dtdc則R(T)=jsg(c)R(T-c)dcxy0x這就是著名的維納-霍夫積分方程。如果輸入是白噪聲，這時x(t)的自相關函數為R(T)二kb(t),R(t-c)二kb(t-c)xx則根據維納-霍夫積分方程可得R(t)二jsg

4、(c)R(Tc)dc二kg(t)xy0x或者R(T)g(T)=-k這樣，只要記錄x(t)、y(t)的值，并計算它們的互相關函數，即可求得脈沖響應函數g(T)。而在系統(tǒng)有正常輸入的情形下，辨識脈沖響應的原理圖如下圖所示。4瞬請I9102-102餾肩毎H尊專芻幽奪:*焦T*坦:H卓排目一中QP巽箸屮早Fo目列炎曲氣呂刖44早pH眞*（習）引妾甜口y召心）人號M渝曰丿（斗）車贛曰寺（i）x音g*（I”Y妙繇助眩占IAI智#ifIAIzz?/(1)6>1u-trxy驟黑衛(wèi)北京工商大學計算機與信息工程學院2014-20151學期六、實驗代碼functionex2clc;clearall;close

5、all;%創(chuàng)建M序列Np=63;%循環(huán)周期delta_T=1;%時鐘節(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;forn=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;%周期數u=repmat(M_XuLie,1,r+1);%將M序列賦給輸入，作為輸入信號%第一

6、步，從u(k)得到x(k),y(k)K=120;T0=1;%采樣時間T1=8.3;T2=6.2;K1=K/(T1*T2);%初始化X(k),Y(k)為0K2=1x(63)=0;y(63)=0fork=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)

7、+T0)*(x(k)-x(k-1)/T0end%獲取沒有白噪聲時候輸出完畢%作圖figure(1);plot(u,'r');holdon;plot(x,'k');plot(y,'b');legend('u(k)','x(k)','y(k)');%第二步，將白噪聲添加入輸出信號%產生白噪聲信號vfangcha=0.5;%隨意指定的方差v=fangcha*randn(1,63*4);%信號疊加，輸出實際信號z(k)z=y+v;figure(2);%打印無白噪聲污染信號plot(y,'b'

8、);holdon;%打印白噪聲信號plot(v,'m');%打印白噪聲污染后的信號plot(z,'k');legend('y(k)','v(k)','z(k)');%計算Rmz(k)fork=1:NpRmz(k)=0;%初始化為0fori=(Np+1):(r+1)*Np)Rmz(k)=Rmz(k)+u(i-k)*z(i);endRmz(k)=Rmz(k)/(r*Np);end%計算cc=-Rmz(Np-1);%計算脈沖響應估計值glg1=Np*(Rmz+c)/(Np+1)*a2*delta_T);%計算理論脈沖g0fork=1:Npg0(k)=K/(T1-T2)*(exp(-k*delta_T/T1)-exp(-k*delta_T/T2);end%計算脈沖響應估計誤差delta_gdelta_g=sqrt(sum(g0-g1).入2)/sum(g0.入2);figure(3);plot(g0,'k');holdon;plot(g1,'r');%axis(0,100,0,10);legend('脈沖響應理論值g0(k)','脈沖響應估計值g1');