白噪声.doc

白噪声白噪声 如果一个零均值、平稳随机过程的谱密度为常数，我们称之为白噪声（由 白色光联想而得） 。白噪声有以下特点特点： 1 0wE(w) 2 ，频谱宽度无限。 2 )( w S 3 ，其中，为 Dirac 函数，即=)()( 2 w R)()( 且 1)(d 4 无记忆性，即 t 时刻的数值与 t 时刻以前的过去值无关，也不影响 t 时 刻以后的将来值。从另一意义上说，即不同时刻的随机信号互不相关。 白噪声的用途用途： 1 作为系统输入时，有，=0,1,2,，即为系统的单位脉冲)()(gRwy 响应。 2 作为被辨识系统输入时，可以激发系统的所有模态，可对系统充分激励； 3 作为被辨识系统输入时，可防止数据病态，保证辨识精度。 4 在辨识过程中，以输出估计误差是否具有白色性来判断辨识方法的优劣， 也可用来判断模型的结构和参数是否合适。 5 产生有色噪声。 白噪声的产生方法产生方法： 1 （0，1）均匀分布白噪声：，i=1,2,3，初值*179 1 ii 取小数 可取为： 350 2 11 2 正态分布白噪声：，),( 2 N 12 1 12)*1( )6()( i ki k 其中 为服从（0，1）均匀分布的白噪声。 例例 2.1 用乘同余法产生随机数用乘同余法产生随机数 选选 A=6,k=8,M=28=256,递推递推 100 次次.用用 MATLAB 编程编程 编程如下： ，=0 0，0 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存储器的数除以256得到小于1的随机数放v1中； v(:,k)=v1; % 将v1中的数（）存放在矩阵存储器v的第k列中，v(:,k) i %表示行不变、列随递推循环次数变化； 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)均匀分布的随机序列 在程序运行结束后，产生的(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.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.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.8438 0.0820 。 编程如下： 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）中； v1=x1/256; %将x1存储器中的数除以256得到小于1的随机数放v1中； v(:,k)=(v1-0.5 )*f; %将v1中的数（）减去0.5再乘以存储器f中的系数，存放 i 在矩阵存储器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)均匀分布的白噪声) 程序运行结果如图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.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.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.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 例例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 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 = 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 through 60 1 1 -1 -1 -1 1 -1 -1 1 1 i1 = 60 最小二乘一次完成算法的最小二乘一次完成算法的MATLAB仿真仿真 考虑仿真对象考虑仿真对象 )()2(5 . 0) 1()2(7 . 0) 1(5 . 1)(kkukukzkzkz 式中是服从正态分布的白噪声 N(0,1),输入信号采用 4 阶 M 序列,幅度为)(k 1。选择如下形式的辨识模型 )()2() 1()2() 1()( 2121 kkubkubkzakzakz 构造 和；数据长度取 L=14；加权阵取;利用 L z L I L 计算参数估计值。 L T LL T L LSz 1 )( LS 设输入信号的取值是从 k=1 到 k=16 的 M 序列。 待辨识参数，观测矩阵和的表达式为： LS L z L ， 2 1 2 1 b b a a LS )16( )4( )3( z z z zL )14()15()14()15( )2()3()2()3( ) 1 ()2() 1 ()2( uuzz uuzz uuzz L 仿真程序流程图： %二阶系统的最小二乘一次完成算法辨识程序。 u=-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1; %系统辨识的输入信号为一个周 期的 M 序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作 为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 赋输入信号初值 u 定义输出观测值的长度并计算系统的输出值 画出输入和输出观测值的图形 给样本矩阵和 zL赋值 L 根据式一次算法计算参数 LS 从中分离出并显示出被辨识参数 a1, a2, b1, b2 LS 停机 最小二乘一次完成算法程序框图 stem(u) %画输入信号 u 的径线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是 1 到 16，步长为 1 plot(i,z) %图形的横坐标是采样时刻 i, 纵坐标是输出观测值 z, 图形 格式为连续曲线 subplot(3,1,3) %画三行一列图形窗口中的第三个图形 stem(z),grid on %画出输出观测值 z 的径线图形，并显示坐标网格 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) 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 %计算并显示 LS %Display Parameters a1=c(1), a2=c(2), b1=c(3),b2=c(4) %从中分离出并显示 a1 、a2、 LS b1、 b2 %End 程序运行结果： u = -1，1，-1，1，1，1，1，-1，-1，-1，1，-1，-1，1，1 z = 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.5185 HL = 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.5250 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.0000 ZL = 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.5185T c = -1.5000，0.7000，1.0000，0.5000T a1 = -1.5000 a2 = 0.7000 b1 = 1.0000 b2 =0.5000 051015 -1 0 1 0246810121416 -10 0 10 0246810121416 -10 0 10 STEM PLOT OF INPUT M SERIES WAVE OF MEASUREMENTS OF OUTPUT STEM PLOT OF MEASUREMENTS OF OUTPUT 从仿真结果表可以看出，由于所用的输出观测值没有任何噪声 成分，所以辨识结果也无任何误差。 最小二乘一次完成算法仿真实例中输入信号和输出观测值 考虑图 3.所示的仿真对象，图中， 是服从 N分布的不)(kv) 1 , 0( 相关随机噪声。且 ， ， )( 1 zG )( )( 1 1 zA zB )( 1 zN )( )( 1 1 zC zD 1)( 5 . 00 . 1)( )(7 . 05 . 11)( 1 211 121 1 1 zD zzzB zCzzazA 选择图所示的辨识模型。仿真对象选择如下的模型结构 )()2() 1()2() 1()( 2121 kvkubkubkzakzakz 其中，是服从正态分布的白噪声 N。输入信号采用 4 位移位)(kv) 1 , 0( 寄存器产生的 M 序 列，幅度为 0.03。按式 )()2(5 . 0) 1()2(7 . 0) 1(5 . 1)(kvkukukzkzkz 构造；加权阵取单位阵；利用式递推公式计算 l(k)、)(kI L 和 P(k)，计算各次参数辨识的相对误差，精度满足要求式适当小 )( k 真 值 -1.5 0.7 1.0 0.5 估计值 -1.5 0.7 1.0 0.5 参 数 a1 a2 b1 b2 表表 最小二乘一次完成算法的辨识结果最小二乘一次完成算法的辨识结果 + e(k) 图 3.最小二乘递推算法辨识实例结构图 y(k) u(k)z(k) v(k) )( 1 zN )( 1 zG 的数(当所有参数估计值变化不大时) 后停机。 下面给出具体程序。 %最小二乘递推算法辨识程序, 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; %第三个移位寄存器的输入是第二个移位寄存器的输出 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 %显示出输入信号径线图并给图形加上网格 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.000000005 c=c0,zeros(4,14); %被辨识参数矩阵的初始值及大小 e=zeros(4,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 列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1*h1*p0*h1+1; %求出 p(k)的值 p0=p1; %给下次用 if e2 c = 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.5020 0.5016 0.5008 0.5002 0.5002 051015 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 Fig.1 Input Signal 051015 -1.5 -1 -0.5 0 0.5 1 1.5 a1 a2 b1 b2 Fig. 2 Parameter Identification with Recursive Least Squares Method 051015 -500 -400 -300 -200 -100 0 100 200 300 Ident