代码白噪声的产生和应用.docx
《代码白噪声的产生和应用.docx》由会员分享，可在线阅读，更多相关《代码白噪声的产生和应用.docx（20页珍藏版）》请在冰豆网上搜索。
代码白噪声的产生和应用
1.乘同余法产生（0,1）均匀分布随机数
%用乘同余法产生（0-1）均匀分布的随机序列v
%也可以产生任意（a,b）均匀分布的随机序列a+v*（b-a）
%初始化
clc;clear;
x0=1;%正奇数
N=200;%递归次数，产生随机数的个数
kk=8;
M=2^kk;
x=zeros（2*M,N）;
%本程序A的取值范围为[3,4*M+1]，用x保存每个A对应的序列，再除M得到均匀分布随机序列
v=zeros（1,N）;%用于暂存一条随机序列
%画出不同的A值对应的随机序列的周期图
fori=1:
2*M
A=2*i+1;%A取奇数
fork=1:
N%乘同余法递推N次
ifk==1
x（i,1）=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是一一对应的
end
i=1:
2*M;
figure
（1）
plot（2*i+1,log2（bb）,'.'）%A=2*i+1
title（'M=2^8时，不同的A值对应的随机序列的周期'）
xlabel（'A'）;ylabel（'log2（周期）'）;
axis（[1,2^10+10,0,7]）
%用矩阵显示出不同的A值对应的随机序列的周期，以便求出A与周期的关系
p=1;TT=zeros（kk-2,8）;%初始化,TT矩阵用于存储周期和对应的A
forj=kk-2:
-1:
0%所有的周期都是2的指数，且最大指数为kk-2
w=find（bb==2^j）;%找出周期为2^j的所有位置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）=2^j;
TT（p,2）=NaN;
p=p+1;
end
display（'不同的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;
fori=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）^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）;
end
end
if（T==0）%即不满足以上几个条件的情况下周期为2
T=2;
end
display（'当前给定A值对应的随机序列的周期：
'）
T
%求随机序列的自相关函数
[r,f]=xcorr（v,'unbiased'）;
figure（3）%画出自相关函数图形
plot（f,r）
title（'随机序列自相关函数'）
输出结果分析：
不同的A值对应的随机序列的周期：
T=
64NaN3511131921…A=4*i-2-（-1）^i
32NaN7923253941…A=8*i-4-（-1）^i*3
16NaN151747497981…A=16*i-8-（-1）^i*7
8NaN31339597159161…A=32*i-16-（-1）^i*15
4NaN6365191193319321…A=64*i-32-（-1）^i*31
2NaN127129255383385511…
1NaN2575137691025…A=256*i+1
A=179时产生的随机序列，周期：
T=64
2.混合同余法产生（0,1）均匀分布随机数
%用混合同余法产生（0-1）均匀分布的随机序列v
%也可以产生任意（a,b）均匀分布的随机序列a+v*（b-a）
%初始化
clc
x0=1;%非负整数
A=2^2+1;%典型值,n在[2,34]内取值,n越大，随机序列的波形越趋于周期三角波
N=1000;%递归次数
kk=8;
M=2^kk;d=10;
x=zeros（2^d,N）;
v=zeros（1,N）;
%画出不同的A值对应的随机序列的周期图
forc=1:
2^d
fork=1:
N%乘同余法递推N次
ifk==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次结束
a=find（v==min（v））;
b（c）=a
（2）-a
（1）;%求出每个c对应的周期
end
figure
（1）
plot（log2（b）,'.'）
title（'M=2^8时，不同的c对应的随机序列的周期'）
xlabel（'c'）;ylabel（'log2（周期）'）;
axis（[1,2^d,0,9]）
%用矩阵显示出不同的c值对应的随机序列的周期，以便求出c与周期的关系
p=1;TT=zeros（kk,8）;%初始化,TT矩阵用于存储周期和对应的c
forj=kk:
-1:
0%所有的周期都是2的指数，且最大指数为kk
w=find（b==2^j）;%找出周期为2^j的所有位置i
%TT每一行第一个数为周期，第二个置NaN作为间隔，其余的为周期对应的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）=2^j;
TT（p,2）=NaN;
p=p+1;
end
display（'M=2^8时，不同的c值对应的随机序列的周期：
'）
TT
%求当前给定c值对应的随机序列
c=1;
k=1:
N;
v=x（c,:
）/M;
figure
（2）
plot（k,v,'r'）;
ylabel（'随机数'）;
title（'（0-1）均匀分布的随机序列'）
%求当前给定c值对应的随机序列的周期
fori=1:
2^（kk-1）;
if（c==2*i-1）
T=2^kk;
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）;
end
end
display（'c值对应的随机序列的周期：
'）
T
%求随机序列的自相关函数
%v=-1+2*v;%产生（-1,1）均匀分布随机序列
[a,b]=xcorr（v,'unbiased'）;
figure（3）%画出自相关函数图形
plot（b,a）
title（'随机序列自相关函数'）
输出结果分析：
M=2^8时，不同的c值对应的随机序列的周期：
T=
256NaN1357911…c=2*i-1
128NaN2610141822…c=4*i-2
64NaN81624324048…c=8*i
32NaN42036526884…c=16*i-12
16NaN124476108140172…c=32*i-20
8NaN2892156220284348…c=64*i-36
4NaN60188316444572700…c=128*i-68
2NaN124380636892…c=256*i-132
1NaN2525087641020…c=256*i-4
当c=1时，得到的随机序列如下，周期：
T=256
3.统计近似抽样法产生正态分布随机数
%产生均值为un,方差为Sigma的正态分布随机数
%资料上的作法是错的
%初始化
clc
N1=600;%产生正态分布随机数的个数
un=0;%正态分布均值为
Sigma=1;%正态分布方差为
x0=1;%正奇数
N2=N1+11;%递归次数，产生[0,1]均匀分布随机数的个数,每相邻12的和记作ksai
kk=10;
M=2^kk;
A=179;%典型值,[0,1]均匀分布随机数周期为2^（kk-2）
x=zeros（1,N2）;
a=zeros（1,N2）;
v=zeros（1,N1）;
%产生随机数
fork=1:
N1
fori=1:
N2
ifi==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）;
end
k=1:
N1;
plot（k,v,'r'）
title（'均值为un,方差为Sigma的正态分布随机数'）
%求随机序列的自相关函数
[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）;%分布的正态性检验
%红线是正态分布累积函数经过变换得到的，
%蓝色的是数据（纵坐标基本上是概率累积，但是经过处理），
%如果数据和正太分布概率累积吻合（几乎在这条直线上），就说明
升级会员 