代碼-白噪聲的產(chǎn)生和應(yīng)用

1.乘同余法產(chǎn)生(0,1)均勻分布隨機(jī)數(shù)%用乘同余法產(chǎn)生(0-1)均勻分布的隨機(jī)序列v%也可以產(chǎn)生任意(a,b)均勻分布的隨機(jī)序列a+v*(b-a)%初始化clc;clear;x0=1;%正奇數(shù)N=200;%遞歸次數(shù)，產(chǎn)生隨機(jī)數(shù)的個(gè)數(shù)kk=8;M=2^kk;x=zeros(2*M,N);%本程序A的取值范圍為[3,4*M+1]，用x保存每個(gè)A對(duì)應(yīng)的序列，再除M得到均勻分布隨機(jī)序列v=zeros(1,N);%用于暫存一條隨機(jī)序列%畫(huà)出不同的A值對(duì)應(yīng)的隨機(jī)序列的周期圖fori=1:2*MA=2*i+1;%A取奇數(shù)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的隨機(jī)數(shù)放v中end%遞推N次結(jié)束aa=find(v==min(v));%找出所有最小值的位置bb(i)=aa(2)-aa(1);%求出每個(gè)A對(duì)應(yīng)的周期,序號(hào)i與A是一一對(duì)應(yīng)的endi=1:2*M;figure(1)plot(2*i+1,log2(bb),'.')%A=2*i+1title('M=2^8時(shí)，不同的A值對(duì)應(yīng)的隨機(jī)序列的周期')xlabel('A');ylabel('log2(周期)');axis([1,2^10+10,0,7])%用矩陣顯示出不同的A值對(duì)應(yīng)的隨機(jī)序列的周期，以便求出A與周期的關(guān)系p=1;TT=zeros(kk-2,8);%初始化,TT矩陣用于存儲(chǔ)周期和對(duì)應(yīng)的Aforj=kk-2:-1:0%所有的周期都是2的指數(shù)，且最大指數(shù)為kk-2w=find(bb==2^j);%找出周期為2^j的所有位置i%TT每一行第一個(gè)數(shù)為周期，第二個(gè)置NaN作為間隔，其余的為周期對(duì)應(yīng)的A值%TT每個(gè)周期只取前6個(gè)值，不足6個(gè)的補(bǔ)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值對(duì)應(yīng)的隨機(jī)序列的周期：')TT%求當(dāng)前給定A值對(duì)應(yīng)的隨機(jī)序列A=179;%典型值k=1:N;v=x((A-1)/2,:)/M;figure(2)plot(k,v,'r');%畫(huà)出隨機(jī)序列圖ylabel('隨機(jī)數(shù)');title('(0-1)均勻分布的隨機(jī)序列')%求當(dāng)前給定A值對(duì)應(yīng)的隨機(jī)序列的周期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)%即不滿(mǎn)足以上幾個(gè)條件的情況下周期為2T=2;enddisplay('當(dāng)前給定A值對(duì)應(yīng)的隨機(jī)序列的周期：')T%求隨機(jī)序列的自相關(guān)函數(shù)[r,f]=xcorr(v,'unbiased');figure(3)%畫(huà)出自相關(guān)函數(shù)圖形plot(f,r)title('隨機(jī)序列自相關(guān)函數(shù)')輸出結(jié)果分析：不同的A值對(duì)應(yīng)的隨機(jī)序列的周期：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時(shí)產(chǎn)生的隨機(jī)序列，周期：T=642.混合同余法產(chǎn)生(0,1)均勻分布隨機(jī)數(shù)%用混合同余法產(chǎn)生(0-1)均勻分布的隨機(jī)序列v%也可以產(chǎn)生任意(a,b)均勻分布的隨機(jī)序列a+v*(b-a)%初始化clcx0=1;%非負(fù)整數(shù)A=2^2+1;%典型值,n在[2,34]內(nèi)取值,n越大，隨機(jī)序列的波形越趨于周期三角波N=1000;%遞歸次數(shù)kk=8;M=2^kk;d=10;x=zeros(2^d,N);v=zeros(1,N);%畫(huà)出不同的A值對(duì)應(yīng)的隨機(jī)序列的周期圖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的隨機(jī)數(shù)放v中end%遞推N次結(jié)束a=find(v==min(v));b(c)=a(2)-a(1);%求出每個(gè)c對(duì)應(yīng)的周期endfigure(1)plot(log2(b),'.')title('M=2^8時(shí)，不同的c對(duì)應(yīng)的隨機(jī)序列的周期')xlabel('c');ylabel('log2(周期)');axis([1,2^d,0,9])%用矩陣顯示出不同的c值對(duì)應(yīng)的隨機(jī)序列的周期，以便求出c與周期的關(guān)系p=1;TT=zeros(kk,8);%初始化,TT矩陣用于存儲(chǔ)周期和對(duì)應(yīng)的cforj=kk:-1:0%所有的周期都是2的指數(shù)，且最大指數(shù)為kkw=find(b==2^j);%找出周期為2^j的所有位置i%TT每一行第一個(gè)數(shù)為周期，第二個(gè)置NaN作為間隔，其余的為周期對(duì)應(yīng)的c值%TT每個(gè)周期只取前6個(gè)值，不足6個(gè)的補(bǔ)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時(shí)，不同的c值對(duì)應(yīng)的隨機(jī)序列的周期：')TT%求當(dāng)前給定c值對(duì)應(yīng)的隨機(jī)序列c=1;k=1:N;v=x(c,:)/M;figure(2)plot(k,v,'r');ylabel('隨機(jī)數(shù)');title('(0-1)均勻分布的隨機(jī)序列')%求當(dāng)前給定c值對(duì)應(yīng)的隨機(jī)序列的周期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值對(duì)應(yīng)的隨機(jī)序列的周期：')T%求隨機(jī)序列的自相關(guān)函數(shù)%v=-1+2*v;%產(chǎn)生(-1,1)均勻分布隨機(jī)序列[a,b]=xcorr(v,'unbiased');figure(3)%畫(huà)出自相關(guān)函數(shù)圖形plot(b,a)title('隨機(jī)序列自相關(guān)函數(shù)')輸出結(jié)果分析：M=2^8時(shí)，不同的c值對(duì)應(yīng)的隨機(jī)序列的周期：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當(dāng)c=1時(shí)，得到的隨機(jī)序列如下，周期：T=2563.統(tǒng)計(jì)近似抽樣法產(chǎn)生正態(tài)分布隨機(jī)數(shù)%產(chǎn)生均值為un,方差為Sigma的正態(tài)分布隨機(jī)數(shù)%資料上的作法是錯(cuò)的%初始化clcN1=600;%產(chǎn)生正態(tài)分布隨機(jī)數(shù)的個(gè)數(shù)un=0;%正態(tài)分布均值為Sigma=1;%正態(tài)分布方差為x0=1;%正奇數(shù)N2=N1+11;%遞歸次數(shù)，產(chǎn)生[0,1]均勻分布隨機(jī)數(shù)的個(gè)數(shù),每相鄰12的和記作ksaikk=10;M=2^kk;A=179;%典型值,[0,1]均勻分布隨機(jī)數(shù)周期為2^(kk-2)x=zeros(1,N2);a=zeros(1,N2);v=zeros(1,N1);%產(chǎn)生隨機(jī)數(shù)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的正態(tài)分布隨機(jī)數(shù)')%求隨機(jī)序列的自相關(guān)函數(shù)[a,b]=xcorr(v,'unbiased');figure(2)%畫(huà)出自相關(guān)函數(shù)圖形plot(b,a)title('隨機(jī)序列自相關(guān)函數(shù)')輸出結(jié)果分析：產(chǎn)生均值為0，方差為1的正態(tài)分布4.變換抽樣法產(chǎn)生正態(tài)分布隨機(jī)數(shù)clc;clear;a=rand(1,10000);%直接用函數(shù)產(chǎn)生(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);%作頻數(shù)直方圖,表示對(duì)c分成20區(qū)間進(jìn)行統(tǒng)計(jì)figure(2);normplot(c);%分布的正態(tài)性檢驗(yàn)%紅線是正態(tài)分布累積函數(shù)經(jīng)過(guò)變換得到的，%藍(lán)色的是數(shù)據(jù)(縱坐標(biāo)基本上是概率累積，但是經(jīng)過(guò)處理)，%如果數(shù)據(jù)和正太分布概率累積吻合（幾乎在這條直線上），就說(shuō)明是正態(tài)分布。[muhat,sigmahat,muci,sigmaci]=normfit(c)%參數(shù)估計(jì)均值,標(biāo)準(zhǔn)差,均值的0.95置信區(qū)間,方差的0.95置信區(qū)間[h,sig,ci]=ttest(c,muhat)%假設(shè)檢驗(yàn)單正態(tài)均值t檢驗(yàn)(方差未知)用統(tǒng)計(jì)量t=sqrt(n)*(mean(x)-u0)/s%h=0，si>0.05,接受原假設(shè),ci為95%置信區(qū)間%求隨機(jī)序列的自相關(guān)函數(shù)[ac,bc]=xcorr(c,'unbiased');figure(3)%畫(huà)出自相關(guān)函數(shù)圖形plot(bc,ac)title('隨機(jī)序列c自相關(guān)函數(shù)')[ad,bd]=xcorr(d,'unbiased');figure(4)%畫(huà)出自相關(guān)函數(shù)圖形plot(bd,ad)title('隨機(jī)序列d自相關(guān)函數(shù)')輸出結(jié)果分析：驗(yàn)證服從正態(tài)分布參數(shù)估計(jì)：均值muhat=0.0021，均值95%置信區(qū)間muci=(-0.0173,0.0216)標(biāo)準(zhǔn)差sigmahat=0.9928，標(biāo)準(zhǔn)差95%置信區(qū)間sigmaci=(0.9793,1.0068)單正態(tài)均值t檢驗(yàn)(方差未知)：h=0，sig=1，均值95%置信區(qū)間ci=(-0.0173,0.0216)h=0，si>0.05,接受原假設(shè)5.M序列clc;clear;%c必須為本原多項(xiàng)式的系數(shù)%c是反饋系數(shù)矩陣，省略c0,從左到右一次是c1到cn%3級(jí)c=[101];c=[01001];%5級(jí)%7級(jí)c=[0010001];%9級(jí)c=[000100001];n=length(c);