1、白噪声白噪声 如果一个零均值、平稳随机过程的谱密度为常数,我们称之为白噪声(由白色光联想而得)。白噪声有以下特点: 1 2 ,频谱宽度无限。,=00,0 3 ,其中,为Dirac函数,即= 且 4 无记忆性,即t时刻的数值与t时刻以前的过去值无关,也不影响t时刻以后的将来值。从另一意义上说,即不同时刻的随机信号互不相关。 白噪声的用途: 1 作为系统输入时,有,=0,1,2,,即为系统的单位脉冲响应。 2 作为被辨识系统输入时,可以激发系统的所有模态,可对系统充分激励; 3 作为被辨识系统输入时,可防止数据病态,保证辨识精度。 4 在辨识过程中,以输出估计误差是否具有白色性来判断辨识方法的优劣
2、,也可用来判断模型的结构和参数是否合适。 5 产生有色噪声。 白噪声的产生方法: 1 (0,1)均匀分布白噪声:,i=1,2,3,初值可取为: 2 正态分布白噪声:, 其中为服从(0,1)均匀分布的白噪声。例2.1 用乘同余法产生随机数选A=6,k=8,M=28=256,递推100次.用MATLAB编程编程如下: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=x1/256; %将x1存储器的数除以
3、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);title(0-1)均匀分布的随机序列) 程序运行结果如图所示 产生的(0-1)均匀分布的随机序列在程序运行结束后,产生的(
4、0,1)均匀分布的随机序列,直接从MATLAB的window界面中copy出来如下(v2中每行存6个随机数): v2 = 0.0234 0.1406 0.8438 0.0820 0.4922 0.9609 0.7852 0.7266 0.3750 0.2578 0.5508 0.3164 0.9023 0.4336 0.6094 0.6680 0.0234 0.1406 0.8438 0.0820 0.4922 0.9609 0.7852 0.7266 0.3750 0.2578 0.5508 0.3164 0.9023 0.4336 0.6094 0.6680 0.0234 0.1406 0
5、.8438 0.0820 0.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.9609 0.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.7266 0.3750 0.2578 0.5508 0.3164 0.9023 0.43
6、36 0.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。1编程如下:A=6; x0=1; M=255; f=2; N=100; %初始化;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)中
7、; 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);title( (-1,+1
8、)均匀分布的白噪声) 程序运行结果如图2.6所示。图2.6 采用MATLAB产生的(-1,+1)均匀分布的白噪声序列 产生的(-1,1)均匀分布的白噪声序列在程序运行结束后,产生的(-1,1)均匀分布的白噪声序列,直接从MATLAB的window界面中copy出来如下(v2中每行存6个随机数):v2 = -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8
9、359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3
10、672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8
11、359例2.3 用移位寄存器产生M序列的MATLAB软件实现 编程如下: 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=0 U(i)=-1; else U(i)=Y4;end endM=U%绘图i1=ik=1:1:i1;plot(k,U,k,U,rx)xlabel(k)ylabel(M序列)title(移位寄存器产生的M序列)
12、程序运行结果如图2.8 所示。图2.8 软件实现的移位寄存器产生的M序列图. 四级移位寄存器产生的M序列M = Columns 1 through 10 -1 1 -1 1 1 1 1 -1 -1 -1 Columns 11 through 20 1 -1 -1 1 1 -1 1 -1 1 1 Columns 21 through 30 1 1 -1 -1 -1 1 -1 -1 1 1 Columns 31 through 40 -1 1 -1 1 1 1 1 -1 -1 -1 Columns 41 through 50 1 -1 -1 1 1 -1 1 -1 1 1 Columns 51 t
13、hrough 60 1 1 -1 -1 -1 1 -1 -1 1 1i1 = 60最小二乘一次完成算法的MATLAB仿真考虑仿真对象 式中是服从正态分布的白噪声N(0,1),输入信号采用4阶M序列,幅度为1。选择如下形式的辨识模型 构造 和;数据长度取L=14;加权阵取;利用计算参数估计值。 设输入信号的取值是从k=1到k=16的M序列。 待辨识参数,观测矩阵和的表达式为: ,仿真程序流程图: %二阶系统的最小二乘一次完成算法辨识程序。u=-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1; %系统辨识的输入信号为一个周期的M序列z=zeros(1,16); %定义输出观
14、测值的长度for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值endsubplot(3,1,1) %画三行一列图形窗口中的第一个图形stem(u) %画输入信号u的径线图形subplot(3,1,2) %画三行一列图形窗口中的第二个图形i=1:1:16; %横坐标范围是1到16,步长为1plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线subplot(3,1,3) %画三行一列图形窗口中的第三个图形stem(z),grid on %画出输出观测值z的径线图形,并显示坐
15、标网格u,z %显示输入信号和输出观测信号%L=14 %数据长度HL=-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12)
16、 u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14) %给样本矩阵HL赋值ZL=z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16) % 给样本矩阵z L赋值%Calculating Parameters c1=HL*HL; c2=inv(c1); c3=HL*ZL; c=c2*c3 %计算并显示%Display Parameters a1=c(1), a2=c(2), b1=c(3),b2=c(4) %从中分离出并显示
17、a1 、a2、 b1、 b2%End程序运行结果: u = -1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1z = 0,0,0.5000,0.2500,0.5250,2.1125, 4.3012,6.4731,6.1988,3.2670,-0.9386, -3.1949,-4.6352,6.2165,-5.5800,-2.5185HL = 0 0 1.0000 -1.0000 -0.5000 0 -1.0000 1.0000 -0.2500 -0.5000 1.0000 -1.0000 -0.5250 -0.2500 1.0000 1.0000 -2.1125 -0.5
18、250 1.0000 1.0000 -4.3012 -2.1125 1.0000 1.0000 -6.4731 -4.3012 -1.0000 1.0000 -6.1988 -6.4731 -1.0000 -1.0000 -3.2670 -6.1988 -1.0000 -1.0000 0.9386 -3.2670 1.0000 -1.0000 3.1949 0.9386 -1.0000 1.0000 4.6352 3.1949 -1.0000 -1.0000 6.2165 4.6352 1.0000 -1.0000 5.5800 6.2165 1.0000 1.0000ZL = 0.5000,
19、0.2500,0.5250,2.1125,4.3012,6.4731,6.1988,3.2670,-0.9386,-3.1949, -4.6352,-6.2165,-5.5800,-2.5185Tc = -1.5000,0.7000,1.0000,0.5000Ta1 = -1.5000a2 = 0.7000b1 = 1.0000b2 =0.5000 从仿真结果表可以看出,由于所用的输出观测值没有任何噪声成分,所以辨识结果也无任何误差。考虑图3.所示的仿真对象,图中, 是服从N分布的不相关随机噪声。且 , ,选择图所示的辨识模型。仿真对象选择如下的模型结构 其中,是服从正态分布的白噪声N。输入信
20、号采用4位移位寄存器产生的M序列,幅度为0.03。按式 构造;加权阵取单位阵;利用式递推公式计算l(k)、和P(k),计算各次参数辨识的相对误差,精度满足要求式适当小的数(当所有参数估计值变化不大时) 后停机。下面给出具体程序。%最小二乘递推算法辨识程序, clear %清理工作间变量L=15; % M序列的周期y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值for i=1:L;%开始循环,长度为L x1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或” x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出 x3=y2;
21、%第三个移位寄存器的输入是第二个移位寄存器的输出 x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出 y(i)=y4; %取出第四个移位寄存器的幅值为0和1的输出信号,即M序列 if y(i)0.5,u(i)=-0.03; %如果M序列的值为1, 辨识的输入信号取“-0.03” else u(i)=0.03; %如果M序列的值为0, 辨识的输入信号取“0.03” end %小循环结束 y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备end %大循环结束,产生输入信号u figure(1); %第一个图形stem(u),grid on %显示出输入信号径线
22、图并给图形加上网格z(2)=0;z(1)=0; %设z的前两个初始值为零for k=3:15; %循环变量从3到15 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %输出采样信号 end%RLS递推最小二乘辨识c0=0.001 0.001 0.001 0.001; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(4,4); %直接给出初始状态P0,即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=c0,zeros(4,14); %被辨识参数矩阵的初始值及大小e=zeros(4,
23、15); %相对误差的初始值及大小for k=3:15; %开始求K h1=-z(k-1),-z(k-2),u(k-1),u(k-2); x=h1*p0*h1+1; x1=inv(x); %开始求l(k) k1=p0*h1*x1;%求出l的值 d1=z(k)-h1*c0; c1=c0+k1*d1; %求被辨识参数c e1=c1-c0; %求参数当前值与上一次的值的差值 e2=e1./c0; %求参数的相对变化 e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1; %新获得的参数作为下一次递推的旧参数 c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最
24、后一列 p1=p0-k1*k1*h1*p0*h1+1; %求出 p(k)的值 p0=p1; %给下次用 if e2c = 0.0010 0 0.0010 -0.4984 -1.2328 -1.4951 -1.4962 -1.4991 -1.4998 -1.4999 0.0010 0 0.0010 0.0010 -0.2350 0.6913 0.6941 0.6990 0.6998 0.6999 0.0010 0 0.2509 1.2497 1.0665 1.0017 1.0020 1.0002 0.9999 0.9998 0.0010 0 -0.2489 0.7500 0.5668 0.502
25、0 0.5016 0.5008 0.5002 0.5002图3.8 最小二乘递推算法的参数辨识仿真-1.5000 -1.5000 -1.5000 -1.4999 -1.49990.7000 0.7000 0.7000 0.7000 -0.70000.9999 0.9999 0.9999 0.9999 0.99990.5000 0.5000 0.5000 0.5000 0.5000e = 0 0 0 -499.4200 1.4734 0.2128 0.0007 0.0020 0.0004 0.0000 0 0 0 0 -235.9916 -3.9416 0.0042 0.0070 0.0012
26、0.0001 0 0 249.8612 3.9816 -0.1466 -0.0607 0.0003 -0.0018 -0.0003 -0.0001 0 0 -249.8612 -4.0136 -0.2443 -0.1143 -0.0007 -0.0016 -0.0012 -0.00010.0001 0.0000 -0.0000 -0.0000 0.00000.0001 -0.0000 0.0000 0.0000 0.00000.0001 0.0000 0.0000 0.0000 -0.0000-0.0004 0.0000 -0.0000 0.0000 -0.0000表3.2 最小二乘递推算法的辨识结果参 数 a1 a2 b1 b2真 值 -1.5 0.7 1.0 0.5 估计值 -1.4999 0.7 0.9999 0.5000仿真结果表明,大约递推到第十步时,参数辨识的结果基本达到稳定状态,即a1=-1.4999, a2= 0.7000, b1=0.9999, b2=0.5000。此时,参数的相对变化量E0.000000005。从整个辨识过程来看,精度的要求直接影响辨识的速度。虽然最终的精度可以达到很小,但开始阶段的相对误差会很大,从图3.8(3)图形中可看出,参数的最大相对误差会达到三维数
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1