功率谱估计方法综述

功率谱估计方法综述：

简介：随机信号的持续时间是无限长的，因此随机信号的总能量是无限的，因而随机过程的任意一个样本寒暑都不满足绝对可积条件，所以其傅里叶变换不存在。尽管随机信号的总能量是无限的，但其平均功率却是有限的，因此，要对随机信号的频域进行分析，应从功率谱出发进行研究才有意义。信号的功率谱密度描述随机信号的功率在频域随频率的分布。功率谱估计(PSD)是用有限长的数据来估计信号的功率谱，即利用给定的N个样本数据估计一个平稳随机信号的功率谱密度。

背景：功率谱估计在实际工程中有重要应用价值，如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域，发挥了重要作用。功率谱估计方法主要分为2大类：非参数化方法（又称经典功率谱估计）和参数化方法（又称现代功率谱估计）。非参数化方法有相关函数法(BT法)、周期图法、平均周期图法、平滑平均周期图法等；而参数化谱估计有R模型法、移动平均模型法(简称MA模型法)、自回归移动平均模型法(简称ARMA模型法)、最大熵谱分析法(AR模型法)、Pisarenko谐波分解法、Prony提取极点法、Prony谱线分解法以及capon最大似然法等，由于涉及许多复杂数学计算，在此未作详细数学推导，以下介绍几种常用的功率谱估计方法

一、非参数化方法（经典法）

经典功率谱估计是将数据工作区外的未知数据假设为零，相当于数据加窗。

1、自相关法

又称相关函数法(BT法)， 根据维纳—辛钦定理：平稳随机过程的自相关函数和功率谱函数是一傅里叶变换对，对于平稳随机信号来说，其相关函数是确定性函数，故其功率谱也是确定的．这样可由平稳随机离散信号的有限个离散值，求出自相关函数，然后作Fourier变换，得到功率谱。由于随机序列{X(n)}的自相关函数R(n)=E[X(n)X(n+m)]定义在离散点m上，设取样间隔为错误！未找到引用源。，可将随机序列的自相关函数用连续时间函数表示为

错误！未找到引用源。

等式两边取傅里叶变换，则随机序列的功率谱密度

错误！未找到引用源。 错误！未找到引用源。 错误！未找到引用源。 BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值错误！未找到引用源。。即 错误！未找到引用源。 其中错误！未找到引用源。可有式得到。

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'); %计算序列的自相关函数，matlab函数 xcorr（求自相关函数）unbiased无偏

CXk=fft(cxn,nfft); %对cxn（即自相关函数）进行快速傅里叶变换，nfft为周期 Pxx=abs(CXk); %对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、周期图法

周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得x(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例2： 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);%直接法计算功率谱密度，xn为功率谱密度信号，window为窗口，nfft为采样点数， fs采样频率 plot(f,10*log10(Pxx));

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

[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法 figure(1)

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

3、平均法

即Bartlett平均周期图法，是将N点的有限长序列x(n)分段求周期图再平均．将长度为N的数据分为L段，先对每段数据用周期图法进行谱估计，然后对L段求平均得到长度为N的数据的功率谱．

平均法可视为周期图法的改进。周期图经过平均后会使它的方差减少，达到一致估计的目的，有一个定理：如果错误！未找到引用源。是不相关的随机变量，且都有个均值错误！未找到引用源。及其方差错误！未找到引用源。，则可以证明它们的算术平均的均值为错误！未找到引用源。。

即：平均法将错误！未找到引用源。的N个数据分成L段（N=ML），若各数据段相互独立，则平方后估计量的方差是原来不分段估计量方差的错误！未找到引用源。。所以当错误！未找到引用源。时，估计量的方差趋于0，达到一致估计的目的。但是，随着分段数L的增加，M点数减少，分辨率减少，使估计变成有偏估计。相反，若L减少，M增加，虽偏差减少，但方差增大。所以，在应用中，必须兼顾分辨率和方差的要求来适当选择M和L的值。 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); figure(2)

plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

Bartlett法方差的改善是以牺牲分辨率为代价的．

4、welch法

Welch法又称修正周期法，Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并在周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。

Welch法谱估计是在上述基础上进行了改进，目的是在保持Bartlett法方差性能的同时，改善其分辨率。其基本原理是先对随机序列分段时，使每一段有部分晕叠，然后对每一段数据用一个合适的窗函数进行平滑处理，最后对各段谱求平均。 Matlab代码 Fs=600; n=0:1/Fs:1;

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

window=boxcar(100);%矩形窗 window1=hamming(100);%海明窗

window2=blackman(100);%blackman窗 noverlap=20;%数据无重叠

range='half';%频率间隔为[O Fs/2]，计算一半的频率 [Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range); [Pxxl,f]=pwelch(xn,window1,noverlap,nfft,Fs,range); [Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range); plot_Pxx=10*log10(Pxx); plot_Pxxl=10*log10(Pxxl); plot_Pxx2=10+log10(Pxx2); figure(1)

plot(f,plot_Pxx); figure(2)

plot(f,plot_Pxxl); figure(3)

plot(f,plot_Pxx2);

函数说明 （pwelch）[Pxx,F] = pwelch(x,window,noverlap,nfft,fs)

x, 为进行功率谱估计的输入有限长序列 window，用于指定采用的窗函数(boxcar, hamming, blackman) noverlap，重叠点数 nfft，设定FFT算法的长度 fs, 采样频率 Pxx，为输出的功率谱估计值 F，为得到的频率点

综述：

从给出一段信号y=cos(2*pi*30*n)+3*cos(2*pi*100*n),利用经典谱估计法的原理，通过不同的谱估计方法，求出信号的功率谱密度函数，并利用GUI界面呈现出不同谱估计方法所得的结果。

BT法：BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值错误！未找到引用源。。

由式错误！未找到引用源。=错误！未找到引用源。估算出错误！未找到引用源。，再对错误！未找到引用源。作FFT变换，得到错误！未找到引用源。。

周期图法：周期图法是根据各态历经随机过程功率谱的定义来进行谱估计的。

获取x(n)后，对x(n)作FFT得到X(W)，再由错误！未找到引用源。得出功率谱。 平均法：平均法可视为周期图法的改进。周期图经过平均后会使它的方差减少，达到一致估计的目的。

获取x(n)后，将x(n)分为10段，对每段用错误！未找到引用源。计算出周期图，对以上10个周期图加以平均得出功率谱。

Welch法：Welch法又称修正周期法，获取x(n)后，先将N个数据分成L段，选择汉明窗w(n)，并用该w(n)依次对每段数据做相应的加权，确定每段的周期图

GM,l?w??1MUMl?1n?M(l?1)2?xnw(n)e?j?n,l?1,2,```,L

1L由GM?w???GM,l?w?得出每段谱估计。

Ll?1

二、参数化方法（现代法）

AR模型功率谱估计（仅举此一种）

AR模型功率谱估计又称为自回归模型，它是一个全极点的模型，要利用AR模型进行功率谱估计须通过levinson_dubin递推算法由正则方程求得AR的参数,在Matlab仿真中可调用Pburg函数直接画出基于burg算法的功率谱估计的曲线图。用周期图法求出的功率谱曲线和burg算法求出AR功率谱曲线。

Matlab代码

用周期图法求出的功率谱曲线和burg算法求出的AR功率谱曲线(p=50) fs=200; n=0:1/fs:1;

xn=cos(2*pi*40*n)+cos(2*pi*41*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); subplot(1,2,1);

plot(f,10+log10(pxx)) xlabel('frequency(hz)');

ylabel('power spectral density(Db/Hz)'); title('periodogram pse estimate'); orderl=50;

range='onesided'; magunits='db'; subplot(1,2,2);

pburg(xn,orderl,nfft,fs,range) 总结

常见谱估计法的比较

通过实验仿真可以直观地看出以下特性：(1)功率谱估计中的相关函数法和周期图法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。(2)平均周期图法和平滑平均周期图法的收敛性较好，曲线平滑，估计的结果方差较小，但是功率谱主瓣较宽，分辨率低。这是由于对随机序列的分段处理引起了长度有限所带来的Gibbs现象而造成的。(3)平滑平均周期图法与平均周期图法相比．谱估值比较平滑，但是分辨率较差。其原因是给每一段序列用适当的窗口函数加权后，在得到平滑的估计结果的同时，使功率谱的主瓣变宽，因此分辨率有所下降。经典功率谱估计的分辨率反比于有效信号的长度，但现代谱估计的分辨率可以不受此限制。这是因为对于给定的N点有限长序列x(n)。虽然其估计出的自相关函数也是有限长的。但是现代谱估计的一些隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度，不象经典谱估计那样受窗函数的影响。因而现代谱的分别率比较高，而且现代谱线要平滑得多。

Matlab中有很多内置的信号处理的专用函数，在此举一下几个常用的信号处理函数：

信号处理函数

conv conv2 fft fft2 ifft ifft2 filter filter2 abs angle unwrap fftshift nextpow2

卷积 2维卷积

快速富里哀变换 2维快速富里哀变换 逆快速富里哀变换 2维逆快速富里哀变换 离散时间滤波器 2维离散时间滤波器 幅值

四个象限的相角

在360°边界清除相角突变 把FFT结果平移到负频率上 2的下一个较高幂次

参考文献：

印勇 随机信号处理教程 北京邮电大学出版社2010 朱冰莲 数字信号处理 电子工业出版社 2011

王春兴 基于MATLAB实现经典功率谱估计[期刊论文] 王春兴 基于MATLAB实现现代功率谱估计[期刊论文]

本文来源：https://www.bwwdw.com/article/b6io.html