经典功率谱估计的分析与研究

1.     引言

信号的频谱分析是研究信号特性的重要手段之一，而信号又有确定性的信号和随机信号的之分，对于确定性的信号，可以用Fourier变换来考察其频谱性质，而对于广义平稳随机信号，由于一般它既不是随机的，又不能满足平方可积，严格来说不能进行Fourier变换，通常是求其功率谱来进行频谱分析。功率谱反映了随机信号各频率成份功率能量的分布情况，可以揭示信号中隐含的周期性及靠得很近的谱峰等有用信息，在许多领域都发挥了重要作用。例如，在语音信号识别，雷达杂波分析，地震勘测信号处理，水声信号处理，系统辨识中非线性系统的识别，物理光学中的透镜干涉，流体力学中的内波分析，太阳黑子活动周期研究等许多领域，发挥了重要的作用，然而，实际应用中的平稳随机信号通常是有限长的，只能根据有限长信号估计原信号的真实功率谱，这就是功率谱估计。

功率谱估计分为经典谱估计和现代谱估计。经典谱估计是将数据工作区外的未知数据假设为零，相当于数据加窗，主要方法有相关法和周期图法；现代谱估计是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱，主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的，应用最广的是AR参数模型。本文主要介绍经典功率谱估计，主要方法有：周期图法，相关函数法。

2．周期图法

    英国科学家牛顿最早给出了“谱”的概念。后来，1822年，法国工程师傅立叶提出了著名的傅立叶谐波分析理论。该理论至今依然是进行信号分析和信号处理的理论基础。

　　傅立叶级数提出后，首先在人们观测自然界中的周期现象时得到应用。19世纪末，Schuster提出用傅立叶级数的幅度平方作为函数中功率的度量，并将其命名为“周期图”。这是经典谱估计的最早提法，这种提法至今仍然被沿用。

周期图法又称直接法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。

周期图是信号功率谱的一个有偏估值；而且，当信号序列的长度增大到无穷时，估值的方差不趋于零。因此，随着所取的信号序列长 度的不同，所得到的周期图也不同，这种现象称为随机起伏。由于随机起伏大，使用周期图不能得到比较稳定的估值。因此需要改进。

为了减小随机起伏，M.S.巴特利特提出平均周期图法，即先把信号序列分为若干段，对每段分别计算其周期图，然后取各个周期图的平均作为功率谱的估值。平均周期图可以减小随机起伏，但是，如果信号序列不是足够长，由于每段序列长度变短，功率谱估值对不同频率成分的分辨能力也随之下降。另一种改进方法是将周期图与一个适当的频域窗函数相褶积，从而对周期图产生平滑作用，以减小随机起伏。加窗处理的结果虽然可以使随机起伏减小，但也会使周期图的分辨能力下降。

韦尔奇提出一种把加窗处理与平均处理结合起来的方法。先把分段的数据乘以窗函数(进行加窗处理)，分别计算其周期图，然后进行平均。韦尔奇方法是较常用的一种计算方法。为了得到较好的功率谱估值，加窗和平均处理均应兼顾减小随机起伏和保证有足够的谱分辨率两个方面。

周期图较差的方差性能促使人们研究另外的分析方法。这就是后面要讲的相关函数法。

3.相关函数法

相关图法是1958年由Blackman与Tukey首先提出的,理论基础是维纳辛钦定理,基本思想是通过改善对相关函数的估计方法,对周期图进行平滑处理以改善周期图谱估计的方差性能,亦称BT法、间接法。

具体的实现方法是先由序列x（n）估计出自相关函数R（n），然后对R（n）进行傅立叶变换，便得到x（n）的功率谱估计。当延迟与数据长度相比很小时，可以有良好的估计精度。

相关函数法又具有如下的缺点：方差很大，有可能是谱估计的质量下降，同时结果也不一定是正值，有可能失去了功率谱的物理意义。

4.MATLAB 仿真

两种功率谱估计方法都可以很容易的在MATLAB中实现仿真。

4.1周期图法的MATLAB仿真

例如，取采样频率为600HZ，N=600，FFT采样点数为512，求随机信号x(n)=cos(2*pi*40*n)+3*cos(2*pi*90*n)的功率谱

具体的仿真程序如下：

Fs=600; %采样频率

n=0:1/Fs:1;%产生含噪声的序列

xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+0.1*randn(size(n));

window=boxcar(length(xn));%矩形窗

Nfft=512;

[Pxx,f]=periodogram(xn,window,Nfft,Fs);%直接法¨

plot(f,10*log10(Pxx));

仿真结果如下：

图1:周期图法的仿真

4.2．相关函数法的MATLAB仿真

取采样频率为500，N=500，FFT采 样点数为512，求随机信号x(n)=cos(2*pi*40*n)+3*cos(2*pi*90*n) 的功率谱。

具体的仿真程序如下：

Fs=500;%采样频率

n=0:1/Fs:1;%产生含噪声的序列

xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));

nfft=512;

cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数

CXk=fft(cxn,nfft);

Pxx=abs(CXk);

index=0:round(nfft/2-1);

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

figure(1)%

plot(k,plot_Pxx);

仿真结果如下：

图2：相关函数法的仿真结果

4.3平均周期图法的MATLAB仿真

对于周期图的功率谱估计，当数据长度N 太大时，谱曲线起伏加剧，若N 太小，谱的分辨率又不好，因此需要改进。两种改进的估计法是平均周期图法和平滑平均周期图法。

Bartlett 法：Bartlett 平均周期图的方法是将N 点的有限长序列x(n)分段求周期图再平均。

Matlab 代码如下：

fs=600;

n=0:1/fs:1;

xn=cos(2*pi*20*n)+3*cos(2*pi*90*n)+randn(size(n));

nfft=512;

window=hamming(nfft); %矩形窗

noverlap=0;%数据无重叠

p=0.9;%置信概率

[Pxx,Pxxc]=psd(xn,nfft,fs,window,noverlap,p);

index=0:round(nfft/2-1);

k=index*fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

plot_Pxxc=10*log10(Pxxc(index+1));

figure(1)

plot(k,plot_Pxx);

所得到的仿真结果如下：

图3.平均周期图法所得到的仿真结果

5．常见的功率谱法的比较

通过实验仿真可以直观地看出以下特性:功率谱估计中的相关函数法和周期图法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。

平均周期图法的收敛性较好，曲线平滑，估计的结果方差较小，但是功率谱主瓣较宽，分辨率低。这是由于对随机序列的分段处理引起了长度有限所带来的Gibbs 现象而造成的。

6.结束语

近十几年来，随着计算机的广泛应用，各种新算法不断提出，谱估计得到了迅速的发展，出现了许多功率谱估计的实现方法。

参考文献

[1]宋宁，关华.经典功率谱估计及其仿真[J].现代电子技术，2008.