实验一 利用相关分析法辨识脉冲响应一、实验目的通过仿真实验掌握利用相关分析法辨识脉冲响应的原理和方法。
二、实验内容下图为本实验的原理框图。
过程传递函数为)(s G ,其中Sec 26T Sec,3812021..,===T K ;)()(k z k u 和分别为过程的输入和输出变量;)(k v 为过程测量白噪声,服从正态分布,均值为零,方差为2v σ,记作),(~)(20v N k v σ;)(k g 0为过程的脉冲响应理论值,)(ˆk g 为过程脉冲响应估计值,)(~k g为过程脉冲响应估计误差。
过程的输入驱动采用M 序列,输出受到白噪声)(k v 的污染。
根据过程的输入和输出数据{})(),(k z k u ,利用相关分析算法根据输出过程的脉冲响应值)(ˆk g ,并与过程脉冲响应理论值)(k g 0比较,得到过程脉冲响应估计误差值)(~g ~三、实验方案设计(1) 采用串联传递函数)(s G 仿真k g =)(ˆ]2T t k /∆21211111T s T s T T K s G //)(++=令211T T KK =,则)(s G 的表达框图为:编程语句可写成:[][][][]};);()();()();()(;/)()(*)(**)(*)(*)(*)(;/)()(*)(***)(*)(**)(*)({);;(;)(;)();/();/();*/(k y k y k x k x k u k u T k x k x T E T T k x E T k y E k y T k u k u T E T K T k u E K T k x E k x k k k y x T T E T T E T T K K =-=-=---+-+--+-=--+-+--+-=++<===-=-==11111111111112521for 0000EXP EXP 002222220011111111202101211(2)白噪声生成● 利用U[0,1]均匀分布的随机数生成正态分布的白噪声),(~)(212106v i i v N k v σξσ⎪⎭⎫⎝⎛-=∑=其中,标准差v σ分别取0,0.1,0.5。
● 编程语句};);.(*};);();;*{);;(;{),;(06FLOAT ,MOD(121for 02521for -=+===++<==++<=ksai Sigma v(k)xi/M ksai ksai M xi xi xi A xi i i i ksai k k k(3)M 序列生成● 用M 序列作为辨识的输入信号,N 序列的循环周期取63126=-=P N ,时钟节拍Sec 1=∆t ,幅度1=a ,逻辑“0”为a ,逻辑“1”为-a ,特征多项式自选,如156⊕⊕=s s s F )(。
● 生成M 序列的结构图● 编程语句};;)()(;)()(};);()({);;(;)()();()()({),;(a k u M a k u M i M i M i i P i for M M M M M k k k -======-=-->====+=++<=then 10ifthen 00if 1000then20if2102521for(4)互相关函数的计算∑++=-=P P N r N i PMz i z k i u rN k R )()()()(111其中,r 为周期数,1+=P N i 表示计算互相关函数所用的数据是从第二个周期开始的,目的是等过程仿真数据进入平稳状态。
(5)计算脉冲响应估计值● 脉冲响应估计值 []c k R ta N N k g Mz P P +∆+=)()()(ˆ21 ● 脉冲响应估计误差 ()∑∑==⎪⎭⎫ ⎝⎛-PPN k N k g k gk g k g 12120)()(ˆ)(=δ四、数据记录理想状态下,即在没有白噪声干扰下的数据输入白噪声标准差sigma(0.5):0 脉冲响应估计误差0.0266输入白噪声标准差sigma(0.5):0.1 脉冲响应估计误差 0.0281输入白噪声标准差sigma(0.5):0.5 脉冲响应估计误差 0.0293输入白噪声标准差sigma(0.5):1 脉冲响应估计误差 0.0391输入白噪声标准差sigma(0.5):2 脉冲响应估计误差 0.0616五、结果分析利用相关分析法分析脉冲响应,得到脉冲响应的估计误差是随着输入白噪声标准差的增加而增大的,带有白噪声污染的输出z,在白噪声标准差为0时与理想输出y是重合的,白噪声的标准差愈小对系统的输出干扰愈小。
六.程序流程源程序清单%利用相关分析法辨识脉冲响应clc;clear all;close all;a=1;Np=63;Ts=1;%采样时间初始化条件%过程仿真参数K=120;T0=1;T1=8.3;T2=6.2;%产生输入u(k)r=4;M = [0 0 0 1 0 1 1]; %-初始状态向量P = 7; %-实际应为6,循环周期Np=2^6-1=63 for k = 1:1:252u(k) = 1-2*M(7); %-取M6(相对M0而言)结果生成M序列 M(1) = xor(M(6),M(7));if M(1) == 2M(1) = 0;endi = P;while i>1M(i) = M(i-1);i = i-1;endend%产生输入数据x(k)%系统仿真方法一,指导书上的K1 = K/(T1*T2);E1 = exp(-T0/T1);E2 = exp(-T0/T2);x(1) = 0;y(1) = 0;for k = 2:1:252x(k) = E1*x(k-1)+T1*K1*(1-E1)*u(k-1)+T1*K1*(T1*(E1-1)+T0)*(u(k)-u(k-1))/T0;y(k) = E2*y(k-1)+T2*(1-E2)*x(k-1)+T2*(T2*(E2-1)+T0)*(x(k)-x(k-1))/T0;end%产生不带噪声的输出数据y(k)%or%sys=tf([120],[8.3*6.2,8.3+6.2,1]);%y1=lsim(sys,u,1:length(u));%输入为M序列时对系统进行仿真%y=y1';sigma=input('输入白噪声标准差sigma(0.5):');v=whitenoise(0,sigma^2,length(y)); %N(0,0.25)高斯白噪声%or%v=wan(sigma);z=y+v;%产生输出数据z(k)ii=1:length(u);plot(ii,u)title('输入u')axis([1 length(u) -1.5 1.5])figure(2)plot(ii,y,'b',ii,z,'r'),title('仿真结果y以及带白噪声的输出z')legend('y','z','Location', 'Best')%--------------------------------------------去掉数据直流分量u = u - mean(u);zc = z(Np+1:4*Np); %-为避开非平稳过程,从第二周期开始采集数据zc = zc - mean(zc);%--------------------------------------------计算互相关函数for k = 1:1:NpRuz(k) = (1/((r-1)*Np))*u(Np+2-k:r*Np+1-k)*zc';end%--------------------------------------------计算脉冲响应估计值和理论值for k = 1:1:Npge(k) = Np/((Np+1)*a^2*Ts)*(Ruz(k)-Ruz(Np));g0(k) = K/(T1-T2)*(exp(-(k-1)*Ts/T1)-exp(-(k-1)*Ts/T2));end%--------------------------------------------计算估计误差deltag = sqrt( (g0-ge)*(g0-ge)'/(g0*g0') );disp('脉冲响应估计误差');disp(deltag);%--------------------------------------------画脉冲响应曲线xk = 0:1:Np-1;figure(3)plot(xk,ge,xk,g0,'r',xk,Ruz,'g');legend('脉冲响应估计值','脉冲响应理论值','互相关函数','Location', 'Best')生成图像曲线M序列的图形曲线白噪声标准差为0时,理想输出y,带干扰的输出z 输入白噪声标准差为0时,脉冲响应理论值与估计值白噪声标准差为0.5时,理想输出y,带干扰的输出z白噪声标准差为0.5时,脉冲响应理论值与估计值七.心得体会为期两周的实验进一步巩固系统辨识方面的专业知识。
理论是实验的基础,实验是对理论的验证和深化。
该实验使我对系统辨识有了更深刻的认识。
同时,提高了我的动手能力和思维能力。
在实验的过程中,我们也遇到了一些错误,程序有任何小的错误都会影响整个程序的运行。
但是通过我们的仔细检查和不断的修正,终于得到了期望的实验结果。
这次实验不仅锻炼了我们的细心,耐心,同时,还让我们意识到,团队精神的重要性。
三人行必有我师。
感谢老师给我们这次锻炼的机会,感谢小组的每一位成员。
电子113杨雨雨(201105102)小组成员:侯婷婷王媛媛杨雨雨郑敏。