功率谱估计 – 小飞侠

对于随机信号，不存在傅里叶变换，其对应的功率谱估计方法主要包括非参数方法（包括韦尔奇）和参数方法（ARMA, Burg）。

周期图方法

fs = 1000;
t = 0:1/fs:1;
xn = sin(2*pi*80*t) + 2*sin(2*pi*140*t) + randn(size(t));
N = length(xn);
y = abs(fft(xn,1024));
Pxx = y.^2/N;
f = [0:N-1]*fs/N;
plot(f(1:N/2+1),Pxx(1:N/2+1));

我们绘制该信号的功率谱图，有

如何降低功率谱估计的方差？

(1). 将长度为N的数据分成若干段，分别求出每一段的功率谱，然后加以评价，甚至允许每段数据有部分重叠； (2). 采用合适的窗函数来消除由矩形窗旁辫带来的谱失真。

我们首先尝试非重叠下的功率谱绘制。

Pxx = (abs(fft(xn(1:256))).^2 + abs(fft(xn(257:512))).^2+...
    abs(fft(xn(513:768))).^2)/(256*3);
figure
plot((0:255)/256*fs,10*log10(Pxx));
title(‘平均周期(非重叠)')
xlabel('Freq'); 
ylabel('PSD(dB)');
grid

然后尝试有一般数据重叠的情况。

Pxx = (abs(fft(xn(1:256))).^2+…
abs(fft(xn(129:384))).^2+...
    abs(fft(xn(257:512))).^2+…
abs(fft(xn(385:640))).^2+...
 abs(fft(xn(513:768))).^2+…
abs(fft(xn(641:896))).^2)/(256*6)
figure
plot((0:255)/256*fs,10*log10(Pxx));
title(‘平均周期图（128个样本重叠）');
xlabel('Freq'); 
ylabel('PSD(dB)');
grid;

由于重叠会使各段之间具有统计相依性，反而会导致方差增大，所以在分段数目与重叠率之间的选择上存在着一个折衷。而且，随着分段的增加，虽然方差减小了，但估计的偏差会变大，因此在使用平均周期图法时，需要在估计期望值偏差和估计方差之间进行权衡。

改进周期图估计方法

非矩形窗方法：在周期图前给每段信号加一个非矩形的数据窗。 (1) 因为窗在两边渐变为0，所以这种方法降低了由于重叠导致的段间统计相依的效应。 (2) 一个非矩形窗可以减小“旁辫效应”，又称为“功率谱泄露”，但同时也使谱峰的宽度增大。 (3) 取合适的窗函数（海明、汉宁或凯塞）和一半段长度的重叠率，可以最有效的降低估计的偏差。

韦尔奇(pwelch)方法

clc,clear
fs = 1000;
t = 0:1/fs:1;
xn = sin(2*pi*80*t) + 2*sin(2*pi*140*t)+randn(size(t));
nfft = 256;
window = hanning(256);
noverlap = 128;
[Pxx, f] = pwelch(xn,window,noverlap,nfft,fs);
figure
plot(f,10*log10(Pxx));
title(‘韦尔奇估计谱(Pxx(f)/fs)');
xlabel('Freq'); ylabel('PSD(dB)');
grid;

结果比前面例子低30dB，这是因为pwelch函数将计算的功率谱结果除以抽样频率fs, 因为fs = 1000Hz,所以10*log10(1000)正好为30dB。

现代谱估计的参数方法

韦尔奇谱估计方法缺点：谱分辨率低。

现代谱估计参数方法：
$x(n)$
是白噪声通过某个模型产生的，不必认为N个以外的样本数据为0.

基本步骤：选择一个正确的模型；用已观测到的样本数据或自相关函数的数据来确定模型的参数；由此模型求出功率谱估计。

对平稳随机信号，3 种常用的线性模型分别是 AR 模型（自回归模型 Auto-Regression Model）、MA模型（滑动平均模型 Moving Average Model）和 ARMA 模型（自回归滑移平均模型 Auto-Regression-Moving Average Model）。

MA 模型：

随机信号
$x(n)$
由当前的激励
$\omega(n)$
和若干次过去的激励
$\omega(n-k)$
线性组合产生：
$x(n)=\sum_{k=0}^{q}b_k \omega(n-k)$

该模型的系统函数是：
$H(z)=\frac{X(z)}{W(z)}=\sum_{k=0}^{q}b_kz^{-k}$

q 表示系统阶数，系统函数只有零点，没有极点，所以该系统一定是稳定的系统，也称为全零点模型，用MA(q)来表示。

AR模型

随机信号
$x(n)$
由本身的若干次过去值
$x(n-k)$
和当前的激励值
$\omega(n)$
线性组合产生：
$x(n)=\omega(n)-\sum_{k=1}^{p}a_kx(n-k)$
.

该模型的系统函数是：
$H(z)=\frac{1}{1+\sum_{k=1}^{p}a_kz^{-k}}$

p是系统阶数，系统函数中只有极点没有零点，也称为全极点模型，系统由于极点的原因，要考虑到系统的稳定性，因而要考虑到极点的分布位置，用AR(p)来表示。

ARMA模型

ARMA模型是AR和MA模型的结合：
$x(n)=\sum_{k=0}^{q}b_k \omega(n-k)-\sum_{k=1}^{p}a_kx(n-k)$
。该模型的系统函数是：
$H(z)=\frac{\sum_{k=0}^{q}b_k z^{-k}}{1+\sum_{k=1}^{p}a_kz^{-k}}$
。

它既有零点又有极点，所以也称为极零点模型，要考虑极零点的分布位置，保证系统的稳定，用ARMA(p,q)表示。在随机信号时域分析中，提出了许多数学模型，用来在已知最大不确定原则下预测将来值，其优点是只需要很少的已知值。但是，它不能用在信号具有确定性的场合。在确定信号的情况下，信号是由确定的数学方程预测的。这点要特别注意。

应用较多的是AR模型，因为建立这种模型的计算工作比较容易。作为数学逼近，3 种模型都可以互相转换。实际中选用哪一种模型就要考虑到节约和计算量，选定模型后，剩下的任务就是用适当的算法估计模型参数
$(a_k, b_k, p, q)$
, 以便用模型对随机信号进行预测。

模型选择：主要考虑原则，模型能够表示谱峰、谱谷能力。对于具有尖峰的谱，应该选用具有极点的模型，如AR和ARMA模型；对于具有平坦的谱峰和深谷的信号，可以选用MA模型；既有极点也有零点的谱应该选用ARMA模型。对于没有模型准确反映其谱特性的信号，可以选用高阶的AR模型近似表示。在选择模型合适基础上，应尽量减少模型的参数。

下图表示的是用MA模型估计二阶AR信号功率谱的例子。

如果能确定信号
$x(n)$
的信号模型，根据信号观测数据求出模型参数，系统函数用
$H(z)$
表示，模型输入白噪声方差为
$\sigma _{\omega}^2$
信号的功率谱可用下式求出：

$P_{xx}(e^{j \omega})=\sigma_{\omega}^2 \left | H(e^{j\omega}) \right|^2$

功率谱估计可分成三个步骤：

(1)选择合适的信号模型：

(2)根据
$x(n)$
有限的观测数据，或者它的有限个自相关函数估计值，估计模型的参数；

(3)计算模型输出功率谱。

这种方法隐藏着数据的外推，不是假设观测区以外的数据为0，它的估计分辨率都比经典谱估计高，故现代谱估计也统称为高分辨率谱估计。

无论采用上述3种模型中的哪一种，有两个关键问题必须解决：一是模型阶数的确定，即模型中p，q的值如何确定；二是在模型阶数确定后，模型参数的确定，即模型中, ak, bk的值如何确定。 Box-Jenkins法是较常见的平稳时间序列建模方法，主要通过自相关函数、偏自相关函数的统计特征为依据来识别模型，其基本步骤有四步，如下所示。 (1)模型识别，即对一个观察序列，用其相关图和偏相关图来从各种模型族中选择一个与其实际过程相吻合的模型结构，找到合适的p和q的值。 (2)参数估计，估计模型中所含自回归和移动平均项的参数。 (3) 模型确认，看所选的模型对数据拟合是否够好。 (4) 预测，利用所选模型对时间序列进行一步或多步的预测。

AR模型参数和自相关函数的关系

根据
$x(n)=w(n)-\sum_{k=1}^{p}a_kx(n-k)$
，对该式两边同时乘以
$x(n-m)$
，然后求均值：
$\small E\left[x(n)x(n-m) \right ]=E \left[ w(n)x(n-m) \right -\sum_{k=1}^{p}a_kx(n-k)x(n-m) ]$
。

因为自相关函数为：
$\small R_{x}(m)=E\left[x(n)x(n+m)\right ]=E\left[x(k-m)x(k)\right ]=E\left[x(n)x(n-m) \right]=R_{xx}(-m)$

原文链接：https://blog.csdn.net/clelouch/article/details/110704763