查看“︁功率譜估計”︁的源代码

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

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

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

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

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

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

其中<math>2T_0</math>是觀察期間。如果平穩[[隨機過程]]的一階和二階[[矩]]([[均值]]和[[自相關函數]])是各態歷經的，那麼
:<math>\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</math>

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

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

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

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

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

其中<math>r_{xx}(m)</math>是<math>x(n)</math>的真實(統計的)自相關序列。因此<math>r_{xx}^'(m)</math>是[[自相關函數]]<math>r_{xx}(m)</math>的無偏差估計。估計<math>r_{xx}^'(m)</math>的[[方差]]近似為:
:<math>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)]</math>

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

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

其偏移為<math>|m|\gamma_{xx}(m)/N</math>，因為其[[均值]]是:
:<math>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)</math>

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

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

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

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

對上兩式的解釋是，估記譜的[[均值]]是窗自相關函數的[[傅立葉轉換]]<math>\tilde{\gamma}_{xx}(m) = (1-\frac{|m|}{N})\gamma_{xx}(m)</math>，其中[[窗函數]]是(三角形的)巴特利特(Bartlett)窗。因此估計譜的[[均值]]是:
:<math>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</math>

其中<math>W_B(f)</math>是巴特利特窗的譜特徵。上式說明了估計譜的均值是真實[[功率譜密度]]<math>\Gamma_{xx}(f)</math>與巴特利特窗[[傅立葉轉換]]<math>W_B(f)</math>的[[旋積]]。結果估計譜的[[均值]]是真實譜的平滑版，受損於有限數據點導致的相同的譜洩漏。注意到估計譜是漸進無偏的，即:
:<math>\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)</math>

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

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

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

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

對於每一段，可計算週期圖:
:<math>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</math>

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

該估計的[[統計]]特性很容易得到。首先[[均值]]是<math>E[P_{xx}^{B}(f)] = \frac{1}{K}\sum_{i=0}^{K-1}E[P_{xx}^{(i)}(f)] = E[P_{xx}^{(i)}(f)]</math>。由式<math>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}</math>和<math>E[P_{xx}(f)] = \sum_{m=-(N-1)}^{N-1} (1-\frac{|m|}{N})\gamma_{xx}(m)e^{-j2\pi fm}</math>及<math>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</math>可得到單個週期圖的[[期望值]]是:
:<math>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)})^2d\alpha</math>

其中<math>W_B(f) = \frac{1}{M}(\frac{sin\pi fM}{sin\pi f})^2</math>是巴特利特窗<math>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} </math>的[[頻率 (物理學)|頻率]]特性。

從<math>E[P_{xx}^{(i)}(f)]</math>的公式注意到，現在真實譜與巴特利特窗的[[頻率 (物理學)|頻率]]特性<math>W_B(f)</math>有關。數據長度從<math>N</math>點減小到<math>M=N/K</math>，導致窗口的譜寬度增加<math>K</math>因子。結果[[頻率 (物理學)|頻率]]分辨率得到減小因子<math>K</math>。分辨率降低的結果使得[[方差]]減小。Bartlett估計的[[方差]]是:
:<math>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)]</math>

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

因此，Bartlett功率譜估計的[[方差]]減小因子為<math>K</math>。

==參考文獻==
*{{Cite book | author = John G. Proakis & Dimitris G. Manolakis(方艳梅、刘永清 等譯) | title = 数字信号处理――原理、算法与应用（第四版） | publisher = 电子工业出版社 | date = 2007 | pages = 709-719 | ISBN = 9787121034961}}