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

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