功率譜估計

自然界出現的許多現象都可以在統計平均意義上很好的表現出來。例如，氣象學中的氣溫與氣壓的變動等，均可以以統計的方式表示為隨機過程。在電阻器和電子設備中生成的熱噪音電壓，也是被抽象為隨機過程模型的物理訊號的例子。由於這些訊號為隨機訊號，我們必須採用一種統計觀點來處理隨機訊號的平均特徵。特別的是隨機訊號的自相關函數很適合用於代表時域中的隨機訊號，並且自相關函數的傅立葉轉換可生成功率譜密度，也可提供時域到頻域的轉換。

基於有限長訊號觀察的功率譜估計
對於有限時間長度的訊號，數據序列的有限記綠長度是對功率譜估計(power spectrum estimation)的主要限制。當處理統計平穩訊號時，數據記錄越長，可從數據提取的訊號估計就越好。另一方面若訊號統計是非平穩的，我們不能選擇任意長度記綠對譜進行估計，在這種情況下我們選擇的數據記綠長度是由訊號統計上的時變速度決定的。最後我們要選擇盡可能短且能解析不同訊號分量譜特徵的數據記綠，所表現的這些訊號分量具有相近空間譜。對基於有限長度數據記綠的古典功率譜估計方法，所面臨的問題之一是我們試圖要估計出的譜會有失真。這一問題無論在確定性訊號的譜計算方面還是在隨機訊號的功率譜估計方面都會遇到，既然很容易觀察到有限長度數據記錄對確定性訊號的效應，我們就先考察確定性訊號的情況然後再討論隨機訊號及其功率譜估計。

能量密度譜計算
考慮有限長度數據序列確定性訊號的譜計算，可參見:功率譜密度。

隨機訊號的自相關和功率譜估計:週期圖
有限能量訊號可進行傅立葉轉換，並在譜域用它們的能量密度譜來表現。另一方面代表為平穩隨機過程的重要類型訊號不具有有限能量，因此不能進行傅立葉轉換。這類訊號具有有限平均功率，因此表現為功率譜密度。如果x(t)是一個平穩隨機過程，它的自相關函數是:
:\gamma_{xx}(\tau) = E(x^*(t)x(t+\tau))

其中E代表統計平均。於是由維納-辛欽定理(Wiener–Khinchin theorem)，平穩隨機過程的功率譜是自相關函數的傅立葉轉換，即:
:\gamma_{xx}(F) = \int_{-\infty}^{\infty} \gamma_{xx}(\tau)e^{-j2\pi F\tau}\mathrm{d}t

實際上我們處理隨機過程的單個實現並從中估計該過程的功率密度譜。由於不知到真實的自相關係數\gamma_{xx}(\tau)，導致我們不能按上式計算傅立葉轉換來得到\Gamma_{xx}(F)。從隨機過程的單個實現，可以計算時間平均自相關函數:
:R_{xx}(\tau) = \frac{1}{2T_0} \int_{-T_0}^{T_0} x^*(t)x(t+\tau)dt

其中2T_0是觀察期間。如果平穩隨機過程的一階和二階矩(均值和自相關函數)是各態歷經的，那麼
:\gamma_{xx}(F) = \lim_{T_0 \to \infty}R_{xx}(\tau) = \lim_{T_0 \to \infty} \frac{1}{2T_0} \int_{-T_0}^{T_0} x^*(t)x(t+\tau)dt

這一關係證實了時間平均自相關函數R_{xx}(\tau)可用做對統計自相關函數\gamma_{xx}(\tau)的估計。更進一步R_{xx}(\tau)的傅立葉轉換提供了對功率密度譜的估計P_{xx}(F),即:
:P_{xx}(F) = \int_{-T_0}^{T_0} R_{xx}(\tau)e^{-j2\pi F\tau}d\tau
::::= \frac{1}{2T_0} \int_{-T_0}^{T_0} [\int_{-T_0}^{T_0} x^*(t)x(t+\tau)dt]e^{-j2\pi F\tau}d\tau
::::= \frac{1}{2T_0} |\int_{-T_0}^{T_0}x(t)e^{-j2\pi Ft}dt|^2

實際功率密度譜是P_{xx}(F)在極限T_0 \rightarrow \infty時的期望值:
:\Gamma_{xx}(F) = \lim_{T_0 \to \infty} E[P_{xx}(F)]
::::= \lim_{T_0 \to \infty} E[\frac{1}{2T_0} |\int_{-T_0}^{T_0}x(t)e^{-j2\pi Ft}dt|^2]

我們將從隨機過程的單個實現樣本考慮功率密度譜估計。假設x_a(t)以F_s>2B取樣，其中B是隨機過程功率密度譜包含的最高頻率。因此通過對x_a(t)取樣，我們得到有限長序列x(n),0\leqslant n\leqslant N-1。從這些樣本我們計算時間平均自相關序列:
:r_{xx}^'(m) = \frac{1}{N-m} \sum_{n=0}^{N-m-1} x^*(n)x(n+m), m = 0,1,...,N-1

並且對於m的負值，我們有r_{xx}^'(m) = [r_{xx}^'(-m)]^*。於是我們計算傅立葉轉換:
:P_{xx}^'(f) = \sum_{m=-N+1}^{N-1} r_{xx}^'(m)e^{-j2\pi fm}

上上式中的歸一化因子N-m導致了均值估計:
:E[r_{xx}^'(m)] = \frac{1}{N-m} \sum_{n=0}^{N-m-1} E[x^*(n)x(n+m)] = r_{xx}(m)

其中r_{xx}(m)是x(n)的真實(統計的)自相關序列。因此r_{xx}^'(m)是自相關函數r_{xx}(m)的無偏差估計。估計r_{xx}^'(m)的方差近似為:
:var[r_{xx}^'(m)]\thickapprox \frac{N}{(N-m)^2}\sum_{n=-\infty}^{\infty}[|\gamma_{xx}(n)|^2+\gamma_{xx}^*(n-m)\gamma_{xx}(n+m)]

這是由Jenkins和Watts於1968年給出的結果，顯然\lim_{N \to \infty} var[r_{xx}^'(m)] = 0有\sum_{n=-\infty}^{\infty}|\gamma(n)|^2 。

因為E[r_{xx}^'(m)] = \gamma_{xx}(m),並且當N \rightarrow \infty時估計的方差收斂於零，所以估計r_{xx}^'(m)是相容的。對於較大值的滯後參數m，特別當m逼近於N時，由r_{xx}^'(m) = \frac{1}{N-m} \sum_{n=0}^{N-m-1} x^(n)x(n+m), m = 0,1,...,N-1給出的估計r_{xx}^'(m)具有較大方差。這是由於很少的數據點數進入大的滯後情況下的估計。作為式r_{xx}^'(m) = \frac{1}{N-m} \sum_{n=0}^{N-m-1} x^(n)x(n+m), m = 0,1,...,N-1的備用方法，我們使用如下的估計:
:r_{xx}(m) = \frac{1}{N} \sum_{n=0}^{N-m-1} x^*(n)x(n+m), 0\leqslant m \leqslant N-1
:r_{xx}(m) = \frac{1}{N} \sum_{n=|m|}^{N-1} x^*(n)x(n+m), m = -1,-2,1,...,1-N

其偏移為|m|\gamma_{xx}(m)/N，因為其均值是:
:E[r_{xx}(m)] = \frac{1}{N} \sum_{n=0}^{N-m-1} E[x^*(n)x(n+m)] = \frac{N-|m|}{N}\gamma_{xx}(m) = (1-\frac{N})\gamma_{xx}(m)

然而該估計具有較小的方差，近似為:
:var[r_{xx}(m)]\thickapprox \frac{1}{N}\sum_{n=-\infty}^{\infty}[|\gamma_{xx}(n)|^2+\gamma_{xx}^*(n-m)\gamma_{xx}(n+m)]

注意到r_{xx}(m)是漸進無偏的，即\lim_{N \to \infty} E(r_xx(m)) = \gamma_{xx}(m)。並且當N \rightarrow \infty時其方差收斂於零。因此估計r_{xx}(m)也是\gamma_{xx}(m)的一致估計。在處理功率譜估計問題時，我們將使用由式r_{xx}(m) = \frac{1}{N} \sum_{n=0}^{N-m-1} x^(n)x(n+m), 0\leqslant m \leqslant N-1和r_{xx}(m) = \frac{1}{N} \sum_{n=|m|}^{N-1} x^(n)x(n+m), m = -1,-2,1,...,1-N給出的估計r_{xx}(m)。相應的功率譜密度是:
:P_{xx}(f) = \sum_{m=-(N-1)}^{N-1} r_{xx}(m)e^{-j2\pi fm}

我們把得到的r_{xx}(m)代入上式，估計P_{xx}(f)可以表示為:
:P_{xx}(f) = \frac{1}{N}|\sum_{n=0}^{N-1}x(n)e^{-j2\pi fn}|^2 = \frac{1}{N}|X(f)|^2

其中X(f)是樣本序列x(n)的傅立葉轉換。這種常見形式的功率密度譜估計稱為週期圖。它最初是由Schuster於1898年引入用來檢測和測量存在於數據中的"隱藏週期"。從式P_{xx}(f) = \sum_{m=-(N-1)}^{N-1} r_{xx}(m)e^{-j2\pi fm}可推出週期圖估計P_{xx}(f)的均值是
:E[P_{xx}(f)] = E[\sum_{m=-(N-1)}^{N-1} r_{xx}(m)e^{-j2\pi fm}] = \sum_{m=-(N-1)}^{N-1} E[r_{xx}(m)]e^{-j2\pi fm}
:E[P_{xx}(f)] = \sum_{m=-(N-1)}^{N-1} (1-\frac{N})\gamma_{xx}(m)e^{-j2\pi fm}

對上兩式的解釋是，估記譜的均值是窗自相關函數的傅立葉轉換\tilde{\gamma}_{xx}(m) = (1-\frac{N})\gamma_{xx}(m)，其中窗函數是(三角形的)巴特利特(Bartlett)窗。因此估計譜的均值是:
:E[P_{xx}(f)]=\sum_{m=-\infty}^{\infty}\tilde{\gamma}_{xx}(m)e^{-j2\pi fm} = \int_{-0.5}^{0.5} \Gamma_{xx}(\alpha)W_B(f-\alpha)d\alpha

其中W_B(f)是巴特利特窗的譜特徵。上式說明了估計譜的均值是真實功率譜密度\Gamma_{xx}(f)與巴特利特窗傅立葉轉換W_B(f)的旋積。結果估計譜的均值是真實譜的平滑版，受損於有限數據點導致的相同的譜洩漏。注意到估計譜是漸進無偏的，即:
:\lim_{N \to \infty} E[\sum_{m=-(N-1)}^{N-1} r_{xx}(m)e^{-j2\pi fm}] = \sum_{m=-\infty)}^{\infty} r_{xx}(m)e^{-j2\pi fm} = \Gamma_{xx}(f)

然而一般來說當N \rightarrow \infty時估計P_{xx}(f)的方差不會衰減到零。例如當數據序列是一個隨機過程時，方差是var[P_{xx}(f)] = \Gamma_{xx}^2(f)[1+(\frac{sin2\pi fN}{Nsin2\pi f})^2],當N \rightarrow \infty時，極限為\lim_{N \to \infty}var[P_{xx}(f)] = \Gamma_{xx}^2(f)。

因此我們認為週期圖不是真實譜密度的一致估計(即不收斂於真正的功率譜密度)。概括來說估計的自相關r_{xx}(m)是真實自相關函數\gamma_{xx}(m)的一致估計。然而它的傅立葉轉換P_{xx}(f)即週期圖不是真實功率譜密度的一致估計。我們注意到P_{xx}(f)是\gamma_{xx}(f)的漸進無偏差估計，但是對於一個有限長序列，從式E[P_{xx}(f)]=\sum_{m=-\infty}^{\infty}\tilde{\gamma}_{xx}(m)e^{-j2\pi fm} = \int_{-0.5}^{0.5} \Gamma_{xx}(\alpha)W_B(f-\alpha)d\alpha得到的P_{xx}(f)均值包含了偏移，說明真實功率譜密度產生了失真。於是估計譜受損於巴特利特窗的平滑效應和具體的洩漏，平滑和洩漏最終限制了我們準確分析緊密譜的能力。

功率譜估計的非參數化方法
非參數方法是由Bartlett(1948)、Blackman和Tukey(1958)及Welch(1967)開發的方法，這些方法沒有假定數據是如何生成的，故被稱作非參數方法。既然這些估計全部基於數據的有限記錄，這些方法的頻率分辨率最好等於長度為N的矩形窗的寬度，也就是在-3dB近似為1/N。這將會更加精確地指定具體方法的頻率分辨率。為了減小譜估計的方差，非參數方法會降低頻率分辨率。

Bartlett方法:平均週期圖
減小週期圖方差的Bartlett方法包含了三個步驟。首先N點序列被劃分為K個不重疊段，每段的長度為M。這樣就生成了K個數據段:
:x_i(n) = x(n+iM), i = 0,1,...,K-1;n = 0,1,...,M-1

對於每一段，可計算週期圖:
:P_{xx}^{(i)}(f) = \frac{1}{M}|\sum_{n=0}^{M-1}x_i(n)e^{-j2\pi fn}|^2, i = 0,1,...,K-1

最後對K段的週期圖進行平均得到Bartlett功率譜估計:
:P_{xx}^{B}(f) = \frac{1}{K}\sum_{i=0}^{K-1}P_{xx}^{(i)}(f)

該估計的統計特性很容易得到。首先均值是E[P_{xx}^{B}(f)] = \frac{1}{K}\sum_{i=0}^{K-1}E[P_{xx}^{(i)}(f)] = E[P_{xx}^{(i)}(f)]。由式E[P_{xx}(f)] = E[\sum_{m=-(N-1)}^{N-1} r_{xx}(m)e^{-j2\pi fm}] = \sum_{m=-(N-1)}^{N-1} E[r_{xx}(m)]e^{-j2\pi fm}和E[P_{xx}(f)] = \sum_{m=-(N-1)}^{N-1} (1-\frac{N})\gamma_{xx}(m)e^{-j2\pi fm}及E[P_{xx}(f)]=\sum_{m=-\infty}^{\infty}\tilde{\gamma}_{xx}(m)e^{-j2\pi fm} = \int_{-0.5}^{0.5} \Gamma_{xx}(\alpha)W_B(f-\alpha)d\alpha可得到單個週期圖的期望值是:
:E[P_{xx}^{(i)}(f)]= \sum_{m=-(M-1)}^{M-1} (1-\frac{M})\gamma_{xx}(m)e^{-j2\pi fm}= \frac{1}{M}\int_{-0.5}^{0.5} \Gamma_{xx}(\alpha)(\frac{sin\pi (f-\alpha)M}{sin\pi (f-\alpha)})^2d\alpha

其中W_B(f) = \frac{1}{M}(\frac{sin\pi fM}{sin\pi f})^2是巴特利特窗w_B(n) =
\begin{cases}
1-\frac{M},& \mbox{if }|m|\leqslant \mbox{ M-1} \\
0,& \mbox{if }|m|\mbox{ is other}
\end{cases} 的頻率特性。

從E[P_{xx}^{(i)}(f)]的公式注意到，現在真實譜與巴特利特窗的頻率特性W_B(f)有關。數據長度從N點減小到M=N/K，導致窗口的譜寬度增加K因子。結果頻率分辨率得到減小因子K。分辨率降低的結果使得方差減小。Bartlett估計的方差是:
:var[P_{xx}^B(f)]= \frac{1}{K^2}\sum_{i=0}^{K-1}var[P_{xx}^{(i)}(f)]= \frac{1}{K}var[P_{xx}^{(i)}(f)]

如果我們利用式var[P_{xx}(f)] = \Gamma_{xx}^2(f)[1+(\frac{sin2\pi fN}{Nsin2\pi f})^2]和上式，得到:
:var[P_{xx}^B(f)]= \frac{1}{K}\Gamma_{xx}^2(f)[1+(\frac{sin2\pi fM}{Msin2\pi f})^2]

因此，Bartlett功率譜估計的方差減小因子為K。

參考文獻
*