程序注释及运行结果

1、第二部分程序注释及运行结果读者须知 :为了便于读者理解， 现将光盘上第一部分可直接在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分别表示 xi和 xi-1 ；x1=mod (x2,M);%将 x2存储器的数除以 M ，取余数放 x1（ xi）中；v1=x

2、1/256;%将 x1存储器的数除以 256得到小于 1的随机数放 v1中；v(:,k)=v1;% 将 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');tktl

3、e('(0-1)均匀分布的随机序列') 程序运行结果如图2.5所示。图2.5 采用 MA TLAB 产生的 (0,1)均匀分布的随机序列图 产生的 (0-1)均匀分布的随机序列在程序运行结束后，产生的 (0 ，1)均匀分布的随机序列， 直接从 MATLAB的 window界面中 copy 出来如下（ v2中每行存6 个随机数）：v2 =0.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.08200.4

4、9220.96090.78520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.08200.49220.96090.7

5、8520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.08200.49220.96090.78520.72660.37500.25780.55080.31640.90230.43360.60940.66800.02340.14060.84380.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

6、;%分别用 x2和x0表示 xi+1和 xi-1 ；x1=mod (x2,M);%取 x2存储器的数除以 M 的余数放 x1（ xi）中；v1=x1/256;%将 x1存储器中的数除以 256得到小于 1的随机数放 v1中；v(:,k)=(v1-0.5 )*f;%将 v1中的数（i）减去 0.5再乘以存储器 f中的系数，存放在矩阵存储器 v的第 k列中， v(:,k) 表示行不变、列随递推循环次数变化；x0=x1;% x i-1= x i；v0=v1;end%递推 100次结束；v2=v%该语句后无； ，实现矩阵存储器 v中随机数放在v2中，且可直接显示在 MA TLAB 的 window 中

7、；k1=k;%grapher%以下是绘图程序；k=1:k1;plot(k,v,k,v,'r');xlabel('k'), ylabel('v');tktle(' (-1,+1)均匀分布的白噪声') 程序运行结果如图2.6所示。图 2.6 采用 MATLAB 产生的 (-1,+1) 均匀分布的白噪声序列 产生的 (-1 ， 1)均匀分布的白噪声序列在程序运行结束后， 产生的 (-1 ，1) 均匀分布的白噪声序列， 直接从 MA TLAB的 window界面中 copy 出来如下（ v2中每行存6 个随机数）：v2 =-0.9531-

8、0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.

9、3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.

10、13280.21880.3359-0.9531-0.71880.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表示 xi+1和 xi-1 ；x1=mod (x2,M);%取 x2存储器的数除以 M 的余数放 x1（ xi）中；v1=x1/256;%将 x1存储器中的数除以 2

11、56得到小于 1的随机数放 v1中；v(:,k)=(v1-0.5 )*f;%将 v1中的数（i ）减去 0.5再乘以存储器 f中的系数，存放在矩阵存储器 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

12、9;);tktle(' (-1,+1)均匀分布的白噪声 ') 程序运行结果如图2.3所示。图2.3 白噪声过程第二章的 M 序列产生程序例2.3 用移位寄存器产生M 序列的MATLAB 软件实现 (见光盘FLch2bzsheg3.m) 编程如下：X1=1;X2=0;X3=1;X4=0;m=60;for i=1:m%1#Y4=X4;Y3=X3;X4=Y3;X3=Y2;% 移位寄存器输入% 置 M 序列总长度Y2=X2;Y1=X1;X2=Y1;Xi 初 T 态（ 0101 ），Yi 为移位寄存器各级输出X1=xor(Y3,Y4);% 异或运算if Y4=0U(i)=-1;elseU

13、(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-11-11111-1-1-1Columns 11 through 201-1-111-11-111Columns 21 through 3011-1-1-11-1-11

14、1Columns 31 through 40-11-11111-1-1-1Columns 41 through 501-1-111-11-111Columns 51 through 6011-1-1-11-1-111i1 =60第五章的递推的极大似然法辨识程序例 5.2 系统模型如图 5.5 所示。试用递推的极大似然法对系统辨识的参数集? v(k)1z 10.2z 21+ e(k)uu(k)z 10.5z 2y(k)+11.2z 10.6 z 2z(k)图 5.5 例 5.2 系统模型v(k) 随机信号，输入信号为幅值为1 的 M 序列或随机信号，要求画出程序流程图，打印出程序（程序中带有注释

15、）和辨识中的参数、误差曲线。解 : 首先解释编程所用的部分字母：由于在MA TLAB 语言中无法用希腊字母描述、无法用上标及下标，故?(k ) 用 o和 o1表示；令 P(k) Pi , (i 0,1) ； 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:1

16、200%产生 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计算 z z(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,

17、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)

18、=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(

19、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',&#

20、39;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(' 误差曲线 ')end

21、subplot(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 所示The parameter idendification of the RML2ra1=-1.2,eta2=0.6em0b1=1.0a

22、rb2=0.5apd1=-1.0-22004006008001000d2=0.2012001.5kr1orre0.500200400600800100012001ktup 0.5ni00200400600800100012005kesionm0odnar-5020040060080010001200k图 5.7 RML 辨识参数曲线第七章的用改进的神经网络MBP算法辨识例 7.1 对具有随机噪声的二阶系统的模型辨识（光盘上编号：FLch7NNeg1）对具有随机噪声的二阶系统的模型辨识，进行标幺化以后系统的参考模型差分方程为y(k)a1y(k - 1)a2 y(k - 2)bu(k- 1)(k)

23、(7.90)式中，a10.3366, a20.6634,b0.68,(k ) 为随机噪声。由于神经网络的输出最大为1，所以，被辨识的系统应先标幺化，识 )结构，神经网络选用3-9-9-1点个数分别为3、9、9、1 个；这里标幺化系数为 5。利用图型，即输入层 i ，隐层 j 包括7.5 正向建模2 级，输出层(并联辨k 的节由于神经网络的最大输出为 1，因此在辨识前应对原系统参考模型标么化处理，辨识结束后再乘以标么化系数才是被辨识系统的辨识结果。1) 编程如下 :%w10ij表示第一隐层权值w1ij ( k2) ， w11ij 表示 w1ij (k 1) ； w120ij 表示第二隐层权值 %

24、 w2ij(k2) ， w121ij表示 w21ij(k 1) ； w20j 表示输出层权值w3ij ( k 2) ， w21j表示 % w3ij(k1) ； q 表示隐层阈值； p 表示输出层阈值；置标幺化系数f1=5 等w10ij=.01 .01 .02; .1 .11 .02; .01 0 .1; .11 .01 .02;.1 .1 .02; .11 .1 .1;.1 .1 .1;0 .1 .1;.1 0 .1;w11ij=.1 .2 .11; .02 .13 .04; .09 .08 .08; .09 .1 .06; .1 .11 .02; .06 0 .1;.1 .1 .1;0 .1

25、 0;.1 .1 .1;w20j=.01;.02;.1;.2;.1;.1;.1;.1;.1;w21j=0; 0.1; .1; .02;0;.1;.1;.1;.1;q0j=.9 .8 .7 .6 .1 .2 .1 .1 .1;q120j=q0j;q11j=.5 .2 .3 .4 .1 .2 .1 .1 .1;q12j=q11j;w121ij=w20j*q0j;w120ij=w20j*q11j;f1=5;q2j=0;% threshold valuep0=.2;k1=1;p1=.3;w=0;xj=1 1 1;% inputserror=0.0001;a1=1 1 1 1;n=1;e1=0;e0=0

26、;e2=0;e3=0;e4=0;yo=0;ya=0;yb=0;y0=0;y1=0;y2=0;y3=0;u=0;u1=0;u2=0.68;u3=.780;u4=u3-u2;k1=1;kn=28;e3=.055; z1=0;z12=0; q123j=0; t2j=0; o12j=0;r=0;r1=0; s=0.1;d2j=0; %+% calculating output of the hidden layerv1=randn(40,1);for m=1:40s1=0.1*v1(m)yn=.3366*y2+.6634*u1+s*s1;y1=y2;y2=yn;yp=yn;u0=u1;u1=u2;ya

27、(m)=yn;fork=1:100% calculating output of the hidden layer(1) for i=1:9x1=w11ij(i,1)*xj(:,1)+w11ij(i,2)*xj(:,2)+w11ij(i,3)*xj(:,3);x=x1+q11j(:,i); o=1/1+exp(-x); o11j(i)=o;end% calculating output of the hidden layer(2)fori=1:9forj=1:9z1=z1+w121ij(i,j)*o11j(:,j);endz=z1+q12j(:,i);o=1/1+exp(-x);o12j(i)

28、=o;end% calculating output of the output layer for i=1:9 yb=yb+w21j(i,:)*o12j(:,i);end yi=yb+p1; y=1/1+exp(-yi);% calculating error value between aim and practice value e0=e1;e1=e2; e2=(yp-y).2/2; e(k)=e2;xj1=e2;xj2=e1;xj3=e0;xj=xj1 xj2 xj3;% revising right value (1)fori=1:9d1=o11j(:,i)*1-o11j(:,i)*

29、d2j*w21j(i,:);%计算第 1 隐层误差反传信号do=o11j(:,i)*d1;qw=q11j(:,i)-q0j(:,i);q2j=q11j(:,i)+.8*do+.4*qw;q3j(:,i)=q2j;for j=1:3dw=w11ij(i,j)-w10ij(i,j);w12ij=w11ij(i,j)+.8*do*xj(j)+.6*dw;w13ij(i,j)=w12ij;endendw10ij=w11ij;w11ij=w13ij;q0j=q11j;q11j=q3j;% revising right value (2)fori=1:9d1=o12j(:,i)*1-o12j(:,i)*d

30、2j*w21j(i,:);%计算第2 隐层误差反传信号do=o12j(:,i)*d1;qw=q12j(:,i)-q120j(:,i);t2j=q12j(:,i)+.8*do+.4*qw;q123j(:,i)=t2j;for j=1:9dw=w121ij(i,j)-w120ij(i,j);w122ij=w121ij(i,j)+.8*do*o11j(j)+.6*dw;w123ij(i,j)=w122ij;endendw120ij=w121ij;w121ij=w123ij;q120j=q12j;q12j=q123j;% revising right value (3)ifm<4,r=0.2;r

31、1=0.0001elser=0.14;r1=0.005;endfori=1:9;d2j=y*(1-y)*(yp-y);%计算输出误差反传信号dw=w21j(i,:)-w20j(i,:);w22j=w21j(i,:)+r*d2j*o12j(i)+.4*dw+r1*e2;w23j(i,:)=w22j;endw20j=w21j;w21j=w23j;ph=p1-p0;p2=p1+.96*(yp-y)+.58*ph+r1*e2;p0=p1;p1=p2;u=y;% k=k+1;ife2<=.005break;elseendendya(m)=yp*f1;e3(m)=e2;ym(m)=y*f1;v(m)

32、=s1;%m=m+1m6=mendw11ij=w13ijw121ij=w123ijw21j=w23jm1=m;% graphersubplot(3,1,1)m=1:m6;plot(m,ya,m,ym,'rx'),title('Identifiedmodelby MBPalgorithm'),xlabel('k'),ylabel('yaandym')legend('ya is system','ym is identified model'); %图标炷endsubplot(3,1,2)m=1:m6;

33、plot(m,e3),xlabel('k'),ylabel('error')endsubplot(3,1,3)m=1:m6;plot(m,v),xlabel('k'),ylabel('random noise')end2) 辨识结果Identified model by MBP algorithm4m 3ya is systemym is identified modelyd2naay100510152025303540x 10-3k3r2orre1005101520253035400.2ke0.1sionm0odna-0.1r-

34、0.25101520253035400k图 7.12MBP 算法对具有随机噪声的二阶系统辨识结果3) 各层权值w11ij =w21j =0.23520.48550.2456-0.0589-0.10090.15850.06900.10110.20930.19890.04940.04780.06070.23640.1214-0.15220.10050.12610.0211-0.1189-0.0157-0.15110.09940.04780.09930.09890.09940.0478-0.00070.0989-0.15060.04780.09930.24890.09940.0478w121ij =

35、0.01650.01820.01430.01040.00230.00330.00230.00230.00230.02480.02980.02130.0129-0.0026-0.0008-0.0026-0.0026-0.00270.14660.16720.12690.08660.00690.01680.00690.00700.00690.30430.34340.26380.18420.02380.04390.02380.02370.02380.15350.17280.13310.09340.01310.02310.01310.01300.01310.14660.16720.12690.08660

36、.00690.01680.00690.00700.00690.14660.16720.12690.08660.00690.01680.00690.00700.00690.14660.16720.12690.08660.00690.01680.00690.00700.00690.14660.16720.12690.08660.00690.01680.00690.00700.0069第七章的多维非线性辨识的MATLAB 程序例 7.2 用神经网络学习二维非线性函数ypcos(2 k1m2 k2 n)sin(2 k2n) （光盘上编号： FLch7NNeg2）% 初始化 :w10ij 表示 wij

37、(k-2) ；w11ij 表示 w ij(k-1) ；w20j 表示 w jk (k-2) ； w21j 表示 wjk (k-1) ； q 和 p 为阈值w10ij=.01 .01 .02; .01 .01 .02; .01 0 .01; .001 .001 .002;.001 0 .002; .0011 .001 .001 ; w11ij=-.1 -.02 .11; -.21 .10 -.19; -.14 .15 -.16; .14 -.13 .17; -.13 .12 .21; -.16 -.23 .13; w20j=.01;.02;.1;.2;.1;.1;w21j=.18;.9;.9;.

38、7;.8;.9;q0j=.5 .8 .4 .6 .1 .2 ;q1j=-.1 .02 .12 .14 -.02 .02;q2j=0;% threshold valuep0=.2;k1=1;p1=.1;w23j=0;0;0;0;0;0;w22j=0;w=0;xj=0.5 0.3 0.2;% inputsya=0 0 0; yp=0; yy=0; m1=0;yam=0;yp1=0; qw=0;yo=0 0 0; ya1=0;error=0.0001;n=1; q=0;e1=0;e0=0;e2=0;e3=0;e4=0;yo=0;ya=0;yb=0;y0=0;y1=0;y2=0;y3=0;u=0;u1

39、=0;u2=0;k1=1;w0=1;kn=28;dj2=0.01;e3=.0055; a1=0.036; a2=0.036; a3=0.08;w0=1;dj2=0.01;%+初始化结束for m=1:36for n=1:36q=pi*0.05*m;p=pi*0.02*n;yn=cos(2*q+2*p)*sin(2*p);yzs=0.1*yn;%xj=0.4+yzs 0.3 0.2; %协调器使网络的输入随样本动态变化 yp=yn;fork=1:6% calculating output of the hidden layerfor i=1:6 x1=w11ij(i,1)*xj(:,1)+w11

40、ij(i,2)*xj(:,2)+w11ij(i,3)*xj(:,3); x=x1+q1j(:,i);o=tanh(x);o1j(i)=o;end% calculating output of the output layer for i=1:3 yb=yb+w21j(i,:)*o1j(:,i);endyi=yb+p1;y=tanh(yi);% calculating error value between aim and practice value e0=e1;e1=e2; e2=(yp-y).2/2;% revising right valuefori=1:6d1=1-o1j(:,i)*d

41、j2*w23j(i,:);% 计算隐层误差反传信号 do=o1j(:,i)*d1;q3j(:,i)=q1j(i);for j=1:3dw=w11ij(i,j)-w10ij(i,j);if e2<0.004, a1=0;a2=0;a3=0; else, a1=0.022; a2=0.03;a3=0.003;end; w12ij=w11ij(i,j)+a1*do*xj(j)+a2*dw+a3*w;w13ij(i,j)=w12ij;endendw10ij=w11ij;w11ij=w13ij;q0j=q1j;q1j=q3j;w=yp-y;w0=w;w=0.36*w0+(yp-y);% 计算积累误

42、差wif e2<0.004, w=0.78*w; endfori=1:6d2j=y*(1-y)*(yp-y);%计算输出层误差反传信号dw=w21j(i,:)-w20j(i,:);w22j=w21j(i,:)+0.132*d2j*o1j(i)+0.26*dw;+0.026*w;w23j(i,:)=w22j;endw20j=w21j;w21j=w23j;ph=p1-p0;p2=p1+.3*(yp-y)+.102*ph;p0=p1;p1=p2;ife2<=0.005break; elseend% 判断误差end%对样本的一次采样值训练k 次结束ypp(m,n)=yn;%存储二维非线性样本yom(m,n)=y;%存储神经网络辨识结果e3(m,n)=e2;%存储训练误差end%变量 n 训练结束m2=mend%变量 m 训练结束w11ij=w13ij %句末无“；