隨機(jī)信號的仿真與分析.doc

實驗十八 隨機(jī)信號的仿真與分析一實驗?zāi)康睦糜嬎銠C(jī)仿真隨機(jī)信號，考察其數(shù)字特征，以此加深對滿足各種分布的隨機(jī)信號的理解。熟悉常用的信號處理仿真軟件平臺：matlab或c/c+語言.二實驗原理 隨機(jī)信號的產(chǎn)生和定義隨機(jī)信號是隨機(jī)變量在時間上推進(jìn)產(chǎn)生的過程量，它同時具有過程性和不確定性。定義如下：給定參量集T與概率空間（, F, P），若對于每個，都有一個定義在（, F, P）上的實隨機(jī)變量X(t)與之對應(yīng)，就稱依賴于參量t的隨機(jī)變量族為一（實）隨機(jī)過程或隨機(jī)信號。 高斯分布隨機(jī)信號統(tǒng)計分布是正態(tài)分布（高斯分布）的隨機(jī)信號為高斯分布隨機(jī)信號。高斯分布的隨機(jī)變量概率密度函數(shù)滿足下式： 均勻分布隨機(jī)信號統(tǒng)計分布是均勻分布的隨機(jī)信號為均勻分布隨機(jī)信號。均勻分布的隨機(jī)變量概率密度函數(shù)滿足下式： 正弦隨機(jī)信號給定具有某種概率分布的振幅隨機(jī)變量A、角頻率隨機(jī)變量與相位隨機(jī)變量，（具體概率分布與特性視應(yīng)用而定），以（時間）參量t建立隨機(jī)變量：。于是，相應(yīng)于某個參量域T的隨機(jī)變量族為正弦隨機(jī)信號（或稱為正弦隨機(jī)過程）。5 貝努里隨機(jī)信號貝努里隨機(jī)變量X(s)基于一個擲幣實驗（s表示基本結(jié)果事件）：1表示s為正面，0表示s不為正面；s不為正面的概率為PX(s)=1=p，s為正面的概率為PX(s)=0=q，其中p+q=1。若無休止地在t=n (n=0, 1, 2, )時刻上，獨立進(jìn)行（相同的）擲幣實驗構(gòu)成無限長的隨機(jī)變量序列：,其中與n和s都有關(guān)，應(yīng)記為X(n,s)，于是，而且有概率：其中, p+q=1。上述的隨機(jī)變量序列：通常被稱為隨機(jī)序列（或隨機(jī)過程），也被稱為（離散）隨機(jī)信號，即貝努里隨機(jī)信號。正式的定義如下：給定某個序列隨機(jī)實驗，觀測某事件B發(fā)生與否，建立事件B的指示函數(shù)，而且，序列隨機(jī)實驗間彼此統(tǒng)計獨立并有相同的概率，于是，是一個（0，1）貝努里隨機(jī)變量，相應(yīng)的隨機(jī)變量序列為（0，1）貝努里隨機(jī)序列（或稱隨機(jī)信號，有時也稱為隨機(jī)過程）。三實驗任務(wù)與要求 用matlab或c/c+語言編程并仿真。 生成滿足幾種概率分布的仿真隨機(jī)信號，自己編寫程序計算幾種概率分布的仿真隨機(jī)信號的特征。具體要求： 隨機(jī)數(shù)的產(chǎn)生與測量：產(chǎn)生的10000個泊松分布隨機(jī)數(shù)，計算它們的均值、方差與概率密度、頻譜、功率譜密度，自相關(guān)函數(shù)，并繪出函數(shù)曲線。確定泊松過程是一個馬爾可夫過程。 產(chǎn)生N(0,3)與N(2,3)的隨機(jī)數(shù)，計算它們的概率密度、頻譜、功率譜密度，自相關(guān)函數(shù)，并繪出函數(shù)曲線。確定它們是否屬于白噪聲。 統(tǒng)計分析：二維正態(tài)分布（X,Y），N(0，1；0，4；0.5)的聯(lián)合概率密度函數(shù)為 其中1=0，1=1，2=0，2=4，=0.5。求 ；） ；）。 計算二維正態(tài)分布（X,Y），N(0，1；0，4；0.5) 的概率密度、頻譜、功率譜密度，自相關(guān)函數(shù)，并繪出函數(shù)曲線。 按要求寫實驗報告。四實驗仿真與分析1隨機(jī)信號的仿真（1）高斯分布隨機(jī)信號的仿真MATLAB程序：N=1000;M=0;SIGMA2=3;SIGMA=sqrt(SIGMA2);X=M+SIGMA.*randn(N,1);t=0:N-1plot(t,X);grid on title(高斯分布的隨機(jī)信號);xlabel(t); 說明：M為均值，SIGMA2為其方差；用randn（N，1）函數(shù)產(chǎn)生一個N個樣本的正態(tài)分布的列矢量；用X=M+SIGMA.*randn(N,1)仿真產(chǎn)生N（M,SIGMA2）的高斯分布的隨機(jī)數(shù)；取t=0:N-1,用plot(t,X)對隨機(jī)數(shù)進(jìn)行仿真繪圖。編程運行后得到的高斯分布隨機(jī)信號的仿真圖如下：(2)均勻分布隨機(jī)信號的仿真程序：a=5;b=8;N=1000; x=b-(b-a).*rand(N,1);t=0:N-1;plot(t,x);title(均勻分布隨機(jī)信號 );說明：用rand（N，1）函數(shù)產(chǎn)生一個N個樣本的（0,1）均勻分布的列矢量；用x=b-(b-a).*rand(N,1)仿真產(chǎn)生（a,b）的均勻分布隨機(jī)數(shù)；取t=0:N-1，用plot(t，x)對隨機(jī)數(shù)進(jìn)行仿真繪圖；編程運行后得到的均勻分布隨機(jī)信號的仿真圖如下：（3）正弦隨機(jī)信號的仿真程序：t=0:pi/50:2*pi;for i=1:100 A=rand(); w=rand()*2; B=rand()*2*pi; y(i)=A*sin(w*t(i)+B);endfigure(1);plot(y,-r);grid on;title(y=Asin(Omegat+Theta);xlabel(t);ylabel(y);編程運行后得到的均勻分布隨機(jī)信號的仿真圖如下：（4）貝努里隨機(jī)信號的仿真：程序：p=0.5;N=50;ind=find(rand(N,1)p);z(1:N)=1;z(ind)=0;stairs(1:N,z(1:N);axis(0 N-1 -0.5 1.5);title(貝努里隨機(jī)信號 );說明：P為某事件B發(fā)生的概率；用rand（N，1）函數(shù)產(chǎn)生一個N個樣本的（0,1）均勻分布的列矢量；find(rand(N,1)p)找出滿足大于p的（0,1）均勻分布的隨機(jī)數(shù)所在的點所在的下標(biāo)，并把那些點仿真為某事件B不發(fā)生，并用z(ind)=0將其置為0；用z(1:N)=1將仿真為某事件B發(fā)生的隨機(jī)數(shù)的點置為1；用階梯函數(shù)stairs(1:N,z(1:N)對其進(jìn)行仿真繪圖；用axis(0 N -0.5 1.5)對坐標(biāo)進(jìn)行控制，是圖像更加美觀。編程運行后得到的貝努里隨機(jī)信號的仿真圖如下：2、幾種概率分布的仿真隨機(jī)信號的特征（1）泊松分布 產(chǎn)生的10000個泊松分布隨機(jī)數(shù)，計算均值和方差Matlab程序如下：y=poissrnd(1,1,10000);%產(chǎn)生泊松隨機(jī)數(shù)10000個，bosong.mjunzhi=mean(y); %求均值disp(均值); %顯示均值disp(junzhi);fangcha=var(y); %求方差disp(方差); %顯示方差 disp(fangcha);figure(1)title(泊松分布隨機(jī)數(shù)10000個);xlabel(x),ylabel(y)plot(y) 圖像如下： 概率密度函數(shù)Matlab程序如下： y=poissrnd(1,1,10000);%產(chǎn)生泊松隨機(jī)數(shù)10000個figure(1)subplot(2,1,1)f,yi=ksdensity(y); %概率密度plot(yi,f,-)grid ontitle(10000個泊松分布隨機(jī)數(shù)的概率密度估計)z=poissfit(y) %泊松分布的參數(shù)估計y=poissrnd(1,1,100); %產(chǎn)生泊松隨機(jī)數(shù)100個，將坐標(biāo)放大figure(1)subplot(2,1,2)f,yi=ksdensity(y); %概率密度plot(yi,f,-)grid ontitle(100個泊松分布隨機(jī)數(shù)的概率密度估計)命令行結(jié)果：均值 0.9993方差 1.0138程序運行后的圖像： 頻譜Matlab 程序如下：y=poissrnd(1,1,10000); %產(chǎn)生泊松隨機(jī)數(shù)10000個pingpu=fft(y) %傅里葉變換求頻譜figure(1)subplot(2,1,1)plot(abs(pingpu)title(10000個泊松分布隨機(jī)數(shù)的頻譜估計)grid onsubplot(2,1,2)y=poissrnd(1,1,100); %產(chǎn)生泊松隨機(jī)數(shù)10000個pingpu=fft(y)plot(abs(pingpu)title(100個泊松分布隨機(jī)數(shù)的頻譜估計)grid on程序運行后的圖像： 功率譜密度Matlab 程序如下：y=poissrnd(1,1,10000); %產(chǎn)生泊松隨機(jī)數(shù)10000個figure(1)window=hann(10000)periodogram(y,window,512,10000)%功率譜密度title(功率譜密度曲線)程序運行后的圖像： 自相關(guān)函數(shù)Matlab 程序如下：y=poissrnd(1,1,10000); %產(chǎn)生泊松隨機(jī)數(shù)10000個R=xcorr(y)figure(1)subplot(2,1,1)plot(R)grid ontitle(10000個泊松隨機(jī)數(shù)的自相關(guān)函數(shù)曲線) y=poissrnd(1,1,100); %產(chǎn)生泊松隨機(jī)數(shù)10000個R=xcorr(y)subplot(2,1,2)plot(R)grid ontitle(100個泊松隨機(jī)數(shù)的自相關(guān)函數(shù)曲線)程序運行后的圖像： 確定泊松過程是馬爾可夫過程證明： 馬爾可夫過程： 如果第n次的取樣試驗只與第n-1次取樣結(jié)果有關(guān)，而與小于n-1次取樣結(jié)果無關(guān)。 設(shè)Xn,n=1,2,為一隨機(jī)序列，其狀態(tài)空間E=a1,a2,aN，若對于任意的n，滿足 PXn=ain|Xn-1=ain-1,Xn-2=ain-2,X1=ai1=PXn=ain|Xn-1=ain-1 則稱Xn為馬爾可夫鏈。 泊松過程具有如下性質(zhì)：1. PX（t0）=0=1；2. 增量N(tn)-N(tn-1)相互獨立3. 增量均勻，即對任意s,t0 PN(s+t)-N(s)=n=PN(t)=n獨立增量過程X1,X2, Xn彼此獨立，并令X(t0)=X(0)=0Xn+1與所有X(ti),in獨立PX(tn+1)=xn+1X(t0)=x0,X(t1)=x1,.,X(tn)=xn= Pxn+Xn+1=xn+1X(tn)=xn=PX(tn+1)=xn+1X(tn)=xn滿足馬爾可夫鏈，所以泊松過程是一個馬爾可夫過程。（2）高斯分布AN(0,3)正態(tài)分布分布 概率密度函數(shù)Matlab 程序如下：N=1024;Fs=10000;m=0;%均值sigma2=3;%方差 sigma=sqrt(sigma2);u=randn(N,1);x=m+sigma.*u;%產(chǎn)生所需的高斯隨機(jī)數(shù)%求高斯隨機(jī)過程的概率密度 t=0:N-1;f,t=ksdensity(x)figure(1);subplot(1,2,1);plot(t,f);axis square;title(概率密度函數(shù)曲線);grid on y=normpdf(t,m,sigma);subplot(1,2,2);plot(t,y);axis square;title(理論概率密度函數(shù)曲線 );grid on 程序運行后的對應(yīng)的圖像： 頻譜Matlab 程序如下：%求高斯隨機(jī)過程的頻譜figure(2);f=(0:N-1)/N*Fs; X=fft(x);%對高斯隨機(jī)信號進(jìn)行傅里葉變換XX=abs(X) %取信號的幅度 plot(f(1:N/2),XX(1:N/2);%畫出頻譜 title(頻譜);xlabel(Frequency/Hz); grid on 程序運行后的對應(yīng)的圖像： 功率譜密度Matlab 程序如下：%利用periodogram函數(shù)計算功率譜 figure(3);nfft=1024;f=(0:N-1)/N*Fs; window=boxcar(length(R); Pxx,f=periodogram(R,window,nfft,Fs) %加窗的周期圖譜估計法plot(f,Pxx); title(功率譜密度);Xlabel(Frequency/Hz) ;grid on 程序運行后的對應(yīng)的圖像： 自相關(guān)函數(shù)Matlab 程序如下：%求高斯隨機(jī)過程的自相關(guān)函數(shù)figure(4);R,lags=xcorr(x,500,unbiased)plot(lags,R);title(自相關(guān)函數(shù));grid on 程序運行后的對應(yīng)的圖像： 分析與確定是否是白噪聲由圖可知：其自相關(guān)函數(shù)只在某點上產(chǎn)生較大的沖擊，沖擊大小約為高斯隨機(jī)過程的方差，在其他點都在零值的附近，近似為零值；在系統(tǒng)的帶內(nèi)，它的功率譜密度基本上是常數(shù)；均值近似為零；在廣義上滿足白噪聲的條件，故其屬于白噪聲。BN(2,3)正態(tài)分布分布 概率密度函數(shù)Matlab程序：N=1024;Fs=10000;m=2;%均值sigma2=3;%方差 sigma=sqrt(sigma2);u=randn(N,1);x=m+sigma.*u;%產(chǎn)生所需的高斯隨機(jī)數(shù)%求高斯隨機(jī)過程的概率密度 t=0:N-1;f,t=ksdensity(x)figure(1);subplot(1,2,1);plot(t,f);axis square;title(概率密度函數(shù)曲線);grid on y=normpdf(t,m,sigma);subplot(1,2,2);plot(t,y);axis square;title(理論概率密度函數(shù)曲線 );grid on 程序運行后的對應(yīng)的圖像： 頻譜Matlab 程序如下：%求高斯隨機(jī)過程的頻譜figure(2);f=(0:N-1)/N*Fs; X=fft(x);%對高斯隨機(jī)信號進(jìn)行傅里葉變換XX=abs(X) %取信號的幅度 plot(f(1:N/2),XX(1:N/2);%畫出頻譜 title(頻譜);xlabel(Frequency/Hz); grid on 程序運行后的對應(yīng)的圖像： 功率譜密度Matlab 程序如下：%利用periodogram函數(shù)計算功率譜 figure(3);nfft=1024;f=(0:N-1)/N*Fs; window=boxcar(length(R); Pxx,f=periodogram(R,window,nfft,Fs) %加窗的周期圖譜估計法plot(f,Pxx); title(功率譜密度);Xlabel(Frequency/Hz) ;grid on 程序運行后的對應(yīng)的圖像： 自相關(guān)函數(shù)Matlab 程序如下：%求高斯隨機(jī)過程的自相關(guān)函數(shù)figure(4);R,lags=xcorr(x,500,unbiased)plot(lags,R);title(自相關(guān)函數(shù));grid on 程序運行后的對應(yīng)的圖像： 分析與確定是否是白噪聲由圖可知：其功率譜為零，不滿足為一個常數(shù)的條件；其自相關(guān)函數(shù)在零點只在某點上產(chǎn)生較大的沖擊，但在其他點都不為零；均值不為零；故其不滿足白噪聲的條件，故其不屬于白噪聲。（3）二維正態(tài)分布相關(guān)計算 邊緣密度函數(shù)Matlab程序：syms x y; f=exp(-2*(x*x-0.5*x*y+0.25*y*y)/3)/( 2*pi*sqrt(3) ); fy=int(f,x,-inf,inf)命令窗口顯示的結(jié)果：fy = 281474976710656/6126469605667549*exp(-1/8*y2)*2(1/2)*3(1/2)*pi(1/2)由計算可知：=281474976710656/6126469605667549*2(1/2)*3(1/2)*pi(1/2)=0.1995因此，經(jīng)過化簡后得到， 的數(shù)學(xué)期望Matlab程序段：x_yf=abs(x+y)*f; Exy=int(int(x_yf,x,-inf,inf),y,-inf,inf)命令窗口顯示的結(jié)果：Exy = 1125899906842624/6126469605667549*42(1/2)*pi(1/2)化簡Exy結(jié)果值，得=0.2111 的概率密度分布Matlab程序段：P=int(int(f,x,-inf,0),y,-inf,0)命令窗口顯示的結(jié)果：P = 140737488355328/18379408817002647*pi*6(1/2)*2(1/2)+70368744177664/6126469605667549*pi*6(1/2)*8(1/2)化簡Exy結(jié)果值，得=0.33333、二維正態(tài)分布隨機(jī)信號的特征 概率密度函數(shù)Matlab程序段：mu1=0,0; Sigma2=1 1; 1 4; % 輸入均值向量和協(xié)方差矩陣X,Y=meshgrid(-4:0.1:4,-4:0.1:4);xy=X(:) Y(:); %產(chǎn)生網(wǎng)格數(shù)據(jù)p=mvnpdf(xy,mu1,Sigma2);P=reshape(p,size(X); %求取聯(lián)合概率密度figure（1）,surf(X,Y,P)title(概率密度)程序運行后的對應(yīng)的圖像： 頻譜Matlab程序段：signal=mvnrnd(mu1,Sigma2,100);signal_frequency=fft(signal);h=abs(signal_frequency);i=(0:length(signal_frequency)-1)*100/length(signal_frequency);figure（2）,plot(i,h),grid on;%繪制頻譜曲線title(幅頻響應(yīng))figure（3）,plot(i,angle(signal_frequency),grid on;title(相頻響應(yīng))程序運行后的對應(yīng)的圖像： 功率譜