所谓频谱分析，又称为功率谱分析或者功率谱密度（Power Spectral Density, PSD）分析，实际就是通过一定方法求解信号的功率power随着频率变化曲线。笔者在这里对目前常用的频谱分析方法做一个总结，并重点介绍目前EEG分析中最常用的频谱分析方法，并给出相应的Matlab程序。
1.频谱分析的方法有哪些？
目前来说，功率谱分析的方法大致可以分为两大类：第一类是经典的功率谱计算方法，第二类是现代功率谱计算方法，如图1所示。
其中第一类经典功率谱分析方法，又可以分为直接法、间接法和改进的直接法。直接法又称之为周期图法，简单地说，其直接利用信号的傅里叶变换系数的幅度平方来计算信号的功率谱。间接法又称为自相关函数法，其先估算出信号的自相关函数，然后对自相关函数求傅里叶变换从而得到信号的功率谱。改进的直接法，是针对直接法存在的缺点改进而来的方法，包括Barlett法、Welch法和Nuttall法。
第二类现代功率谱计算方法，又可以分为基于参数建模的功率谱计算和基于非参数建模的功率谱计算。基于参数建模的功率谱计算方法又分为基于AR模型、MA模型、ARMA模型等方法；基于非参数建模的功率谱计算方法主要基于矩阵特征分解的功率谱估计，主要包括基于MUSIC算法的功率谱估计和基于特征向量的功率谱估计。
本文，笔者主要对经典功率谱分析方法中的直接法（周期图法）以及在EEG频谱分析中最常用的改进直接法（Welch法）进行详细的介绍，并给出相应的Matlab程序。
图12.直接法计算PSD
直接法又称之为周期图法，是由Schuster于1899年提出，其方法很简单：对于长度为N的序列x(n)，其傅里叶变换为X(k),那么x(n)在每个频率（或k）处的功率可以表示为
即傅里叶变换系数的幅值平方除以数据长度N（注意：上述公式是Power或者功率的计算方法，对于功率谱密度PSD来说，还需要除以采样率！）。
根据直接法求解PSD的定义，可以直接通过调用Matlab中的fft函数（fft函数是计算信号的傅里叶变换）进行计算；
此外，Matlab中有专门的函数periodogram实现直接法的PSD计算。
例1：按照直接法计算PSD的定义，利用Matlab中的fft函数直接计算信号y=5sin(2pif1t)+3cos(2pif2*t)+δ（其中f1=20Hz，f2=35Hz，δ为一随机噪声）的PSD。
结果如图2所示，Matlab程序可以在公众号后台输入“PSDcode”进行下载，下载后可以直接在Matlab中运行出以下结果。
Matlab中有专门的函数periodogram，也可以实现直接法的PSD计算，关于其用法，这里笔者就不再赘述，各位可以直接在Matlab中help一下这个函数即可，其使用方法很简单。
图2
3.直接法计算PSD
上述直接法计算PSD，虽然可以直接FFT，计算速度快，但是频率分辨率比较低。因此，研究者提出了改进的直接法来计算PSD。**其中Welch方法就是其中的一种。Welch方法的思路是：**先把长度为N的信号分成L段，每段数据长为M，则N=LM；然后把窗函数w加到每段数据上，求出每段数据的功率谱；最后对每段数据的功率谱进行平均，得到整个信号的功率谱。各个数据段之间可以有重叠，窗函数w可以选择如Hanning、Hamming等任意一种窗口。
Welch方法是EEG的PSD计算中最常用的一种方法，在Matlab中直接采用pwelch函数可以实现Welch方法的功率谱估计，其一般调用格式如下：
[Pxx,F] = pwelch(X,WINDOW,NOVERLAP,NFFT,Fs)
[Pxx,F] = pwelch(X,WINDOW,NOVERLAP,NFFT,Fs,REQRANGE)
关于函数中的参数含义，各位可以在Matlab命令窗口中输入“help pwelch”即可查看其详细用法，这里笔者不再赘述。
例2：采用Welch方法计算信号y=5sin(2pif1t)+3cos(2pif2*t)+δ（其中f1=20Hz，f2=35Hz，δ为一随机噪声）的PSD。
计算结果如图3所示，与图2利用直接法求得的PSD基本一样；Matlab程序可以在公众号后台输入“PSDcode”进行下载，下载后可以直接在Matlab中运行出以下结果。
图3
4.总结
本文首先对目前进行PSD计算的不同方法进行了总结和简单介绍，重点详细介绍了如何利用直接法和改进的直接法（Welch法）来计算信号的PSD，并给出了Matlab程序。其中Welch法是目前计算EEG的PSD最常用的方法之一，理解和学会使用Welch法进行频率分析对于我们做EEG研究来说至关重要。
