功率谱密度相关方法的MATLAB实现

1.基本方法

周期图法是直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法。假定有限长随机信号序列为x(n)。它的Fourier变换和功率谱密度估计存在下面的关系：

式中，N为随机信号序列x(n)的长度。在离散的频率点f=kΔf,有：

其中，FFT[x(n)]为对序列x(n)的Fourier变换，由于FFT[x(n)]的周期为N，求得的功率谱估计以N为周期，因此这种方法称为周期图法。下面用例子说明如何采用这种方法进行功率谱

用有限长样本序列的Fourier变换来表示随机序列的功率谱，只是一种估计或近似，不可避免存在误差。为了减少误差，使功率谱估计更加平滑，可采用分段平均周期图法（Bartlett法）、加窗平均周期图法（Welch 法）等方法加以改进。

2. 分段平均周期图法（Bartlett法）

将信号序列x(n)，n=0,1,…,N-1，分成互不重叠的P个小段，每小段由m个采样值，则P*m=N。对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列x(n)的功率谱估计。

平均周期图法还可以对信号x(n)进行重叠分段，如按2:1重叠分段，即前一段信号和后一段信号有一半是重叠的。对每一小段信号序列进行功率谱估计，然后再取平均值作为整个序列x(n)的功率谱估计。这两种方法都称为平均周期图法，一般后者比前者好。程序运行结果为图9-5，上图采用不重叠分段法的功率谱估计，下图为2:1重叠分段的功率谱估计，可见后者估计曲线较为平滑。与上例比较，平均周期图法功率谱估计具有明显效果（涨落曲线靠近0dB）。

3.加窗平均周期图法

加窗平均周期图法是对分段平均周期图法的改进。在信号序列x(n)分段后，用非矩形窗口对每一小段信号序列进行预处理，再采用前述分段平均周期图法进行整个信号序列x(n)的功率谱估计。由窗函数的基本知识（第7章）可知，采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。

其中上图采用无重叠数据分段的加窗平均周期图法进行功率谱估计，而下图采用重叠数据分段的加窗平均周期图法进行功率谱估计，显然后者是更佳的，信号谱峰加宽，而噪声谱均在

0dB附近，更为平坦（注意采用无重叠数据分段噪声的最大的下降分贝数大于5dB，而重叠数据分段周期图法噪声的最大下降分贝数小于5dB）。

4. Welch法估计及其MATLAB函数

Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch 法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计(PSD如上例中的下半部分的求法)和两个信号序列的互功率谱估计（CSD）。

MATLAB信号处理工具箱函数提供了专门的函数PSD和CSD自动实现Welch法估计，而不需要自己编程。

（1）函数psd利用Welch法估计一个信号自功率谱密度，函数调用格式为：

[Pxx[,f]]=psd(x[,Nfft,Fs,window,Noverlap,’dflag’])

式中，x为信号序列；Nfft为采用的FFT长度。这一值决定了功率谱估计速度，当Nfft采用2的幂时，程序采用快速算法；Fs为采样频率；Window定义窗函数和x分段序列的长度。窗函数长度必须小于或等于Nfft，否则会给出错误信息；Noverlap为分段序列重叠的采样点数（长度），它应小于Nfft;dflag为去除信号趋势分量的选择项：’linear’,去除线性趋势分量，’mean’去除均值分量，’none’不做去除趋势处理。Pxx为信号x的自功率谱密度估计。f为返回的频率向量，它和Pxx对应，并且有相同长度。

在psd函数调用格式中，缺省值为：Nfft=min(256,length(x)),Fs=2Hz,

window=hanning(Nfft),noverlap=0. 若x是实序列，函数psd仅计算频率为正的功率注意程序前半部分中频率向量f的创建方法。它与函数psd的输出Pxx长度的关系如下：若x为实序列，当Nfft为奇数时，f=(0:(Nfft+1)/2-1)/Nfft;当Nfft为偶数时，

f=(0:Nfft/2)/Nfft。

函数还有一种缺省返回值的调用格式，用于直接绘制信号序列x的功率谱估计曲线。函数还可以计算带有置信区间的功率谱估计，调用格式为：

[Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,Noverlap,p)

式中，p为置信区间，0<=p<=1。

由此可知，滤波器输入白噪声序列的输出信号的功率谱或自相关可以确定滤波器的频率特性。

（2）函数csd利用welch法估计两个信号的互功率谱密度，函数调用格式为： [Pxy[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,’dflag’)

[Pxy,Pxyc[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,p)

这里，x,y为两个信号序列；Pxy为x,y的互功率谱估计；其他参数的意义同自功率谱函数psd。

可以看到，两个白噪声信号的互功率谱（上图）杂乱无章，看不出周期成分，大部分功率谱在-5dB以下。然而白噪声与带有噪声的周期信号的功率谱在其周期（频率为1000Hz）处有一峰值，清楚地表明了周期信号的周期或频率。因此，利用未知信号与白噪声信号的互功率谱也可以检测未知信号中所含有的频率成分。

5 多窗口法

多窗口法(Multitaper method,简称MTM法)利用多个正交窗口（Tapers）获得各自独立的近似功率谱估计，然后综合这些估计得到一个序列的功率谱估计。相对于普通的周期图法，这种功率谱估计具有更大的自由度，并在估计精度和估计波动方面均有较好的效果。普通的功率谱估计只利用单一窗口，因此在序列始端和末端均会丢失相关信息，而且无法找回。而MTM法估计增加窗口使得丢失的信息尽量减少。

MTM法简单地采用一个参数：时间带宽积(Time-bandwidth product)NW，这个参数用以定义计算功率谱所用窗的数目，为2*NW-1。NW越大，功率谱计算次数越多，时间域分辨率越高，而频率域分辨率降低，使得功率谱估计的波动减小。随着NW增大，每次估计中谱泄漏增多，总功率谱估计的偏差增大。对于每一个数据组，通常有一个最优的NW使得在估计偏差和估计波动两方面求得折中，这需要在程序中反复调试来获得。

MATLAB信号处理工具箱中函数PMTM就是采用MTM法估计功率谱密度。函数调用格式为：

[Pxx[,f]]=pmtm(x[,nw,Nfft,Fs])

式中，x为信号序列；nw为时间带宽积，缺省值为4。通常可取2，5/2，3，7/2；Nfft为FFT长度；Fs为采样频率。

上面的函数还可以通过无返回值而绘出置信区间，如pmtm(x,nw,Nfft,Fs,’option’,p)绘制带置信区间的功率谱密度估计曲线，0<=p<=1。

6 最大熵法（Maxmum entropy method， MEM法）

如上所述，周期图法功率谱估计需要对信号序列“截断”或加窗处理，其结果是使估计的功率谱密度为信号序列真实谱和窗谱的卷积，导致误差的产生。