代码-白噪声的产生和应用

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*MA=2*i+1;%A取奇数fork=1:N%乘同余法递推N次ifk==1x(i,1)=mod(A*x0,M);elsex(i,k)=mod(A*x(i,k-1),M);endv(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=2^8时，不同的A值对应的随机序列的周期')xlabel('A');ylabel('log2(周期)');axis([1,2^10+10,0,7])%用矩阵显示出不同的A值对应的随机序列的周期，以便求出A与周期的关系p=1;TT=zeros(kk-2,8);%初始化,TT矩阵用于存储周期和对应的Aforj=kk-2:-1:0%所有的周期都是2的指数，且最大指数为kk-2w=find(bb==2^j);%找出周期为2^j的所有位置i%TT每一行第一个数为周期，第二个置NaN作为间隔，其余的为周期对应的A值%TT每个周期只取前6个值，不足6个的补0if(length(w)<6)TT(p,3:2+length(w))=2*w(1:length(w))+1;elseTT(p,3:8)=2*w(1:6)+1;endTT(p,1)=2^j;TT(p,2)=NaN;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;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)%周期为1T=2^(kk-8);endendif(T==0)%即不满足以上几个条件的情况下周期为2T=2;enddisplay('当前给定A值对应的随机序列的周期：')T%求随机序列的自相关函数[r,f]=xcorr(v,'unbiased');figure(3)%画出自相关函数图形plot(f,r)title('随机序列自相关函数')输出结果分析：不同的A值对应的随机序列的周期：T=64NaN3511131921…A=4*i-2-(-1)^i32NaN7923253941…A=8*i-4-(-1)^i*316NaN151747497981…A=16*i-8-(-1)^i*78NaN31339597159161…A=32*i-16-(-1)^i*154NaN6365191193319321…A=64*i-32-(-1)^i*312NaN127129255383385511…1NaN2575137691025…A=256*i+1A=179时产生的随机序列，周期：T=642.混合同余法产生(0,1)均匀分布随机数%用混合同余法产生(0-1)均匀分布的随机序列v%也可以产生任意(a,b)均匀分布的随机序列a+v*(b-a)%初始化clcx0=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^dfork=1:N%乘同余法递推N次ifk==1x(c,k)=mod(A*x0+c,M);elsex(c,k)=mod(A*x(c,k-1)+c,M);endv(k)=x(c,k)/M;%将x除以M得到小于1的随机数放v中end%递推N次结束a=find(v==min(v));b(c)=a(2)-a(1);%求出每个c对应的周期endfigure(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矩阵用于存储周期和对应的cforj=kk:-1:0%所有的周期都是2的指数，且最大指数为kkw=find(b==2^j);%找出周期为2^j的所有位置i%TT每一行第一个数为周期，第二个置NaN作为间隔，其余的为周期对应的c值%TT每个周期只取前6个值，不足6个的补0if(length(w)<6)TT(p,3:2+length(w))=w(1:length(w));elseTT(p,3:8)=w(1:6);endTT(p,1)=2^j;TT(p,2)=NaN;p=p+1;enddisplay('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);endenddisplay('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-1128NaN2610141822…c=4*i-264NaN81624324048…c=8*i32NaN42036526884…c=16*i-1216NaN124476108140172…c=32*i-208NaN2892156220284348…c=64*i-364NaN60188316444572700…c=128*i-682NaN124380636892…c=256*i-1321NaN2525087641020…c=256*i-4当c=1时，得到的随机序列如下，周期：T=2563.统计近似抽样法产生正态分布随机数%产生均值为un,方差为Sigma的正态分布随机数%资料上的作法是错的%初始化clcN1=600;%产生正态分布随机数的个数un=0;%正态分布均值为Sigma=1;%正态分布方差为x0=1;%正奇数N2=N1+11;%递归次数，产生[0,1]均匀分布随机数的个数,每相邻12的和记作ksaikk=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:N1fori=1:N2ifi==1x(1)=mod(A*x0,M);elsex(i)=mod(A*x(i-1),M);enda(i)=x(i)/M;endksai=sum(a(k:k+11));%每相邻12的和记作ksaiv(k)=un+Sigma*(ksai-6);endk=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);%分布的正态性检验%红线是正态分布累积函数经过变换得到的，%蓝色的是数据(纵坐标基本上是概率累积，但是经过处理)，%如果数据和正太分布概率累积吻合（几乎在这条直线上），就说明是正态分布。[muhat,sigmahat,muci,sigmaci]=normfit(c)%参数估计均值,标准差,均值的0.95置信区间,方差的0.95置信区间[h,sig,ci]=ttest(c,muhat)%假设检验单正态均值t检验(方差未知)用统计量t=sqrt(n)*(mean(x)-u0)/s%h=0，si>0.05,接受原假设,ci为95%置信区间%求随机序列的自相关函数[ac,bc]=xcorr(c,'unbiased');figure(3)%画出自相关函数图形plot(bc,ac)title('随机序列c自相关函数')[ad,bd]=xcorr(d,'unbiased');figure(4)%画出自相关函数图形plot(bd,ad)title('随机序列d自相关函数')输出结果分析：验证服从正态分布参数估计：均值muhat=0.0021，均值95%置信区间muci=(-0.0173,0.0216)标准差sigmahat=0.9928，标准差95%置信区间sigmaci=(0.9793,1.0068)单正态均值t检验(方差未知)：h=0，sig=1，均值95%置信区间ci=(-0.0173,0.0216)h=0，si>0.05,接受原假设5.M序列clc;clear;%c必须为本原多项式的系数%c是反馈系数矩阵，省略c0,从左到右一次是c1到cn%3级c=[101];c=[01001];%5级%7级c=[0010001];%9级c=[000100001];n=length(c);