《系统辨识大作业(共26页)270.pdf》由会员分享，可在线阅读，更多相关《系统辨识大作业(共26页)270.pdf（27页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、一、问题(wnt)描述 考虑(kol)仿真(fn zhn)对象(duxing)：式中，u(k)和 z(k)是输入输出数据，v(k)是零均值、方差为 1 的不相关的随机噪声；u(k)采用与 e(k)不相关的随机序列。1.设计实验，产生输入输出数据；2.使用基本最小二乘法估计参数；3.考虑其他适用于有色噪声的辨识方法估计参数；4.模型验证。二、问题分析 对于单输入单输出系统（Single Input Single Output,SISO），如图 1所示，将待辨识的系统看作“灰箱”，它只考虑系统的输入输出特性，而不强调系统的内部机理。图 1 中，输入 u(k)和输出 z(k)是可以测量的，是系统模型
2、，用来描述系统的输入输出特性，y(k)是系统的实际输出。是噪声模型，v(k)是均值为零的不相关随机噪声，e(k)是有色噪声。图 1 SISO 系统的“灰箱”结构 对于 SISO 随机系统，被辨识模型为：其相应的差分方程为 若考虑被辨识系统或观测信息中含有噪声，被辨识模型(mxng)可改写为 式中，z(k)为系统(xtng)输出量的第 k 次观测(gunc)值；y(k)为系统(xtng)输出量的第 k 次真值，y(k-1)为系统输出量的第 k-1 次真值，以此类推；u(k)为系统的第 k 个输入值，u(k-1)为系统的第 k-1 个输入值；v(k)为均值为 0 的不相关随机噪声。1.数据生成 本
3、部分需要生成系统的输入输出数据以及噪声数据。1）白噪声的生成 辨识数据通常包含有噪声，如果该噪声相关性较弱或者强度很小，可近似看作白噪声。本次实验问题中，即服从标准正态分布，可以将噪声看作为服从正态分布的白噪声过程，在Matlab 中可以由 randn 函数生成。2）输入数据的生成 伪随机二进制序列（Pseudo Random Binary Sequence,PRBS）是广泛应用的一种伪随机序列，所谓“二进制”是指序列中每个随机变量只有“0”或“1”两种逻辑状态。伪随机二进制序列可由多级线性反馈移位寄存器组成的随机信号发生器产生，其中具有最长循环周期的线性移位寄存器序列是伪随机二进制序列最常见
4、的一种形式，简称 M 序列（Maximal Length Sequence）。M 序列由于具有近似白噪声的性质，而且工程上易于实现，能够保证较好的系统辨识精度，是普遍采用的一种辨识用输入信号。图 2 线性反馈移位寄存器产生伪随机二进制序列结构图 以一个(y)4 级线性反馈移位(y wi)寄存器产生伪随机二进制序列为例，如图 2 所示。假设(jish)4 个移位(y wi)寄存器输出的初态非全零，移位寄存器的工作原理是：一个移位脉冲来到后，每级移位寄存器的输出移到下一级移位寄存器作为输入，最末一级移位寄存器的输出即为伪随机二进制序列。3）输出数据的生成 根据给定的SISO系统，可以求出z(k)的
5、表达式：其理想系数值为.可以根据生成的白噪声序列和输入序列，以及必要的 0 初始值，带入表达式即可得到采样输出数据。2.差分模型阶检验 在实际场景中，辨识模型的阶数和纯时延往往是未知的，在很多情况下仅仅依靠猜测。在模型的阶数和纯时延不确定时，设系统模型为 其中 n 为模型的阶数，模型的阶估计可以采用多种方法，本实验采用比较简单易行的损失函数检验法。定义预报误差（噪声方差的估计值）的平方和为损失函数，即 当 n 从小增大时，应随之减小，当n 增大到某一值时，NJ应近似白噪声过程。采用以下的检验原则：在n-1 这一点，NJ最后一次出现陡峭的下降，此后就近似地保持不变或只有微小的下降，则取。3.参数
6、辨识模型 在系统辨识和参数估计领域中，最广泛的估计方法时最小二乘法和极大似然估计法。最小二乘法作为一种最基本的估计方法应用极为广泛，其他的大多数估计算法都与最小二乘法有关。它既可用于动态系统，也可用于静态系统；可用于线性系统，也可用于非线性系统；可用于离线估计，也可用于在线估计。在随机环境下利用(lyng)最小二乘法时，无须知道观测数据的概率统计信息，而这种方法获得的结果，却有相当好的统计性质。最小二乘参数估计方法来源于数理统计的回归分析，它能提供一个在最小方差意义上与实验数据最好的你和的模型(mxng)。该估计在一定条件下有最佳的统计特性，即它们是一致的，无偏的和有效的。最小二乘法时一个经典
7、的方法，概念简明，适应范围广，在一些情况下，可得到与极大似然法一样好的统计效果，它能很方便地与其他辨识算法建立关系。（1）递推最小二乘算法(sun f)当获得一批数据后，可一次求得相应的参数估计值，这样处理问题的方法称为一次完成算法或批处理算法。它在理论研究方面有许多方便之处，但当矩阵的维数增加时，矩阵求逆运算的计算量会急剧增加，将给计算机的计算速度和存储量带来负担，而且不适合在线辨识，无法跟踪参数随时间变化的情况。为了减少计算量，减少数据在计算机中所占的存储量，也为了实时地辨识出动态系统的特性，在用最小二乘法(chngf)进行参数估计时，把它转化成参数递推的估计。参数递推估计是指被辨识的系统
8、，每取得一次新的测量数据后，就在前一次估计结果的基础上，利用新引入的测量数据对前一次估计的结果进行修正，从而递推地得出新的参数估计值。这样，随着新测量数据的引入，一次接一次地进行参数估计，直到估计值达到满意的精确程度为止。最小二乘递推算法的基本思想可以概括为：当前的估计值=上次估计值+修正项 即新的估计值()k是在旧的估计值(1)k的基础上，利用新的观测数据对旧的估计值进行修正而成的。考 虑如下(r xi)模型：，其中(q zh ng)，分别是系统(xtng)的输入和输出；是均值(jn zh)为零，方差为一的不相关白噪声。且满足：令 则使用加权最小二乘参数估计递推算法（Recursive We
9、ighted Least Squares,RWLS）有：图 3 最小二乘递推辨识参数估计过程中信息的变换 可以看出，取的时候，加权最小二乘估计就退化成了最小二乘参数估计的递推算法（Recursive Least Squares,RLS）。加权参数可以在范围内选择，如果选，所有的采样数据都是等同加权的，如果,则表示对新近获得的数据给予充分大的加权因子，而削弱历史观测数据的影响。实际计算时，需要首先确定初始参数和。一般说来选取(xunq)对于这样的系统，使用最小二乘法参数估计的递推算法进行(jnxng)辨识可以得到无偏估计，但是如果噪声模型必须用表示时，此时就无法得到无偏估计了，因为该方法没有把噪
10、声模型(mxng)考虑进去。使用递推最小二乘法进行系统(xtng)辨识的过程如图 4 所示。（2）增广最小二乘法 当噪声均值为 0 时，最小二乘参数估计算法为无偏估计；当噪声的均值不为 0 时，最小二乘参数估计算法为有偏估计。为了解决最小二乘参数估计的有偏性，将噪声模型的辨识同时考虑进去，因此称为增广最小二乘法。该算法可以看成是对一般最小二乘参数辨识算法的简单推广或扩充，因此又称为扩充最小二乘算法。考虑如下模型：，其中()u k，()z k 分别是系统的输入和输出；()v k是均值为零，方差为一的不相关白噪声。且满足：在模型阶次已经确定的情况下，令 由于()v k是白噪声，所以利用增广最小二乘
11、法这一形式上的变换，即可获得参数的无偏估计。不过由于数据向量中包含着不可测的噪声量，因此要用相应的估计值代替。当；当，即可得递推公式如下：一般说来选取(xunq)初值 增广(zn un)最小二乘递推算法扩充了最小二乘法的参数向量和数据(shj)向量()h k的维数，把噪声模型的辨识同时考虑进去。最小二乘法的许多(xdu)结论对它都是适用的，但最小二乘法只能获得模型的参数估计，对于有色噪声，也就是噪声模型必须用1()()C zv k表示时，只能采用增广最小二乘递推算法，方可获得无偏估计。这是最小二乘参数估计的递推算法不可代替的。增广最小二乘算法的流程如图 5 所示。图 4 递推最小二乘参数辨识流
12、程图 图 5 增广最小二乘算法流程图 三、Matlab 仿真(fn zhn)及运行结果 在本次(bn c)实验中，取值为 1，则系统(xtng)辨识对象为:1.差分(ch fn)模型的阶检验 在实际应用场景中，模型的阶数和纯时延往往是未知的，在此种情况下，一般采用如下方法来确定模型的阶数：固定一个阶数n，当 n 从小增大时，定义损失函数NJ应随之减小，当 n 增大到某一值时，NJ应近似白噪声过程。在n-1 这一点，NJ最后一次出现陡峭的下降，此后就近似地保持不变或只有微小的下降，则取 nn。使用 Matlab 进行模型阶数确定代码如下：%噪声方差的估计值和损失函数(hnsh)法定阶%clear
13、 close all;L=100;y1=1;y2=1;y3=1;y4=0;for i=1:L+5;x1=xor(y3,y4);%第一个移位(y wi)寄存器的输入信号 x2=y1;%第二个移位(y wi)寄存器的输入信号 x3=y2;%第三个移位寄存器的输入(shr)信号 x4=y3;%第四个移位寄存器的输入信号 y(i)=y4;%第四个移位寄存器的输出，即 M 序列，幅值为“0”或“1”if y(i)0.5;u(i)=-1;%M 序列的值为 1 时，辨识的输入取-1 else u(i)=1;%M 序列为 0 时，输入取 1 end y1=x1;y2=x2;y3=x3;y4=x4;%移位寄存器
14、的输出序列 end v=randn(L+5,1);%-ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)-%ksi(2)=0;ksi(1)=0;%设 ksi 的前两个初始值为 0 lamda=1;for k=3:L+5;ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2);end z(3)=0;z(2)=0;z(1)=0;%取 z 的前三个初始值为零 for k=4:L+5;z(k)=-0.9*z(k-1)-0.15*z(k-2)-0.02*z(k-3)+0.7*u(k-1)-1.5*u(k-2)+ksi(k);%理想辨识输出(shch
15、)采样信号 end%模型(mxng)阶次 n=1 for i=1:L H1(i,1)=z(i);H1(i,2)=u(i);end estimate1=inv(H1*H1)*H1*z(2:L+1);D1=(z(2:L+1)-H1*estimate1)*(z(2:L+1)-H1*estimate1)/L;%噪声(zoshng)方差的估计值 AIC1=L*log(D1)+4*1;%模型(mxng)阶次 n=2 for i=1:L H2(i,1)=z(i+1);H2(i,2)=z(i);H2(i,3)=u(i+1);H2(i,4)=u(i);end estimate2=inv(H2*H2)*H2*z(
16、3:L+2);D2=(z(3:L+2)-H2*estimate2)*(z(3:L+2)-H2*estimate2)/L;%噪声方差的估计值 AIC2=L*log(D2)+4*2;%模型阶次 n=3 for i=1:L H3(i,1)=z(i+2);H3(i,2)=z(i+1);H3(i,3)=z(i);H3(i,4)=u(i+2);H3(i,5)=u(i+1);H3(i,6)=u(i);end estimate3=inv(H3*H3)*H3*z(4:L+3);D3=(z(4:L+3)-H3*estimate3)*(z(4:L+3)-H3*estimate3)/L;%噪声(zoshng)方差的估
17、计值 AIC3=L*log(D3)+4*3;%模型(mxng)阶次 n=4 for i=1:L H4(i,1)=z(i+3);H4(i,2)=z(i+2);H4(i,3)=z(i+1);H4(i,4)=z(i);H4(i,5)=u(i+3);H4(i,6)=u(i+2);H4(i,7)=u(i+1);H4(i,8)=u(i);end estimate4=inv(H4*H4)*H4*z(5:L+4);D4=(z(5:L+4)-H4*estimate4)*(z(5:L+4)-H4*estimate4)/L;%噪声(zoshng)方差的估计值 AIC4=L*log(D4)+4*4;%模型(mxng)
18、阶次 n=5 for i=1:L H5(i,1)=z(i+4);H5(i,2)=z(i+3);H5(i,3)=z(i+2);H5(i,4)=z(i+1);H5(i,5)=z(i);H5(i,6)=u(i+4);H5(i,7)=u(i+3);H5(i,8)=u(i+2);H5(i,9)=u(i+1);H5(i,10)=u(i);end estimate5=inv(H5*H5)*H5*z(6:L+5);D5=(z(6:L+5)-H5*estimate5)*(z(6:L+5)-H5*estimate5)/L;%噪声(zoshng)方差的估计值 AIC5=L*log(D5)+4*5;plot(1:5,
19、AIC1 AIC2 AIC3 AIC4 AIC5)grid on title(损失(snsh)函数法判断模型阶次)xlabel(阶次)ylabel(损失(snsh)函数)程序运行结果(ji gu)如图 6 所示：图 6 模型阶次的判定 由图 6 可知，当系统阶次从 1 到 5 变化时，在 n=3 及以后，损失函数只有微小的下降，在 n=4 以后又出现上升的趋势，所以取 n=3 作为模型的阶次。2.参数辨识（1）递推最小二乘法 根据前文所述递推最小二乘法算法流程，编写递推最小二乘法程序如下：%递推最小二乘辨识(bin sh)clc clear all%采用 M 序列(xli)方法产生输入数据 L
20、=90;%90 个数据(shj)namenda=1;y1=1;y2=1;y3=1;y4=0;u_f=2;for i=1:L x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;if y(i)0.5 u(i)=-u_f;else u(i)=u_f;end y1=x1;y2=x2;y3=x3;y4=x4;end v=1*randn(L,1);eb(2)=0;eb(1)=0;for k=3:L eb(k)=-eb(k-1)-0.41*eb(k-2)+namenda*v(k);end z(3)=0;z(2)=0;z(1)=0;for k=4:L;z(k)=-0.9*z(k-
21、1)-0.15*z(k-2)-0.02*z(k-3)+0.7*u(k-1)-1.5*u(k-2)+eb(k);end figure(1)i=1:30;subplot(2,1,1)stem(u(i);title(输入(shr)M 序列)grid on subplot(2,1,2)stem(z(i);title(输出(shch)序列)grid on;%c0=0.001 0.001 0.001 0.001 0.001;p0=103*eye(5,5);E=0.0000005;c=c0,zeros(5,L-1);e=zeros(5,L);%for k=4:L h1=-z(k-1),-z(k-2),-z(
22、k-3),u(k-1),u(k-2);k1=p0*h1*inv(h1*p0*h1+1);c1=c0+k1*(z(k)-h1*c0);e1=c1-c0;e2=e1./c0;e(:,k)=e2;c0=c1;c(:,k)=c1;p1(:,:,k)=p0-k1*k1*h1*p0*h1+1;p0=p1(:,:,k);if abs(e2)0.5;u(i)=-1;%M 序列的值为 1 时，辨识的输入为-1 else u(i)=1;%M 序列的值为 0 时，辨识的输入为 1 end y1=x1;y2=x2;y3=x3;y4=x4;%移位寄存器的输出序列,为下一次的输入信号做准备 end figure(1);s
23、ubplot(3,1,1);stem(u),grid on%画出 M 序列输入信号径线图并给图形加上网格 title(输入信号（M 序列）)%-随机噪声的产生-%v=randn(1,L);%产生一组 L=60 个正态分布的噪声 subplot(3,1,2);plot(v),grid on%画出随机噪声信号 title(随机噪声信号 e(k)%-实际(shj)噪声的产生-%-ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)-%ksi(2)=0;ksi(1)=0;%设 ksi 的前两个(lin)初始值为0 lamda=1;for k=3:L;ksi(k)=lamda
24、*v(k)-ksi(k-1)-0.41*ksi(k-2);end subplot(3,1,3)plot(ksi),grid on title(实际(shj)噪声信号xi(k)%-辨识(bin sh)过程-%z(3)=0;z(2)=0;z(1)=0;zm(3)=0;zm(2)=0;zm(1)=0;%输出采样，不考虑噪声时系统输出，不考虑噪声时模型输出，模型输出的初值 c0=0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001;%直接给出辨识被辨识参数的初始值，即一个充分小的实向量 p0=106*eye(8,8);%直接给出初始向量 P0，即一个充分大的实
25、数单位矩阵 c=c0,zeros(8,L-1);%被辨识参数参数矩阵的初始值及大小 e=zeros(8,L);%相对误差的初始值及大小 for k=4:L;%开始求 K z(k)=-0.9*z(k-1)-0.15*z(k-2)-0.02*z(k-3)+0.7*u(k-1)-1.5*u(k-2)+ksi(k);%系统在 M 序列下的输出采样信号 h1=-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),v(k),ksi(k-1),ksi(k-2);%为求 K（k）做准备 x=h1*p0*h1+1;x1=inv(x);k1=p0*h1*x1;%K d1=z(k)-h1*c0
26、;c1=c0+k1*d1;%辨识参数 c e1=c1-c0;e2=e1./c0;%求参数误差的相对变化 e(:,k)=e2;c0=c1;%给下一次用 c(:,k)=c1;%把递推出的辨识参数 c 的列向量加入辨识参数矩阵 p1=p0-k1*k1*h1*p0*h1+1;%find p(k)p0=p1;%给下次(xi c)的计算使用 end%整个(zhngg)循环结束 c,e,%显示被辨识参数(cnsh)及参数收敛情况%分离(fnl)变量 a1=c(1,:);a2=c(2,:);a3=c(3,:);b1=c(4,:);b2=c(5,:);%分离出 a1,a2,a3,b1,b2 d1=c(6,:);
27、d2=c(7,:);d3=c(8,:);%分离出的 d1,d2,d3%分离变量的收敛情况 ea1=e(1,:);ea2=e(2,:);ea3=e(3,:);eb1=e(4,:);eb2=e(5,:);%分离出 a1,a2,b1,b2 的收敛情况 ed1=e(6,:);ed2=e(7,:);ed3=e(8,:);%分离出 d1,d2,d3 的收敛情况 figure(2);i=1:L;plot(i,a1,r,i,a2,r:,i,a3,r+,i,b1,b,i,b2,b:,i,d1,g,i,d2,g:,i,d3,g+)grid on legend(a1,a2,a3,b1,b2,c1,c2,c3)tit
28、le(辨识曲线)figure(3);i=1:L;plot(i,ea1,r,i,ea2,r:,i,ea3,r+,i,eb1,b,i,eb2,b:,i,ed1,g,i,ed2,g:,i,ed3,g+)%画出各个参数的收敛情况 grid on legend(ea1,ea2,ea3,eb1,eb2,ec1,ec2,ec3)title(辨识误差曲线)figure(4);i=1:L;plot(i,z(i),g),grid on;title(辨识系统的输出响应)%画出被辨识系统的输出响应 程序运行结果如图 10 至图 13 所示。图 10 输入数据(shj)及噪声 图 11 辨识(bin sh)曲线 图 1
29、2辨识误差(wch)曲线 图 13辨识的输出(shch)响应 使用增广最小二乘算法进行系统辨识(bin sh)的结果如表 2 所示。表 2 增广(zn un)最小二乘算法的辨识(bin sh)结果 参数 真值 0.9 0.15 0.02 0.7-1.5 1.0-1.0-0.4 估计值 0.9 0.15 0.02 0.7-1.5 1.0-1.0-0.4 从上表中可以看出(kn ch)，对于含有有色噪声的系统，采用增广最小二乘法进行系统辨识可以取得无偏估计。四、模型与误差分析 对系统进行辨识后，需要对辨识所得的系统进行进一步的模型验证，以观察所得的模型是否是对当前观测数据的最佳选择。对模型进行验证
30、，可以采用事先生成的输入数据，分别输入到已知模型和经过系统辨识所得的模型中，检验所得模型是否存在问题。如果检验失败，可能存在的问题是：辨识所用的一组数据包含的信息部足或所选的模型类不合适。对本实验中的已知对象及经过系统辨识所得的模型，选取一组输入信号为：比较原参考对象的输出以及辨识模型的输出。具体代码如下：%模型检验%clear all close all u=1,-1,1,-1,1,1,-1,1,1,-1,1,1,-1,-1,-1,1;%系统输入信号 L=20 z=zeros(1,L);x=zeros(1,L);y=zeros(1,L);v=randn(1,L);%-实际噪声的产生-%-ks
31、i(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)-%ksi(2)=0;ksi(1)=0;%设 ksi 的前两个初始值为0 lamda=1;for k=3:L;ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2);end for k=4:16 z(k)=-0.9*z(k-1)-0.15*z(k-2)-0.02*z(k-3)+0.7*u(k-1)-1.5*u(k-2)+ksi(k)%真实(zhnsh)模型的输出 y(k)=-1.2317*y(k-1)-0.4735*y(k-2)-0.0600*y(k-3)+0.8534*u(k-1)-1.23
32、73*u(k-2)%递推最小二乘辨识模型(mxng)输出 x(k)=-0.9*x(k-1)-0.15*x(k-2)-0.02*x(k-3)+0.7*u(k-1)-1.5*u(k-2)+ksi(k)%增广最小二乘辨识模型(mxng)的输出 end figure(1);subplot(3,1,1)stem(u),grid on title(系统(xtng)输入信号)subplot(3,1,2)i=1:1:L;plot(i,z,r,i,y,b),grid on legend(真实模型输出,递推最小二乘辨识模型)title(真实模型输出和递推最小二乘辨识模型输出值)subplot(3,1,3)i=1:
33、1:L;plot(i,z,r,i,x,g),grid on legend(真实模型输出,增广最小二乘辨识模型)title(真实模型输出和增广最小二乘辨识模型输出值)figure(2);subplot(2,1,1)i=1:1:L;plot(i,z-y),grid on;title(真实模型与递推最小二乘辨识模型的差值)subplot(2,1,2)i=1:1:L plot(i,z-x),grid on;title(真实模型与增广最小二乘辨识模型的差值)程序运行结果如图 14 和图 15 所示。图 14 不同辨识(bin sh)模型输出和真实系统输出对比 图 15 真实模型与辨识(bin sh)模型
34、输出信号的差值 由上图并结合(jih)表 1 和表 2，分析可得 表 1 递推最小二乘算法(sun f)的辨识结果 参数 真值 0.9 0.15 0.02 0.7-1.5 估计值 0.9367 0.1591-0.0138 0.6275-1.3558 表 2 增广(zn un)最小二乘算法的辨识(bin sh)结果 参数 真值 0.9 0.15 0.02 0.7-1.5 1.0-1.0-0.4 估计值 0.9 0.15 0.02 0.7-1.5 1.0-1.0-0.4 辨识给定的含有色噪声的系统时，最小二乘递推辨识算法的参数辨识结果是有偏差的，无法达到稳定(wndng)值,但是这个误差是在-0.
35、5,0.5范围内，还是可以接受的。而与最小二乘递推辨识算法相比，增广最小二乘递推辨识算法考虑了噪声模型，并且在递推到第九步时辨识结果就达到了稳定值，具有辨识速度快、辨识结果精确的特点。因此，对于含有有色噪声的系统，使用增广二乘辨识算法能够取得较好的效果。五、实验总结 通过本次仿真实验，使我加深了对系统辨识算法的认识和了解，对于最小二乘参数辨识这种系统辨识领域中最基本最常用的方法有了较好的掌握，并通过 MATLAB 仿真熟悉了系统辨识的一般内容和步骤，对系统辨识的基本过程有了一个较为清晰的认识。最后，王立琦老师这一学期给我们所带的系统辨识这门课程使我受益良多，对我的科研工作有很大的启发借鉴意义。在此十分感谢王老师的指导和帮助。内容摘要 （1）问题描述 考虑仿真对象：式中，u(k)和 z(k)是输入输出数据，v(k)是零均值、方差为 1 的不相关的随机噪声（2）H5(i,6)=u(i+4)（3）H5(i,7)=u(i+3)（4）p0=103*eye(5,5)（5）c1=c0+k1*(z(k)-h1*c0)（6）%移位寄存器产生的 M 序列周期（取 L=60）y1=1（7）%系统输入信号 L=20 z=zeros(1,L)（8）在此十分感谢王老师的指导和帮助