數(shù)字信號處理實驗三

1、隨機信號的功率譜估計方法數(shù)字信號處理()實驗三隨機信號的功率譜估計方法班級：碩1081班姓名：呂海洋學(xué)號：31113130191 實驗?zāi)康?利用自相關(guān)函數(shù)法和周期圖法實現(xiàn)隨機信號的功率譜估計。2觀察數(shù)據(jù)長度、自相關(guān)序列長度、信噪比、窗函數(shù)、平均次數(shù)等對譜估計的分辨率、穩(wěn)定性、主瓣寬度和旁瓣效應(yīng)的影響。3學(xué)習(xí)使用FFT 提高譜估計的運算速度。4體會線性譜估計方法的優(yōu)缺點。2 實驗原理與方法2.1 自相關(guān)函數(shù)法假設(shè)我們已知隨機信號的長的自相關(guān)序列，利用自相關(guān)函數(shù)法可以得到的功率譜估計： (3-4)利用窗函數(shù)，上式又可表達為 (3-5)其中，為矩形窗函數(shù)，定義為 (3-6)因此，實際上是真正功率譜與

2、窗函數(shù)傅立葉變換的卷積。矩形窗函數(shù)不僅降低了譜估計的分辨率，而且使譜估計產(chǎn)生了旁瓣。為了降低旁瓣影響，可以采用具有較小旁瓣的窗函數(shù)，如Hamming 窗，它定義為： (3-7)這種窗函數(shù)可以有效的抑制旁瓣，但是，此時主瓣寬度增大，從而降低了譜估計的分辨率，這種主瓣和旁瓣之間的矛盾在線性譜估計方法中是無法解決的。2.2 周期圖方法假設(shè)已知隨機信號的個樣本，利用周期圖方法，信號的功率譜估計為 (3-8)利用上述方法得到的譜估計方差與信號的功率譜平方成正比，為了減小它的方差，可以將信號序列進行分段處理，然后再求各分段結(jié)果的平均，這就是平均周期圖方法，即Bartlett方法。（1）Bartlett平均

3、周期圖方法將一個隨機序列分成段，每段長度為，各段之間互不重迭，因而，可以想到，第i段的信號序列可表示為 (3-9)對于每一段的周期圖又可寫成 (3-10)于是，功率譜估計定義為 (3-11)因此，對于固定的記錄長度來講，分段數(shù)增大可使譜估計的方差減小，但是由于L的減小，相應(yīng)的功率譜主瓣增寬，譜分辨率降低，顯然，方差和分辨率也是矛盾的。除了分辨率降低以外，分段處理還會引起序列的長度有限所帶來的旁瓣效應(yīng)。為減小這種影響，最有效的辦法是給分段序列用適當?shù)拇昂瘮?shù)加權(quán)，可以得到較平滑的譜估計，當然，相應(yīng)的分辨率也有所下降。（2）平滑平均周期圖方法這時一種改進的Bartlett周期圖方法，它特別適用于FF

4、T直接計算功率譜估值。將長度為N 的平穩(wěn)隨機信號序列x(n)分成K段，每段長度為L，即L=N/K。但這里在計算周期圖之前，先用窗函數(shù)給每段序列加權(quán)，K個修正的周期圖定義為 (3-12)其中表示窗函數(shù)序列的能量， (3-13)在這種情況下，功率譜估計可按下面表達式給出： (3-14)本實驗主要是利用自相關(guān)函數(shù)法和周期圖方法對下面受噪聲干擾的正弦信號進行譜估計： (3-15)其中為正弦個數(shù),，和分別為第i個正弦信號的數(shù)字頻率、相位和幅度，隨機的分布在之間，為零均值方差等于的復(fù)高斯白噪聲。3 實驗內(nèi)容及步驟1 仔細閱讀有關(guān)線性譜估計的內(nèi)容，按照圖3.1給出的框圖編制自相關(guān)函數(shù)法譜估計程序。運行程序，

5、輸入，選擇矩形窗。觀察譜峰位置是否正確（注意：由于窗效應(yīng)可能引起譜估計的非正定）。如圖3.1所示，譜峰位置大約為，誤差范圍內(nèi)是精確的。圖3.1 M=10，矩形窗功率譜估計的譜峰位置2 觀察并記錄參數(shù)變化對譜估計性能的影響。（1） 改變，其他輸入同步驟1，觀察功率譜估計的主瓣寬度和旁瓣大小隨自相關(guān)序列長度的變化情況。圖3.2 M=5，矩形窗功率譜估計的譜峰位置從圖3.2可以看出，自相關(guān)序列長度M減小時，會引起主瓣和旁瓣寬度增大，主瓣幅度下降，分辨率下降，旁瓣震蕩減弱，旁瓣個數(shù)減少。（2） 選擇窗函數(shù)為Hamming窗，其他輸入同步驟1，觀察不同窗函數(shù)對譜估計性能影響。圖3.3 M=10，Hamm

6、ing窗功率譜估計的譜峰位置圖3.4 M=10，hamming窗和矩形窗功率譜估計對比從圖3.3和3.4可以看出，與矩形窗相比，加 Hamming 窗，主瓣寬度增加一倍，分辨率降低，但旁瓣大大得到抑制?？梢?，改變或者優(yōu)化窗的形式，也就是使得窗更加光滑，可以降低旁瓣，但這往往是以加寬主瓣為代價。（3） 改變 ，其它輸入同步驟1，觀察初始相位的變化對譜估計性能的影響。圖3.5 M=10，相位變化時矩形窗功率譜估計對比由圖3.5可以看出，改變初始相位對譜估計沒有影響。這是因為，在沒有噪聲的情況下，共軛運算就消除了初始相位對譜估計的影響。而在有噪聲的情況下，初始相位的選擇會對譜估計產(chǎn)生影響。（4） 改

7、變 ，其它輸入同步驟 1，觀察信噪比變化對譜估計性能的影響。圖3.6 M=10，加噪聲矩形窗功率譜估計對比從上圖可以看出，增加噪聲功率，降低信噪比，使得譜估計過程中旁瓣的震蕩幅度增加，但是可以看到噪聲對于信號的功率譜估計的影響不大。（5） 改變 N=10， 其它輸入同步驟 1，結(jié)合（4）的內(nèi)容，觀察數(shù)據(jù)長度及信噪比對譜估計性能的影響。圖3.7 數(shù)據(jù)長度對矩形窗功率譜估計的影響從上圖可以看出，數(shù)據(jù)長度變小，譜估計的性能下降，這是由于較短的數(shù)據(jù)長度，在求自相關(guān)的時候，m接近N-1的幾個點處，可用數(shù)據(jù)點過少，最后僅剩一個點，會導(dǎo)致求出的自相關(guān)不準確，大大影響了功率譜的估計。另一方面，噪聲方差不為零，

8、降低了信噪比，估計性能也隨之受到影響。但是，可以看到即使噪聲方差不為零，但只要數(shù)據(jù)的長度較長，仍能獲得較為滿意的功率譜估計結(jié)果，主瓣和旁瓣之間的差異明顯，旁瓣沒有影響到譜峰數(shù)目及位置的確定；而數(shù)據(jù)長度過短的時候，這一點就不能保證，可以看到，會出現(xiàn)多個峰值，甚至噪聲的譜峰高過了有用信號的譜峰。3 運行程序，輸入，選擇矩形窗，調(diào)整自相關(guān)序列長度M，使得兩個正弦頻率分量臨界分辨出來，記錄此時M值，并繪出功率譜圖。同樣，在加Hamming窗的情況下，記錄使兩個正弦頻率分量臨界分開的M值，并繪出此時的功率譜圖。（1） 矩形窗：圖3.8 M=7、8的矩形窗功率譜估計矩形窗條件下，從圖中可以看出，隨著窗長的

9、增加，頻譜的分辨率越高，兩個主瓣越明顯，旁瓣的震蕩頻率越高，可以使兩個正弦頻率分量臨界分辨出來的M 值為M=8。（2） Hamming窗圖3.9 M=9、10的Hamming窗功率譜估計Hamming窗條件下，從圖中可以看出，隨著窗長的增加，頻譜的分辨率越高，兩個主瓣越明顯，旁瓣的震蕩頻率越高，M=9的時候只有非常微弱的峰值，而可以使兩個正弦頻率分量臨界清晰分辨出來的M值為M=10。4 運行自相關(guān)函數(shù)法譜估計程序，輸入，選擇矩形窗，觀察利用自相關(guān)函數(shù)法得到白噪聲信號譜估計。改變 M=3，20，觀察M的變化對白噪聲譜估計的影響。圖3.10 M對白噪聲矩形窗功率譜估計的影響從圖3.10可以看出，窗

10、寬越寬，白噪聲功率譜震蕩越明顯，這是因為，矩形窗函數(shù)的傅里葉變換為sinc函數(shù)，在時域窗寬越寬，頻域分辨率就越高，即sinc主瓣越窄，和噪聲的功率譜做卷積后，就能夠使噪聲的每個頻點的功率譜幅值越精確，體現(xiàn)出其隨機性越明顯。5 根據(jù)圖 3.1 給出的框圖，編制周期圖法譜估計程序，子程序 FFT 可以直接調(diào)用Matlab中的函數(shù)。運行程序，選擇窗函數(shù)為矩形窗，輸入 ，觀察譜峰位置是否正確，并與步驟1結(jié)果比較。從下圖中可以看出，與利用矩形窗的自相關(guān)函數(shù)法相比，周期圖法主瓣幅度比較低，并且主瓣分辨率也相對較差，但是周期圖法具有旁瓣能量小，振蕩小的優(yōu)點，所以周期圖法更為穩(wěn)定。圖3.11 不同方法功率譜估

11、計的比較6利用周期圖法重復(fù)2，3，4 的內(nèi)容，這里相當于自相關(guān)函數(shù)法中的，觀察周期圖法譜估計和自相關(guān)函數(shù)法譜估計在分辨率和穩(wěn)定性方面的差別。（1） 改變 L=5，其它輸入同步驟 5，觀察功率譜估計的主瓣寬度和旁瓣大小隨序列分段長度的變化情況。圖3.11 不同M對周期圖法功率譜估計的影響從上圖中可以看出，分段長度的減小使得主瓣的幅度下降，并且變寬，從而使得分辨率下降，但是同時旁瓣的振蕩減小。（2） 選擇窗函數(shù)為Hamming窗，其他輸入同步驟5，觀察不同窗函數(shù)對譜估計性能影響。圖3.12 不同窗函數(shù)對周期圖法功率譜估計的影響從圖中可以看到，加hamming窗后主瓣幅度下降，寬度變寬，分辨率下降，

12、但是對旁瓣的抑制有所改進。（3） 改變，其它輸入同步驟5，觀察初始相位的變化對譜估計性能的影響。由下圖可見，初始相位的變化對譜估計性能并沒有影響。，噪聲為0時，對做FFT，取模，所以消除了相位影響。圖3.13 信號相位對周期圖法功率譜估計的影響（4） 改變 ，其它輸入同步驟5，觀察信噪比變化對譜估計性能的影響。圖3.14 信噪比對周期圖法功率譜估計的影響從圖3.14可以看出，增加噪聲功率，降低信噪比，主瓣的幅度變大，這是因為疊加上了噪聲功率。在噪聲的影響下，其峰值稍有所偏移，譜估計產(chǎn)生了一定的誤差。（5） 改變 N=10， 其它輸入同步驟 1，結(jié)合（4）的內(nèi)容，觀察數(shù)據(jù)長度及信噪比對譜估計性能

13、的影響。圖3.14 數(shù)據(jù)長度及信噪比對周期圖法功率譜估計的影響從上圖可以看出，數(shù)據(jù)長度變小，譜估計的性能主要是對譜峰位置估計有影響，而對于譜的分辨率和振蕩性影響不大，這是因為周期圖法每段數(shù)據(jù)都使用了L=10點數(shù)據(jù)做FFT，保證了每一段數(shù)據(jù)的功率譜估計都是較為準確的。（6） 運行程序，輸入，選擇矩形窗，調(diào)整L，使得兩個正弦頻率分量臨界分辨出來，記錄此時L值，并繪出功率譜圖。同樣，在加Hamming窗的情況下，記錄使兩個正弦頻率分量臨界分開的L值，并繪出此時的功率譜圖。1）矩形窗：圖3.15 L=8、9的矩形窗周期圖法功率譜估計矩形窗條件下，從圖中可以看出，隨著窗長的增加，頻譜的分辨率越高，兩個主

14、瓣越明顯，旁瓣的震蕩頻率越高，可以使兩個正弦頻率分量臨界分辨出來的L 值為L=9。2）Hamming窗：圖3.16 L=9、10的Hamming窗周期圖法功率譜估計Hamming窗條件下，從圖中可以看出，隨著窗長的增加，頻譜的分辨率越高，兩個主瓣越明顯，旁瓣的震蕩頻率越高，可以使兩個正弦頻率分量臨界清晰分辨出來的L值為L=10。（7） 輸入，選擇矩形窗，運行周期圖方法譜估計程序。改變 L=3，20，觀察L的變化對白噪聲譜估計的影響。圖3.17 L的變化對周期圖法白噪聲譜估計的影響從圖3.10可以看出，窗寬越寬，白噪聲功率譜震蕩越明顯，這是因為，矩形窗函數(shù)的傅里葉變換為sinc函數(shù)，在時域窗寬越

15、寬，頻域分辨率就越高，即sinc主瓣越窄，和噪聲的功率譜做卷積后，就能夠使噪聲的每個頻點的功率譜幅值越精確，體現(xiàn)出其隨機性越明顯。4 思考題1、證明式（34）可以表示為其中表示取實部。證明：為實信號，2、 證明：由于是實信號所以以及把帶入上式得 其中，為矩形窗。令 則 5 主要實驗結(jié)論本次實驗使用了兩種經(jīng)典的方法進行信號的功率譜估計，兩種方法各有優(yōu)缺點，對比如下。1、 自相關(guān)函數(shù)法分辨率相對較高，但是隨著相關(guān)函數(shù)序列長度的減小，主瓣幅度下降，分辨率下降，但是旁瓣振蕩減弱，旁瓣個數(shù)減少。在數(shù)據(jù)長度明顯減少時，特別是在N與M接近的時候，自相關(guān)函數(shù)的兩邊序列的樣本數(shù)很少，估計很不準確，在有噪聲干擾的

16、情形下，基本上很難有效估計出功率譜峰位置，因此，穩(wěn)定性較差。2、 周期圖法主瓣幅度比較低，并且主瓣分辨率也相對較差，與自相關(guān)函數(shù)法相似，在L值（與自相關(guān)函數(shù)法M相對應(yīng)）減小的時候，周期圖法的譜估計主瓣變寬，分辨率下降。但是周期圖法具有旁瓣能量小，振蕩小的優(yōu)點，特別是在數(shù)據(jù)長度N明顯減小的時候，即便有噪聲的影響，依然能夠較為有效地估計出功率譜峰的位置，并且不會因為噪聲的影響而錯誤的產(chǎn)生出有效信號中本來沒有的譜線，而自相關(guān)函數(shù)法則不具備這一優(yōu)點，因此相比較而言，周期圖法更有較好的穩(wěn)定性。3、 優(yōu)化加窗用的窗函數(shù)，或者改進的平滑平均周期圖法都是利用較為光滑的窗函數(shù)代替矩形窗，減小了截斷產(chǎn)生的吉布斯效

17、應(yīng)帶來的振蕩性，有效抑制旁瓣，但需要看到，抑制旁瓣的結(jié)果是以主瓣加寬分辨率下降為代價的，因此，在需要對功率譜譜線位置精確估計時，還是希望主瓣越窄越好，因此，使用窗函數(shù)時可以輔以增加自相關(guān)序列長度（或者周期圖方法中的增大L值），當然，性能的提高也意味著計算開銷的增大。6 附錄（主要源程序）6.1 自相關(guān)函數(shù)法clear allclcN=100;M=10;Ns=1; a1=1;omiga1=0.6*pi;phi1=0; a2=1;omiga2=0.8*pi;phi2=0; a=a1 a2;omiga=omiga1 omiga2;phi=phi1 phi2; sigma_2=0;% variance

18、 of the noise %generate the signals=zeros(1,N);for n=1:N for i=1:Ns s(n)=s(n)+a(i)*exp(1i*(omiga(i)*(n-1)+phi(i); endendw=sqrt(sigma_2)*randn(1,N)+1i*sqrt(sigma_2)*randn(1,N);% complex gauss white noisex=s+w;% compute the Autocorrelation function% rectangular windowc,lags=xcorr(x,M-1,unbiased) % Add

19、 rectangular window whose half length is M% Add Hamming window whose half length is MW_Haming = hamming(2*M-1);c_Hamming=c.*W_Haming.;% half of Autocorrelation functionc1=c(M:M+M-1);c1_Hamming=c_Hamming(M:M+M-1);% compute the PSDx_omiga=0:2*pi/128:2*pi*127/128;P_Rect1=zeros(1,128);P_Hamming1=zeros(1

20、,128);for kk=1:1:128 for m=-M+1:1:M-1 P_Rect1(kk)=P_Rect1(kk)+c(m+M)*exp(-1i*x_omiga(kk)*m); P_Hamming1(kk)=P_Hamming1(kk)+c_Hamming(m+M)*exp(-1i*x_omiga(kk)*m); endend% P_Rect1=2*real(fft(c1,128)-c1(1)% P_Hamming1=2*real(fft(c1_Hamming,128)-c1_Hamming(1)6.2 周期圖法clear allclc% input parameters:N=100;L=3;K=floor(N/L);Ns=0;a1=1;omiga1=0.6*pi;phi1=0; a2=1;omiga2=0.8*pi;phi2=0; a=a1 a2;omiga=omiga1 omiga2;phi=p