一. 概述

随机信号的功率谱估计一直是常见的任务, 此文将一个统一的滤波视角, 罗列一下周期图谱估计(傅立叶变换), 自回归(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 的思想, 简单却十分强大.

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

最后

以上就是义气鞋子为你收集整理的几种经典功率谱密度估计方法的记录(从滤波的视角出发)一. 概述二. 统一的角度: 滤波三. 几种经典的功率谱密度估计方法四. 后记的全部内容，希望文章能够帮你解决几种经典功率谱密度估计方法的记录(从滤波的视角出发)一. 概述二. 统一的角度: 滤波三. 几种经典的功率谱密度估计方法四. 后记所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。