功率谱估计

自然界出现的许多现象都可以在统计平均意义上很好的表现出来。例如，气象学中的气温与气压的变动等，均可以以统计的方式表示为随机过程。在电阻器和电子设备中生成的热噪音电压，也是被抽象为随机过程模型的物理讯号的例子。由于这些讯号为随机讯号，我们必须采用一种统计观点来处理随机讯号的平均特征。特别的是随机讯号的自相关函数很适合用于代表时域中的随机讯号，并且自相关函数的傅立叶转换可生成功率谱密度，也可提供时域到频域的转换。

基于有限长讯号观察的功率谱估计

对于有限时间长度的讯号，数据序列的有限记绿长度是对功率谱估计(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估计的方差是: