實驗一隨機序列的產生及數字特征估計實驗目的1、學習和掌握隨機數的產生方法。2、實現隨機序列的數字特征估計。實驗原理1、隨機數的產生隨機數指的是各種不同分布隨機變量的抽樣序列（樣本值序列）。進行隨機信號仿真分析時，需要模擬產生各種分布的隨機數。在計算機仿真時，通常利用數學方法產生隨機數，這種隨機數稱為偽隨機數。偽隨機數是按照一定的計算公式產生的，這個公式稱為隨機數發生器。偽隨機數本質上不是隨機的，而且存在周期性，但是如果計算公式選擇適當，所產生的數據看似隨機的，與真正的隨機數具有相近的統計特性，可以作為隨機數使用。(0，1)均勻分布隨機數是最最基本、最簡單的隨機數。(0，1)均勻分布指的是在[0，1]區間上的均勻分布，即U(0，1)。實際應用中有許多現成的隨機數發生器可以用于產生(0，1)均勻分布隨機數，通常采用的方法為線性同余法，公式如下：序列為產生的(0，1)均勻分布隨機數。下面給出了上式的3組常用參數：1、，周期；2、（IBM隨機數發生器）周期；3、（ran0）周期；由均勻分布隨機數，可以利用反函數構造出任意分布的隨機數。定理1.1若隨機變量X具有連續分布函數FX(x)，而R為(0,1)均勻分布隨機變量，則有由這一定理可知，分布函數為FX(x)的隨機數可以由(0,1)均勻分布隨機數按上式進行變換得到。2、MATLAB中產生隨機序列的函數（1）(0，1)均勻分布的隨機序列函數：rand用法：x=rand(m,n)功能：產生m×n的均勻分布隨機數矩陣。（2）正態分布的隨機序列函數：randn用法：x=randn(m,n)功能：產生m×n的標準正態分布隨機數矩陣。如果要產生服從分布的隨機序列，則可以由標準正態隨機序列產生。（3）其他分布的隨機序列MATLAB上還提供了其他多種分布的隨機數的產生函數，下表列出了部分函數。MATLAB中產生隨機數的一些函數3、隨機序列的數字特征估計對于遍歷過程，可以通過隨機序列的一條樣本函數來獲得該過程的統計特性。這里我們假定隨機序列X(n)為遍歷過程，樣本函數為x(n)，其中n=0,1,2，…，N-1。那么，X(n)的均值、方差和自相關函數的估計為利用MATLAB的統計分析函數可以分析隨機序列的數字特征。（1）均值函數函數：mean用法：m=mean(x)功能：返回按上面第一式估計X(n)的均值，其中x為樣本序列x(n)。（2）方差函數函數：var用法：sigma2=var(x)功能：返回按上面第二式估計X(n)的方差，其中x為樣本序列x(n)，這一估計為無偏估計。（3）互相關函數函數：xcorr用法：c=xcorr(x,y)c=xcorr(x)c=xcorr(x,y,'opition')c=xcorr(x,'opition')功能：xcorr(x,y)計算X(n)與Y(n)的互相關，xcorr(x)計算X(n)的自相關。option選項可以設定為：'biased'有偏估計，即'unbiased'無偏估計，即按上面第三式估計。'coeff'm=0時的相關函數值歸一化為1。'none'不做歸一化處理。實驗內容1、采用線性同余法產生均勻分布隨機數1000個，計算該序列均值和方差與理論值之間的誤差大小。改變樣本個數重新計算。代碼及結果：>>Num=input('Num=');>>N=2^31;>>k=2^16+3;>>Y=zeros(1,num);>>X=zeros(1,num);>>Y(1)=1;>>fori=2:num>>Y(i)=mod(k*Y(i-1),N);>>end>>X=Y/N;>>a=0;>>b=1;>>m0=(a+b)/2;>>sigma0=(b-a)^2/12;>>m=mean(X);>>sigma=var(X);>>delta_m=abs(m-m0);>>delta_sigma=abs(sigma-sigma0);>>plot(X,'k');>>xlabel('n');>>ylabel('X(n)');>>delta_m>>delta_sigma>>axistightNum=1000delta_=0.0110delta_sigma=0.0011Num=5000delta_m=2.6620e-04delta_sigma=0.0020可以看出，樣本值取得越大，實際值與理論值越接近，相差越小。2、參數為的指數分布的分布函數為利用反函數法產生參數為0.5的指數分布隨機數1000個，測試其方差和相關函數。實驗代碼及結果：>>R=rand(1,1000);>>lambda=0.5;>>X=-log(1-R)/lambda;>>DX=var(X);>>[Rm,m]=xcorr(X);>>subplot(211);>>plot(X,'k');xlabel('n');ylabel('X(n)');axistight;>>subplot(212);>>plot(m,Rm,'k');xlabel('m');ylabel('R(m)');axistight;DX=4.1201理論上方差的值為1/(0.5^2)=4，實際值為4.1201，因為取樣個數有限，導致存在一定偏差。但大體相近。3、產生一組N(1,4)分布的高斯隨機數（1000個樣本），估計該序列的均值、方差和相關函數。實驗代碼及結果：>>X=normrnd(1,2,[1,1000]);>>Mx=mean(X);Dx=var(X);>>[Rm,m]=xcorr(X);>>subplot(211);>>plot(X,'k');xlabel('n');ylabel('X(n)');axistight;>>subplot(212);>>plot(m,Rm,'k');xlabel('m');ylabel('R(m)');axistight;Mx=0.9937Dx=3.8938理論上的均值為1，方差為4。而在實驗中得到的均值為0.9937，方差為3.8938。考慮到取樣點有限，誤差可以接受，理論值和實驗值基本相同。實驗二隨機過程的模擬與數字特征一、實驗目的1、學習利用MATLAB模擬產生隨機過程的方法。2、熟悉和掌握特征估計的基本方法及其MATLAB實現。二、實驗原理1、正態分布白噪聲序列的產生MATLAB提供了許多產生各種分布白噪聲序列的函數，其中產生正態分布白噪聲序列的函數為randn。函數：randn用法：x=randn(m,n)功能：產生m×n的標準正態分布隨機數矩陣。如果要產生服從分布的隨機序列，則可以由標準正態隨機序列產生。如果N(0,1)，則。2、相關函數估計MATLAB提供了函數xcorr用于自相關函數的估計。函數：xcorr用法：c=xcorr(x,y)c=xcorr(x)c=xcorr(x,y,'opition')c=xcorr(x,'opition')功能：xcorr(x,y)計算X(n)與Y(n)的互相關，xcorr(x)計算X(n)的自相關。option選項可以設定為：'biased'有偏估計。'unbiased'無偏估計。'coeff'm=0時的相關函數值歸一化為1。'none'不做歸一化處理。3、功率譜估計MATLAB函數periodogram實現了周期圖法的功率譜估計。函數：periodogram用法：[Pxx,w]=periodogram(x)[Pxx,w]=periodogram(x,window)[Pxx,w]=periodogram(x,window,nfft)[Pxx,f]=periodogram(x,window,nfft,fs)periodogram(...)功能：實現周期圖法的功率譜估計。其中：Pxx為輸出的功率譜估計值；f為頻率向量；w為歸一化的頻率向量；window代表窗函數，這種用法對數據進行了加窗，對數據加窗是為了減少功率譜估計中因為數據截斷產生的截斷誤差，下圖列出了產生常用窗函數的MATLAB函數。nfft設定FFT算法的長度；fs表示采樣頻率；三、實驗內容1、按如下模型產生一組隨機序列其中是均值為1，方差為4的正態分布白噪聲序列。估計過程的自相關函數和功率譜。實驗代碼及結果：>>y0=randn(1,500);%產生一長度為500的隨機序列>>y=1+2*y0;>>x(1)=y(1);>>n=500;>>fori=2:1:n>>x(i)=0.8*x(i-1)+y(i);%按題目要求產生隨機序列>>x(n)=0.8x(n-1)+w(n)>>end>>subplot(311);>>plot(x);>>title('x(n)');>>subplot(312);>>c=xcorr(x);%用xcorr函數求x(n)的自相關函數>>plot(c);>>title('R(n)');>>p=periodogram(x);%用periodogram函數求功率譜密度>>subplot(313);>>plot(p);>>title('S(w)');得到長度為500的樣本序列分布、自相關函數及功率譜如下：2、設信號為其中，為正態分布白噪聲序列，試在N=256和N=1024點時，分別產生隨機序列x(n)，畫出x(n)的波形并估計x(n)的相關函數和功率譜。實驗代碼及結果：（1）N=256時>>N=256;>>w=randn(1,N);%用randn函數產生一個長度為256的正態分布白噪聲序列>>n=1:1:N;>>f1=0.05;>>f2=0.12;>>x=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);%產生題目所給信號>>R=xcorr(x);%求x(n)的自相關函數>>p=periodogram(x);%求x的功率譜>>subplot(311);>>plot(x);title('x(n)');>>subplot(312);>>plot(R);title('R(n)');>>subplot(313);>>plot(p);title('S(w)');得到長度為256的樣本序列分布、自相關函數及功率譜：（2）N=1024時>>N=1024;>>w=randn(1,N);%用randn函數產生一個長度為256的正態分布白噪聲序列>>n=1:1:N;>>f1=0.05;>>f2=0.12;>>x=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);%產生題目所給信號>>R=xcorr(x);%求x(n)的自相關函數>>p=periodogram(x);%求x的功率譜>>subplot(311);>>plot(x);title('x(n)');>>subplot(312);>>plot(R);title('R(n)');>>subplot(313);>>plot(p);title('S(w)');得到長度為1024的樣本序列、序列的自相關函數、序列的功率譜如下：實驗三隨機過程通過線性系統的分析一、實驗目的1、理解和分析白噪聲通過線性系統后輸出的特性。2、學習和掌握隨機過程通過線性系統后的特性，驗證隨機過程的正態化問題。二、實驗原理1、白噪聲通過線性系統設連續線性系統的傳遞函數為H()或H(s)，輸入白噪聲的功率譜密度為SX()=N0/2，那么系統輸出的功率譜密度為SY(ω)=|H(ω)|2?N02 輸出自相關函數為RY(τ)=N04π-∞∞|H(ω)|2輸出相關系數為γyτ=輸出相關時間為τ0=0∞γy輸出平均功率為EY2(τ)=N02π0∞|上述式子表明，若輸入端是具有均勻譜的白噪聲，則輸出端隨機信號的功率譜主要由系統的幅頻特性|H(ω)|決定，不再是常數。2、等效噪聲帶寬在實際中，常常用一個理想系統等效代替實際系統的H(ω)，因此引入了等效噪聲帶寬的概念，他被定義為理想系統的帶寬。等效的原則是，理想系統與實際系統在同一白噪聲的激勵下，兩個系統的輸出平均功率相等，理想系統的增益等于實際系統的最大增益。實際系統的等效噪聲帶寬為?ωe=1|H(ω)|max20∞|H(ω)|或?ωe=12j|H(ω)|max3、線性系統輸出端隨機過程的概率分布（1）正態隨機過程通過線性系統若線性系統輸入為正態過程，則該系統輸出仍為正態過程。（2）隨機過程的正態化隨機過程的正態化指的是，非正態隨機過程通過線性系統后變換為正態過程。任意分布的白噪聲通過線性系統后輸出是服從正態分布的；寬帶噪聲通過窄帶系統，輸出近似服從正態分布。三、實驗內容1、仿真一個平均功率為1的白噪聲帶通系統，白噪聲為高斯分布，帶通系統的兩個截止頻率分別為3kHz和4kHz，估計輸出的自相關函數和功率譜密度函數。（假設采樣頻率為10kHz）實驗代碼及結果：>>Fs=10000;%采樣頻率10KHz>>x=randn(1000,1);%產生隨機序列>>figure(1);>>subplot(3,1,1);>>plot(x);gridon;%x的分布>>xlabel('t');>>subplot(3,1,2);>>x_corr=xcorr(x,'unbiased');%自相關函數>>plot(x_corr);>>gridon;>>subplot(3,1,3);>>[Pxx,w]=periodogram(x);%功率譜密度>>x_Px=Pxx;plot(x_Px);>>gridon;>>figure(2);>>subplot(2,1,1);>>xlabel('f/Hz');>>[x_pdf,x1]=ksdensity(x);%平滑密度分布函數估計>>plot(x1,x_pdf);>>gridon;>>subplot(2,1,2);>>f=(0:999)/1000*Fs;>>X=fft(x);%求DFT>>mag=abs(X);>>plot(f(1:1000/2),mag(1:1000/2));%幅頻特性>>gridon;>>xlabel('f/Hz');>>figure(3);>>subplot(3,1,1);>>[b,a]=ellip(10,0.5,50,[3000,4000]*2/Fs);%構造帶通濾波器>>[H,w]=freqz(b,a);>>plot(w*Fs/(2*pi),abs(H));%畫出濾波器頻率特性>>set(gcf,'color','white')>>xlabel('f/Hz');>>ylabel('H(w)');>>gridon;>>subplot(3,1,2);>>y=filter(b,a,x);%通過帶通濾波器>>[y_pdf,y1]=ksdensity(y);%繪出通過后的概率分布>>plot(y1,y_pdf);>>gridon;>>y_corr=xcorr(y,'unbiased');%自相關函數>>subplot(3,1,3);>>plot(y_corr);>>gridon;>>figure(4);>>Y=fft(y);%DFT>>magY=abs(Y);>>subplot(2,1,1);>>plot(f(1:1000/2),magY(1:1000/2));%幅頻特性>>gridon;>>xlabel('f/Hz');>>subplot(2,1,2);>>nfft=1024;>>index=0:round(nfft/2-1);%下標>>ky=index.*Fs./nfft;>>window=boxcar(length(y_corr));%矩形窗>>[Pyy,fy]=periodogram(y_corr,window,nfft,Fs);%求功率譜密度>>y_Py=Pyy(index+1);>>plot(ky,y_Py);>>gridon;得到：2、設白噪聲通過下圖所示的RC電路，分析輸出的統計特性。（1）試推導系統輸出的功率譜密度、相關函數、相關時間和系統的等效噪聲帶寬。（2）采用MATLAB模擬正態分布白噪聲通過上述RC電路，觀察輸入和輸出的噪聲波形以及輸出噪聲的概率密度。（3）模擬產生均勻分布的白噪聲通過上述RC電路，觀察輸入和輸出的噪聲波形以及輸出噪聲的概率密度。（4）改變RC電路的參數（電路的RC值），重做（2）和（3），與之前的結果進行比較。（1）由圖中所示電路，根據電路分析的相關知識，可推導出輸出功率譜密度為：相關函數為：相關時間為：等效噪聲帶寬為：（2）實驗代碼及結果：>>R=100;>>C=0.01;>>b=1/(R*C);>>n=1:1:500;>>h=b*exp(-n*b);>>x=randn(1,1000);>>y=conv(x,h);>>[fyy1]=ksdensity(y)>>subplot(3,1,1);>>plot(x);>>title('x(n)');>>subplot(3,1,2);>>plot(y);>>title('y(n)');>>subplot(3,1,3);>>plot(fy);>>title('fy');（3）實驗代碼及結果：>>R=100;>>C=0.01;>>b=1/(R*C);>>n=1:1:500;>>h=b*exp(-n*b);>>x=rand(1,1000);>>y=conv(x,h);>>[fyy1]=ksdensity(y);>>subplot(3,1,1);>>plot(x);>>title('x(n)');>>subplot(3,1,2);>>plot(y);>>title('y(n)');>>subplot(3,1,3);>>plot(fy);>>title('fy');實驗結果：（4）改變RC值R=200，C=0.01;正態分布：均勻分布：R=10，C=0.01;正態分布:均勻分布:由圖可得，系統相關時間與系統帶寬成反比。另外，由輸入輸出波形可以看出，正態隨機過程通過一個線性系統后，輸出仍服從正態分布。而對于任意分布的白噪聲，通過一個線性系統后，輸出也服從正態分布。實驗四窄帶隨機過程的產生及其性能測試一、實驗目的1、基于隨機過程的萊斯表達式產生窄帶隨機過程。2、掌握窄帶隨機過程的特性，包括均值（數學期望）、方差、相關函數及功率譜密度等。二、實驗原理1.窄帶隨機過程的萊斯表達式任何一個實平穩窄帶隨機過程X(t)都可以表示為上式稱為萊斯表達式，根據上式可以模擬產生窄帶隨機過程，具體過程下圖所示。2.窄帶隨機過程包絡與相位的概率密度包絡的概率密度為fA相位的概率密度為fφφ3.窄帶隨機過程包絡平方的概率密度包絡平方的概率密度為fU三、實驗內容1、按上圖所示結構框圖，基于隨機過程的萊斯表達式，用MATLAB產生一滿足條件的窄帶隨機過程。實驗代碼及結果：>>n=1:1:1000;>>h=exp(-n);>>c1=randn(1,1000);>>a=conv(c1,h);>>c2=randn(1,1000);>>b=conv(c2,h);>>fc=10000;>>x=zeros(1,1000);>>fori=1:1000>>x(i)=a(i)*cos(2*pi*fc*i)-b(i)*sin(2*pi*fc*i);>>end>>plot(x);>>title('窄帶隨機過程');實驗結果：2、畫出該隨機過程的若干次實現，觀察其形狀。代碼同上，得到圖形：3、編寫MATLAB程序計算該隨機過程的均值函數