隨機(jī)信號(hào)分析實(shí)驗(yàn)課件

蒙特卡洛方法：也稱(chēng)為統(tǒng)計(jì)試驗(yàn)方法，它是采用統(tǒng)計(jì)的抽樣理論來(lái)近似求解數(shù)學(xué)問(wèn)題或物理問(wèn)題，它即可以求解概率問(wèn)題，也可以求解非概率問(wèn)題，蒙特卡洛方法是系統(tǒng)模擬的重要方法。蒙特卡洛方法：1用一個(gè)例子來(lái)說(shuō)明蒙特卡洛的基本思想：用一個(gè)例子來(lái)說(shuō)明蒙特卡洛的基本思想：2蒙特卡洛模擬的基本步驟建立合適的概率模型進(jìn)行多次重復(fù)試驗(yàn)對(duì)重復(fù)試驗(yàn)結(jié)果進(jìn)行統(tǒng)計(jì)分析（估計(jì)頻率、均值等）、分析精度重復(fù)試驗(yàn)的次數(shù)稱(chēng)為蒙特卡洛仿真次數(shù)，試驗(yàn)次數(shù)越多，精度越高蒙特卡洛方法可以求解復(fù)雜系統(tǒng)的計(jì)算問(wèn)題，如雷達(dá)檢測(cè)系統(tǒng)的檢測(cè)概率蒙特卡洛模擬的基本步驟建立合適的概率模型進(jìn)行多次重復(fù)試驗(yàn)對(duì)重3二、隨機(jī)序列的產(chǎn)生1、均勻隨機(jī)數(shù)的產(chǎn)生利用MATLAB函數(shù)x=rand(m,n)x=rand(100,1)二、隨機(jī)序列的產(chǎn)生1、均勻隨機(jī)數(shù)的產(chǎn)生利用MATLAB函數(shù)x42、任意分布隨機(jī)數(shù)的產(chǎn)生反函數(shù)法變換法（1）反函數(shù)法定理：如果隨機(jī)變量X具有連續(xù)分布函數(shù)FX(x)，而r=是(0,1)上均勻分布的隨機(jī)變量，則X=Fx-1(r)由此等式,根據(jù)(0,1)隨機(jī)序列可以產(chǎn)生服從分布fX(x)的隨機(jī)序列xi2、任意分布隨機(jī)數(shù)的產(chǎn)生反函數(shù)法（1）反函數(shù)法定理：如果隨機(jī)5舉例：指數(shù)分布隨機(jī)數(shù)的產(chǎn)生或舉例：指數(shù)分布隨機(jī)數(shù)的產(chǎn)生或6%指數(shù)分布隨機(jī)數(shù)的產(chǎn)生N=200;r=rand(N,1);l=0.1;x=-log(r)/l;plot(x);%指數(shù)分布隨機(jī)數(shù)的產(chǎn)生7瑞利分布：瑞利分布：8%產(chǎn)生瑞利分布隨機(jī)數(shù)N=500;sigma=1;r=rand(N,1);x=sigma*sqrt(-2*log(r));subplot(2,1,1);plot(x);y=ksdensity(x)subplot(2,1,2);plot(y);%產(chǎn)生瑞利分布隨機(jī)數(shù)9韋泊分布雷達(dá)地雜波或海浪雜波服從該分布%產(chǎn)生韋泊分布隨機(jī)數(shù)N=500;b=1;a=1.2;r=rand(N,1);x=b*(-log(r)).^(1/a);subplot(2,1,1);plot(x);y=ksdensity(x)subplot(2,1,2);plot(y);韋泊分布雷達(dá)地雜波或海浪雜波服從該分布%產(chǎn)生韋泊分布隨機(jī)數(shù)10(2)變換法N(m,2)的正態(tài)隨機(jī)數(shù)的產(chǎn)生(2)變換法N(m,2)的正態(tài)隨機(jī)數(shù)的產(chǎn)生113MATLAB的隨機(jī)數(shù)生成函數(shù)1）獨(dú)立同分布白噪聲序列的產(chǎn)生（1）(0,1)均勻分布的白噪聲序列rand()用法：x=rand(m,n)功能：產(chǎn)生mn的均勻分布隨機(jī)數(shù)矩陣，例如，x=rand(100,1)，產(chǎn)生一個(gè)100個(gè)樣本的均勻分布白噪聲列矢量。3MATLAB的隨機(jī)數(shù)生成函數(shù)1）獨(dú)立同分布白噪聲序12（2）正態(tài)分布白噪聲序列randn()用法：x=randn(m,n)功能：產(chǎn)生mn的標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)矩陣，例如，x=randn(100,1），產(chǎn)生一個(gè)100個(gè)樣本的正態(tài)分布白噪聲列矢量。如果要產(chǎn)生服從N(,2)分布的隨機(jī)矢量，則可以通過(guò)標(biāo)準(zhǔn)正態(tài)隨機(jī)矢量來(lái)產(chǎn)生，MATLAB的語(yǔ)句為:x=+.*randn(100,1)。（2）正態(tài)分布白噪聲序列randn()13（3）韋伯分布白噪聲序列weibrnd()用法：x=weibrnd(A,B,m,n);功能：產(chǎn)生mn的韋伯分布隨機(jī)數(shù)矩陣，其中A、B是韋伯分布的兩個(gè)參數(shù)。例如，x=weibrnd(1,1.5,100,1)，產(chǎn)生一個(gè)100個(gè)樣本的韋分布白噪聲列矢量，韋伯分布參數(shù)a=1，b=1.5。其他分布的隨機(jī)數(shù)產(chǎn)生函數(shù)還有瑞利分布、伽瑪分布、指數(shù)分布等，在此不一一列舉。（3）韋伯分布白噪聲序列weibrnd()144、相關(guān)正態(tài)隨機(jī)矢量的產(chǎn)生產(chǎn)生N維正態(tài)隨機(jī)矢量，要求服從如下概率密度其中K為協(xié)方差矩陣是對(duì)稱(chēng)正定矩陣4、相關(guān)正態(tài)隨機(jī)矢量的產(chǎn)生產(chǎn)生N維正態(tài)隨機(jī)矢量，要求服從如下15基本方法是先產(chǎn)生零均值、單位方差，且各個(gè)分量相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)隨機(jī)矢量U，然后做變換X=AU+M基本方法是先產(chǎn)生零均值、單位方差，且各個(gè)分量相互獨(dú)立的標(biāo)準(zhǔn)正16其中A由協(xié)方差矩陣K確定A可以用矩陣分解函數(shù)得到Chol()其中A由協(xié)方差矩陣K確定A可以用矩陣分解函數(shù)得到Chol()175相關(guān)正態(tài)隨機(jī)序列的產(chǎn)生----已知相關(guān)函數(shù)產(chǎn)生一個(gè)正態(tài)隨機(jī)序列，要求相關(guān)函數(shù)滿(mǎn)足a<1產(chǎn)生公式5相關(guān)正態(tài)隨機(jī)序列的產(chǎn)生----已知相關(guān)函數(shù)產(chǎn)生一個(gè)正態(tài)隨18a=0.8;sigma=2;N=500;u=randn(N,1);x(1)=sigma*u(1)/sqrt(1-a^2);fori=2:Nx(i)=a*x(i-1)+sigma*u(i);endplot(x);MATLAB程序a=0.8;MATLAB程序19如果要產(chǎn)生任意形式的相關(guān)函數(shù)的相關(guān)正態(tài)隨機(jī)序列根據(jù)相關(guān)函數(shù)確定協(xié)方差矩陣對(duì)協(xié)方差矩陣進(jìn)行矩陣分解產(chǎn)生N維標(biāo)準(zhǔn)正態(tài)隨機(jī)矢量做變換X=AU+M在進(jìn)行矩陣分解時(shí)可以利用MATLAB的Cholesky矩陣分解函數(shù)chol()，利用chol()函數(shù)可以直接得到A矩陣如果要產(chǎn)生任意形式的相關(guān)函數(shù)的相關(guān)正態(tài)隨機(jī)序列根據(jù)相關(guān)函數(shù)確20隨機(jī)信號(hào)分析的MATLAB函數(shù)一、特征估計(jì)對(duì)于各態(tài)歷經(jīng)過(guò)程，我們可以通過(guò)對(duì)隨機(jī)序列的一條樣本函數(shù)來(lái)獲得該過(guò)程的統(tǒng)計(jì)特性，利用MATLAB的統(tǒng)計(jì)分析函數(shù)我們可以分析隨機(jī)序列的統(tǒng)計(jì)特性。在以下的介紹中，我們假定隨機(jī)序列X(n)和Y(n)是各態(tài)歷經(jīng)過(guò)程，他們的樣本分別為x(n)和y(n)，其中n=0,1,2,…N-1。隨機(jī)信號(hào)分析的MATLAB函數(shù)一、特征估計(jì)對(duì)于各態(tài)歷經(jīng)過(guò)程，211均值函數(shù)mean()用法：m=mean(x)功能：返回X(n)按估計(jì)的均值，其中x為樣本序列x(n)(n=1,2,…,N-1)構(gòu)成的數(shù)據(jù)矢量。2方差函數(shù)var()用法：sigma2=var(x)功能：返回X(n)按估計(jì)的方差，這一估計(jì)是無(wú)偏估計(jì)。在實(shí)際中也經(jīng)常采用下式估計(jì)方差，1均值函數(shù)mean()2方差函數(shù)var()22互相關(guān)函數(shù)的估計(jì)3互相關(guān)函數(shù)估計(jì)xcorrc=xcorr(x,y)c=xcorr(x)c=xcorr(x,y,'option')c=xcorr(x,'option')xcorr(x,y)計(jì)算X與Y的互相關(guān)，矢量X表示序列x(n)，矢量Y表示序列y(n)。xcorr(x)計(jì)算X的自相關(guān)。option選項(xiàng)是：互相關(guān)函數(shù)的估計(jì)3互相關(guān)函數(shù)估計(jì)xcorr23'biased''unbiased''coeff''none'Normalizesthesequencesotheautocorrelationsatzerolagareidentically1.0tousetheraw,unscaledcross-correlations(default)'biased'Normalizesthesequen24功率譜估計(jì)功率譜估計(jì)25隨機(jī)信號(hào)分析實(shí)驗(yàn)課件264概率密度估計(jì)概率密度的估計(jì)有兩個(gè)函數(shù)：ksdensity()，hist()ksdensity()函數(shù)直接估計(jì)隨機(jī)序列概率密度的估計(jì)，它的用法是：[f,xi]=ksdensity(x)它的功能是估計(jì)用矢量x表示的隨機(jī)序列在xi處的概率密度f(wàn)。也可以指定xi，估計(jì)對(duì)應(yīng)點(diǎn)的概率密度值，用法為：f=ksdensity(x,xi)4概率密度估計(jì)概率密度的估計(jì)有兩個(gè)函數(shù)：ksdensity27MATLAB程序如下：a=0.8;sigma=2;N=200;u=randn(N,1);x(1)=sigma*u(1)/sqrt(1-a^2);fori=2:Nx(i)=a*x(i-1)+sigma*u(i);end[f,xi]=ksdensity(x);plot(xi,f);xlabel('x');ylabel('f(x)');axis([-151500.13]);MATLAB程序如下：28直方圖hist(),他的用法為hist(y,x)，他的功能是畫(huà)出用矢量y表示的隨機(jī)序列的直方圖，參數(shù)x表示計(jì)算直方圖劃分的單元，也是用矢量表示。例產(chǎn)生一組隨機(jī)序列，并畫(huà)出他的直方圖。MATLAB程序如下：x=-2.9:0.1:2.9;y=normrnd(0,1,1000,1);hist(y,x);以上程序產(chǎn)生1000個(gè)標(biāo)準(zhǔn)正態(tài)隨機(jī)數(shù)，畫(huà)出的直方圖如圖所示。直方圖hist(),他的用法為hist(y,x)，他的功能是29隨機(jī)信號(hào)分析實(shí)驗(yàn)課件30蒙特卡洛方法：也稱(chēng)為統(tǒng)計(jì)試驗(yàn)方法，它是采用統(tǒng)計(jì)的抽樣理論來(lái)近似求解數(shù)學(xué)問(wèn)題或物理問(wèn)題，它即可以求解概率問(wèn)題，也可以求解非概率問(wèn)題，蒙特卡洛方法是系統(tǒng)模擬的重要方法。蒙特卡洛方法：31用一個(gè)例子來(lái)說(shuō)明蒙特卡洛的基本思想：用一個(gè)例子來(lái)說(shuō)明蒙特卡洛的基本思想：32蒙特卡洛模擬的基本步驟建立合適的概率模型進(jìn)行多次重復(fù)試驗(yàn)對(duì)重復(fù)試驗(yàn)結(jié)果進(jìn)行統(tǒng)計(jì)分析（估計(jì)頻率、均值等）、分析精度重復(fù)試驗(yàn)的次數(shù)稱(chēng)為蒙特卡洛仿真次數(shù)，試驗(yàn)次數(shù)越多，精度越高蒙特卡洛方法可以求解復(fù)雜系統(tǒng)的計(jì)算問(wèn)題，如雷達(dá)檢測(cè)系統(tǒng)的檢測(cè)概率蒙特卡洛模擬的基本步驟建立合適的概率模型進(jìn)行多次重復(fù)試驗(yàn)對(duì)重33二、隨機(jī)序列的產(chǎn)生1、均勻隨機(jī)數(shù)的產(chǎn)生利用MATLAB函數(shù)x=rand(m,n)x=rand(100,1)二、隨機(jī)序列的產(chǎn)生1、均勻隨機(jī)數(shù)的產(chǎn)生利用MATLAB函數(shù)x342、任意分布隨機(jī)數(shù)的產(chǎn)生反函數(shù)法變換法（1）反函數(shù)法定理：如果隨機(jī)變量X具有連續(xù)分布函數(shù)FX(x)，而r=是(0,1)上均勻分布的隨機(jī)變量，則X=Fx-1(r)由此等式,根據(jù)(0,1)隨機(jī)序列可以產(chǎn)生服從分布fX(x)的隨機(jī)序列xi2、任意分布隨機(jī)數(shù)的產(chǎn)生反函數(shù)法（1）反函數(shù)法定理：如果隨機(jī)35舉例：指數(shù)分布隨機(jī)數(shù)的產(chǎn)生或舉例：指數(shù)分布隨機(jī)數(shù)的產(chǎn)生或36%指數(shù)分布隨機(jī)數(shù)的產(chǎn)生N=200;r=rand(N,1);l=0.1;x=-log(r)/l;plot(x);%指數(shù)分布隨機(jī)數(shù)的產(chǎn)生37瑞利分布：瑞利分布：38%產(chǎn)生瑞利分布隨機(jī)數(shù)N=500;sigma=1;r=rand(N,1);x=sigma*sqrt(-2*log(r));subplot(2,1,1);plot(x);y=ksdensity(x)subplot(2,1,2);plot(y);%產(chǎn)生瑞利分布隨機(jī)數(shù)39韋泊分布雷達(dá)地雜波或海浪雜波服從該分布%產(chǎn)生韋泊分布隨機(jī)數(shù)N=500;b=1;a=1.2;r=rand(N,1);x=b*(-log(r)).^(1/a);subplot(2,1,1);plot(x);y=ksdensity(x)subplot(2,1,2);plot(y);韋泊分布雷達(dá)地雜波或海浪雜波服從該分布%產(chǎn)生韋泊分布隨機(jī)數(shù)40(2)變換法N(m,2)的正態(tài)隨機(jī)數(shù)的產(chǎn)生(2)變換法N(m,2)的正態(tài)隨機(jī)數(shù)的產(chǎn)生413MATLAB的隨機(jī)數(shù)生成函數(shù)1）獨(dú)立同分布白噪聲序列的產(chǎn)生（1）(0,1)均勻分布的白噪聲序列rand()用法：x=rand(m,n)功能：產(chǎn)生mn的均勻分布隨機(jī)數(shù)矩陣，例如，x=rand(100,1)，產(chǎn)生一個(gè)100個(gè)樣本的均勻分布白噪聲列矢量。3MATLAB的隨機(jī)數(shù)生成函數(shù)1）獨(dú)立同分布白噪聲序42（2）正態(tài)分布白噪聲序列randn()用法：x=randn(m,n)功能：產(chǎn)生mn的標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)矩陣，例如，x=randn(100,1），產(chǎn)生一個(gè)100個(gè)樣本的正態(tài)分布白噪聲列矢量。如果要產(chǎn)生服從N(,2)分布的隨機(jī)矢量，則可以通過(guò)標(biāo)準(zhǔn)正態(tài)隨機(jī)矢量來(lái)產(chǎn)生，MATLAB的語(yǔ)句為:x=+.*randn(100,1)。（2）正態(tài)分布白噪聲序列randn()43（3）韋伯分布白噪聲序列weibrnd()用法：x=weibrnd(A,B,m,n);功能：產(chǎn)生mn的韋伯分布隨機(jī)數(shù)矩陣，其中A、B是韋伯分布的兩個(gè)參數(shù)。例如，x=weibrnd(1,1.5,100,1)，產(chǎn)生一個(gè)100個(gè)樣本的韋分布白噪聲列矢量，韋伯分布參數(shù)a=1，b=1.5。其他分布的隨機(jī)數(shù)產(chǎn)生函數(shù)還有瑞利分布、伽瑪分布、指數(shù)分布等，在此不一一列舉。（3）韋伯分布白噪聲序列weibrnd()444、相關(guān)正態(tài)隨機(jī)矢量的產(chǎn)生產(chǎn)生N維正態(tài)隨機(jī)矢量，要求服從如下概率密度其中K為協(xié)方差矩陣是對(duì)稱(chēng)正定矩陣4、相關(guān)正態(tài)隨機(jī)矢量的產(chǎn)生產(chǎn)生N維正態(tài)隨機(jī)矢量，要求服從如下45基本方法是先產(chǎn)生零均值、單位方差，且各個(gè)分量相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)隨機(jī)矢量U，然后做變換X=AU+M基本方法是先產(chǎn)生零均值、單位方差，且各個(gè)分量相互獨(dú)立的標(biāo)準(zhǔn)正46其中A由協(xié)方差矩陣K確定A可以用矩陣分解函數(shù)得到Chol()其中A由協(xié)方差矩陣K確定A可以用矩陣分解函數(shù)得到Chol()475相關(guān)正態(tài)隨機(jī)序列的產(chǎn)生----已知相關(guān)函數(shù)產(chǎn)生一個(gè)正態(tài)隨機(jī)序列，要求相關(guān)函數(shù)滿(mǎn)足a<1產(chǎn)生公式5相關(guān)正態(tài)隨機(jī)序列的產(chǎn)生----已知相關(guān)函數(shù)產(chǎn)生一個(gè)正態(tài)隨48a=0.8;sigma=2;N=500;u=randn(N,1);x(1)=sigma*u(1)/sqrt(1-a^2);fori=2:Nx(i)=a*x(i-1)+sigma*u(i);endplot(x);MATLAB程序a=0.8;MATLAB程序49如果要產(chǎn)生任意形式的相關(guān)函數(shù)的相關(guān)正態(tài)隨機(jī)序列根據(jù)相關(guān)函數(shù)確定協(xié)方差矩陣對(duì)協(xié)方差矩陣進(jìn)行矩陣分解產(chǎn)生N維標(biāo)準(zhǔn)正態(tài)隨機(jī)矢量做變換X=AU+M在進(jìn)行矩陣分解時(shí)可以利用MATLAB的Cholesky矩陣分解函數(shù)chol()，利用chol()函數(shù)可以直接得到A矩陣如果要產(chǎn)生任意形式的相關(guān)函數(shù)的相關(guān)正態(tài)隨機(jī)序列根據(jù)相關(guān)函數(shù)確50隨機(jī)信號(hào)分析的MATLAB函數(shù)一、特征估計(jì)對(duì)于各態(tài)歷經(jīng)過(guò)程，我們可以通過(guò)對(duì)隨機(jī)序列的一條樣本函數(shù)來(lái)獲得該過(guò)程的統(tǒng)計(jì)特性，利用MATLAB的統(tǒng)計(jì)分析函數(shù)我們可以分析隨機(jī)序列的統(tǒng)計(jì)特性。在以下的介紹中，我們假定隨機(jī)序列X(n)和Y(n)是各態(tài)歷經(jīng)過(guò)程，他們的樣本分別為x(n)和y(n)，其中n=0,1,2,…N-1。隨機(jī)信號(hào)分析的MATLAB函數(shù)一、特征估計(jì)對(duì)于各態(tài)歷經(jīng)過(guò)程，511均值函數(shù)mean()用法：m=mean(x)功能：返回X(n)按估計(jì)的均值，其中x為樣本序列x(n)(n=1,2,…,N-1)構(gòu)成的數(shù)據(jù)矢量。2方差函數(shù)var()用法：sigma2=var(x)功能：返回X(n)按估計(jì)的方差，這一估計(jì)是無(wú)偏估計(jì)。在實(shí)際中也經(jīng)常采用下式估計(jì)方差，1均值函數(shù)mean()2方差函數(shù)var()52互相關(guān)函數(shù)的估計(jì)3互相關(guān)函數(shù)估計(jì)xcorrc=xcorr(x,y)c=xcorr(x)c=xcorr(x,y,'option')c=xcorr(x,'option')xcorr(x,y)計(jì)算X與Y的互相關(guān)，矢量X表示序列x(n)，矢量Y表示序列y(n)。xcorr(x)計(jì)算X的自相關(guān)。option選項(xiàng)是：互相關(guān)函數(shù)的估計(jì)3互相關(guān)函數(shù)估計(jì)xcorr53'biased''unbiased''coeff''none'Normalizesthesequencesotheautocorrelationsatzerolagareidentically1.0tousetheraw,unscaledcross-correlations