代码白噪声的产生和应用.docx

1、代码白噪声的产生和应用1.乘同余法产生(0,1)均匀分布随机数%用乘同余法产生(0-1)均匀分布的随机序列v%也可以产生任意(a,b)均匀分布的随机序列a+v*(b-a)%初始化clc;clear;x0=1;%正奇数N=200;%递归次数，产生随机数的个数kk=8;M=2kk;x=zeros(2*M,N);%本程序A的取值范围为3,4*M+1，用x保存每个A对应的序列，再除M得到均匀分布随机序列v=zeros(1,N);%用于暂存一条随机序列%画出不同的A值对应的随机序列的周期图for i=1:2*M A=2*i+1;%A取奇数 for k=1:N %乘同余法递推N次 if k=1 x(i,1

2、)=mod(A*x0,M); else x(i,k)=mod(A*x(i,k-1),M); end v(k)=x(i,k)/M; %将x除以M得到小于1的随机数放v中 end %递推N次结束 aa=find(v=min(v);%找出所有最小值的位置 bb(i)=aa(2)-aa(1);%求出每个A对应的周期,序号i与A是一一对应的endi=1:2*M;figure(1)plot(2*i+1,log2(bb),.)%A=2*i+1title(M=28时，不同的A值对应的随机序列的周期)xlabel(A);ylabel(log2(周期);axis(1,210+10,0,7)%用矩阵显示出不同的A值

3、对应的随机序列的周期，以便求出A与周期的关系p=1;TT=zeros(kk-2,8);%初始化,TT矩阵用于存储周期和对应的Afor j=kk-2:-1:0%所有的周期都是2的指数，且最大指数为kk-2 w=find(bb=2j);%找出周期为2j的所有位置i %TT每一行第一个数为周期，第二个置NaN作为间隔，其余的为周期对应的A值 %TT每个周期只取前6个值，不足6个的补0 if(length(w)6) TT(p,3:2+length(w)=2*w(1:length(w)+1; else TT(p,3:8)=2*w(1:6)+1; end TT(p,1)=2j; TT(p,2)=NaN;

4、p=p+1;enddisplay(不同的A值对应的随机序列的周期：)TT%求当前给定A值对应的随机序列A=179;%典型值k=1:N ;v=x(A-1)/2,:)/M;figure(2)plot(k,v,r);%画出随机序列图ylabel(随机数);title(0-1)均匀分布的随机序列)%求当前给定A值对应的随机序列的周期T=0;for i=1:2*M; if(A=4*i-2-(-1)i) T=2(kk-2); elseif(A=8*i-4-(-1)i*3) T=2(kk-3); elseif(A=16*i-8-(-1)i*7) T=2(kk-4); elseif(A=32*i-16-(-1

5、)i*15) T=2(kk-5); elseif(A=64*i-32-(-1)i*31) T=2(kk-6); elseif(A=256*i+1)%周期为1 T=2(kk-8); endendif(T=0)%即不满足以上几个条件的情况下周期为2 T=2;enddisplay(当前给定A值对应的随机序列的周期：)T%求随机序列的自相关函数r,f=xcorr(v,unbiased);figure(3)%画出自相关函数图形plot(f,r)title(随机序列自相关函数)输出结果分析：不同的A值对应的随机序列的周期：T =64 NaN 3 5 11 13 19 21 A=4*i-2-(-1)i32

6、NaN 7 9 23 25 39 41 A=8*i-4-(-1)i*316 NaN 15 17 47 49 79 81 A=16*i-8-(-1)i*78 NaN 31 33 95 97 159 161 A=32*i-16-(-1)i*154 NaN 63 65 191 193 319 321 A=64*i-32-(-1)i*312 NaN 127 129 255 383 385 511 1 NaN 257 513 769 1025 A=256*i+1A=179时产生的随机序列，周期：T =642.混合同余法产生(0,1)均匀分布随机数%用混合同余法产生(0-1)均匀分布的随机序列v%也可以产

7、生任意(a,b)均匀分布的随机序列a+v*(b-a)%初始化clcx0=1;%非负整数A=22+1;%典型值,n在2,34内取值,n越大，随机序列的波形越趋于周期三角波N=1000;%递归次数kk=8;M=2kk;d=10;x=zeros(2d,N);v=zeros(1,N);%画出不同的A值对应的随机序列的周期图for c=1:2d for k=1:N %乘同余法递推N次 if k=1 x(c,k)=mod(A*x0+c,M); else x(c,k)=mod(A*x(c,k-1)+c,M); end v(k)=x(c,k)/M; %将x除以M得到小于1的随机数放v中 end %递推N次结束

8、 a=find(v=min(v); b(c)=a(2)-a(1);%求出每个c对应的周期endfigure(1)plot(log2(b),.)title(M=28时，不同的c对应的随机序列的周期)xlabel(c);ylabel(log2(周期);axis(1,2d,0,9)%用矩阵显示出不同的c值对应的随机序列的周期，以便求出c与周期的关系p=1;TT=zeros(kk,8);%初始化,TT矩阵用于存储周期和对应的cfor j=kk:-1:0%所有的周期都是2的指数，且最大指数为kk w=find(b=2j);%找出周期为2j的所有位置i %TT每一行第一个数为周期，第二个置NaN作为间隔，

9、其余的为周期对应的c值 %TT每个周期只取前6个值，不足6个的补0 if(length(w)6) TT(p,3:2+length(w)=w(1:length(w); else TT(p,3:8)=w(1:6); end TT(p,1)=2j; TT(p,2)=NaN; p=p+1;enddisplay(M=28时，不同的c值对应的随机序列的周期：)TT%求当前给定c值对应的随机序列c=1;k=1:N ;v=x(c,:)/M;figure(2)plot(k,v,r);ylabel(随机数);title(0-1)均匀分布的随机序列)%求当前给定c值对应的随机序列的周期for i=1:2(kk-1)

10、; if(c=2*i-1) T=2kk; elseif(c=4*i-2) T=2(kk-1); elseif(c=8*i) T=2(kk-2); elseif(c=16*i-12) T=2(kk-3); elseif(c=32*i-20) T=2(kk-4); elseif(c=64*i-36) T=2(kk-5); elseif(c=128*i-68) T=2(kk-6); elseif(c=256*i-132) T=2(kk-7); elseif(c=256*i-4) T=2(kk-8); endenddisplay(c值对应的随机序列的周期：)T%求随机序列的自相关函数%v=-1+2*v

11、;%产生(-1,1)均匀分布随机序列a,b=xcorr(v,unbiased);figure(3)%画出自相关函数图形plot(b,a)title(随机序列自相关函数)输出结果分析：M=28时，不同的c值对应的随机序列的周期：T =256 NaN 1 3 5 7 9 11 c=2*i-1128 NaN 2 6 10 14 18 22 c=4*i-264 NaN 8 16 24 32 40 48 c=8*i32 NaN 4 20 36 52 68 84 c=16*i-1216 NaN 12 44 76 108 140 172c=32*i-208 NaN 28 92 156 220 284 348

12、c=64*i-364 NaN 60 188 316 444 572 700c=128*i-682 NaN 124 380 636 892 c=256*i-1321 NaN 252 508 764 1020 c=256*i-4当c=1时，得到的随机序列如下，周期：T =2563.统计近似抽样法产生正态分布随机数%产生均值为un,方差为Sigma的正态分布随机数%资料上的作法是错的%初始化clcN1=600;%产生正态分布随机数的个数un=0;%正态分布均值为Sigma=1;%正态分布方差为x0=1;%正奇数N2=N1+11;%递归次数，产生0,1均匀分布随机数的个数,每相邻12的和记作ksaik

13、k=10;M=2kk;A=179;%典型值,0,1均匀分布随机数周期为2(kk-2)x=zeros(1,N2);a=zeros(1,N2);v=zeros(1,N1);%产生随机数for k=1:N1 for i=1:N2 if i=1 x(1)=mod(A*x0,M); else x(i)=mod(A*x(i-1),M); end a(i)=x(i)/M; end ksai=sum(a(k:k+11);%每相邻12的和记作ksai v(k)=un+Sigma*(ksai-6);endk=1:N1;plot(k,v,r)title(均值为un,方差为Sigma的正态分布随机数)%求随机序列的自

14、相关函数a,b=xcorr(v,unbiased);figure(2)%画出自相关函数图形plot(b,a)title(随机序列自相关函数)输出结果分析：产生均值为0，方差为1的正态分布4.变换抽样法产生正态分布随机数clc;clear;a=rand(1,10000);%直接用函数产生(0,1)均匀分布b=rand(1,10000);c=sqrt(-2*log(a).*cos(2*pi*b);d=sqrt(-2*log(a).*sin(2*pi*b);figure(1)hist(c,20);%作频数直方图 ,表示对c分成20区间进行统计figure(2);normplot(c);%分布的正态性检验 %红线是正态分布累积函数经过变换得到的，%蓝色的是数据(纵坐标基本上是概率累积，但是经过处理)，%如果数据和正太分布概率累积吻合（几乎在这条直线上），就说明