功率谱估计性能分析及Matlab仿真

功率谱估计性能分析及Matlab 仿真

1 引言

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

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

2 经典功率谱估计

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

2.1 周期图法( Periodogram )

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

取平稳随机信号()x n 的有限个观察值(0),(1),...,(1)x x x n -，求出其傅里叶变换

1

0()()N j j n N n X e x n e ωω---==∑

然后进行谱估计

21()()j N S X e N

ωω-= 周期图法应用比较广泛，主要是由于它与序列的频谱有直接的对应关系，并且可以采用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 周期图法 4096N =

图2 周期图法 128N =

说明：

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

图3 复序列 (),[0,127]x n n ∈

(2) 从仿真图可以清晰看到，1f 和2f 不能完全分开，仅在波形的顶部能看出

是两个频率分量；此外，当数据长度N 太大时（图1），谱曲线呈现较大的起伏；当数据长度N 太小时（图2），谱的分辨率又不好。据此，周期图法不满足一致性估计条件。

2.2 自相关法( BT 法)

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

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

101()()()N m x n R m x n x n m N --==+∑

然后在(,)M M -内对()x R m 做傅里叶变换，得到功率谱

()()M j n x m M S R m e ωω-=-=

∑

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

clear all;

load test x;

N=4096;

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

Mlag=64;

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;

M=

图4 自相关法不加窗64

M=

图5 自相关法不加窗32

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

说明：

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

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

2.3 Welch 法

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

2

110

1()()L M i j m i n S x n e MUL ωω--===∑∑ 其中

10()M n U n ω-==∑