[转载]功率谱密度相关方法的MATLAB实现

原文地址：功率谱密度相关方法的MATLAB实现作者：xd_fly

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法）

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

最大熵功率谱估计的目的是最大限度地保留截断后丢失的“窗口”以外信号的信息，使估计谱的熵最大。主要方法是以已知的自相关序列rxx(0),rxx(1),…,rxx(p)为基础，外推自相关序列rxx(p+1),rxx(p+2),…，保证信息熵最大。

最大熵功率谱估计法假定随机过程是平稳高斯过程，可以证明，随机信号的最大熵谱与AR自回归（全极点滤波器）模型谱是等价的。

MATLAB信号处理工具箱提供最大熵功率谱估计函数pmem，其调用格式为：

[Pxx,f,a]=pmem(x,p,Nfft,Fs,’xcorr’)

式中，x为输入信号序列或输入相关矩阵；p为全极点滤波器阶次；a为全极点滤波器模型系数向量；’xcorr’是把x认为是相关矩阵。

 

比较最大熵功率谱估计(MEM)和改进的平均周期图功率谱估计，可见，MEM法估计的功率谱曲线较光滑。在这一方法中，MEM法选定全极点滤波器的阶数取得越大，能够获得的窗口外的信息越多，但计算量也越大，需要根据情况折中考虑。

 

7 多信号分类法

MATLAB信号处理工具箱还提供另一种功率谱估计函数pmusic。该函数执行多信号分类法（multiple signal classification， Music法)。将数据自相关矩阵看成由信号自相关矩阵和噪声自相关矩阵两部分组成，即数据自相关矩阵R包含有两个子空间信息：信号子空间和噪声子空间。这样，矩阵特征值向量（Eigen vector）也可分为两个子空间：信号子空间和噪声子空间。为了求得功率谱估计，函数pmusic计算信号子空间和噪声子空间的特征值向量函数，使得在周期信号频率处函数值最大，功率谱估计出现峰值，而在其他频率处函数值最小。其调用格式为：

[Pxx,f,a]=pmusic(x,p[[,thresh],Nfft,Fs,window,Noverlap])

式中,x为输入信号的向量或矩阵；p为信号子空间维数；thresh为阈值，其他参数的意义与函数psd相同。

 

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

%%

%分段平均周期图法（Bartlett法）

%运用信号不重叠分段估计功率谱 

Nsec=256;n=0:sigLength-1;t=n/Fs;  %数据点数，分段间隔，时间序列

pxx1=abs(fft(y(1:256),Nsec).^2)/Nsec;   %第一段功率谱

pxx2=abs(fft(y(257:512),Nsec).^2)/Nsec; %第二段功率谱

pxx3=abs(fft(y(515:768),Nsec).^2)/Nsec; %第三段功率谱

pxx4=abs(fft(y(769:1024),Nsec).^2)/Nsec; %第四段功率谱

Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);   %平均得到整个序列功率谱

f=(0:length(Pxx)-1)*Fs/length(Pxx);     %给出功率谱对应的频率

%%plot(f(1:Nsec/2),Pxx(1:Nsec/2));       %绘制功率谱曲线

figure,plot(f(1:Nsec),Pxx(1:Nsec));       %绘制功率谱曲线

xlabel('频率/Hz');ylabel('功率谱 /dB');

title('平均周期图(无重叠) N=4*256');

grid on

%%

%运用信号重叠分段估计功率谱 

pxx1=abs(fft(y(1:256),Nsec).^2)/Nsec;   %第一段功率谱

pxx2=abs(fft(y(129:384),Nsec).^2)/Nsec; %第二段功率谱

pxx3=abs(fft(y(257:512),Nsec).^2)/Nsec; %第三段功率谱

pxx4=abs(fft(y(385:640),Nsec).^2)/Nsec; %第四段功率谱

pxx5=abs(fft(y(513:768),Nsec).^2)/Nsec; %第五段功率谱

pxx6=abs(fft(y(641:896),Nsec).^2)/Nsec; %第六段功率谱

pxx7=abs(fft(y(769:1024),Nsec).^2)/Nsec; %第七段功率谱

Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7); %功率谱平均并转化为dB

f=(0:length(Pxx)-1)*Fs/length(Pxx);  %频率序列

figure,plot(f(1:Nsec/2),Pxx(1:Nsec/2));   %绘制功率谱曲线

xlabel('频率/Hz');  ylabel('功率谱/dB');

title('平均周期图(重叠一半) N=1024');

grid on

 

%%

Nsec=256;n=0:sigLength-1;t=n/Fs; %数据长度，分段数据长度、时间序列

w=hanning(256);   %采用的窗口数据

%采用不重叠加窗方法的功率谱估计

pxx1=abs(fft(w.*y(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方

pxx2=abs(fft(w.*y(257:512),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方

pxx3=abs(fft(w.*y(513:768),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方

pxx4=abs(fft(w.*y(769:1024),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方

Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);    %求得平均功率谱，转换为dB

f=(0:length(Pxx)-1)*Fs/length(Pxx);  %求得频率序列

figure

subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线

xlabel('频率/Hz');ylabel('功率谱/dB');

title('加窗平均周期图(无重叠) N=4*256');

grid on

%采用重叠加窗方法的功率谱估计

pxx1=abs(fft(w.*y(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方

pxx2=abs(fft(w.*y(129:384),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方

pxx3=abs(fft(w.*y(257:512),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方

pxx4=abs(fft(w.*y(385:640),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方

pxx5=abs(fft(w.*y(513:768),Nsec).^2)/norm(w)^2; %第五段加窗振幅谱平方

pxx6=abs(fft(w.*y(641:896),Nsec).^2)/norm(w)^2; %第六段加窗振幅谱平方

pxx7=abs(fft(w.*y(769:1024),Nsec).^2)/norm(w)^2; %第七段加窗振幅谱平方

Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7);%平均功率谱转换为dB

f=(0:length(Pxx)-1)*Fs/length(Pxx); %频率序列

subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线

xlabel('频率/Hz');ylabel('功率谱/dB');

title('加窗平均周期图(重叠一半)N=1024');

grid on

%%

%4分段平均周期图法（hanning窗）

Nsec=256;

n=0:sigLength-1;

t=n/Fs;

w=hanning(256);

Pxx1=abs(fft(w.*y(1:256),Nsec).^2)/Nsec;

Pxx2=abs(fft(w.*y(257:512),Nsec).^2)/Nsec;

Pxx3=abs(fft(w.*y(513:768),Nsec).^2)/Nsec;

Pxx4=abs(fft(w.*y(769:1024),Nsec).^2)/Nsec;

Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);

f=(0:length(Pxx)-1)*Fs/length(Pxx);

figure

subplot(2,1,1)

plot(f,Pxx);

xlabel('频率/Hz');ylabel('功率谱/dB');

title('Averaged Modified Periodogram (none overlap) N=4*256');

grid

%4分段（2：1重叠）平均周期图法（hanning窗）

Nsec=256;

n=0:sigLength-1;

t=n/Fs;

w=hanning(256);

Pxx1=abs(fft(w.*y(1:256),Nsec).^2)/Nsec;

Pxx2=abs(fft(w.*y(129:384),Nsec).^2)/Nsec;

Pxx3=abs(fft(w.*y(257:512),Nsec).^2)/Nsec;

Pxx4=abs(fft(w.*y(385:640),Nsec).^2)/Nsec;

Pxx5=abs(fft(w.*y(513:768),Nsec).^2)/Nsec;

Pxx6=abs(fft(w.*y(641:896),Nsec).^2)/Nsec;

Pxx7=abs(fft(w.*y(769:1024),Nsec).^2)/Nsec;

Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);

f=(0:length(Pxx)-1)*Fs/length(Pxx);

subplot(2,1,2);plot(f,Pxx);

xlabel('频率/Hz');ylabel('Power Spectrum (dB)');

title('Averaged Modified Periodogram (half overlap) N=1024');

grid

%%

%PSD_WELCH方法

%采样频率

Nfft=256;n=0:sigLength-1;t=n/Fs; %数据长度、时间序列

window=hanning(256);   %选用的窗口

noverlap=128;   %分段序列重叠的采样点数（长度）

dflag='none';  %不做趋势处理

[Pxx,Pxxc,f]=psd(y,Nfft,Fs,window,noverlap,0.95);   %功率谱估计,并以0.95的置信度给出置信区间，无返回值是绘制出置信区间

figure

plot(f,10*log10(Pxx));  %绘制功率谱

xlabel('频率/Hz');ylabel('功率谱/dB');

title('PSD—Welch方法'); grid on

    

 

%%

%最大熵法（MEM法）

Nfft=256;n=0:sigLength-1;t=n/Fs;   %数据长度、分段长度和时间序列

window=hanning(256);   %采用窗口

[Pxx1,f]=pmem(x,20,Nfft,Fs);   %采用最大熵法，采用滤波器阶数14，估计功率谱

figure,subplot(2,1,1),plot(f,10*log10(Pxx1));   %绘制功率谱

xlabel('频率/Hz');ylabel('功率谱/dB');

title('最大熵法 Order=20原始信号功率谱');grid on

[Pxx1,f]=pmem(y0,20,Nfft,Fs);   %采用最大熵法，采用滤波器阶数14，估计功率谱

subplot(2,1,2),plot(f,10*log10(Pxx1));   %绘制功率谱

xlabel('频率/Hz');ylabel('功率谱/dB');axis([0 4000 -20 0])

title('最大熵法 Order=20滤波后的信号功率谱');grid on

%%

%%功率谱密度

%PSD_WELCH方法

Nfft=512;n=0:sigLength-1;t=n/Fs; %数据长度、时间序列

window=hanning(256);   %选用的窗口

noverlap=128;   %分段序列重叠的采样点数（长度）

dflag='none';  %不做趋势处理

[Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,noverlap,0.95);   %功率谱估计,并以0.95的置信度给出置信区间，无返回值是绘制出置信区间

figure;subplot(211);

plot(f,10*log10(Pxx));  %绘制功率谱

xlabel('频率/Hz');ylabel('功率谱/dB');

grid on;title('PSD—Welch方法的原始信号功率谱')

subplot(212)

[Pxx,Pxxc,f]=psd(y0,Nfft,Fs,window,noverlap,0.95);   %功率谱估计,并以0.95的置信度给出置信区间，无返回值是绘制出置信区间

plot(f,10*log10(Pxx));  %绘制功率谱

xlabel('频率/Hz');ylabel('功率谱/dB');axis([0 4000 -30 0])

grid on;title('PSD—Welch方法的滤波后的信号功率谱')

%%

%用多窗口法(MTM)

n=0:sigLength-1;t=n/Fs;  %数据长度、分段数据长度，时间序列

[Pxx1,f]=pmtm(x,2,Nfft,Fs); %用多窗口法(NW=4)估计功率谱

figure;subplot(2,1,1),plot(f,10*log10(Pxx1));  %绘制功率谱

xlabel('频率/Hz');ylabel('功率谱/dB');

title('多窗口法(MTM) nw=2原始信号功率谱');grid on

[Pxx,f]=pmtm(y0,2,Nfft,Fs); %用多窗口法(NW=2)估计功率谱

subplot(2,1,2),plot(f,10*log10(Pxx));  %绘制功率谱

xlabel('频率/Hz');ylabel('功率谱/dB');axis([0 4000 -80 -20])

title('多窗口法(MTM) nw=2滤波后的信号功率谱');grid on

    

%%

%采用Welch方法估计功率谱

noverlap=128;    %重叠数据

dflag='none'; Nfft=1024;

figure;subplot(2,1,1)

psd(x,Nfft,Fs,window,noverlap,dflag);  %采用Welch方法估计功率谱

xlabel('频率/Hz');ylabel('功率谱/dB')

title('Welch方法原始信号功率谱');grid on

subplot(2,1,2)

psd(y0,Nfft,Fs,window,noverlap,dflag);  %采用Welch方法估计功率谱

xlabel('频率/Hz');ylabel('功率谱/dB'),axis([0 4000 -30 0])

title('Welch方法滤波后的信号功率谱');grid on

pmusic(y0,[7,1.1],Nfft,Fs,32,16);   %采用多信号分类法估计功率谱

xlabel('频率/Hz'); ylabel('功率谱/dB')

title('通过MUSIC法估计的伪谱')