數(shù)字信號(hào)處理實(shí)驗(yàn)三

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

17、應(yīng)帶來的振蕩性，有效抑制旁瓣，但需要看到，抑制旁瓣的結(jié)果是以主瓣加寬分辨率下降為代價(jià)的，因此，在需要對(duì)功率譜譜線位置精確估計(jì)時(shí)，還是希望主瓣越窄越好，因此，使用窗函數(shù)時(shí)可以輔以增加自相關(guān)序列長(zhǎng)度（或者周期圖方法中的增大L值），當(dāng)然，性能的提高也意味著計(jì)算開銷的增大。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