相关分析法辩识离散线性系统的脉冲响应特性.doc

北京工业大学电控学院 朱熙 S201002153相关分析法辩识离散线性系统的脉冲响应特性朱熙S201002153指导老师：余春暄2010年11月13日原理介绍辩识系统如图1所示。待辨识系统图1 辨识系统设其在静态工作点附近作局部线性化所得动态模型为：。为了进行系统冲激响应特性的辩识，在对象输入端施加电平幅值V为1的M序列信号（其波形见图2）作为动态测试信号。本次实验的基本任务就是在输入、输出的实际观测值和的基础上，用批量或递推算法估计出脉冲响应序列。0123456789101112131415图2 M序列在相关法中，采用伪随机信号M序列对系统模型进行在线辩识。由于M序列的自相关函数很接近于周期性的函数，可近似地当作白噪声信号叠加在系统输入信号上，按其与输出信号间的互相关函数，可得到系统的脉冲响应曲线。设对象输入为具有平稳遍历性的随机信号，其自相关函数为。当系统初始过渡过程结束后，输入和输出之间的互相关函数为。其中，表示统计平均，表示时间平均。由维纳-霍甫方程有。当为伪随机信号M序列时，自相关函数的离散值为。式中，、和分别为M序列的电平值、基本电平时间和离散周期长度，为时差的离散值。相应的功率谱密度函数为。式中，。若将M序列的自相关函数近似为并取其周期大于对象的抽样时间，即则有。这样就可以得到脉冲响应曲线。批量算法（Batch algorithm）递推算法（Recursive algorithm）M序列、L序列测试信号的生成方法M序列又称最大长度伪随机二进制序列，它实际上是一种特殊的周期序列，其电平值只取+V和-V（或者说其逻辑电平值只取1或0）。M序列的当前值可以由过去的值进行异或运算得到。选取恰当的，便可以产生最长周期为(2n-1)的M序列。很容易由计算机程序来产生M序列。如取n=4，=3，便可按依次推出，最后得出图2所示的M序列。L序列又称逆重复M序列，可由M序列和序列进行异或运算得到：。初始化生成数据u(k-d) , Y(k) , k=1-r-d,-1,0,Np-1k1 Np + r + d生成新的数据 ， 构造h, 向量，并送入h1向量，等待循环计算k-k1iTp?打印结果k k +1图3 递推算法流程图初始化生成数据 u(k-d) , Y(k) , k=1-r-d,-1,0,Np-1k1 Np + r + d生成新的数据 k-k1iTp?打印结果k k +1图4 批量算法流程图程序代码及实验结果clear all;close all;clc;a1=-1;a2=0.5;b1=1;b2=0.5;U0=0.1;Y0=0.2;%系统参数d=0;%系统延时Ts=0.02;%系统稳定时间% M 序列的产生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;Y1=X1;X4=Y3;X3=Y2;X2=Y1;X1=xor(Y3,Y4);U(i)=Y4;if Y4=0M(i)=-a;elseM(i)=a;endend%draw figurei=1:1:m_length;plot(i,M)所产生M序列如下图%L序列的产生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;elseL(i)=a;endendi=1:1:m_length;plot(i,L)所产生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)所产生白噪声如下图%formulationy(1)=0;y(2)=0;yp(1)=0;yp(2)=0;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)所产生图像为%批量法U=M;%或U=L;mid1=ones(Np);mid2=eye(Np);mid=mid1 + mid2;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)所产生批量法结果如下图%递推算法U=M;%或U=L;m=10;% m次h=zeros(Np,1);for i=1:1:mfor j=1:1:NpUm(j)=U(i+1-j+Np+1);endhm1=h+(1/(i+1)*(2*Np*Um*y(Np