(完整版)功率谱估计性能分析及Matlab仿真

1 引言

随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱[1]。

信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类：经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。这是不符合实际情况的，因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg法等。

2 经典功率谱估计

经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的个样本数据估计出其功率谱[1]。

2.1 周期图法( Periodogram )

Schuster首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。

取平稳随机信号的有限个观察值，求出其傅里叶变换

然后进行谱估计

周期图法应用比较广泛，主要是由于它与序列的频谱有直接的对应关系，并且可以采用FFT快速算法来计算。但是，这种方法需要对无限长的平稳随机序列进行截断，相当于对其加矩形窗，使之成为有限长数据。同时，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，从而产生频谱泄露，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，使得谱分辨率较低[1]。

该方法基于Matlab实现的程序：

clear all;

load test x;

N=4096;

Fn=-0.5:1/N:0.5-1/N;

px=fft(x,N);

pmax=max(px);%归一化

px=px/pmax;

px=10*log10(px+0.000001);

plot(Fn,fftshift(px));grid on;

图1 周期图法 

图2 周期图法 

说明：

(1) 本报告仿真中所采用的用于功率谱估计的数据文件来自参考文献[3]的test.dat。该数据为128点复序列（图3），由复数噪声加上四个复正弦组成。其归一化频率分别是：。

图3 复序列 

(2) 从仿真图可以清晰看到，和不能完全分开，仅在波形的顶部能看出是两个频率分量；此外，当数据长度太大时（图1），谱曲线呈现较大的起伏；当数据长度太小时（图2），谱的分辨率又不好。据此，周期图法不满足一致性估计条件。

2.2 自相关法( BT法)

自相关法的理论基础是维纳—辛钦定理。1958年Blackman和Tukey给出了这一方法的具体实现。

对于平稳随机信号来说，其自相关函数是确定性函数，故其功率谱也是确定的。这样可由平稳随机离散信号的有限个离散值求出自相关函数

然后在内对做傅里叶变换，得到功率谱

该方法基于Matlab实现的程序：

clear all;

load test x;

N=4096;

Fn=-0.5:1/N:0.5-1/N;

Mlag=;

rx=xcorr(x,Mlag,'unbiased');

px=fft(rx,N);

pmax=max(px);%归一化

px=px/pmax;

px=10*log10(px+0.000001);

plot(Fn,fftshift(px));

grid on;

图4 自相关法不加窗 

图5 自相关法不加窗 

图6 自相关法使用汉明窗( Hamming )

说明：

(1) 该方法先由序列估计出自相关函数，然后对进行傅里叶变换，便得到的功率谱估计。当延迟与数据长度之比很小时，可以有良好的估计精度。

(2) 图4是用自相关法（BT法）求出的功率谱，没有加窗；图5也是用自相关法（BT法）求出的功率谱，，没有加窗；图6同样是采用自相关法求出的功率谱，，使用了汉明窗。显然，自相关函数的延迟越小，谱变得越平滑。

2.3 Welch法

该方法的基本原理是在对随机序列分段时，使每一段有部分重叠，然后对每一段数据用一个合适的窗函数进行平滑处理，最后对各段谱求平均。这样可得功率谱

其中

    这里为窗函数。

该方法基于Matlab实现的程序：

clear all;

load test x;

N=4096;

Fn=-0.5:1/N:0.5-1/N;

xpsd=pwelch(x,hamming(33),16,N,'whole');

mmax=max(xpsd);%归一化

xpsd=xpsd/mmax;

xpsd=10*log10(xpsd+0.000001);

plot(Fn,fftshift(xpsd));

grid on;

图7 Welch法 不叠合 使用汉明窗( Hamming )

图8 Welch法 叠合16点 使用汉明窗( Hamming )

图9 Welch法 叠合16点 使用矩形窗( Boxcar )

图10 Welch法 叠合16点 使用布莱克曼窗( Blackman )

说明：

(1) 因为Welch法允许各段数据交叠，所以数据段数会增加，使方差得到更大的改善，但是数据的交叠减小了每一段数据的不相关性，使方差的减小不会达到理论程度。另外，采用合适的窗函数可以减少信号的频谱泄露，同时也可以增加谱峰的宽度，从而提高分辨率。

(2) 图7是利用Welch法求出的周期图，共分四段，每段32点，没有叠合，使用了汉明窗；图8也是利用Welch法求出的周期图，共分四段，每段32点，使用了汉明窗；图9是利用Welch法求出的周期图，共分四段，每段32点，交叠数为16，且使用了矩形窗；图10是利用Welch法求出的周期图，共分四段，每段32点，交叠数为16，使用了布莱克曼窗。从图中可以看出，由矩形窗处理的谱估计的主瓣宽度最窄，分辨率最好，但是其起伏性较大，所以其方差特性最差。由汉明窗和布莱克曼窗得到的谱估计的主瓣宽度最宽，因此其分辨率相对较差，但其旁瓣较小，大大改善了由矩形窗处理的谱估计旁瓣较大所产生的谱失真。因此，选择不同的窗函数其主瓣宽度不一样，造成谱估计的分辨率也不相同。

2.4 经典功率谱估计的性能比较

由以上的Matlab仿真图形和相关结果分析，我们得到了经典谱估计算法性能的直观比较：

(1) 周期图法得到的功率谱分辨率最高，但是方差性能最差，功率谱起伏剧烈，容易出现虚假谱峰。

(2) 自相关法（BT法）由于使用了平滑窗对周期图法估计的功率谱进行了平滑，因此方差性能较好，功率谱比周期图法估计的要平滑，但其分辨率比周期图法低。

(3) Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的，估计的功率谱也最为平滑，但这是以分辨率的下降及偏差的增大为代价的。

综合上述讨论，我们可以对经典谱估计的算法作大致的总结[3]：

(1) 功率谱估计，不论是直接法还是间接法都可以用FFT快速计算，且物理概念明确，因而仍是目前较常用的谱估计方法。

(2) 谱的分辨率较低，它正比于，是所使用的数据长度。

(3) 方差性能不好，不是真实功率谱的一致估计，且增大时，功率谱起伏加剧。

(4) 周期图的平滑和平均是和窗函数的使用紧密关联的，平滑和平均主要是用来改善周期图的方差性能，但往往又减小了分辨率且增加了偏差，没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善，因此使用窗函数只是改进估计质量的一个技巧问题，并不能从根本上解决问题。

3 现代功率谱估计

由前一章的讨论我们可知，经典功率谱估计方法的方差性能较差，分辨率较低。而现代谱估计技术的目标都是旨在努力改善谱估计的分辨率。参数模型法是现代谱估计的主要内容，参数模型主要分为AR模型、MA模型和ARMA模型。由于AR模型具有一系列好的性能，因此是被研究最多并获得广泛应用的一种模型。本报告中现代功率谱估计的仿真基于的是AR模型。

3.1 自相关法

假定观察到得数据为，而对于无法观察到得区间（即），的样本假定为0，观测数据区间之外的数据为0，在均方误差意义下使得数据的预测误差功率最小。由于自相关矩阵是Toeplitz矩阵，而且又为正定的，故可利用Levinson-Durbin递归算法高效求解，得到AR模型参数。

该方法基于Matlab实现的程序：

clear all;

load test x;

N=4096;

fn=-0.5:1/N:0.5-1/N;

xpsd=pyulear(x,20,N);

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

plot(fn,fftshift(xpsd));

grid on;

图11 自相关法 

图12 自相关法 

图13 自相关法 

说明：

(1) 图11、12和13是用自相关法求出的AR谱曲线，阶次分别等于10，20和30。可以看出，在阶次较低时（图11），分辨率和检测能力均不好。当时，和处的两个正弦刚刚可以分开，在和处的两个正弦也可以检出。因此必须通过提高阶次来达到分辨出间隔较小的频率点的效果。

(2) AR模型的自相关法等效于对前向预测的误差序列前后加窗，加窗的结果是使得自相关法的分辨率降低。数据越短，分辨率越不好。

3.2 协方差法

协方差法与自相关法的区别主要在于预测误差功率求和式的上下限取得不同。由于协方差法对于观察区间外样本并未假定为0，故预测误差功率表达式中的总是落在观察区间中，为此预测误差功率的求和上下限必须在之间。但由此得到的自相关矩阵是对此的半正定矩阵，且不具有Toeplitz性质，故不能采用Levinson-Durbin递归算法求解，因此得到的AR模型可能不稳定。

该方法基于Matlab实现的程序：

clear all;

load test x;

N=4096;tn=-0.5:1/N:0.5-1/N;

xpsd=pcov(x,10,N);

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

plot(tn,fftshift(xpsd));

grid on;

图14 协方差法 

说明：

(1) 可以看到谱图在信号源频率处：谱线狭窄突出，其他处谱线起伏较为剧烈。

(2) 采用协方差法对信号进行建模，能够较好地反映出信号真正的模型。

3.3修正的协方差法

AR谱估计的协方差算法基于的是最小化前向预测误差。而修正的协方差算法基于的是最小化前向和后向预测误差。这样使得它的误差功率的计算是在相对于协方差法多一倍的数据点上进行，这在观察数据长度很短的情况下，是非常有利的，但这要求信号在正反两个方向上呈现相同的特性。此外，由此得到的自相关矩阵不具有Toeplitz性质，故其正则方程不能采用Levinson-Durbin递归算法求解。

该方法基于Matlab实现的程序：

clear all;

load test x;

N=4096;

fn=-0.5:1/N:0.5-1/N;

xpsd=pmcov(x,10,N);

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

plot(fn,fftshift(xpsd));grid on;

图15 修正的协方差法 

图16 修正的协方差法

说明：

(1) 修正的协方差法较协方差法而言，谱估计图大致相同，但前者在信号源频率处的谱峰更加突出、尖锐，易于辨别。

(2) 通过图15和16的对比，我们可发现，在阶次较高的情况下能得到非常满意的结果。

3.4 Burg法

Burg算法是较早提出的建立在数据基础上的AR系数求解的有效算法。它基于最小化前向后向预测误差的同时满足Levinson-Durbin递归。对比其它的AR估计方法，Burg法避免了对自相关函数的计算，改而直接估计反射系数。

对于短数据的估计，Burg法求出的AR功率谱密度估计非常逼近于真值。另外，它能确保产生一个稳定的AR模型，并且能高效计算。

Burg法由于具有上述优点，所以分辨率比自相关法高，但对于混有白噪声的正弦信号，有时可能会出现谱线现象。

该方法基于Matlab实现的程序：

clear all;

N=4096;

fn=-0.5:1/N:0.5-1/N;

xpsd=pburg(x,10,N);

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

plot(fn,fftshift(xpsd));

grid on;

图17 Burg法 

图18 Burg法 

说明：

(1) 图17和18是采用Burg算法对数据文件所做的功率谱估计，在阶次的情况下，得到了非常满意的效果。

(2) 通过图13和17的对比，我们可明显看出Burg法的分辨率比自相关法高，在阶次较低的情况下（图17）也能较好的分辨出间隔小的频率点。

3.5 MUSIC法

刚才所讨论的自相关法、协方差法和Burg法都是基于参数建模的功率谱估计，而基于非参数建模的功率谱估计也是现代功率谱估计的重要内容。该方法是基于自相关矩阵的特征分析或者特征值分解的功率谱估计，它将相关矩阵的特征向量空间分解为信号子空间和噪声子空间，由此衍生出EV( Eigenvector )算法与MUSIC( Multiple Signal Classification )算法的信号功率谱估计[4]。其中EV谱估计与MUSIC算法谱估计都是基于噪声子空间的功率谱估计。这类方法对线谱（正弦信号的谱）最合适，对检测混有白噪声的正弦信号很有效，特别是低信噪比的情况。

MUSIC估计由下面方程给出

此处是复正弦信号向量。

该方法基于Matlab实现的程序：

clear all;

load test x;

N=4096;

fn=-0.5:1/N:0.5-1/N;

xpsd=pmusic(x',10,N);

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

for i=1:N

    xxpsd(i)=xpsd(N+1-i);

end

plot(fn,fftshift(xxpsd));

grid on;

图19 MUSIC法 

说明：

(1) 从谱图上看，MUSIC算法的谱图较为圆滑，在信号源频率处，谱线峰值明显，但峰带频率跨度较大。

(2) 通过13和19的对比，我们可明显看出MUSIC法的分辨率要好于AR模型的自相关法，在阶次较低的情况下也能较好的分辨出间隔小的频率点。

4 总结

本报告首先讨论了经典功率谱估计中三种主要方法的定义、算法及估计性能，其次讨论了现代功率谱估计中的主要内容——参数模型法谱估计，最后简要介绍了基于相关阵特征分解的谱估计算法。对所有的谱估计算法，报告均给出了Matlab程序、仿真结果及相应的分析。

参 考 文 献

[1]  王春兴．基于MATLAB实现经典功率谱估计 [J] ．曲阜师范大学学报，2011,37(2) ．

[2]  宋宁，关华．经典功率谱估计及其仿真[J]．现代电子技术，2008,31(11)．

[3]  胡广书．数字信号处理[M]．北京：清华大学出版社，2003．

[4]  王福杰，潘宏霞．MATLAB中几种功率谱估计函数的分析比较与选择[J]．计算机科学与技术，2009,27(6)．