功率谱分析笔记-------脑电相关

1：功率谱分析的方法介绍

功率谱分析的方法大致可以分为两大类：第一类是经典的功率谱计算方法，第二类是现代功率谱计算方法，如图1所示。其中第一类经典功率谱分析方法，又可以分为直接法、间接法和改进的直接法。直接法又称之为周期图法，简单地说，其直接利用信号的傅里叶变换系数的幅度平方来计算信号的功率谱。间接法又称为自相关函数法，其先估算出信号的自相关函数，然后对自相关函数求傅里叶变换从而得到信号的功率谱。改进的直接法，是针对直接法存在的缺点改进而来的方法，包括Barlett法、Welch法和Nuttall法。第二类现代功率谱计算方法，又可以分为基于参数建模的功率谱计算和基于非参数建模的功率谱计算。基于参数建模的功率谱计算方法又分为基于AR模型、MA模型、ARMA模型等方法；基于非参数建模的功率谱计算方法主要基于矩阵特征分解的功率谱估计，主要包括基于MUSIC算法的功率谱估计和基于特征向量的功率谱估计。

 图1

2:周期图法

定义：

取平稳随机信号 x(n)的有限个观察点 x(0)、x(1)、…x(n-1)，则X(n)的傅立叶变换为:

由于序列x(n)的离散傅里叶变换  具有周期性，因而这种功率谱 也具有周期性，称为周期图：

周期图是信号功率谱的一个有偏估值。

x(n)在每个频率（或k）处的功率即傅里叶变换的系数的幅值平方除以数据长度N。根据直接法求解PSD的定义（注：PSD是上述公式再除以信号的采样率fs），可以直接通过调用Matlab中的fft函数（fft函数是计算信号的傅里叶变换）进行计算；此外，Matlab中有专门的函数periodogram实现直接法的PSD计算。

周期图作为频率的曲线有非常大的方差。在数学上可以推导出，周期图的方差即周期图的平方值，与采样点数目无关。这意味着无法通过使用更多的数据采样来降低周期图的方差。因此，周期图通常表现出非常高的方差，这使得信号频谱峰值难以被清楚观察和精确定位。对于需要平滑频谱估计的脑电应用而言，周期图并不适合。摘自《脑电信号处理与特征提取 第5章》

matlab调用方法：

示例：

3：Welch方法

定义：

周期图法估计出的功率谱性能并不好，当是数据长度过大时，谱的曲线起伏剧烈，而当数据长度过小时，谱的分辨率又不好，因此在其基础上又产生了一系列的改进方法，共同的原则是将周期图法进行平滑，使得估计方差减小，其中以Welch方法为代表。

Welch方法是一种最为广泛的经典功率谱估计方法，该法在周期图法的基础上主要进行了两处改进：其一是将采样数据进行分段，并允许各段数据有交叠的部分；其二是每段数据加窗处理时，窗函数不局限于矩形窗，也可使用汉宁窗或海明窗等，从而可以改善谱失真的现象。

简单来理解就是：分段的每一段用周期图法得到功率谱密度，然后加权平均（笔者个人理解：分段的每一段数据进行加窗（数据同窗函数相乘）后，利用周期图法计算加窗后数据的功率谱，然后将所有段的功率谱相加平均）。

Welch方法被广泛用于估计脑电图的频谱。不过窗函数的选择对脑电频谱估计结果的影响并不大，所以这个问题在脑电分析领域很少被讨论。摘自《脑电信号处理与特征提取 第5章》

matlab调用方法：

示例：

4：频谱估计方法的比较

非参数方法在脑电信号处理中的使用更加广泛，因为它的优点与脑电信号的特征相匹配。首先，如果频谱本身较光滑，非参数方法则更适当、更合理。其次，因参数方法对噪声水平非常敏感，故当噪声大时，非参数方法则更加准确。最后，如果信号的长度足够，非参数方法可以得到较准确的频谱估计结果。由于脑电信号通常具有宽且平滑的频谱，包含大量噪声，并且数据量足够（例如，时长1s且采样率为1000Hz的脑电数据就有1000个采样点)，因此非参数方法是比较适合进行脑电分析的。当然，如果预先指定的模型准确，参数方法，如基于自回归模型的频谱估计，也可以获得满意的结果。图5.7显示并比较了本章介绍的几种常见的非参数和参数方法在估计一个脑电信号频谱时的结果。可以清楚地看出，如果为自回归模型选择适当的模型阶数（在本例中，P=20)，则基于自回归模型的频谱估计与Welch法的估计就会非常相似。

 

 

5：计算功率谱函数（matlab）

目前matlab可以计算功率谱估计的函数有：pcov, pmcov, pyulear, pmtm, pmusic, peig, pwelch, pburg ，periodogram, arburg, prony。

Matlab调用方法：

Fs=1000;  

n=0:1/Fs:1;  

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

nfft=1024;  

window=boxcar(100); %矩形窗  

window1=hamming(100); % Harmming窗  

noverlap=20; %数据无重叠  

range='onesided'; %频率间隔为[0 Fs/2]，只计算一半的频率  

  

[Pxx1,f]=pwelch(xn,[],noverlap,nfft,Fs,range);     %pwelch  

[Pxx2,f]=pwelch(xn,[],noverlap,nfft,Fs,range);     %pwelch  

[Pxx3,f]=pcov(xn,20,nfft,Fs,range);                %1.pcov  

[Pxx4,f]=pmcov(xn,35,nfft,Fs,range);               %2.pmcov  

[Pxx5,f]=pyulear(xn,35,nfft,Fs,range);             %3.pyulear  

[Pxx6,f]=pmtm(xn,5,nfft,Fs,range);                 %4.pmtm   

[Pxx7,f]=pmusic(xn,6,nfft,Fs,range);               %5.pmusic  

[Pxx8,f]=peig(xn,6,nfft,Fs,range);                 %6.peig  

[Pxx9,f]=pburg(xn,30,nfft,Fs,range)                %7.pburg  

[Pxx10,f]=periodogram(xn,[],nfft,Fs,range);        %8.periodogram  

  

%[Pxx2,f]=arburg(xn,nfft);                         %9.arburg  

%[Pxx2,f]=prony(xn,500,20);                        %10.prony  

参考：

做EEG频谱分析，看这一篇文章就够了！ - 知乎

matlab在信号与图像处理中的应用第7章new

故障诊断作业2_叨叨Echo的博客-CSDN博客_periodogram pwelch