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

相关分析法辩识离散线性系统的脉冲响应特性 朱熙 S201002153 指导老师：余春暄 2010年11月13日 原理介绍 辩识系统如图1所示。 图1 辨识系统 设其在静态工作点附近作局部线性化所得动态模型为： ， 式中参数建议值为 ，，，。 为了进行系统冲激响应特性的辩识，在对象输入端施加电平幅值V为1的M序列信号（其波形见图2）作为动态测试信号。本次实验的基本任务就是在输入、输出的实际观测值和的基础上，用批量或递推算法估计出脉冲响应序列。 图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序列和序列进行异或运算得到： 。 程序代码及实验结果 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_length Y4=X4;Y3=X3;Y2=X2;Y1=X1; X4=Y3;X3=Y2;X2=Y1; X1=xor(Y3,Y4); U(i)=Y4; if Y4==0 M(i)=-a; else M(i)=a; end end %draw figure i=1:1:m_length; plot(i,M) 所产生M序列如下图 %L序列的产生 for i=1:1:m_length; tt=(-1)^i; if tt==-1 t_l(i)=0; else t_l(i)=1; end l_L(i)=xor(U(i),t_l(i)); if l_L(i)==0; L(i)=-a; else L(i)=a; end end i=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_length noise=noise +0.5*rands(1,m_length); end noise=noise/12; %figure i=1:1:m_length; plot(i,noise) 所产生白噪声如下图 %formulation y(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_length y(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); end for k =3:1:Np yp(k)=b1*impulse(k-1-d)+b2*impulse(k-2-d)-a1*yp(k-1)-a2*yp(k-2); end %draw figure i=1:1:m_length; plot(i,Y) %figure i=1:1:Np;%m_length plot(i,yp(1:Np)) 所产生图像为 %批量法 U=M;%或U=L; mid1=ones(Np);mid2=eye(Np); mid=mid1 + mi