第二部分程序注释及运行结果读者须知:为了便于读者理解，现将光盘上第一部分可直接在MATLAB6.I 下运行的MATLAB程序的编号和书本上的内容对应如下，每个程序题目括号内的file.m是对应书本上的内容在光盘上第一部分的程序编号。

第二章的随机序列产生程序例2.1 用乘同余法产生随机数(见光盘FLch2sjxleg1.m)①编程如下：A=6; N=100; %初始化；x0=1; M=255;for k=1:N %乘同余法递推100次；x2=A*x0; %x2和x0分别表示x i和x i-1；x1=mod (x2,M); %将x2存储器的数除以M，取余数放x1（x i）中；v1=x1/256; %将x1存储器的数除以256得到小于1的随机数放v1中；）存放在矩阵存储器v的第k列中，v(:,k) v(:,k)=v1; % 将v1中的数（i%表示行不变、列随递推循环次数变化；x0=x1; %x i-1= x i；v0=v1;end %递推100次结束；v2=v %该语句末无‘；’，实现矩阵存储器v中随机数放在v2中，%且可直接显示在MA TLAB的window中；k1=k;%grapher %以下是绘图程序；k=1:k1;plot(k,v,k,v,'r');xlabel('k'), ylabel('v');tktle('(0-1)均匀分布的随机序列')②程序运行结果如图2.5所示。

图2.5 采用MA TLAB产生的(0,1)均匀分布的随机序列图③产生的(0-1)均匀分布的随机序列在程序运行结束后，产生的(0，1)均匀分布的随机序列，直接从MATLAB的window 界面中copy出来如下（v2中每行存6个随机数）：v2 =0.0234 0.1406 0.8438 0.0820 0.4922 0.96090.7852 0.7266 0.3750 0.2578 0.5508 0.31640.9023 0.4336 0.6094 0.6680 0.0234 0.14060.8438 0.0820 0.4922 0.9609 0.7852 0.72660.3750 0.2578 0.5508 0.3164 0.9023 0.43360.6094 0.6680 0.0234 0.1406 0.8438 0.08200.4922 0.9609 0.7852 0.7266 0.3750 0.25780.5508 0.3164 0.9023 0.4336 0.6094 0.66800.0234 0.1406 0.8438 0.0820 0.4922 0.96090.7852 0.7266 0.3750 0.2578 0.5508 0.31640.9023 0.4336 0.6094 0.6680 0.0234 0.14060.8438 0.0820 0.4922 0.9609 0.7852 0.72660.3750 0.2578 0.5508 0.3164 0.9023 0.43360.6094 0.6680 0.0234 0.1406 0.8438 0.08200.4922 0.9609 0.7852 0.7266 0.3750 0.25780.5508 0.3164 0.9023 0.4336 0.6094 0.66800.0234 0.1406 0.8438 0.0820第二章的白噪声产生程序例2.2 用乘同余法产生(见光盘FLch2bzsheg2.m)①编程如下：A=6; x0=1; M=255; f=2; N=100；%初始化；x0=1; M=255;for k=1: N %乘同余法递推100次；x2=A*x0; %分别用x2和x0表示x i+1和x i-1；x1=mod (x2,M); %取x2存储器的数除以M的余数放x1（x i）中；v1=x1/256; %将x1存储器中的数除以256得到小于1的随机数放v1中；）减去0.5再乘以存储器f中的系数，存放v(:,k)=(v1-0.5 )*f; %将v1中的数（i在矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环次数变化；x0=x1; % x i-1= x i；v0=v1;end %递推100次结束；v2=v %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MA TLAB的window中；k1=k;%grapher %以下是绘图程序；k=1:k1;plot(k,v,k,v,'r');xlabel('k'), ylabel('v');tktle(' (-1,+1)均匀分布的白噪声')②程序运行结果如图2.6所示。

图2.6 采用MA TLAB产生的(-1,+1)均匀分布的白噪声序列③产生的(-1，1)均匀分布的白噪声序列在程序运行结束后，产生的(-1，1)均匀分布的白噪声序列，直接从MA TLAB的window 界面中copy出来如下（v2中每行存6个随机数）：v2 =-0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.92190.5703 0.4531 -0.2500 -0.4844 0.1016 -0.36720.8047 -0.1328 0.2188 0.3359 -0.9531 -0.71880.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531-0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.13280.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359-0.0156 0.9219 0.5703 0.4531 -0.2500 -0.48440.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359-0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.92190.5703 0.4531 -0.2500 -0.4844 0.1016 -0.36720.8047 -0.1328 0.2188 0.3359 -0.9531 -0.71880.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531-0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.13280.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359-0.0156 0.9219 0.5703 0.4531 -0.2500 -0.48440.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359-0.9531 -0.7188 0.6875 -0.8359*另外，书中图2.3白噪声的产生如下：显然，只要在例2.2程序的初始化部分中给N=300，f=6，运行程序就可以得到如图2.3所示的（-3，3）的白噪声过程.①编程如下：A=6; x0=1; M=255; f=6; N=300；%初始化；x0=1; M=255;for k=1: N %乘同余法递推100次；x2=A*x0; %分别用x2和x0表示x i+1和x i-1；x1=mod (x2,M); %取x2存储器的数除以M的余数放x1（x i）中；v1=x1/256; %将x1存储器中的数除以256得到小于1的随机数放v1中；）减去0.5再乘以存储器f中的系数，存放v(:,k)=(v1-0.5 )*f; %将v1中的数（i在矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环次数变化；x0=x1; % x i-1= x i；v0=v1;end %递推100次结束；v2=v %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MA TLAB的window中；k1=k;%grapher %以下是绘图程序；k=1:k1;plot(k,v,k,v,'r');xlabel('k'), ylabel('v');tktle(' (-1,+1)均匀分布的白噪声')②程序运行结果如图2.3所示。

图2.3 白噪声过程第二章的M序列产生程序例2.3 用移位寄存器产生M序列的MATLAB软件实现(见光盘FLch2bzsheg3.m)①编程如下：X1=1;X2=0;X3=1;X4=0; %移位寄存器输入Xi初T态（0101），Yi为移位寄存器各级输出m=60; %置M序列总长度for i=1:m %1#Y4=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%绘图i1=ik=1:1:i1;plot(k,U,k,U,'rx')xlabel('k')ylabel('M序列')title('移位寄存器产生的M序列')②程序运行结果如图2.8 所示。

图2.8 软件实现的移位寄存器产生的M序列图.③'四级移位寄存器产生的M序列M =Columns 1 through 10-1 1 -1 1 1 1 1 -1 -1 -1 Columns 11 through 201 -1 -1 1 1 -1 1 -1 1 1 Columns 21 through 301 1 -1 -1 -1 1 -1 -1 1 1Columns 31 through 40-1 1 -1 1 1 1 1 -1 -1 -1 Columns 41 through 501 -1 -1 1 1 -1 1 -1 1 1 Columns 51 through 601 1 -1 -1 -1 1 -1 -1 1 1 i1 =60第五章的递推的极大似然法辨识程序例5.2系统模型如图5.5所示。

试用递推的极大似然法对系统辨识的参数集?=θv(k)随机信号，输入信号为幅值为的M 序列或随机信号，要求画出程序流程图，打印出程序（程序中带有注释）和辨识中的参数、误差曲线。

解: 首先解释编程所用的部分字母：由于在MA TLAB 语言中无法用希腊字母描述、无法用上标及下标，故)(ˆk θ用 ‘o ’和‘o1’表示；令)1,0(,)(==i Pi k P ；K k K =)(；产生M 序列时，a(i),b(i),c(i),d(i)表示四级移位寄存器的第1，2，3，4级寄存器的输出; ① 编程如下(光盘上该程序：FLch5RMLeg2.m ，可在MATLAB6.I 下直接运行):编程如下：clear %清零a(1)=1;b(1)=0;c(1)=1;d(1)=0;u(1)=d(1);z(1)=0;z(2)=0; %初始化for i=2:1200 %产生m 序列u(i)a(i)=xor(c(i-1),d(i-1));b(i)=a(i-1);c(i)=b(i-1);d(i)=c(i-1);u(i)=d(i);endu; %若取去‘；’可以在程序运行中观测到m 序列v=randn(1200,1); %产生正态分布随机数V=0; %计算噪声方差for i=1:1200V=V+v(i)*v(i);endV1=V/1200;for k=3:1200 %根据v和u计算zz(k)=1.2*z(k-1)-0.6*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2); endo1=0.001*ones(6,1);p0=eye(6,6); %赋初值zf(1)=0.1;zf(2)=0.1;vf(2)=0.1;vf(1)=0.1;uf(2)=0.1;uf(1)=0.1;%迭代计算参数值和误差值for k=3:1200h=[-z(k-1);-z(k-2);u(k-1);u(k-2);v(k-1);v(k-2)];hf=h;K=p0*hf*inv(hf'*p0*hf+1);p=[eye(6,6)-K*hf']*p0;v(k)=z(k)-h'*o1;o=o1+K*v(k) ;p0=p;o1=o;a1(k)=o(1);a2(k)=o(2);b1(k)=o(3);b2(k)=o(4);d1(k)=o(5);d2(k)=o(6);e1(k)=abs(a1(k)+1.2);e2(k)=abs(a2(k)-0.6);e3(k)=abs(b1(k)-1.0);e4(k)=abs(b2(k)-0.5);e5(k)=abs(d1(k)+1.0);e6(k)=abs(d2(k)-0.2);zf(k)=z(k)-d1(k)*zf(k-1)-d2(k)*zf(k-2);uf(k)=u(k)-d1(k)*uf(k-1)-d2(k)*uf(k-2);vf(k)=v(k)-d1(k)*vf(k-1)-d2(k)*vf(k-2);hf=[-zf(k-1);-zf(k-2);uf(k-1);uf(k-2);vf(k-1);vf(k-2)];endo1 %若取去‘；’可以在程序运行中观测到参数V1%绘图subplot(4,1,1)k=1:1200;plot(k,a1,'k:',k,a2,'b',k,b1,'r',k,b2,'m:',k,d1,'g',k,d2,'k');xlabel('k')ylabel('parameter')legend('a1=-1.2,','a2=0.6','b1=1.0','b2=0.5','d1=-1.0','d2=0.2'); %图标炷title('The parameter idendification of the RML');endsubplot(4,1,2)k=1:1200;plot(k,e1,'k',k,e2,'b',k,e3,'r',k,e4,'m',k,e5,'g',k,e6,'k');xlabel('k')ylabel('error')%title('误差曲线')endsubplot(4,1,3)k=1:1200;plot(k,u);xlabel('k')ylabel('input')%title('系统输入信号')endsubplot(4,1,4)k=1:1200;plot(k,v);xlabel('k')ylabel('random noise')%title('系统所加的随机噪声')end②程序运行结果如图5.7 所示p a r a m e t e rThe parameter idendification of the RMLe r r ori n p u tkr a n d o m n o i s e图5.7 RML 辨识参数曲线第七章的用改进的神经网络MBP 算法辨识例7.1 对具有随机噪声的二阶系统的模型辨识（光盘上编号：FLch7NNeg1）对具有随机噪声的二阶系统的模型辨识，进行标幺化以后系统的参考模型差分方程为)(1)-bu(k 2)-y(k 1)-y(k a )(21k a k y ξ+++= (7.90)式中，)(,68.0,6634.0,3366.0a 21k b a ξ===为随机噪声。