几种经典功率谱密度估计方法的记录(从滤波的视角出发)

一. 概述

随机信号的功率谱估计一直是常见的任务, 此文将一个统一的滤波视角, 罗列一下周期图谱估计(傅立叶变换), 自回归(AR)模型, 最小方差无失真响应(MVDR)法, 正交子空间(MUSIC)方法的公式, 看一下它们是如何对自相关矩阵的特征向量, 特征向量进行加权的.
本文不涉及原理的推导, 具体原理和统计特性分析可以参考: <<现代数字信号处理及其应用 (何子述.等)>>.

1.1 前提!

一般来说, 只针对广义平稳随机过程; 由于样本统计特性都是用时间平均计算的, 故还需要满足各态历经性.

二. 统一的角度: 滤波

由于不同的专业或课本都有不同的说法, 此处为了统一, 就采用课本中的滤波视角, 即对于一段有限的数字信号向量 $\vec{\mathbf{x}}(n)$, 想用一个权向量(列向量) $\vec{\mathbf{w}}$ 对其滤波(相乘), 来提取出某一个频率的(相对)功率谱密度值: P = E { ∣ w ⃗ H x ⃗ ( n ) ∣ 2 } P = \mathbb{E}\{ | \vec{\boldsymbol w}^H \vec{\boldsymbol x}(n) |^2\} P=E{∣w Hx (n)∣2}. 对于不同的频率成分, 只要另 $\vec{\mathbf{w}}\to \vec{\mathbf{w}}(\omega )$ 即可, 则功率谱密度为:
P ( ω ) = E { ∣ w ⃗ ( ω ) H x ⃗ ( n ) ∣ 2 } = w ⃗ ( ω ) H R   w ⃗ ( ω ) . P(\omega) = \mathbb{E}\{ | \vec{\boldsymbol w}(\omega)^H \vec{\boldsymbol x}(n) |^2\} = \vec{\boldsymbol w}(\omega)^H \boldsymbol{R}\ \vec{\boldsymbol w}(\omega). P(ω)=E{∣w (ω)Hx (n)∣2}=w (ω)HR w (ω).

2.1 一些符号定义

上式中粗体为矩阵或向量, 非粗体一般为标量,

• $\vec{\mathbf{x}}(n)=[x(n)x(n-1)x(n-2)\cdots x(n-M+1){]}^{\top }$ 为 $n$ 时刻下采样到的列向量, 长为 $M$;
• R = E { x ⃗ ( n ) x ⃗ ( n ) H } \boldsymbol{R}=\mathbb E \{\vec{\boldsymbol x}(n)\vec{\boldsymbol x}(n)^H\} R=E{x (n)x (n)H} 是它的自相关矩阵(平稳), 同时它具有: Hermitian 正定, Toeplitz 的性质. 如果对他做一个谱分解/酉相似对角化/特征值分解, 并把特征值从大到小排列, 则有
  $$\mathbf{R}=\mathbf{Q}\mathbf{\Lambda }{\mathbf{Q}}^H=\sum _{k=1}^M{\lambda }_k{\vec{\mathbf{q}}}_k{\vec{\mathbf{q}}}_k^H$$

另外, 有一种很常见的 $\vec{\mathbf{w}}(\omega )$, 即标准频率列向量 $\vec{\mathbf{a}}(\omega )=[1e^{j\omega }{e}^{j2\omega }{e}^{j3\omega }\cdots e^{j(M-1)\omega }{]}^{\top }$.

基于此, 标准 DTFT 变换可以写为:
$$\mathrm{D}\mathrm{T}\mathrm{F}\mathrm{T}[\vec{\mathbf{x}}]=X(\omega )=\sum _{n=-\infty }^{+\infty }x(n)e^{-j\omega n}=\vec{\mathbf{a}}(\omega {)}^H\vec{\mathbf{x}},$$
其中 $\vec{\mathbf{x}}=\vec{\mathbf{x}}(M-1)=[x(M-1)x(M-1)\cdots x(1)x(0){]}^{\top }$, 只有时间 $0\sim M-1$ 内有值, 其余为 $0$.

三. 几种经典的功率谱密度估计方法

3.1 经典法–周期图谱估计(直接法)

基本思想

即对序列进行傅立叶变换并取模方. 这种方式相当于认为给这一段随机信号外插 $0$, 显然精度不会很好, 方差特性也不大行.

有一个和这个很接近的方法是 BT 法 (间接法), 它是将序列先求自相关函数后(此时可以取中间一部分), 再进行傅立叶变换, 根据维纳辛钦定理, 自相关函数的傅立叶变换即为功率谱密度. 两种方法在 BT 保留全部自相关函数时完全等价. 更详细分析可以看课本.

权向量

没错, 此法的权向量就是:
$$\vec{\mathbf{w}}(\omega )=\vec{\mathbf{a}}(\omega )={\left[1e^{j\omega }{e}^{j2\omega }{e}^{j3\omega }\cdots e^{j(M-1)\omega }\right]}^{\top }.$$

功率谱密度

P ( ω ) = E { ∣ a ⃗ ( ω ) H x ⃗ ( n ) ∣ 2 } = E { ∣ D T F T [ x ⃗ ( n ) ] ∣ 2 } = a ⃗ ( ω ) H R   a ⃗ ( ω ) = ∑ k = 1 M λ k ∣ a ⃗ ( ω ) H q ⃗ k ∣ 2 P(\omega) = \mathbb{E}\{ | \vec{\boldsymbol a}(\omega)^H \vec{\boldsymbol x}(n) |^2\} =\mathbb{E}\{ |\mathrm{DTFT}[\vec{\boldsymbol x}(n)]|^2 \}\\=\vec{\boldsymbol a}(\omega)^H \boldsymbol{R}\ \vec{\boldsymbol a}(\omega)=\sum_{k=1}^{M} \lambda_k|\vec{\boldsymbol a}(\omega)^H \vec{\boldsymbol q}_k|^2 P(ω)=E{∣a (ω)Hx (n)∣2}=E{∣DTFT[x (n)]∣2}=a (ω)HR a (ω)=k=1∑M​λk​∣a (ω)Hq ​k​∣2
此处的 $n$ 其实取的是 $M-1$, 因为基本都把第一个采样点为时间起点. 周期图法求某个频率 $\omega$ 处的功率密度的做法其实是, 把一个固定的频率向量向自相关矩阵的特征向量上投影, 再把内积的模方用特征值加权后相加.

PS: 虽然实际计算时, 若只有一个采样样本, 是没有取期望的操作的; 但为了分析统计特性, 假设有无穷多个 $M$ 长的样本(外插为 $0+0j$) 做平均, 即取期望. 以下算法一样, 实际情况下的自相关矩阵都要用数据进行估计, 而公式中都假设为统计情况下的矩阵.

3.2 基于自回归(AR)模型的功率谱估计

基本思想

上一种方式是完全没有任何先验, 直接根据一次观测进行估计. 这个方法是引入了先验知识. 首先, 从 Wold 分解定理出发, 我们认为一个随机过程通常为具有连续谱的离散时间规则随机过程 (即没有若干较强的复正弦信号); 进而, 如果还满足 Paley-Wiener 条件, 则它可以由一个白噪声激励一个最小相位 LTI 系统产生. 利用信号系统的知识, 此时就可以把白噪声(输入)和此随机过程(输出)用差分方程联系起来, 即可以写成一个 ARMA 模型, 此方法只取一个 AR, 此时系统函数是一个全极点模型.
此模型由于可以用公式无限递推, 所以相当于对一段随机过程样本能做很长很长的外插, 而不是暴力补零.

权向量

可以从功率谱密度结论凑得,
$$\vec{\mathbf{w}}(\omega )=\frac{{\mathbf{R}}_{p+1}^{-1/2}\cdot {\vec{\mathbf{e}}}_1}{\vec{\mathbf{a}}(\omega {)}^H\cdot {\mathbf{R}}_{p+1}^{-1}\cdot \sigma {\vec{\mathbf{e}}}_1}$$
其中, $p<<M$ 为远小于样本点数的 AR 模型阶数, ${\vec{\mathbf{e}}}_1=[1000\cdots 0{]}^{\top }$ 为仅第一个元素为 $1$ 的列向量, $\sigma$ 为白噪声的功率开平方根.

功率谱密度

$$\begin{array}{rl}P(\omega )& =\frac{{\sigma }^2}{|1+{\sum }_{k=1}^p{a}_k{e}^{-j\omega k}{|}^2}\stackrel{{\vec{\mathbf{\theta }}}_{p+1}=[1a_1{a}_2\cdots a_p{]}^{\top }}{=}\frac{{\sigma }^2}{|\vec{\mathbf{a}}(\omega {)}^H{\vec{\mathbf{\theta }}}_{p+1}{|}^2}\\ & \underset{{\mathbf{R}}_{p+1}^{-1}={\sum }_{k=1}^{p+1}{\lambda }_k{\vec{\mathbf{q}}}_k{\vec{\mathbf{q}}}_k^H}{\overset{{\vec{\mathbf{\theta }}}_{p+1}=({\mathbf{R}}_{p+1}^{\ast }{)}^{-1}\cdot \sigma {\vec{\mathbf{e}}}_1}{=}}\frac{{\sigma }^2}{|\vec{\mathbf{a}}(\omega {)}^H{\mathbf{R}}_{p+1}^{-1}\cdot \sigma {\vec{\mathbf{e}}}_1{|}^2}\\ & =\frac{1}{|{\sum }_{k=1}^{p+1}\frac{\sigma }{{\lambda }_k}{q}_{k1}^{\ast }\cdot \vec{\mathbf{a}}(\omega {)}^H{\vec{\mathbf{q}}}_k{|}^2}\end{array}$$
其中, $q_{k1}$ 表示 ${\vec{\mathbf{q}}}_k$ 的第 $1$ 个元素. 式子中的额外标注均遵循课本表示, 如 ${\vec{\mathbf{\theta }}}_{p+1}$ 可由 Yule-Walker 方程求得. 我们能看到, 这种方式的 PSD 和 $p+1$ 阶自相关矩阵特征向量有关.

3.3 最小方差无失真响应(MVDR)

基本思想

这个算法在很多地方都有, 如波束成形, 空域滤波等等. 它的基本思想恰如其名. 首先, MVDR 引入的先验假设是输入随机过程为 $K$ 个复正弦信号加上白噪声:
$$u(n)=\sum _{k=1}^K{\alpha }_k{\mathrm{e}}^{\mathrm{j}{\omega }_{k}n}+v(n)$$
它这样来选取某个感兴趣频率下的权向量 $\vec{\mathbf{w}}({\omega }_1)$: 使这个复正弦信号 $e^{j{\omega }_{1}n}$ 能无失真通过滤波器, 同时最小化此时的平均输出功率, 也就抑制了其他频率成分和噪声, 最后写成一个简单的凸优化问题来求解. 它的倒数是 $0\sim M$ 阶 AR 方法倒数的和, 有点调和平均的意思, 所以它的分辨率可能比 AR 方法稍低.

权向量

可以用拉格朗日乘子法求得:
$$\vec{\mathbf{w}}(\omega )=\frac{{\mathbf{R}}^{-1}\vec{\mathbf{a}}(\omega )}{\vec{\mathbf{a}}(\omega {)}^H\cdot {\mathbf{R}}^{-1}\vec{\mathbf{a}}(\omega )}$$

(相对)功率谱密度

$$P(\omega )=\vec{\mathbf{w}}(\omega {)}^H\mathbf{R}\vec{\mathbf{w}}(\omega )=\frac{1}{\vec{\mathbf{a}}(\omega {)}^H\cdot {\mathbf{R}}^{-1}\vec{\mathbf{a}}(\omega )}=\frac{1}{{\sum }_{k=1}^M\frac{1}{{\lambda }_k}|\vec{\mathbf{a}}(\omega {)}^H{\vec{\mathbf{q}}}_k{|}^2}$$

3.4 基于特征子空间的多重信号分类(MUSIC)算法

基本思想

这是一种超分辨率的算法, 它同样也是基于信号为 $K$ 个复正弦+白噪声的假设. 它将自相关矩阵的特征值分解结果(对角阵中的特征值从大到小排)中的酉矩阵 $\mathbf{Q}$ 的前面 $K$ 列视为信号子空间, 包含了噪声+信号, 对应的特征值比较大; 而剩下的列构成了噪声子空间 $G$, 仅仅包含噪声, 它们对应的特征值比较小, 就是白噪声的方差. 由于正交性, 如果含有频率 ${\omega }_1$, 则 $\vec{\mathbf{a}}({\omega }_1)$ 会和 $G$ 正交, 所以 $\vec{\mathbf{a}}({\omega }_1)$ 和 $G$ 的所有列的内积的模方和就会很小, 取个倒数就会很大. 就这样构造一个展现频率成分强度的函数, 来展示此样本的 PSD 趋势, 不是真的功率谱.

方便看懂书上推导的结论:
$\mathrm{R}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(AP{A}^H)\subseteq \mathrm{R}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(A)$, 又因二者秩相同都为信号数 $K$, 所以二者的值域空间相同, 即 $\mathrm{R}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(AP{A}^H)=\mathrm{R}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(A)$.
由于 $AP{A}^H$ 为纯信号, 和噪声独立, 所以它和噪声子空间正交, 也就是 $A$ 与噪声子空间正交, 也就是 $A$ 中的每一列和 $G$ 的任意一列都正交.

权向量

可以凑出来:
$$\vec{\mathbf{w}}(\omega )=\frac{\mathbf{G}{\mathbf{\Lambda }}_N^{-1/2}{\mathbf{G}}^H\vec{\mathbf{a}}(\omega )}{\vec{\mathbf{a}}(\omega {)}^H\mathbf{G}{\mathbf{G}}^H\vec{\mathbf{a}}(\omega )}$$
其中, ${\mathbf{\Lambda }}_N\in {\mathbb{R}}_{+}^{(M-K)\times (M-K)}$ 为噪声子空间对应的特征值构成的对角阵, 也就是自相关矩阵特征值构成的, 从大到小排好的对角阵 $\mathbf{\Lambda }$ 的右下角一块儿.

功率谱密度(伪谱, 假的)

$$\tilde{P}(\omega )=\frac{1}{\vec{\mathbf{a}}(\omega {)}^H\mathbf{G}{\mathbf{G}}^H\vec{\mathbf{a}}(\omega )}=\frac{1}{{\sum }_{k=K+1}^M|\vec{\mathbf{a}}(\omega {)}^H{\vec{\mathbf{q}}}_k{|}^2}$$
可以看出, 它直接找出了噪声部分的特征向量, 抓住了主要矛盾.

四. 后记

功率谱估计是数字信号处理的第一站, 除了知道一些基本算法及其原理外, 我们需要从其中归纳出共性的东西:

• 随机信号按时间采样得到的 $M$ 维随机向量(零均值), 从本文看得到, 它经常和自相关矩阵的特征向量有关; 进一步, 它在自相关矩阵特征值分解得到的酉矩阵 $Q$ 张成的空间中, 即它是各个归一化特征向量的加权和, 这就是所谓的 Karhunen-Loeve 展开.
• 信号在高维下是在一个流形(Manifold)上, 所以我们不能仅仅认为多采样几个点只是多了几个点, 能提高传统估计方法的性能; 而是从高维空间去认识信号, 就像 MUSIC 的思想, 简单却十分强大.

另外, 机器学习中也常常认为数据分布于一个流形上, 用高维思维去看待问题会很有帮助.