- 微信
- 微博
  
  分享文章到微博
- 复制链接
  
  复制链接到剪贴板

信号功率谱估计

aqhs 发表于 2022/06/21 09:03:03 2022/06/21

【摘要】信号的功率谱密度用于描述信号的功率随频率的变化情况，简称为功率谱。这里介绍几种功率谱估计常用方法，包括周期图法，Welch方法，基于自相关函数傅里叶变换的间接方法，自回归(AR)模型方法。为了减少数据截断而起泄露的估计误差，给出了几种数据加窗函数：矩形窗、Parzen窗、Hamming窗、Hanning窗、Welch窗、ExactBlackman窗、三角窗。最后用实例说明了功率谱估计。

信号功率谱估计

信号的功率谱密度用于描述信号的功率随频率的变化情况，简称为功率谱。这里介绍几种信号功率谱估计(Power Spectrum Density Estimate)常用方法，包括利用快速傅里叶变换的周期图法(periodogram)，长信号分段叠加计算的Welch方法，基于自相关函数傅里叶变换的间接方法(BT法)，自回归(AR)模型方法。为了减少由于数据截断而引起泄露所带来的功率谱估计误差，给出了几种数据加窗函数：矩形窗、Parzen窗、Hamming窗、Hanning窗、Welch窗、ExactBlackman窗、三角窗。最后用实例说明了功率谱估计。

(1)用FFT算法做功率谱估计

(1.1)周期图法

设等间隔采样的信号序列为x(i), i = 0 , 1 , … , N-1, 它的付利叶变换为X(k), k = 0, 1, …, N-1, 则功率谱的周期图(periodogram)估计Ｐ(k)为：

其中k对应于信号的零频率和正频率。如果信号的采样率为Fs，则k=0对应零频率，k=N/2代表Fs/2。

(1.2)窗函数

为了减少由于数据截断而引起泄露所带来的误差，一般在时域上对数据加窗，用户应根据信号的不同性质和不同处理目的来选择不同的窗函数。设信号为x(i)，窗函数为W(i)，那么数据加窗就是原数据乘以窗函数，即x’(i)=x(i)w(i) ， i = 0, 1, ……, N-1。

(1.3)韦尔奇(Welch)方法估计功率谱

1. 先将长度为N的信号序列x(i)分段，每段长度为L，相邻两段之间有M点数据重叠。

2. 将每段数据在时域上乘以窗函数。

3. 将加窗后的数据用FFT算法做周期图功率谱估计。

4. 将各段功率谱做平均，结果作为整个信号的功率谱估计。

韦尔奇(Welch)又称加权交叠平均法，是应用较广的一种方法。

(2)极大熵功率谱估计

(2.1)时间序列AR模型

平稳高斯随机序列x(n)可以假定为高斯白噪声w(n)通过线性系统产生，即：

写成时域表达式为：

x(t) + a(1)x(t-1) + … + a(p)x(t-p) = w(t)

上式称为时间序列AR(p)模型。

(2.2)极大熵AR功率谱估计

由随机信号的自回归(AR)模型计算功率谱：

通过估计AR模型的阶次和参数，然后由AR模型的传递函数计算出功率谱密度，这种方法称为极大熵谱估计方法。

(2.3)AR模型定阶

1 FPE准则

FPE(Final Prediction Error)准则是由模型的预报误差来判别AR模型的阶次是否恰当，其原理是如果数据所符合的真实模型是AR(p)，而我们用AR(n)(n>p或n<p)来拟合，那么导致预报误差的方差增大。

最终预报误差准则(FPE)定义为：

上式中N为数据的长度，R(i)为x(t)的自相关函数。

我们从低阶到高阶建立AR模型，得到估计参数并计算出FPE值，当FPE值达到最小时对应的阶次为AR模型的阶次。

2 AIC准则

AIC准则是由日本统计学家Akaike于1973年提出的，全称是最小化信息量准则(Akaike Information Criterion)。AIC准则定义为：

这里M是最大可能的阶数。

3 BIC准则

当样本容量N很大时，在AIC准则中拟合误差提供的信息就要受到样本容量的放大。BIC(Bayesian Information Criterion)贝叶斯信息准则是Schwartz在1978年根据Bayes理论提出的判别准则，BIC准则函数定义如下：

阶数p的选择为：

BIC(p) = min(BIC(n)), n=0, 1 ,.., M

M是最大可能的阶数。

(2.4)AR参数的Burg算法求解

设自回归模型AR(p)的系数为： a(i), i=1,2, …, p，则AR参数BURG递阶估计为：

当k=p时, 求得的系数φpj, j=1,2, …, p就是AR模型参数估计值, σ2p是残差方差。

由上述参数估计过程可以看出，可利用σ2k和AIC或BIC准则自动选定AR模型的阶数。

(2.5)AR参数的Yule-Walker方程求解

p阶AR过程x(t)满足差分方程：

x(t) + a(1)x(t-1)+ … +a(p)x(t-p)=w(t)

其中a(1), a(2), …, a(p)为AR系数，w(t)为纯随机过程。

在上式两边同时乘以x(t-k), k>1并取期望得：

r(k) + a(1)r(k-1) + a(2)r(k-2) + … + a(p)r(k-p)=0

上式中r(n)为x(t)的自相关函数，且r(-n)=r(n)。由此可得出下述线性方程组：

此方程称为Yule-Walker方程。

(2.6)AR参数Levinson-Durbin递推算法求解

直接求解Yule-Walker方程中的AR参数需要作矩阵求逆运算，当阶数增大时，运算量也显著增大。Levinson-Durbin算法对Yule-Walker方程的求解提供了一个高效的解法，其算法可描述如下：

(2.7)AR参数的正交变换求解

设AR(p)模型

正交变换法与Burg算法一样都适合短序列信号，而且还具有以下优点：较小的频率偏移；没有谱线分裂；较高的频率分辨率。

(3)自相关函数估计功率谱

利用信号的自相关函数估计计算功率谱，是功率谱估计的间接法，又称BT法。基于维纳-辛钦定理，1958年Blackman和Tukey给出了这一方法的具体实现，即先由信号x(i), i=0,1, …, N-1，估计出自相关函数r(i)，然后对r(i)求傅里叶变换便得到了x(i)的功率谱估计。

k = -m, -m+1, …, 0, …, m，且r(k)=r(-k)。

(4)功率谱幅度的归一化

功率谱体现功率随频率的变化情况，通常主要关注频谱幅度的相对变化，而对单个幅度值并不在意。但是有时也要关注各频率点上的幅度，或者频率区域的面积值，此时就要对幅度进行归一化。假设信号的均值为零，用信号的方差作为信号的总功率，计算各谱线上的功率分布值如下：

q(i)，i=0,1, …, m-1是归一化前的谱估计结果，p(i)是归一化后结果，D是信号的方差，Fs是采样率，频率轴坐标范围0-Fs/2，那么：

上述单位是每1Hz的功率，为了让p(i)变成每毫赫兹的功率则p(i)=p(i)/1000。

参考文献

[1] 胡广书编著. 数字信号处理—理论、算法与实现[M]. 北京清华大学出版社，1997年8月第1版，p324-341.
[2] 杨位钦，顾岚编著. 时间序列分析与动态数据建模[M]. 北京理工大学出版社，1986年第1版.
[3] 林代庸编著. 信号理论和应用[M]. 北京电子工业出版社，1990年3月第1版.

(5)信号功率谱估计举例

举例说明使用多道信号分析软件进行功率谱估计。

例1：一段0.5s的语音信号，采样率为16000Hz，分别采用周期图法、自相关法和AR模型法估计其功率谱密度。采用FFT法的Welch分段叠加周期图法，菜单操作：《分析》→《通用分析》→《功率谱估计》，计算结果如下。