實驗隨機信號的功率譜估計方法

1、隨機信號的功率譜估計方法1 / 22一、 實驗目的1、 利用自相關函數(shù)法和周期圖法實現(xiàn)對隨機信號的功率譜估計 2、 觀察數(shù)據(jù)長度、自相關序列長度、信噪比、窗函數(shù)、平均次 數(shù)等譜估計的分辨率、穩(wěn)定性、主瓣寬度和旁瓣效應的影響。3、 學習使用FFT提高譜估計的運算速度。4、 體會非參數(shù)化功率譜估計方法的優(yōu)缺點。二、 實驗原理假設信號x(n)為平穩(wěn)隨機過程，其自相關函數(shù)定義為 mEx*(n)x(n+m) (3-1)其中E表示取數(shù)學期望，*表示共軛運算。根據(jù)定義，x(n)的功率譜密度與自相關序列存在下面關系： (3-2) (3-3)但是，實際中我們很難得到準確的自相關序列，只能通過隨機信號的一段樣本序

2、列來估計信號的自相關序列，進而得到信號的功率譜估計。目前，常用的線性譜估計方法有兩種，即相關函數(shù)法和周期圖方法，本實驗對這兩種方法分別予以討論。1 自相關函數(shù)法假設我們已知隨機信號x(n)的M長的自相關序列，利用自相關函數(shù)法可以得到x(n)的功率譜估計： (3-4)利用窗函數(shù)，上式又可表達為 (3-5)其中，為矩形窗函數(shù)，定義為 (3-6)因此，實際上是真正功率譜與窗函數(shù)傅立葉變換的卷積。矩形窗函數(shù)不僅降低了譜估計的分辨率，而且使譜估計產(chǎn)生了旁瓣。為了降低旁瓣影響，可以采用具有較小旁瓣的窗函數(shù)，如Hamming 窗，它定義為 (3-7)這種窗函數(shù)可以有效的抑制旁瓣，但是，此時主瓣寬度增大，從而

3、降低了譜估計的分辨率，這種主瓣和旁瓣之間的矛盾在線性譜估計方法中是無法解決的。2 周期圖方法假設已知隨機信號的N個樣本，利用周期圖方法，信號x(n)的功率譜估計為 (3-8)利用上述方法得到的譜估計方差與信號的功率譜平方成正比，為了減小它的方差，可以將信號序列進行分段處理，然后再求各分段結果的平均，這就是平均周期圖方法，即Bartlett方法。（1）Bartlett 平均周期圖方法將一個隨機序列 (0nN）分成K段，每段長度為L，各段之間互不重迭，因而N=KL，可以想到，第i段的信號序列可表示為 (3-9)對于每一段的周期圖又可寫成, (3-10)于是，功率譜估計定義為 (3-11)因此，對于

4、固定的記錄長度來講，分段數(shù)K增大可使譜估計的方差減小，但是由于L的減小，相應的功率譜主瓣增寬，譜分辨率降低，顯然，方差和分辨率也是矛盾的。除了分辨率降低以外，分段處理還會引起序列的長度有限所帶來的旁瓣效應。為減小這種影響，最有效的辦法是給分段序列用適當?shù)拇昂瘮?shù)加權，可以得到較平滑的譜估計，當然，相應的分辨率也有所下降。（2）平滑平均周期圖方法這時一種改進的Bartlett周期圖方法，它特別適用于FFT直接計算功率譜估值。將長度為N的平穩(wěn)隨機信號序列x(n)分成K段，每段長度為L，即L=N/K。但這里在計算周期圖之前，先用窗函數(shù)給每段序列加權，K個修正的周期圖定義為, (3-12)其中U表示窗函

5、數(shù)序列的能量， (3-13)在這種情況下，功率譜估計可按下面表達式給出： (3-14)本實驗主要是利用自相關函數(shù)法和周期圖方法對下面受噪聲干擾擾的正弦信號進行譜估計： (3-15)其中NS為正弦個數(shù)，i,i和ai分別為第i個正弦信號的數(shù)字頻率、相位和幅度，i隨機的分布在（0, 2）之間，w(n)為零均值方差等于w2的復高斯白噪聲。三、實驗內容和步驟1 仔細閱讀有關線性譜估計的內容，根據(jù)給出的框圖3-1編制自相關函數(shù)法譜估計的程序。運行程序，輸入N=100,M=10, 1=0.6,1=0,a1=1,w2=0選擇矩形窗。觀察譜峰位置是否正確（注意：由于窗效應可能引起譜估計的非正定）。如下圖所示，圖

6、（1） 由上圖（1）的仿真結果看出，譜峰的位置處于0.6處，故估計較精確。開始輸入?yún)?shù)：數(shù)據(jù)長度N，自相關函數(shù)個數(shù)M，平均次數(shù)K信號產(chǎn)生：輸入正弦個數(shù)Ns，每個正弦信號的數(shù)字頻率、相位和幅度，白噪聲信號的方差w2，按照公式（3.16）產(chǎn)生復正弦加白噪聲信號的N個采樣結束自相關函數(shù)法由N個x(n)估計出自相關序列(M長),并對此自相關序列加矩形窗或Hamming窗，利用公式(3.5)計算0，2）之間的128個功率譜抽樣點周期圖方法輸入 FFT 點數(shù) ，按照公式(3.11 )、(3.12 )和(3.13)計算NF點功率譜圖3-12 觀察并記錄參數(shù)變化對譜估計性能的影響。（1）改變M=5，其它輸入同

7、步驟1，觀察功率譜估計的主瓣寬度和旁瓣大小隨自相關序列長度的變化情況。 圖（2）上圖（2）為M=5是，矩形窗函數(shù)的功率譜估計的仿真圖。比較M=5與M=10時的圖形可以看出：M越小，主瓣寬度越大，分辨率越低，譜峰高度越低低。（2）選擇窗函數(shù)為Hamming窗，其它輸入同步驟1，觀察不同的窗函數(shù)對譜估計性能的影響。圖（3）比較圖（1）圖（3）可以看出，圖（3）的主瓣寬度增加了近一倍，旁瓣相應也減少了許多，同時發(fā)現(xiàn)Hamming窗的分辨率降低了。（3）改變，其它輸入同步驟1，觀察初始相位的變化對譜估計性能的影響。圖（4）比較圖（1）、（4）可以看出，說明相位變化對譜估計的性能沒有影響。（4）改變w2

8、=1，其它輸入同步驟1，觀察信噪比變化對譜估計性能的影響。（5）改變N=10，其它輸入同步驟1，結合(4)的內容，觀察數(shù)據(jù)長度及信噪比對譜估計性能的影響。 圖（5） 圖（6）如上圖（5）為w2=1時，加矩形窗普估計的仿真圖，對比圖（1）可以看出，當增大高斯白噪聲的方差時，對譜估計的峰值沒有產(chǎn)生影響，但對旁瓣的寬度和幅值有影響。如上圖（6）所示，在N較大時，改變信噪比對譜估計影響比較小，而當N較小時，改變信噪比對譜估計的影響較大。有圖（6）可以看出此時頻譜峰值已偏離了0.6，同時旁瓣的幅值增大了許多。3運行程序，輸入，選擇矩形窗，調整自相關序列長度M，使得兩個正弦頻率分量臨界分辨出來，紀錄此時的

9、M值，并繪制此時的功率譜圖。同樣，在加Hamming窗的情況下，記錄使兩個正弦頻率分量臨界分開的M值，并繪制此時的功率譜圖。 圖（7） 圖（8）圖（7）、（8）分別是M=7和M=8時，加矩形窗的普估計仿真圖，可以看出M=8時可以使兩個正弦頻率分量臨界分辨出來。圖（9） 圖（10）圖（9）、（10）分別是M=9和M=10時，加Hamming窗的普估計仿真圖，可以看出M=10時可以使兩個正弦頻率分量臨界分辨出來。4運行自相關函數(shù)法譜估計程序，輸入，選擇矩形窗，觀察利用自相關函數(shù)法得到的白噪聲信號譜估計。改變M=3，20，觀察M的變化對白噪聲譜估計的影響。如下圖（11）、（12）、（13）分別為M=

10、3、10、20時，加矩形窗的普估計的仿真圖，可以看出，M越大，譜估計的各個分量越越清晰。圖（11）圖（12） 圖（13）5. 根據(jù)框圖，編制周期圖法譜估計程序。自程序FFT可以直接調用Matlab中的函數(shù)。運行程序，輸入。選擇窗函數(shù)為矩形窗，觀察譜峰位置是否正確，并與步驟1結果比較。圖（14）對比圖（1），上圖中的譜峰位置同樣也出現(xiàn)在0.6處，譜的幅值明顯降低，主瓣寬度增加，旁瓣數(shù)量減少。6利用周期圖方法重復步驟2、3、4的內容，這里L=N/K相當于自相關函數(shù)中的M，觀察周期圖法譜估計和自相關函數(shù)法譜估計在分辨率和穩(wěn)定性方面的差別。1）、圖（15） 圖（16）如圖（15）、（16）分別為N=1

11、00，K=20時，加矩形窗和漢明窗后的普估計圖。N/K減小，周期圖法譜估計的主瓣寬度越大，分辨率降低，幅值降低，而且旁瓣數(shù)量減少。這和自相關函數(shù)法中減小M時，譜估計出現(xiàn)的規(guī)律是一致的。加漢明窗時的普估計，對比于圖（16）和圖（3）可以看出，周期圖法得到的譜的主瓣寬度加倍，旁瓣減小許多。周期圖法譜估計和自相關法譜估計出現(xiàn)的規(guī)律相似。2）、改變，其它輸入同步驟1，觀察初始相位的變化對譜估計性能的影響。圖（17）對比圖（14）以看出，改變初始相位的值對譜估計沒影響。改變w2=1，其它輸入同步驟1，觀察信噪比變化對譜估計性能的影響。圖（18）如上圖所示，增大高斯白噪聲的方差時，對譜估計的峰值沒有影響，

12、但卻使旁瓣增大。3）、改變N=50，其它輸入同步驟1，結合(4)的內容，觀察數(shù)據(jù)長度及信噪比對譜估計性能的影響。 圖（19）由圖可知，當N減小，同時增大噪聲方差時，矩形窗的周期圖譜估計的主瓣寬度增加。4）、運行程序，輸入，選擇矩形窗，調整自相關序列長度M，使得兩個正弦頻率分量臨界分辨出來，紀錄此時的M值，并繪制此時的功率譜圖。同樣，在加Hamming窗的情況下，記錄使兩個正弦頻率分量臨界分開的M值，并繪制此時的功率譜圖。固定K=10，分別取N=30和40，運行程序可以得到下圖（20）、（21）所示的加矩形窗普估計的圖形。 圖（20） 圖（21）圖（20）、（21）分別是N/K=3、4時，加矩形

13、窗的普估計仿真圖，可以看出N/K=4時可以使兩個正弦頻率分量臨界分辨出來。只不過此時譜峰的位置已不在0.6和0.8處了，普估計不準確了。5）、固定K=10，分別取N=80和90，運行程序可以得到下圖（22）、（23）所示的加漢明窗普估計的圖形。 圖（22） 圖（23） 圖（22）、（23）分別是N/K=8、9時，加漢明窗時的普估計仿真圖，可以看出N/K=9時可以使兩個正弦頻率分量臨界分辨出來。此時估計仍發(fā)生了錯誤。6）、運行自相關函數(shù)法譜估計程序，輸入，選擇矩形窗，觀察利用自相關函數(shù)法得到的白噪聲信號譜估計。改變M=3，20，觀察M的變化對白噪聲譜估計的影響。固定N=100，分別取K=50，1

14、0，5，即N/K=2、10、20，運行程序，對應的有圖（24）、（25）、（26）的仿真結果。圖（24） 圖（25） 圖（26）對比以上三圖可以看出，N/K的值越大，譜估計的各個分量就越清晰，也就越容易分辨了。四、實驗主要結論1、使用不同的窗函數(shù)譜估計的質量是不一樣的，矩形窗的主瓣較窄，分辨率較好，但方差較大，噪聲水平較高；而Hamming窗的主瓣較寬，分辨率較低，但方差較小，噪聲水平較低。2、周期圖法在分辨率和穩(wěn)定性方面是優(yōu)于自相關函數(shù)法的譜估計：增大 M(或 N/K)值，會使主瓣寬度減小，即分辨率增加；反之，M 越小，主瓣越寬，分辨率越小。另一方面，對于周期圖法，仍存在著頻率分辨率低、方差

15、性能不好的問題，原因是譜估計時對數(shù)據(jù)加窗截斷，用有限個數(shù)據(jù)或其自相關函數(shù)來估計無限個數(shù)據(jù)的功率譜，當數(shù)據(jù)很短時，這個問題更為突出。3、初始相位對譜估計沒有影響。五、思考題1. 證明：式（3-4）可以表示為：其中，Real·表示取實部。證明：w=m=-M+1M-1me-jwm =m=0M-1me-jwm+m=-M+10me-jwm-(0) =m=0M-1me-jwm+m=0M-1-mejwm-(0)又自相關函數(shù)為偶函數(shù)，則上式可化簡為： w=m=0M-1me-jwm+n=0M-1mejwm-(0)=m=0M-1m（e-jwm+ejwm）-(0)即：2、 已知實信號的 N 個樣本，可以定

16、義的估計如下：Aw=n=-N+1N-1me-jwm其中 試證明：證明：Aw=m=-N+1N-1me-jwm =m=-N+1N-11Nn=0N-m-1x*(n)xn+me-jwm令 m+n=k,上式為：=1Nk=0n-1xk n=0N-1x*ne-jwk-n =1Nk=0n-1x(k)e-jwk· n=0N-1xnejwn*即有：附錄一clearclcN=100;M=10;Ns=1;w1=0.6*pi;pha1=0;a1=1;sigma_w2=0;L=256;wn=sqrt(sigma_w2)*(randn(N,1); xn=zeros(N,1);for n=1:Nxn(n)=a1*e

17、xp(j*(w1*n+pha1)+wn(n); end r_xx=xcorr(xn,'unbiased');win_p=2*M-1; W_rect=rectwin(win_p);rect1=r_xx(N-M+1:N+M-1).*W_rect;rect2=fft(rect1,L); rect3=abs(rect2(1:L/2);figure(1),clfstem(0:L/2-1)/(L/2),rect3,'b')%歸一化 title('矩形窗功率譜估計');W_han=hamming(win_p); han1=r_xx(N-M+1:N+M-1).*

18、 W_han;han2=fft(han1,L);han3=abs(han2(1:L/2);figure(2),clfstem(0:L/2-1)/(L/2),han3,'b')title('hamming窗功率譜估計');附錄二clearclcN=100;M=10;Ns=2;w1=0.6*pi;w2=0.8*pi;pha1=0;pha2=0;a1=1;a2=1;sigma_w2=0;L=256; wn=sqrt(sigma_w2)*(randn(N,1);xn=zeros(N,1);for n=1:N xn(n)=a1*exp(j*(w1*n+pha1)+a2*e

19、xp(j*(w2*n+pha2)+wn(n); end r_xx=xcorr(xn,'unbiased');win_p=2*M-1;W_rect=rectwin(win_p);rect1=r_xx(N-M+1:N+M-1).*W_rect;rect2=fft(rect1,L); rect3=abs(rect2(1:L/2);figure(1),clfstem(0:L/2-1)/(L/2),rect3,'b')%歸一化 title('矩形窗功率譜估計');W_han=hamming(win_p); han1=r_xx(N-M+1:N+M-1).* W_han;han2=fft(ham1,L);han3=abs(ham2(1:L/2);figure(2),clfstem(0:L/2-1)/(L/2),han3,'b')title('hamming窗功率譜估計');附錄三clearclcN=100;Nf=128;K=30;Ns=1;w1=0.6*pi;pha1=0;a1=1;si