程序注释及运行结果.docx
《程序注释及运行结果.docx》由会员分享，可在线阅读，更多相关《程序注释及运行结果.docx（53页珍藏版）》请在冰豆网上搜索。
程序注释及运行结果
第二部分程序注释及运行结果
读者须知:
为了便于读者理解，现将光盘上第一部分可直接在MATLAB6.I下运行的MATLAB程序的编号和书本上的内容对应如下，每个程序题目括号内的file.m是对应书本上的内容在光盘上第一部分的程序编号。
第二章的随机序列产生程序
例2.1用乘同余法产生随机数（见光盘FLch2sjxleg1.m）
①编程如下：
A=6;N=100;%初始化；
x0=1;M=255;
fork=1:
N%乘同余法递推100次；
x2=A*x0;%x2和x0分别表示xi和xi-1；
x1=mod（x2,M）;%将x2存储器的数除以M，取余数放x1（xi）中；
v1=x1/256;%将x1存储器的数除以256得到小于1的随机数放v1中；
v（:
k）=v1;%将v1中的数（
）存放在矩阵存储器v的第k列中，v（:
k）
%表示行不变、列随递推循环次数变化；
x0=x1;%xi-1=xi；
v0=v1;
end%递推100次结束；
v2=v%该语句末无‘；’，实现矩阵存储器v中随机数放在v2中，%且可直接显示在MATLAB的window中；
k1=k;
%grapher%以下是绘图程序；
k=1:
k1;
plot（k,v,k,v,'r'）;
xlabel（'k'）,ylabel（'v'）;tktle（'（0-1）均匀分布的随机序列'）
②程序运行结果如图2.5所示。
图2.5采用MATLAB产生的（0,1）均匀分布的随机序列图
③产生的（0-1）均匀分布的随机序列
在程序运行结束后，产生的（0，1）均匀分布的随机序列，直接从MATLAB的window界面中copy出来如下（v2中每行存6个随机数）：
v2=
0.02340.14060.84380.08200.49220.9609
0.78520.72660.37500.25780.55080.3164
0.90230.43360.60940.66800.02340.1406
0.84380.08200.49220.96090.78520.7266
0.37500.25780.55080.31640.90230.4336
0.60940.66800.02340.14060.84380.0820
0.49220.96090.78520.72660.37500.2578
0.55080.31640.90230.43360.60940.6680
0.02340.14060.84380.08200.49220.9609
0.78520.72660.37500.25780.55080.3164
0.90230.43360.60940.66800.02340.1406
0.84380.08200.49220.96090.78520.7266
0.37500.25780.55080.31640.90230.4336
0.60940.66800.02340.14060.84380.0820
0.49220.96090.78520.72660.37500.2578
0.55080.31640.90230.43360.60940.6680
0.02340.14060.84380.0820
第二章的白噪声产生程序
例2.2用乘同余法产生（见光盘FLch2bzsheg2.m）
1编程如下：
A=6;x0=1;M=255;f=2;N=100；%初始化；
x0=1;M=255;
fork=1:
N%乘同余法递推100次；
x2=A*x0;%分别用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中的数（
）减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v（:
k）表示行不变、列随递推循环次数变化；
x0=x1;%xi-1=xi；
v0=v1;
end%递推100次结束；
v2=v%该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MATLAB的window中；
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）均匀分布的白噪声序列，直接从MATLAB的window界面中copy出来如下（v2中每行存6个随机数）：
v2=
-0.9531-0.71880.6875-0.8359-0.01560.9219
0.57030.4531-0.2500-0.48440.1016-0.3672
0.8047-0.13280.21880.3359-0.9531-0.7188
0.6875-0.8359-0.01560.92190.57030.4531
-0.2500-0.48440.1016-0.36720.8047-0.1328
0.21880.3359-0.9531-0.71880.6875-0.8359
-0.01560.92190.57030.4531-0.2500-0.4844
0.1016-0.36720.8047-0.13280.21880.3359
-0.9531-0.71880.6875-0.8359-0.01560.9219
0.57030.4531-0.2500-0.48440.1016-0.3672
0.8047-0.13280.21880.3359-0.9531-0.7188
0.6875-0.8359-0.01560.92190.57030.4531
-0.2500-0.48440.1016-0.36720.8047-0.1328
0.21880.3359-0.9531-0.71880.6875-0.8359
-0.01560.92190.57030.4531-0.2500-0.4844
0.1016-0.36720.8047-0.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;
fork=1:
N%乘同余法递推100次；
x2=A*x0;%分别用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中的数（
）减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v（:
k）表示行不变、列随递推循环次数变化；
x0=x1;%xi-1=xi；
v0=v1;
end%递推100次结束；
v2=v%该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MATLAB的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序列总长度
fori=1:
m%1#
Y4=X4;Y3=X3;Y2=X2;Y1=X1;
X4=Y3;X3=Y2;X2=Y1;
X1=xor（Y3,Y4）;%异或运算
ifY4==0
U（i）=-1;
else
U（i）=Y4;
end
end
M=U
%绘图
i1=i
k=1:
1:
i1;
plot（k,U,k,U,'rx'）
xlabel（'k'）
ylabel（'M序列'）
title（'移位寄存器产生的M序列'）
②程序运行结果如图2.8所示。
图2.8软件实现的移位寄存器产生的M序列图
.
③'四级移位寄存器产生的M序列
M=
Columns1through10
-11-11111-1-1-1
Columns11through20
1-1-111-11-111
Columns21through30
11-1-1-11-1-111
Columns31through40
-11-11111-1-1-1
Columns41through50
1-1-111-11-111
Columns51through60
11-1-1-11-1-111
i1=
60
第五章的递推的极大似然法辨识程序
例5.2系统模型如图5.5所示。
试用递推的极大似然法对系统辨识的参数集
u
v（k）随机信号，输入信号为幅值为
的M序列或随机信号，要求画出程序流程图，打印出程序（程序中带有注释）和辨识中的参数、误差曲线。
解:
首先解释编程所用的部分字母：
由于在MATLAB语言中无法用希腊字母描述、无法用上标及下标，故
用‘o’和‘o1’表示；令
；
；产生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;%初始化
fori=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）;
end
u;%若取去‘；’可以在程序运行中观测到m序列
v=randn（1200,1）;%产生正态分布随机数
V=0;%计算噪声方差
fori=1:
1200
V=V+v（i）*v（i）;
end
V1=V/1200;
fork=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）;
end
o1=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;
%迭代计算参数值和误差值
fork=3:
1200
h=[-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）];
end
o1%若取去‘；’可以在程序运行中观测到参数
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（'TheparameteridendificationoftheRML'）;
end
subplot（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
subplot（4,1,3）
k=1:
1200;
plot（k,u）;
xlabel（'k'）
ylabel（'input'）
%title（'系统输入信号'）
end
subplot（4,1,4）
k=1:
1200;
plot（k,v）;
xlabel（'k'）
ylabel（'randomnoise'）
%title（'系统所加的随机噪声'）
end
2程序运行结果如图5.7所示
图5.7RML辨识参数曲线
第七章的用改进的神经网络MBP算法辨识
例7.1对具有随机噪声的二阶系统的模型辨识（光盘上编号：
FLch7NNeg1）
对具有随机噪声的二阶系统的模型辨识，进行标幺化以后系统的参考模型差分方程为
（7.90）
式中，
为随机噪声。
由于神经网络的输出最大为1，所以，被辨识的系统应先标幺化，这里标幺化系数为5。
利用图7.5正向建模（并联辨识）结构，神经网络选用3-9-9-1型，即输入层i，隐层j包括2级，输出层k的节点个数分别为3、9、9、1个；
由于神经网络的最大输出为1，因此在辨识前应对原系统参考模型标么化处理，辨识结束后再乘以标么化系数才是被辨识系统的辨识结果。
1）编程如下:
%w10ij表示第一隐层权值
，w11ij表示
；w120ij表示第二隐层权值%
，w121ij表示
；w20j表示输出层权值
，w21j表示%
；q表示隐层阈值；p表示输出层阈值；置标幺化系数f1=5等
w10ij=[.01.01.02;.1.11.02;.010.1;.11.01.02;.1.1.02;.11.1.1;.1.1.1;0.1.1;.10.1];
w11ij=[.1.2.11;.02.13.04;.09.08.08;.09.1.06;.1.11.02;.060.1;.1.1.1;0.10;.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;%thresholdvalue
p0=.2;
k1=1;
p1=.3;
w=0;
xj=[111];%inputs
error=0.0001;
a1=[1111];
n=1;
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=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;
%+++++++++++++++++++++++++++++++++%calculatingoutputofthehiddenlayer
v1=randn（40,1）;
form=1:
40
s1=0.1*v1（m）
yn=.3366*y2+.6634*u1+s*s1;
y1=y2;
y2=yn;
yp=yn;
u0=u1;
u1=u2;
ya（m）=yn;
fork=1:
100
%calculatingoutputofthehiddenlayer
（1）
fori=1:
9
x1=[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
%calculatingoutputofthehiddenlayer
（2）
fori=1:
9
forj=1:
9
z1=z1+w121ij（i,j）*o11j（:
j）;
end
z=z1+q12j（:
i）;
o=1/[1+exp（-x）];
o12j（i）=o;
end
%calculatingoutputoftheoutputlayer
fori=1:
9
yb=yb+w21j（i,:
）*o12j（:
i）;
end
yi=yb+p1;
y=1/[1+exp（-yi）];
%calculatingerrorvaluebetweenaimandpracticevalue
e0=e1;
e1=e2;
e2=[（yp-y）.^2]/2;
e（k）=e2;
xj1=e2;
xj2=e1;
xj3=e0;
xj=[xj1xj2xj3];
%revisingrightvalue
（1）
fori=1:
9
d1=o11j（:
i）*[1-o11j（:
i）]*d2j*w21j（i,:
）;%计算第1隐层误差反传信号
do=o11j（:
i）*d1;
qw=q11j（:
i）-q0j（:
i）;
q2j=q11j（:
i）+.8*do+.4*qw;
q3j（:
i）=q2j;
forj=1:
3
dw=w11ij（i,j）-w10ij（i,j）;
w12ij=w11ij（i,j）+.8*do*xj（j）+.6*dw;
w13ij（i,j）=w12ij;
end
end
w10ij=w11ij;
w11ij=w13ij;
q0j=q11j;
q11j=q3j;
%revisingrightvalue
（2）
fori=1:
9
d1=o12j（:
i）*[1-o12j（:
i）]*d2j*w21j（i,:
）;%计算第2隐层误差反传信号
do=o12j（:
i）*d1;
qw=q12j（:
i）-q120j（:
i）;
t2j=q12j（:
i）+.8*do+.4*qw;
q123j（:
i）=t2j;
forj=1:
9
dw=w121ij（i,j）-w120ij（i,j）;
w122ij=w121ij（i,j）+.8*do*o11j（j）+.6*dw;
w123ij（i,j）=w122ij;
end
end
w120ij=w121ij;
w121ij=w123ij;
q120j=q12j;
q12j=q123j;
%revisingrightvalue（3）
ifm<4,r=0.2;r1=0.0001;
else
r=0.14;r1=0.005;
end
fori=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;
end
w20j=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;
else
end
end
ya（m）=yp*f1;
e3（m）=e2;
ym（m）=y*f1;
v（m）=s1;
%m=m+1
m6=m
end
w11ij=w13ij
w121ij=w123ij
w21j=w23j
m1=m;
%grapher
subplot（3,1,1）
m=1:
m6;
plot（m,ya,m,ym,'rx'）,title（'IdentifiedmodelbyMBPalgorithm'）,xlabel（'k'）,ylabel（'yaandym'）
legend（'yaissystem','ymisidentifiedmodel'）;%图标炷
end
subplot（3,1,2）
m=1:
m
升级会员 