北理工隨機(jī)信號分析實(shí)驗(yàn)

實(shí)驗(yàn)一隨機(jī)序列的產(chǎn)生及數(shù)字特征估計(jì)實(shí)驗(yàn)?zāi)康?、學(xué)習(xí)和掌握隨機(jī)數(shù)的產(chǎn)生方法。2、實(shí)現(xiàn)隨機(jī)序列的數(shù)字特征估計(jì)。實(shí)驗(yàn)原理1、隨機(jī)數(shù)的產(chǎn)生隨機(jī)數(shù)指的是各種不同分布隨機(jī)變量的抽樣序列（樣本值序列）。進(jìn)行隨機(jī)信號仿真分析時(shí)，需要模擬產(chǎn)生各種分布的隨機(jī)數(shù)。在計(jì)算機(jī)仿真時(shí)，通常利用數(shù)學(xué)方法產(chǎn)生隨機(jī)數(shù)，這種隨機(jī)數(shù)稱為偽隨機(jī)數(shù)。偽隨機(jī)數(shù)是按照一定的計(jì)算公式產(chǎn)生的，這個(gè)公式稱為隨機(jī)數(shù)發(fā)生器。偽隨機(jī)數(shù)本質(zhì)上不是隨機(jī)的，而且存在周期性，但是如果計(jì)算公式選擇適當(dāng)，所產(chǎn)生的數(shù)據(jù)看似隨機(jī)的，與真正的隨機(jī)數(shù)具有相近的統(tǒng)計(jì)特性，可以作為隨機(jī)數(shù)使用。(0，1)均勻分布隨機(jī)數(shù)是最最基本、最簡單的隨機(jī)數(shù)。(0，1)均勻分布指的是在[0，1]區(qū)間上的均勻分布，即U(0，1)。實(shí)際應(yīng)用中有許多現(xiàn)成的隨機(jī)數(shù)發(fā)生器可以用于產(chǎn)生(0，1)均勻分布隨機(jī)數(shù)，通常采用的方法為線性同余法，公式如下：序列為產(chǎn)生的(0，1)均勻分布隨機(jī)數(shù)。下面給出了上式的3組常用參數(shù)：1、，周期；2、（IBM隨機(jī)數(shù)發(fā)生器）周期；3、（ran0）周期；由均勻分布隨機(jī)數(shù)，可以利用反函數(shù)構(gòu)造出任意分布的隨機(jī)數(shù)。定理1.1若隨機(jī)變量X具有連續(xù)分布函數(shù)FX(x)，而R為(0,1)均勻分布隨機(jī)變量，則有由這一定理可知，分布函數(shù)為FX(x)的隨機(jī)數(shù)可以由(0,1)均勻分布隨機(jī)數(shù)按上式進(jìn)行變換得到。2、MATLAB中產(chǎn)生隨機(jī)序列的函數(shù)（1）(0，1)均勻分布的隨機(jī)序列函數(shù)：rand用法：x=rand(m,n)功能：產(chǎn)生m×n的均勻分布隨機(jī)數(shù)矩陣。（2）正態(tài)分布的隨機(jī)序列函數(shù)：randn用法：x=randn(m,n)功能：產(chǎn)生m×n的標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)矩陣。如果要產(chǎn)生服從分布的隨機(jī)序列，則可以由標(biāo)準(zhǔn)正態(tài)隨機(jī)序列產(chǎn)生。（3）其他分布的隨機(jī)序列MATLAB上還提供了其他多種分布的隨機(jī)數(shù)的產(chǎn)生函數(shù)，下表列出了部分函數(shù)。MATLAB中產(chǎn)生隨機(jī)數(shù)的一些函數(shù)3、隨機(jī)序列的數(shù)字特征估計(jì)對于遍歷過程，可以通過隨機(jī)序列的一條樣本函數(shù)來獲得該過程的統(tǒng)計(jì)特性。這里我們假定隨機(jī)序列X(n)為遍歷過程，樣本函數(shù)為x(n)，其中n=0,1,2，…，N-1。那么，X(n)的均值、方差和自相關(guān)函數(shù)的估計(jì)為利用MATLAB的統(tǒng)計(jì)分析函數(shù)可以分析隨機(jī)序列的數(shù)字特征。（1）均值函數(shù)函數(shù)：mean用法：m=mean(x)功能：返回按上面第一式估計(jì)X(n)的均值，其中x為樣本序列x(n)。（2）方差函數(shù)函數(shù)：var用法：sigma2=var(x)功能：返回按上面第二式估計(jì)X(n)的方差，其中x為樣本序列x(n)，這一估計(jì)為無偏估計(jì)。（3）互相關(guān)函數(shù)函數(shù)：xcorr用法：c=xcorr(x,y)c=xcorr(x)c=xcorr(x,y,'opition')c=xcorr(x,'opition')功能：xcorr(x,y)計(jì)算X(n)與Y(n)的互相關(guān)，xcorr(x)計(jì)算X(n)的自相關(guān)。option選項(xiàng)可以設(shè)定為：'biased'有偏估計(jì)，即'unbiased'無偏估計(jì)，即按上面第三式估計(jì)。'coeff'm=0時(shí)的相關(guān)函數(shù)值歸一化為1。'none'不做歸一化處理。實(shí)驗(yàn)內(nèi)容1、采用線性同余法產(chǎn)生均勻分布隨機(jī)數(shù)1000個(gè)，計(jì)算該序列均值和方差與理論值之間的誤差大小。改變樣本個(gè)數(shù)重新計(jì)算。代碼及結(jié)果：>>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可以看出，樣本值取得越大，實(shí)際值與理論值越接近，相差越小。2、參數(shù)為的指數(shù)分布的分布函數(shù)為利用反函數(shù)法產(chǎn)生參數(shù)為0.5的指數(shù)分布隨機(jī)數(shù)1000個(gè)，測試其方差和相關(guān)函數(shù)。實(shí)驗(yàn)代碼及結(jié)果：>>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，實(shí)際值為4.1201，因?yàn)槿觽€(gè)數(shù)有限，導(dǎo)致存在一定偏差。但大體相近。3、產(chǎn)生一組N(1,4)分布的高斯隨機(jī)數(shù)（1000個(gè)樣本），估計(jì)該序列的均值、方差和相關(guān)函數(shù)。實(shí)驗(yàn)代碼及結(jié)果：>>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。而在實(shí)驗(yàn)中得到的均值為0.9937，方差為3.8938。考慮到取樣點(diǎn)有限，誤差可以接受，理論值和實(shí)驗(yàn)值基本相同。實(shí)驗(yàn)二隨機(jī)過程的模擬與數(shù)字特征一、實(shí)驗(yàn)?zāi)康?、學(xué)習(xí)利用MATLAB模擬產(chǎn)生隨機(jī)過程的方法。2、熟悉和掌握特征估計(jì)的基本方法及其MATLAB實(shí)現(xiàn)。二、實(shí)驗(yàn)原理1、正態(tài)分布白噪聲序列的產(chǎn)生MATLAB提供了許多產(chǎn)生各種分布白噪聲序列的函數(shù)，其中產(chǎn)生正態(tài)分布白噪聲序列的函數(shù)為randn。函數(shù)：randn用法：x=randn(m,n)功能：產(chǎn)生m×n的標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)矩陣。如果要產(chǎn)生服從分布的隨機(jī)序列，則可以由標(biāo)準(zhǔn)正態(tài)隨機(jī)序列產(chǎn)生。如果N(0,1)，則。2、相關(guān)函數(shù)估計(jì)MATLAB提供了函數(shù)xcorr用于自相關(guān)函數(shù)的估計(jì)。函數(shù)：xcorr用法：c=xcorr(x,y)c=xcorr(x)c=xcorr(x,y,'opition')c=xcorr(x,'opition')功能：xcorr(x,y)計(jì)算X(n)與Y(n)的互相關(guān)，xcorr(x)計(jì)算X(n)的自相關(guān)。option選項(xiàng)可以設(shè)定為：'biased'有偏估計(jì)。'unbiased'無偏估計(jì)。'coeff'm=0時(shí)的相關(guān)函數(shù)值歸一化為1。'none'不做歸一化處理。3、功率譜估計(jì)MATLAB函數(shù)periodogram實(shí)現(xiàn)了周期圖法的功率譜估計(jì)。函數(shù)：periodogram用法：[Pxx,w]=periodogram(x)[Pxx,w]=periodogram(x,window)[Pxx,w]=periodogram(x,window,nfft)[Pxx,f]=periodogram(x,window,nfft,fs)periodogram(...)功能：實(shí)現(xiàn)周期圖法的功率譜估計(jì)。其中：Pxx為輸出的功率譜估計(jì)值；f為頻率向量；w為歸一化的頻率向量；window代表窗函數(shù)，這種用法對數(shù)據(jù)進(jìn)行了加窗，對數(shù)據(jù)加窗是為了減少功率譜估計(jì)中因?yàn)閿?shù)據(jù)截?cái)喈a(chǎn)生的截?cái)嗾`差，下圖列出了產(chǎn)生常用窗函數(shù)的MATLAB函數(shù)。nfft設(shè)定FFT算法的長度；fs表示采樣頻率；三、實(shí)驗(yàn)內(nèi)容1、按如下模型產(chǎn)生一組隨機(jī)序列其中是均值為1，方差為4的正態(tài)分布白噪聲序列。估計(jì)過程的自相關(guān)函數(shù)和功率譜。實(shí)驗(yàn)代碼及結(jié)果：>>y0=randn(1,500);%產(chǎn)生一長度為500的隨機(jī)序列>>y=1+2*y0;>>x(1)=y(1);>>n=500;>>fori=2:1:n>>x(i)=0.8*x(i-1)+y(i);%按題目要求產(chǎn)生隨機(jī)序列>>x(n)=0.8x(n-1)+w(n)>>end>>subplot(311);>>plot(x);>>title('x(n)');>>subplot(312);>>c=xcorr(x);%用xcorr函數(shù)求x(n)的自相關(guān)函數(shù)>>plot(c);>>title('R(n)');>>p=periodogram(x);%用periodogram函數(shù)求功率譜密度>>subplot(313);>>plot(p);>>title('S(w)');得到長度為500的樣本序列分布、自相關(guān)函數(shù)及功率譜如下：2、設(shè)信號為其中，為正態(tài)分布白噪聲序列，試在N=256和N=1024點(diǎn)時(shí)，分別產(chǎn)生隨機(jī)序列x(n)，畫出x(n)的波形并估計(jì)x(n)的相關(guān)函數(shù)和功率譜。實(shí)驗(yàn)代碼及結(jié)果：（1）N=256時(shí)>>N=256;>>w=randn(1,N);%用randn函數(shù)產(chǎn)生一個(gè)長度為256的正態(tài)分布白噪聲序列>>n=1:1:N;>>f1=0.05;>>f2=0.12;>>x=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);%產(chǎn)生題目所給信號>>R=xcorr(x);%求x(n)的自相關(guān)函數(shù)>>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的樣本序列分布、自相關(guān)函數(shù)及功率譜：（2）N=1024時(shí)>>N=1024;>>w=randn(1,N);%用randn函數(shù)產(chǎn)生一個(gè)長度為256的正態(tài)分布白噪聲序列>>n=1:1:N;>>f1=0.05;>>f2=0.12;>>x=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);%產(chǎn)生題目所給信號>>R=xcorr(x);%求x(n)的自相關(guān)函數(shù)>>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的樣本序列、序列的自相關(guān)函數(shù)、序列的功率譜如下：實(shí)驗(yàn)三隨機(jī)過程通過線性系統(tǒng)的分析一、實(shí)驗(yàn)?zāi)康?、理解和分析白噪聲通過線性系統(tǒng)后輸出的特性。2、學(xué)習(xí)和掌握隨機(jī)過程通過線性系統(tǒng)后的特性，驗(yàn)證隨機(jī)過程的正態(tài)化問題。二、實(shí)驗(yàn)原理1、白噪聲通過線性系統(tǒng)設(shè)連續(xù)線性系統(tǒng)的傳遞函數(shù)為H()或H(s)，輸入白噪聲的功率譜密度為SX()=N0/2，那么系統(tǒng)輸出的功率譜密度為SY(ω)=|H(ω)|2?N02 輸出自相關(guān)函數(shù)為RY(τ)=N04π-∞∞|H(ω)|2輸出相關(guān)系數(shù)為γyτ=輸出相關(guān)時(shí)間為τ0=0∞γy輸出平均功率為EY2(τ)=N02π0∞|上述式子表明，若輸入端是具有均勻譜的白噪聲，則輸出端隨機(jī)信號的功率譜主要由系統(tǒng)的幅頻特性|H(ω)|決定，不再是常數(shù)。2、等效噪聲帶寬在實(shí)際中，常常用一個(gè)理想系統(tǒng)等效代替實(shí)際系統(tǒng)的H(ω)，因此引入了等效噪聲帶寬的概念，他被定義為理想系統(tǒng)的帶寬。等效的原則是，理想系統(tǒng)與實(shí)際系統(tǒng)在同一白噪聲的激勵(lì)下，兩個(gè)系統(tǒng)的輸出平均功率相等，理想系統(tǒng)的增益等于實(shí)際系統(tǒng)的最大增益。實(shí)際系統(tǒng)的等效噪聲帶寬為?ωe=1|H(ω)|max20∞|H(ω)|或?ωe=12j|H(ω)|max3、線性系統(tǒng)輸出端隨機(jī)過程的概率分布（1）正態(tài)隨機(jī)過程通過線性系統(tǒng)若線性系統(tǒng)輸入為正態(tài)過程，則該系統(tǒng)輸出仍為正態(tài)過程。（2）隨機(jī)過程的正態(tài)化隨機(jī)過程的正態(tài)化指的是，非正態(tài)隨機(jī)過程通過線性系統(tǒng)后變換為正態(tài)過程。任意分布的白噪聲通過線性系統(tǒng)后輸出是服從正態(tài)分布的；寬帶噪聲通過窄帶系統(tǒng)，輸出近似服從正態(tài)分布。三、實(shí)驗(yàn)內(nèi)容1、仿真一個(gè)平均功率為1的白噪聲帶通系統(tǒng)，白噪聲為高斯分布，帶通系統(tǒng)的兩個(gè)截止頻率分別為3kHz和4kHz，估計(jì)輸出的自相關(guān)函數(shù)和功率譜密度函數(shù)。（假設(shè)采樣頻率為10kHz）實(shí)驗(yàn)代碼及結(jié)果：>>Fs=10000;%采樣頻率10KHz>>x=randn(1000,1);%產(chǎn)生隨機(jī)序列>>figure(1);>>subplot(3,1,1);>>plot(x);gridon;%x的分布>>xlabel('t');>>subplot(3,1,2);>>x_corr=xcorr(x,'unbiased');%自相關(guān)函數(shù)>>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);%平滑密度分布函數(shù)估計(jì)>>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);%構(gòu)造帶通濾波器>>[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');%自相關(guān)函數(shù)>>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);%下標(biāo)>>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、設(shè)白噪聲通過下圖所示的RC電路，分析輸出的統(tǒng)計(jì)特性。（1）試推導(dǎo)系統(tǒng)輸出的功率譜密度、相關(guān)函數(shù)、相關(guān)時(shí)間和系統(tǒng)的等效噪聲帶寬。（2）采用MATLAB模擬正態(tài)分布白噪聲通過上述RC電路，觀察輸入和輸出的噪聲波形以及輸出噪聲的概率密度。（3）模擬產(chǎn)生均勻分布的白噪聲通過上述RC電路，觀察輸入和輸出的噪聲波形以及輸出噪聲的概率密度。（4）改變RC電路的參數(shù)（電路的RC值），重做（2）和（3），與之前的結(jié)果進(jìn)行比較。（1）由圖中所示電路，根據(jù)電路分析的相關(guān)知識(shí)，可推導(dǎo)出輸出功率譜密度為：相關(guān)函數(shù)為：相關(guān)時(shí)間為：等效噪聲帶寬為：（2）實(shí)驗(yàn)代碼及結(jié)果：>>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）實(shí)驗(yàn)代碼及結(jié)果：>>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');實(shí)驗(yàn)結(jié)果：（4）改變RC值R=200，C=0.01;正態(tài)分布：均勻分布：R=10，C=0.01;正態(tài)分布:均勻分布:由圖可得，系統(tǒng)相關(guān)時(shí)間與系統(tǒng)帶寬成反比。另外，由輸入輸出波形可以看出，正態(tài)隨機(jī)過程通過一個(gè)線性系統(tǒng)后，輸出仍服從正態(tài)分布。而對于任意分布的白噪聲，通過一個(gè)線性系統(tǒng)后，輸出也服從正態(tài)分布。實(shí)驗(yàn)四窄帶隨機(jī)過程的產(chǎn)生及其性能測試一、實(shí)驗(yàn)?zāi)康?、基于隨機(jī)過程的萊斯表達(dá)式產(chǎn)生窄帶隨機(jī)過程。2、掌握窄帶隨機(jī)過程的特性，包括均值（數(shù)學(xué)期望）、方差、相關(guān)函數(shù)及功率譜密度等。二、實(shí)驗(yàn)原理1.窄帶隨機(jī)過程的萊斯表達(dá)式任何一個(gè)實(shí)平穩(wěn)窄帶隨機(jī)過程X(t)都可以表示為上式稱為萊斯表達(dá)式，根據(jù)上式可以模擬產(chǎn)生窄帶隨機(jī)過程，具體過程下圖所示。2.窄帶隨機(jī)過程包絡(luò)與相位的概率密度包絡(luò)的概率密度為fA相位的概率密度為fφφ3.窄帶隨機(jī)過程包絡(luò)平方的概率密度包絡(luò)平方的概率密度為fU三、實(shí)驗(yàn)內(nèi)容1、按上圖所示結(jié)構(gòu)框圖，基于隨機(jī)過程的萊斯表達(dá)式，用MATLAB產(chǎn)生一滿足條件的窄帶隨機(jī)過程。實(shí)驗(yàn)代碼及結(jié)果：>>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('窄帶隨機(jī)過程');實(shí)驗(yàn)結(jié)果：2、畫出該隨機(jī)過程的若干次實(shí)現(xiàn)，觀察其形狀。代碼同上，得到圖形：3、編寫MATLAB程序計(jì)算該隨機(jī)過程的均值函數(shù)