实验一利用相关函数辨识脉冲响应

1、北京工商大学系统辨识课程实验报告（2014-2015 1 学期）课程名称 ：系统辨识题 目：利用相关分析法辨识脉冲响应专业班级 ：控制工程学生姓名： 指导教师：文U文U成 绩：2015年1月18日、实验目的通过仿真实验掌握利用相关分析法辨识脉冲响应的原理和方法。、实验内容图1为本实验的原理框图。过程传递函数为G(s),其中K 120, 8.3Sec,T2 6.2Sec； u(k)和z(k)分别为过程的输入和输出变量；v(k)为过程测量白噪声，服从正态分布，均值为零，方差为v2，记作v(k)N(o, v2)；?g0(k)为过程的脉冲响应理论值，g(k)为过程脉冲响应估计值，(k)为过程脉 冲响应

2、估计误差。过程的输入驱动采用M序列，输出受到白噪声v(k)的污染。根据过程的输入Q和输出数据u(k),z(k)，利用相关分析算法根据输出过程的脉冲响应值 g(k)，并 与过程脉冲响应理论值 go(k)比较，得到过程脉冲响应估计误差值 g(k)，当k 时,图1相关分析法辨识脉冲响应原理框图三、实验要求进行方案设计，模拟过程传递函数，获得输出数据，用 时列作为辨识的输 入信号，噪声采用标准正态分布的白噪声， 计算互相关函数，不同 值的脉冲响 应估计值、脉冲响应理论值和脉冲响应估计误差， 计算信噪比，画出实验流程图， 用MATLA编程实现。四、实验原理1、采用串联传递函数G(s)仿真G(s)T2 s

3、 1 /T1 s 1/T22、来描述g( c)则 y(t) 0 g( )x(t )d上式两端同乘x(t)，进而取时间均值，有1 T 1 TTim T 0 y(t)x(t )dt 0 g( )Tm 。x(t )x(t )dtd 则Rxy( )0 g( )Rx()d这就是著名的维纳霍夫积分方程。如果输入是白噪声，这时x(t)的自相关函数为Rx( ) k ( ),Rx() k ()Rxy()k这样，只要记录x(t)y(t)的值，并计算它们的互相关函数，即可求得脉则根据维纳霍夫积分方程可得Rxy( )0 g()Rx()dkg()或者g()冲响应函数g( T )。而在系统有正常输入的情形下，辨识脉冲响应

4、的原理图如下图所示止席输人止席输人五、实验框图止席输人幵殆L产生21丿产少扩磁21序句作为 牛俞A»u(k>aR ;Wx( k)求徉v5)产土广"I聚声I 口憾弹f%v(町吾 力11 彳Ef 3lJz(k)信1邑0占古i十-L吴delta_gTV 釜1( 结瑕 )六、实验代码function ex2clc;clear all ;close all ;%创建M序列Np=63; %循环周期delta_T = 1;%时钟节拍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; %周期数u=repmat(M_XuLie,1,r+1);%将M序列赋给输入，作为输入信号%第一步，从 u(k) 得到 x(k),y(k)K = 120;T0 = 1; % 采样时间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%获取没有白噪声时候输出完毕 %作图figure(1);

7、plot(u, 'r' ); hold on;plot(x, 'k' ); plot(y, 'b' );legend( 'u(k)' , 'x(k)' , 'y(k)' ); %第二步，将白噪声添加入输出信号 %产生白噪声信号 v fangcha = 0.5; %随意指定的方差 v = fangcha * randn(1,63*4);%信号叠加 , 输出实际信号 z(k) z = y + v;figure(2); %打印无白噪声污染信号 plot(y, 'b' );hold on;

8、%打印白噪声信号 plot(v, 'm' );%打印白噪声污染后的信号 plot(z, 'k' );legend( 'y(k)' , 'v(k)' , 'z(k)' );%计算 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);%十算脉冲响应估十值 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%十算脉冲响应估计误差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('脉冲响应理论值g0(k)','脉冲响应估计值g1');七、实验结果1、输入u(k),中间输入x(k)，无干扰输入(k)2、白噪声标准差为1.5时，理想输出y(k)，带干扰的输出z(k)，干扰v(k)3、输入白噪声标准差为1.5,周期数r为3时，脉冲响应理论值与估计值脉冲晌电理论倩010203040606070脉冲响应估计误差：0.0