功率譜估計：BT ，週期圖，Bartlett ，AR ，MVDR，APES，MUSIC

二：

週期圖法

經典譜估計中的週期圖法是用得較多且最具代表性的方法。我們先複習基本的週期圖法，接著針對它譜解析度比較低等缺點，利用 Matlab 實現了幾種改進的週期圖法。

1、基本的週期圖法

基本的週期圖法可以提高計算效率，不需要計算自相關函式，但譜解析度較低。基本週期圖的基本原理是對觀測到的資料直接進行

傅立葉變換

，然後取模的平方就是功率譜。取平穩隨機訊號 x（n）的有限個觀察點 x（0）、x（1）、…x（n-1），則傅立葉變換

matlab實現程式碼如下

clc；

clear；

t = 0：0。001：0。6；

y = sin（2*pi*50*t）+sin（2*pi*120*t）；

x = y + randn（size（t））；

N=length（x）；%觀察資料點數

nfft=1024；%

傅立葉變換

點數

Xk=fft（x，nfft）；

Pxx=abs（Xk）。^2/N；

index=0：round（nfft/2-1）；

k=index*N/nfft；% 頻率（Hz）

plot_Pxx=10*log10（Pxx（index+1））；%功率譜 dB

plot（k，plot_Pxx）；

xlabel（‘Hz’）；

ylabel（‘功率譜 dB’）；

title（‘用基本週期圖法（直接計算）計算功率譜’）；

grid on；

2、Welch 法

P。O。

韋爾奇

提出一種把加窗處理與平均處理結合起來的方法。先把分段的資料乘以窗函式（進行加窗處理），分別計算其週期圖，然後進行平均，得到的功率譜值為：

為了得到較好的功率譜估值，加窗和平均處理均應兼顧減小隨機起伏和保證有足夠的譜解析度兩個方面。

matlab程式

clc；

clear；

t = 0：0。001：0。6；

y = sin（2*pi*50*t）+sin（2*pi*120*t）；%訊號

x = y+randn（size（t））；

N=length（x）；%觀察資料點數

n=1：1：N；

%引數設定

Fs=1000；%每百年取樣 100 次

nfft=512；%傅立葉變換點數

%window=hamming（25）；

window=boxcar（100）；

noverlap=20

range=‘half’；

%計算序列的 psd

［Pxx，f］=pwelch（x，window，noverlap，nfft，Fs，range）；

plot_Pxx=10*log10（Pxx）；%功率譜 dB

plot（f，plot_Pxx）

xlabel（‘Hz’）；

ylabel（‘功率譜 dB’）；

title（‘用 welch 法計算的功率譜’）；

grid on；

3、

多視窗的週期圖法

由於普通的週期圖功率譜估計只利用單一視窗，因此在序列始端和末端均會丟失相關資訊。多視窗法利用多個

正交視窗

獲得各自獨立的近似功率譜估計，然後綜合這些估計，最終得到平穩訊號序列的功率譜估計。

Matlab 程式如下：

%訊號

t = 0：0。001：0。6；

y = sin（2*pi*50*t）+sin（2*pi*120*t）；

x = y+randn（size（t））；

N=length（x）；%觀察資料點數

n=1：1：N；

%引數設定

Fs=1000；%每百年取樣 100 次

nfft=512；%傅立葉變換點數

%window=hamming（25）；

window=boxcar（100）；

noverlap=20

range=‘half’；

%計算序列的 psd

nw=3。5；p=0。99；

［Pxx，Pxxc，f］=pmtm（x，nw，nfft，Fs，p）；

plot_Pxx=10*log10（Pxx）；%功率譜 dB

plot（f，plot_Pxx）

xlabel（‘Hz’）；

ylabel（‘功率譜 dB’）；

title（‘用多視窗法計算的功率譜’）；

grid on；

總結

透過實驗可以看出，基本的週期圖法譜解析度較低，會有很大的失真。 Welch 法可以減小隨機起伏，收斂性較好，曲線平滑，估計的結果方差較小，但是功率譜主瓣較寬，解析度低。但是，如果訊號序列不是足夠長，由於每段序列長度變短，功率譜估值對不同頻率成分的分辨能力也隨之下降。多視窗的週期圖法功率譜比 Welch 法的功率譜譜峰窄一些，譜峰高度也有所增大，因此多視窗的週期圖法的功率譜表現得更加準確。普通週期圖法的功率譜估計利用單一視窗，所以在序列的始端和末端均會丟失相關資訊，而多視窗的週期圖法以增加視窗來充分利用這些相關資訊。