實驗一：隨機過程的模擬與特征估計

1、.實驗一：隨機過程的模擬與特征估計一、實驗?zāi)康牧私怆S機過程特征估計的基本概念和方法，學(xué)會運用MATLAB軟件產(chǎn)生各種隨機過程，對隨機過程的特征進行估計，并通過實驗了解不同估計方法所估計出來的結(jié)果之間的差異。二、實驗原理（ 1）高斯白噪聲的產(chǎn)生利用 MATLAB函數(shù) randn 產(chǎn)生（ 2）自相關(guān)函數(shù)的估計1 Nm 1m) x(n)對有偏估計x(n?N n 0Rx (m)?1N m 1xnmxn 對無偏估計Rx ( m)N mn 0MATLAB自帶的函數(shù)為 xcorr()，闡述 xcorr 的用法R=xcorr(x,y)或 R=xcorr(x,y,option )用來求序列 x(n) 與 y(n

2、) 的互相關(guān)函數(shù)R=xcorr(x) 或 R=xcorr(x, option ) 用來求序列 x(n) 的自相關(guān)函數(shù)option 選項是 :biased 有偏估計， unbiased 無偏估計，coeff m=0 的相關(guān)函數(shù)值歸一化為 1 none不作歸一化處理（ 3）功率譜的估計?12( )X ( )利用周期圖方法估計功率譜， GxN提示： MATLAB自帶的函數(shù)為 periodogram() ，闡述 periodogram() 的用法 ;闡述其它譜估計方法的用法。Pxx,w=periodgram(x)Pxx 為對應(yīng)頻率 w的功率譜密度值。Pxx,w=periodgram(x,window)

3、window =boxcar(n) 矩形窗（ Rectangle Window ）.window =triang(n)三角窗（ Triangular Window ）window =hanning(n) 漢寧窗（ Hanning Window）window =hamming(n)海明窗（ Hamming Window）window =blackman(n) 布拉克曼窗（ Blackman Window）window=kaiser(n,beta)愷撒窗（ Kaiser Window ）Window代表與 x 等長度的窗序列，對數(shù)據(jù)進行加窗。其它譜估計方法：相關(guān)函數(shù)法 (BT 法)該方法先由序列

4、x(n) 估計出自相關(guān)函數(shù) R(n), 然后對 R(n) 進行傅立葉變換 , 便得到 x(n)的功率譜估計。r=xccor(x);R=fft(r);Pxx=abs(R);平均周期圖法和平滑平均周期圖法對于周期圖的功率譜估計 , 當(dāng)數(shù)據(jù)長度 N太大時 , 譜曲線起伏加劇 , 若 N太小, 譜的分辨率又不好 , 因此需要改進。兩種改進的估計法是平均周期圖法和平滑平均周期圖法。 Bartlett 法:Bartlett 平均周期圖的方法是將 N 點的有限長序列 x(n) 分段求周期圖再平均。 Welch 法:Welch 法對 Bartlett 法進行了兩方面的修正 , 一是選擇適當(dāng)?shù)拇昂瘮?shù) w(n),

5、 并在周期圖計算前直接加進去 , 加窗的優(yōu)點是無論什么樣的窗函數(shù)均可使譜估計非負。 二是在分段時 , 可使各段之間有重疊 , 這樣會使方差減小。現(xiàn)代譜估計 - 非參數(shù)法MTM法 pmtm：使用正交窗口來截取獲得相互獨立的改進周期圖法功率譜估計，然后再把這些估計結(jié)果結(jié)合得到最終的估計。隨著NW的增大，窗的個數(shù)增多，會有更多的譜估計，從而譜估計的方差得到減小， 但同時帶來譜泄露的增大， 而且正的譜估計的結(jié)果將會有更大的偏差。MUSIC法 pmusic：基于矩陣特征分解的譜估計非參數(shù)方法， 它把相關(guān)數(shù)據(jù)矩陣中的信息分類，把信息分配到信號的子空間或噪聲的子空間。 它適合于普遍情況下的正弦信號參數(shù)估計的

6、方法，是多信號分類法的簡稱。特征向量法peig ：也是一種基于矩陣特征分解的譜估計非參數(shù)方法，它主要適用于混有噪.聲的正弦信號的功率譜估計，此方法利用相關(guān)矩陣的特征值來對MUSIC法公式中的求和進行加權(quán)得到的。Matlab 函數(shù)（ 4）均值的估計?1 N1()mxN n1x nMATLAB自帶的函數(shù)為 mean()（ 5）方差的估計21N1 x(n)x 2?xN n 1MATLAB自帶的函數(shù)為 var()(6) AR(1) 模型的理論自相關(guān)函數(shù)和理論功率譜對于 AR(1) 模型 X (n)aX (n 1) W ( n) ，自相關(guān)函數(shù)為 RX (m)2 a|m|，其功率a212譜為 GX( )j

7、 ) 2。(1 ae三、實驗容（帶* 為選作）.1. 相關(guān)高斯隨機序列的產(chǎn)生按如下模型產(chǎn)生一組隨機序列x(n)ax(n1)w(n) ，其中 w(n) 為均值為 1，方差為4 的正態(tài)分布白噪聲序列。（ 1）產(chǎn)生并畫出 a=0.8 和 a=0.2 的 x(n) 的波形；（ 2）估計 x(n) 的均值和方差；（ 3）估計 x(n) 的自相關(guān)函數(shù)，并畫出相關(guān)函數(shù)的圖形當(dāng) a=0.8 時程序：clc,clearall ;a=0.8;w =1+sqrt(4)* randn(1,1000);%w（ n） 均值為1 方差為 4 的正態(tài)分布白噪聲x(1)=w(1)/(sqrt(1-a2);%x(1) 初值條件f

8、orn=2:1000x(n)=a*x(n-1) +w(n);endsubplot(2,1,1);plot (x);%畫出 x （ n）的波形title('x（ n）波形圖 ' );axis(0 1000 -10 20);%設(shè)定坐標(biāo)區(qū)間mean=mean(x) %估計 x（ n）的均值var=var(x)%估計 x（ n）的方差R=xcorr(x,'coeff' ); %估計 x（ n）自相關(guān)函數(shù) , 歸一化subplot(2,1,2);plot(R);%畫出自相關(guān)函數(shù)title('x(n)自相關(guān)函數(shù)' );axis(0 2000 -0.3 1);

9、結(jié)果：估計 x（ n）的均值及方差mean =4.6889var =10.8210.X（ n）的波形及其估計自相關(guān)函數(shù)圖形：x（ n） 波 形 圖20151050-5-1001002003004005006007008009001000x(n)自 相 關(guān) 函 數(shù)0.80.60.40.20-0.20200400600800100012001400160018002000當(dāng) a=0.2 時程序：clc,clearall ;a=0.2;w =1+sqrt(4)* randn(1,1000);%w（ n）均值為1 方差為 4 的正態(tài)分布白噪聲x(1)=w(1)/(sqrt(1-a2);%x(1) 初值

10、條件forn=2:1000x(n)=a*x(n-1) +w(n);endsubplot(2,1,1);plot (x);%畫出 x （ n）的波形title('x （ n）波形圖 ' );axis(0 1000 -6 10);%設(shè)定坐標(biāo)區(qū)間mean=mean(x) %估計 x（ n）的均值var=var(x)%估計 x（ n）的方差R=xcorr(x,'coeff'); %估計 x（ n）自相關(guān)函數(shù), 歸一化subplot(2,1,2);plot(R);%畫出自相關(guān)函數(shù)title('x(n)自相關(guān)函數(shù)' );axis(0 2000 -0.3 1)

11、.結(jié)果：估計 x（ n）的均值及方差mean =1.2322var =3.9102X（ n）的波形及其估計自相關(guān)函數(shù)圖形：x（ n） 波 形 圖1050-501002003004005006007008009001000x(n) 自 相 關(guān) 函 數(shù)0.80.60.40.20-0.20200400600800100012001400160018002000對比 a=0.8 和 a=0.2 時的情況，當(dāng) a>0 時，當(dāng) a 越大， x（ n）變化幅度越大，變化越快，自相關(guān)函數(shù)下降越緩慢。2. 兩個具有不同頻率的正弦信號的識別設(shè)信號為x( n)sin(2f1n)2cos(2f2n)w(n) ，

12、 n1,2,K , N ，其中 w(n) 為零.均值正態(tài)白噪聲，方差為2 。（ 1）假定21 ，針對 f10.05 , f20.08 和 f1 0.05, f2 0.20 兩種情況，使用周期圖 periodogram()的方法估計功率譜。取 N=1000，由于本題中 f<0.5,故周期圖函數(shù)使用默認采樣頻率Fs=1Hz程序：clc,clearall;f1=0.05;f2=0.08;w =sqrt(1)* randn(1,1000);%w（ n） 均值為 0方差為 1 的正態(tài)分布白噪聲n=1:1000;x(n)=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);su

13、bplot(2,1,1);periodogram(x);title('f1=0.05,f2=0.08時功率譜 ' )f1=0.05;f2=0.20;y(n)=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);subplot(2,1,2);periodogram(y);title('f1=0.05,f2=0.20時功率譜 ' ).f1=0.05,f2=0.08 時 功 率 譜40)elpm20as/dar/B0d(ycn-20euqerf/r-40ewoP-6000.10.20.30.40.50.60.70.80.91Normalized

14、Frequency(rad/sample)f1=0.05,f2=0.20時功率譜40)elpm20as/dar/B0d(ycne-20uqerf/r-40ewoP-6000.10.20.30.40.50.60.70.80.91Normalized Frequency (rad/sample)f 與圖像上角頻率對應(yīng)關(guān)系f2Fs從兩幅圖中可以看出，在 f1和 f2 對應(yīng)的頻率上功率能量明顯較大，第一幅圖w1=0.1pi ，w2=0.16pi ，第二幅圖 w1=0.1pi ， w2=0.2pi ，且 w2 對應(yīng)能量幅度大于 w1.（ 2）假 定 f1 0.05 ,f2 0.08 ， 針 對21 和2

15、4兩種情況，用周期圖periodogram() 的方法估計功率譜程序clc,clearall ;f1=0.05;f2=0.08;w1 =sqrt(1)* randn(1,1000);%w（ n） 均值為 0 方差為 1 的正態(tài)分布白噪聲w2 =sqrt(4)* randn(1,1000);%w（ n） 均值為 0 方差為 4 的正態(tài)分布白噪聲n=1:1000;x(n)=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w1(n);subplot(2,1,1);periodogram(x);title(' 噪聲方差為1 時功率譜 ' ).y(n)=sin(2*pi*

16、f1*n)+2*cos(2*pi*f2*n)+w2(n);subplot(2,1,2);periodogram(y);title(' 噪聲方差為4 時功率譜 ' )噪聲方差為1時功率譜)30el20pmas/10dar/B0d(yc-10neuq-20erf/rew -30oP-4000.10.20.30.40.50.60.70.80.91Normalized Frequency (rad/sample)噪聲方差為4時功率譜30)elpm20as/dar10/Bd(y 0c neuq-10efr/rew-20oP-3000.10.20.30.40.50.60.70.80.91N

17、ormalized Frequency (rad/sample)對比兩幅圖像，可以看出，噪聲方差越大能量越大。2*(3)假定 f10.05 , f20.08 ,4 , 選用不同的譜估計方法進行估計，并進行比較。clc,clearall ;f1=0.05;f2=0.08;w =sqrt(4)* randn(1,1000);%w（ n）均值為 0 方差為 4 的正態(tài)分布白噪聲n=1:1000;x(n)=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);%自相關(guān)法R=xcorr(x);subplot(5,1,1);Pw=fft(R/2000);.f=2*(0:length(P

18、w)-1)/length(Pw);plot(f,10*log10(abs(Pw);axis(0 1 -50 50 );title(' 理論 x（ n）功率譜 ' );gridon;title('自相關(guān)法 ' )subplot(5,1,2);periodogram(x);xlabel(' ');ylabel(' ');title('周期圖法 ' )subplot(5,1,3);window=hann(1000);periodogram(x,window);xlabel(' ');title('

19、 漢寧周期圖法 ' )subplot(5,1,4)pwelch(x,512);xlabel(' ');ylabel(' ');title(' 平滑平均周期圖法')subplot(5,1,5)pyulear(x,999);xlabel(' ');ylabel(' ');title('Yule-Walker方法').自相關(guān)法500-500.10.20.30.40.50.60.70.80.910周期圖法500)-50e00.10.20.30.40.50.60.70.80.91lpmas漢寧周期圖法

20、/da50r/Bd(0ycneuq-50erf00.10.20.30.40.50.60.70.80.91/rewo平滑平均周期圖法P500-500.10.20.30.40.50.60.70.80.91020Yule-Walker 方 法0-200.10.20.30.40.50.60.70.80.910周期圖法比自相關(guān)方產(chǎn)生圖像更為平滑， 漢寧窗周期圖法產(chǎn)生的波形比無窗函數(shù)周期法減少了大量高次諧波，而平滑平均周期圖法比漢寧窗周期圖法產(chǎn)生波形更為平滑。從圖中看出， f1 和 f2對應(yīng)頻點能量明顯高。3. 理論值與估計值的對比分析設(shè)有 AR(1) 模型，X (n)0.8 X (n1)W ( n) ，

21、W(n)是零均值正態(tài)白噪聲，方差為4。用 MATLAB模擬產(chǎn)生 X(n) 的 500 個樣本，并估計它的均值和方差；畫出 X(n) 的理論的自相關(guān)函數(shù)和功率譜；估計 X(n) 的自相關(guān)函數(shù)和功率譜。.clc,clearall ;holdoff ;a=-0.8;w =1+sqrt(4)* randn(1,500);%w（ n） 均值為1 方差為4 的正態(tài)分布白噪聲x(1)=w(1)/(sqrt(1-a2);%x(1) 初值條件forn=2:500x(n)=a*x(n-1) +w(n);endsubplot(111)plot (x);%畫出 x （ n）的波形title('x （ n）波形

22、圖 ' );axis(0 500 -20 20);%設(shè)定坐標(biāo)區(qū)間mean=mean(x) %估計 x（ n）的均值var=var(x)%估計 x（ n）的方差fori=1:500;R(i)=4*(ai)/(1-a2);end;forj=1:999ifj<500R1(j)=R(500-j);elseifj>500R1(j)=R(j-499);endendendR2=xcorr(x);%估計 x（ n）自相關(guān)函數(shù), 歸一化figure;subplot(2,1,1);i=1:999;plot(i-500,R1(i);%畫出自相關(guān)函數(shù)title(' 理論 x(n) 自相關(guān)函

23、數(shù) ' );subplot(2,1,2);plot(R2);%畫出自相關(guān)函數(shù)title(' 估計 x(n) 自相關(guān)函數(shù) ' );figure;subplot(2,1,1);Pw=fft(R/50000);f=2*(0:length(Pw)-1)/length(Pw);plot(f,10*log10(abs(Pw);title(' 理論 x（ n）功率譜 ' );subplot(2,1,2);periodogram(x,'twosided');title(' 估計 x（ n）功率譜 ' );.axis(0 2 -40 40

24、)mean =0.5657var =11.5339x（ n） 波 形 圖20151050-5-10-15-20050100150200250300350400450500.理論 x(n) 自 相 關(guān) 函 數(shù)1050-5-10-400-300-200-1000100200300400500-500估計 x(n) 自 相 關(guān) 函 數(shù)6000400020000-2000-4000-600001002003004005006007008009001000理 論 x （ n） 功 率 譜-30-35-40-450.20.40.60.811.21.41.61.820)el估 計 x （ n） 功 率 譜p

25、m40as/dar20/Bd(y0cneuq-20erf/re-40w0.20.40.60.811.21.41.61.82o0PNormalized Frequency (rad/sample)可以看出估計的自相關(guān)函數(shù)與理論值大致吻合， 但高頻上有些噪聲誤差。 但隨著樣本數(shù)的增多，估計曲線接近理論值。.若采樣點數(shù)為5000理 論 x(n) 自 相 關(guān) 函 數(shù)10.50-0.5-1-4900-4800-4700-4600-4500-4400-4300-4200-4100-4000-5000x 104估 計 x(n) 自 相 關(guān) 函 數(shù)1050-5100020003000400050006000700080009000100000理 論 x（ n） 功 率 譜-40-45-50-550.20.40.60.811.21.41.61.820)估 計 x（ n） 功 率 譜el40pmas/d20ar/Bd(0ycneuq-20erf/rew-40o0.20.40.60.811.21.41.61.82P0Normalized Frequency (rad/sample)4. 隨機信號通過線性系統(tǒng)分析考慮圖示系統(tǒng).wnxn+0.9-0.1z-1z-1假定 w 為正態(tài)分布的隨機序列和均勻分布的