功率譜估計

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

基於有限長訊號觀察的功率譜估計

對於有限時間長度的訊號，數據序列的有限記綠長度是對功率譜估計(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}<\infty$。

因為$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 {|m|}{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_{x}x(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 {|m|}{N}})\gamma _{xx}(m)e^{-j2\pi fm}$

對上兩式的解釋是，估記譜的均值是窗自相關函數的傅立葉轉換${\tilde {\gamma }}_{xx}(m)=(1-{\frac {|m|}{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 {|m|}{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|}{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 )}})^{2}d\alpha$

其中$W_{B}(f)={\frac {1}{M}}({\frac {sin\pi fM}{sin\pi f}})^{2}$是巴特利特窗$w_{B}(n)={\begin{cases}1-{\frac {|m|}{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$。

參考文獻

• John G. Proakis & Dimitris G. Manolakis(方艳梅、刘永清 等譯). 数字信号处理――原理、算法与应用（第四版）. 电子工业出版社. 2007: 709–719. ISBN 9787121034961.