现代数字信号处理——功率谱和信号频率估计

最新推荐文章于 2024-08-21 20:05:03 发布

程勇uestc

最新推荐文章于 2024-08-21 20:05:03 发布

阅读量7.3k

点赞数 7

分类专栏：现代数字信号处理文章标签：功率谱估计频率估计现代数字信号处理

本文链接：https://blog.csdn.net/qq_40438165/article/details/108044914

版权

现代数字信号处理专栏收录该内容

7 篇文章

订阅专栏

如何使用随机过程 $u (n) 的 N$ 个观测数据 $u_N(0),u_N(1),\cdots,u_N(N-1)$ 估计出随机过程的功率谱 $S (w)$ 。估计方法分为三个流派：

1.经典功率谱估计
2.参数模型法估计
3.基于相关矩阵特征分解的信号频率估计

一、经典功率谱估计

$\quad$ 基于传统傅里叶变换的思想，有BT法和周期图法及其相关改进。

1、BT法

$\quad$ 已知N个观测值 $u_N(n)$ ，则 $u (n)$ 的自相关函数估计值为 $\hat{r}(m)=\frac{1}{N}\sum_{n=0}^{N-1}u_N(n)n^*_N(n-m),|M|\le N-1$ 根据维纳=辛钦定理可以得到功率谱估计值： $\hat{S}_{BT}(w)=\sum_{m=-M}^{M}\hat{r}(m)e^{-jwm}$ $\quad$ BT法需要先计算自相关函数估计值，再以此为依据估计功率谱密度，称为间接法估计。显然这是一个 $n^2$ 级别的算法，如果观测点数目较多则难以计算，因此使用FFT加速计算过程。
$\quad$ 接下来，我们来看看自相关函数的估计性能，即估计值与真实值的差距，已知 $\hat{r}(m)=\frac{1}{N}\sum_{n=0}^{N-1}u_N(n)n^*_N(n-m)\\r(m)=lim_{N\rightarrow \infty}\frac{1}{N}\sum_{n=0}^{N-1}u_N(n)n^*_N(n-m)\\ 故而，E\{\hat{r}(m)\}=\frac{N-|m|}{N}r(m)$

1.固定 $∣ m ∣$ ，当 $N\rightarrow \infty$ 时， $E\{\hat{r}(m)\}=r(m)$ ， $\hat{r}(m)$ 是对 $r (m)$ 的渐进无偏估计
2.固定 $N$ ， $∣ m ∣$ 越接近 $N$ 时估计的偏差越大

$\quad$ 另一种估计是 $\hat{r}(m)=\frac{1}{N-|m|}\sum_{n=0}^{N-1}u_N(n)n^*_N(n-m)$ ，此时 $E\{\hat{r}(m)\}=r(m)$ ，是无偏估计。

2、周期图法

$\quad$ $\hat{r}(m)=\frac{1}{N}\sum_{n=0}^{N-1}u_N(n)n^*_N(n-m)=\frac{1}{N}u_N(m)u^*_N(-m)\\两边取傅里叶变换，得到\hat{S}_{PER}=\frac{1}{N}|U_N(w)|^2,U_N(w)=\sum_{n=1}^{N-1}u_N(n)e^{-jwn}$

$\quad$ 周期法和BT法异同：

当 $M = N - 1$ 时二者相同时一样
当 $M < N - 1$ 时BT法相当于是对长度为 $2 N - 1$ 的 $\hat{r}(m)$ 进行了截断处理，即加了一个矩形窗，所以BT法相当于是对周期法的平滑。

3、Bartlett法

$\quad$ 周期图法估计出的谱性能不好，当观测数据太长时谱的曲线起伏加剧；数据太短时分辨率不好（观测点越多分辨率越高）。Bartlett法的思想是：将N个观测点分为L段，每段长度为 $M = N / L$ 。对于每段数据先使用周期图法求其功率，再对每段功率谱结果平均，得到平均周期图。
$\quad$ Bartlett法功率谱估计频率分辨率下降为原来的 $1 / L$ ，方差也减小为原来的 $1 / L$ 。因此，Bartlett估计的功率谱比周期图法估计的功率谱更加平滑。

4、Welch法

$\quad$ Welch法又称修正平均周期法，是对Bartlett法的改进。Welch法也是对N个观测点分段，与Bartlett不同的是，Welch法允许各段的信号点有所交叠，通常取相邻两端信号交叠的一半。若每段信号长度为 $M$ ，信号被划分为 $L=\frac{N-M/2}{M/2}$ 段。将每段信号与窗函数 $w (n)$ 相乘后得到每段的功率谱估计，最后加权得到功率谱估计。
$\quad$ Welch法运行分段数据样本的重叠，于是可以得到更多的周期图估计，从而进一步减小功率谱密度的方差。通过窗函数的加权，可以减小样本段之间的相关性。
在这里插入图片描述

二、参数模型

$\quad$ 参数模型功率谱估计的基本思想是，认为随机过程 $u (n)$ 就是白噪声通过LTI离散时间系统得到的响应，利用观测样本估计出模型参数，也就得到了功率谱 $S(w)=|H(w)|^2\sigma_v^2$ 。
在这里插入图片描述
$\quad$ 规则过程的参数模型有三种：

AR参数模型
MA参数模型
ARMA参数模型

1、AR（P）参数模型

$\quad$ AR模型的输入 $v (n)$ 是零均值，方差为 $\sigma^2$ 的白噪声，输出为 $u (n)$ ，满足如下差分方程： $u(n)=-\sum_{k=1}^p a_ku(n-k)+v(n)$ 。阶数为 $P$ ，可以推得：

当 $m > 0$ ， $r_{uv}(m)=0$
当 $m = 0$ ， $r_{uv}(0)=\sigma^2$

$\quad$ 进而，可以推得AR模型的正则方程，也叫Yule-Walker方程：
在这里插入图片描述
$\quad$ 上述式子等价于：

$\quad$ 这样就可以求出参数 $a_k$ 和 $u (n)$ 功率谱：

$\quad$ 最终算法模型如下：

$\quad$ 因为矩阵求逆gao时间复杂度很高，一般使用AR参数模型的Levinson-Durbin迭代算法来进行求解。利用低阶推高阶。
$\quad$ AR模型阶数的选择对信号功率谱估计的质量有很大的影响，阶数越高分辨率更高。

AR模型性能

$\quad$ BT法分辨率随信号长度增加而提高，而AR模型避免了窗函数的影响，因此可以得到高的谱分辨率，且与真实值差距很小。AR模型是否稳定取决于 $u (n)$ 。对于高斯信号，最大熵谱估计与AR模型等价。缺点是质量受阶数p的影响，阶次太低，谱太平滑，反映不出谱峰。阶次选得过大可能会产生虚
假的谱峰。

2、MA模型

$\quad$ MA(q)为q阶模型，方程为 $u(n)=v(n)+\sum_{k=1}^qb_kv(n-k)$ ，可以计算得到 $r_u(m)=\sigma^2\sum_{k=0}^{q-m}b_k^*b_{k+m},m=0,1,\cdots,q$ 。这是一个非线性方程，求解复杂，因此可以用高阶AR模型来近似MA模型。
$\quad$ 无穷阶AR模型为 $H_{\infty}=\frac{1}{1+\sum_{k=1}^\infty a_kz^{-k}}$ ，等价于q阶MA模型为 $H_q(z)=1+\sum_{k=1}^qb_kz^{-k}=\frac{1}{1+\sum_{k=1}^\infty a_kz^{-k}}$ ，因此 $(1+\sum_{k=1}^qb_kz^{-k})(1+\sum_{m=1}^\infty a_mz^{-m})=1$ 。
$\quad$ 可以根据 $u (n)$ 估计AR模型的参数 $a_m$ ，根据 $e(m)=\hat{a}_m^{(p)}+\sum_{k=1}^qb_k\hat{a}_{m-k}^{(p)}\\\hat{a}_m^{(p)}=e(m)-\sum_{k=1}^qb_k\hat{a}_{m-k}^{(p)}$ 可以通过将 $e (m)$ 看作AR模型的输入，计算 $a_m$ 的相关函数来估计 $b_k$ 的值。
$\quad$ 根据正则方程，令 $m = 0$ 可以得到输入白噪声方差的估计为 $\hat{\sigma}^2=\frac{r_u(0)}{\sum_{k=0}^q|\hat{b}_k|^2}$ 。总结一下，MA模型求解过程如下：
在这里插入图片描述

3、ARMA模型

$\quad$ ARMA(p,q)模型综合AR和MA模型，表达式为 $u(n)=-\sum_{k=1}^pa_ku(n-k)+\sum_{k=1}^qb_kv(n-k)$
在这里插入图片描述

三、频率估计

1、MVDR

$\quad$ 三个常用的梯度：

$\Delta J(w)=2\frac{\partial}{\partial w^*}(c^Hw)=0$
$\Delta J(w)=2\frac{\partial}{\partial w^*}(w^Hc)=2c$
$\Delta J(w)=2\frac{\partial}{\partial w^*}(w^HRw)=2Rw$

$\quad$ M抽头的FIR滤波器如下：
在这里插入图片描述
$\quad$ 输入为随机过程 $x (n)$ ，输出为 $y(n)=\sum_{i=0}^{M-1}w_i^*x(n-i)=w^Hx(n)=x^T(n)w^*$ ，输出的平均功率为 $P=E\{|y(n)|^2\}=E\{w^Hx(n)x^H(n)w\}=w^HRw$ ，其中 $R=E\{x(n)x^H(n)\}$ 。
$\quad$ 假设滤波器输入信号 $x(n)=\sum_{k=1}^K\alpha_ke^{jw_kn}+v(n)=\sum_{i=1}^K\alpha(w_i)e^{jw_in}+v(n)$ 。设感兴趣的期望信号是角频率为 $w_1$ 的复正弦信号，则选择滤波器权向量 $w$ 的原则是，使复正弦信号 $e^{jw_1n}$ 无失真地通过滤波器，而尽量抑制其余频率的信号和噪声。因此，滤波器权重 $w$ 满足：

$w^Ha(w_1)=1$
$P=w^HRw$ 最小化

$\quad$ 利用拉格朗日乘子法，可以求得最优权向量为 $\hat{w}_{MVDR}=\frac{\hat{R}^{-1}a(w)}{a^H(w)\hat{R}^{-1}a(w)}$ 滤波器最小输出功率正相关于 $P_{MVDR}(w)=\frac{1}{a^H(w)\hat{R}^{-1}a(w)}$ ，这也是MVDR谱估计。
在这里插入图片描述
$\quad$ MVDR谱分辨率随着抽头数M 的增大而提高；MVDR不是信号的功率谱，但是它描述了信号功率的相对强度。

2、相关矩阵特征分解

信号子空间和噪声子空间

$\quad$ $x(n)=\sum_{k=1}^K\alpha_ke^{jw_kn}+v(n)=As(n)+v(n)$ ， $R=E\{x(n)x^H(n)\}\\=APA^H+E\{v(n)x^v(n)\}\\=APA^H+\sigma_v^2I \\ 其中，P=E\{s(n)x^s(n)\}=diag\{|\alpha_1|^2,|\alpha_1|^2,\cdots,|\alpha_k|^2\}$

$\quad$ 矩阵 $APA^H$ 共有k个非零特征值，对其进行特征值分解，可以得到 $APA^H=\sum_{i=1}^k\lambda_iu_iu_i^H$ 。
信号子空间 $E_s$ ：由特征向量 $u_1,\cdots,u_k$ 生成的子空间， $E_s=span\{u_1,\cdots,u_k\}$
噪声子空间 $E_N$ ：由特征向量 $u_{K+1},\cdots,u_{M}$ 生成的子空间， $E_N=span\{u_{K+1},\cdots,u_M\}$

MUSIC算法

$\quad$ 信号频率向量 $a(w_k)$ 与噪声子空间的特征向量正交，有 $\sum_{i=K+1}^M|a^H(w_k)u_i|^2$ ，用噪声子空间的向量构成矩阵 $G=[u_{K+1},\cdots,u_M]$ ，可得 $a^H(w_k)GG^Ha(w_k)=0\\ \hat{P}_{MUSIC}=\frac{1}{a^H(w_k)GG^Ha(w_k)}$ 信号的角频率由几个峰值确定。
在这里插入图片描述

Root-MUSIC

$\quad$ 定义向量 $a(z)=[1,z^{-1},\cdots,z^{-(M-1)}]^T$ ，求出 $a^T(z^{-1})GG^Ha(z)=0$ 的根。

3、Pisarenko谐波提取方法

$\quad$ $x(n)=\sum_{k=1}^K\alpha_ke^{jw_kn}+v(n)$ ，求解 $f(z)=u_0+u_1z+\cdots+u_Kz^K=0$ ，得到K个根即为信号频率的估计。
在这里插入图片描述
$\quad$ Pisarenko谐波分解方法还可以估计出各频率信号的功率和噪声功率，噪声功率估计为 $\sigma_v^2=r(0)-\sum_{k=1}^K|\alpha_k|^2$ 。