1、实验一 辩识离散线性系统脉冲响应特性的相关分析法实验目的（1） 掌握利用相关分析法辩识离散线性系统冲激响应特性的基本步骤和使用要点；（2） 了解递推和批量两种算法如何具体编程；（3） 熟悉估计算法中各种参数选择对辩识效果的影响。实验原理待辩识系统的基本情况和辩识工作的要求待辩识系统的运行情况如图 1 所示，设其在静态工作点（U 0，Y 0）附近作局部线性化所得动态模型为y(k)+a1y(k-1)+a2y(k-2)=b1u(k-1-d)+b2u(k-2-d) (1)式中参数建议值为 。5.0,5.0,2ba图 1 待辩识系统的基本运行情况为了进行系统冲激响应特性的辩识，在对象输入端施加电平幅值
2、V=1 的 M 序列信号（其波形见图 2）作为动态测试信号 u(k)。本次实验的基本任务就是在输入、输出的实际观测值 Y(k)=Y0+y(k) + (k) 和 u(k)的基础上，用批量或递推算法估计出脉冲响应序列h (k)。u(k)图 2 M 序列信号的波形M 序列测试信号的生成方法M 序列又称最大长度伪随机二进制序列，它实际上是一种特殊的周期序列，其电平值 u(k)只取 +V 和-V ( 或者说其逻辑电平值 只取 1 或 0)。)(ku待辩识系统U0u(k)Y(k)=Y0+y(k) (k)Y(k)=Y0+y(k) + (k) U0+ u(k)0 1 2 3 4 5 1 6 5 1 7 1 8
3、 1 9 1 10 1 11 1 12 1 13 1 14 1 15 1 )()()(0nkkkuuM 序列的当前值可以由过去的值进行异或运算得到。选取恰当的 k0，便可以产生最长周期为 Np=2n-1 的 M 序列。很容易由计算机程序来产生 M 序列。如取n=4，k 0=3， ，便可按1)3()2()1()(u)4()3()( kkkuu依次推出 ，最后得出图 2 所示的 M 序列。,)5(4M 序列的循环长度为 NP=2 -1=15,由此得到 N=4.NM=1 1 1 1 1 1 1 1 1 1 1 1 1 1 -1,寄存器初始状态为（1111） 。基本原理和基本关系式在相关法中，采用伪随
4、机信号 M 序列对系统模型进行在线辩识。由于 M 序列的自相关函数很接近于周期性的 函数，可近似地当作白噪声信号叠加在系统输入信号上，按其与输出信号间的互相关函数，可得到系统的脉冲响应曲线。设对象输入 u 为具有平稳遍历性的随机信号，其自相关函数为(2)(r u(t) E(t) (Ru当系统初始过渡过程结束后，输入和输出之间的互相关函数为(3)(r (t)y u(t)y ( uyuy式中 E表示统计平均， 表示时间平均。由维纳 霍甫方程有(4) pTuu dvdv00uy )(rh1)(rh)(r当 u 为伪随机信号 M 序列时，自相关函数的离散值为(5)式中,V、和 N p分别为 M 序列的
5、电平值、基本电平时间和离散周期长度， 为时差的离散值。相应的功率谱密度函数为(6)(2)sin()( 020 iSjiu 式中otherwisNVNlrl pppluu , ,2,|)()2 (7)210VSpN若将 M 序列的自相关函数近似为(8)ipuyiThr)()(0并取的周期大于对象的整定时间，即(9)spNT则有(10)puy ThSr0),()(0这样就可以得到脉冲响应曲线(11)puy,)(0批量算法和递推算法1） 批量算法（Batch algorithm）：(12)10 )(|)(pNkluyuy kylurl(13)()1() 2uyuyprVh2) 递推算法（Recurs
6、ive algorithm）：(14) 1)()1()1(1)1()1(0 )(2 mNmmpmNm hyuVhh pp 辩识对象运行情况的仿真取系统对象模型 中初值为)2()1()2()1()( dkdkkkk ubyay 。将 M 序列作为 ，即可求 。设工作点0)()1(y)(u,3,)(（U 0,Y0）分别为 0.1 和 0.2，则对象的实际输入、输出 U(k)和 Y(k)分别为)(0)(kkU(18)()()(yY其中正态白噪声序列 (k)根据中心极限定理采用 12 个在-0.5,0,5中均匀分布的随机数叠加的方法生成，其方差可由遭受幅度因子调节。在待辩识对象在大范围内是非线性的情况
7、下，通过辩识得出的是对象在工作点附近的局部线性化动态模型，因此还应去处工作点 Y0 ，以得出输出端动态分量的实测值。但工作点无法单独测量，只能根据最新得到的推算：)(43(kpYN个(19)10)()(0Nikk或递推计算如下：(20)()(1)(0)( NkNkkYY式中 Y0(k)的初值 Y0(0)可设为 0。于是，输出端动态分量的实测值可近似估算为(21)()()(koky实验步骤批量算法如图 3：初始化生成数据 u(k-d) , Y(k) , k=1-r-d,-1,0,Np-1k1 N p + r + d生成新的数据 )()(,kkYuk-k1iTp?打印结果k k +1图 3 批量算
8、法流程图10 )(|)()( pNkluyuy kylurl 0)2uyyprrVh 第一步：产生 M 序列，如图所示NP=15,N=4第二步：求系统输出，编程结果如图所示y(k)+a1y(k-1)+a2y(k-2)=b1u(k-1-d)+b2u(k-2-d) 。5.0,5.0,ba第三步：利用相关法求脉冲响应，如图所示 10 )(|)(pNkluyuy kylurl)()1() 2uyuyprVh递推算法如图 3：第一、 二步同前第三步：脉冲响应如图所示初始化生成数据 u(k-d) , Y(k) , k=1-r-d,-1,0,Np-1k1 N p + r + d生成新的数据 ， 构造 h,
9、)()(,kkYu向量，并送入 h1 向量，等待循环计算)(ku 1)()1()(1)()1(0 )(2 mNmpmNm hyuVhhpp k-k1iTp?打印结果k k +1 1)()1()1(1)1()1(0 )(2 mNmmpmNm hyuVhh pp 如下图所示，相关分析法很好地辨识了系统的脉冲响应附录程序：（1） 批量算法x1=1;x2=1;x3=1;x4=1;m=15;for i=1:my4=x4;y3=x3;y2=x2;y1=x1;x4=y3;x3=y2;x2=y1;x1=xor(y3,y4);if y4=0u(i)=-1;elseu(i)=y4;endendm=u%graphe
10、ri1=i;k=1:1:i1;subplot(3,1,1)plot(k,u,k,u,rx)xlabel(k)ylabel(m 序列)title(移位寄存器产生的 m 序列 )a1=-1;a2=0.5;b1=1;b2=0.5;Y(1)=0;Y(2)=0;for K=3:1:16Y(K)=Y(K-1)-0.5*Y(K-2)+u(K-1)+0.5*u(K-2);endsubplot(3,1,2)k=1:1:i1+1;plot(k,Y,k,Y,rx)xlabel(k)ylabel(输出)title(系统输出 )H=1/16;p=ones(15);q=eye(15);R=p+q;for j=1:15fo
11、r i=1:16-jM(i,i+j-1)=u(j);endendfor t=1:15Z(t,1)=Y(1,t);endg=H*R*M*Zsubplot(3,1,3)k=1:1:i1;plot(k,g,k,g,rx)xlabel(k)ylabel(输出)title(脉冲响应估计值)for r=1:15z(r)=m(r)*g(r);end%subplot(4,1,4)%k=1:1:i1;%plot(k,z,k,z,rx)%xlabel(k)%ylabel(输出)%title(辨识输出 )（2）递推算法x1=1;x2=1;x3=1;x4=1;m=15;for i=1:my4=x4;y3=x3;y2=
12、x2;y1=x1;x4=y3;x3=y2;x2=y1;x1=xor(y3,y4);if y4=0u(i)=-1;elseu(i)=y4;endendm=u%grapheri1=i;k=1:1:i1;subplot(3,1,1)plot(k,u,k,u,rx)xlabel(k)ylabel(m 序列)title(移位寄存器产生的 m 序列 )a1=-1;a2=0.5;b1=1;b2=0.5;Y(1)=0;Y(2)=0;for K=3:1:16Y(K)=Y(K-1)-0.5*Y(K-2)+u(K-1)+0.5*u(K-2);endsubplot(3,1,2)k=1:1:i1+1;plot(k,Y,
13、k,Y,rx)xlabel(k)ylabel(输出)title(系统输出 )p=ones(15);q=eye(15);R=p+q;for t=1:15Z(t,1)=Y(1,t);endfor j=1:15for i=1:16-jM(i,i+j-1)=u(j);endendg(1)=0;for i=2:15W=M; m=W(i,:);g(i)=i/(i+1)*g(i-1)+1/(i+1)*m*Z;endsubplot(3,1,3)k=1:1:i1;plot(k,g,k,g,rx)xlabel(k)ylabel(输出)title(脉冲响应估计值)for r=1:15z(r)=m(r)*g(r);end%subplot(4,1,4)%k=1:1:i1;%plot(k,z,k,z,rx)%xlabel(k)%ylabel(输出)%title(辨识输出 )
