代碼-白噪聲的產生和應用

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.統(tǒng)計近似抽樣法產生正態(tài)分布隨機數%產生均值為un,方差為Sigma的正態(tài)分布隨機數%資料上的作法是錯的%初始化clcN1=600;%產生正態(tài)分布隨機數的個數un=0;%正態(tài)分布均值為Sigma=1;%正態(tài)分布方差為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的正態(tài)分布隨機數')%求隨機序列的自相關函數[a,b]=xcorr(v,'unbiased');figure(2)%畫出自相關函數圖形plot(b,a)title('隨機序列自相關函數')輸出結果分析：產生均值為0，方差為1的正態(tài)分布4.變換抽樣法產生正態(tài)分布隨機數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區(qū)間進行統(tǒng)計figure(2);normplot(c);%分布的正態(tài)性檢驗%紅線是正態(tài)分布累積函數經過變換得到的，%藍色的是數據(縱坐標基本上是概率累積，但是經過處理)，%如果數據和正太分布概率累積吻合（幾乎在這條直線上），就說明是正態(tài)分布。[muhat,sigmahat,muci,sigmaci]=normfit(c)%參數估計均值,標準差,均值的0.95置信區(qū)間,方差的0.95置信區(qū)間[h,sig,ci]=ttest(c,muhat)%假設檢驗單正態(tài)均值t檢驗(方差未知)用統(tǒng)計量t=sqrt(n)*(mean(x)-u0)/s%h=0，si>0.05,接受原假設,ci為95%置信區(qū)間%求隨機序列的自相關函數[ac,bc]=xcorr(c,'unbiased');figure(3)%畫出自相關函數圖形plot(bc,ac)title('隨機序列c自相關函數')[ad,bd]=xcorr(d,'unbiased');figure(4)%畫出自相關函數圖形plot(bd,ad)title('隨機序列d自相關函數')輸出結果分析：驗證服從正態(tài)分布參數估計：均值muhat=0.0021，均值95%置信區(qū)間muci=(-0.0173,0.0216)標準差sigmahat=0.9928，標準差95%置信區(qū)間sigmaci=(0.9793,1.0068)單正態(tài)均值t檢驗(方差未知)：h=0，sig=1，均值95%置信區(qū)間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);