Python中MNE库进行PSD分析（计算不同频率区间的累加和）

使用的代码和数据：https://download.csdn.net/download/zhoudapeng01/12545345

在做脑波数据分析的时候，免不了需要进行频率域的数据分析，功率谱密度是常用的一种分析方法，在MNE库中有psd算法的实现。开始使用之前推荐个博客，https://zhuanlan.zhihu.com/p/49328001 讲解了什么是PSD。有一点需要注意，PSD通常根据频率分辨率做归一化也就是其分母为采样点数。

对应的官方实例：https://mne.tools/stable/auto_tutorials/time-freq/plot_sensors_time_frequency.html#sphx-glr-auto-tutorials-time-freq-plot-sensors-time-frequency-py

本文是以BCI竞赛的数据为例，文章中的数据和代码如下：

https://download.csdn.net/download/zhoudapeng01/12545345

数据是挑选有效epoch后保存为fif格式的，有关保存epoch格式的数据详情见博客https://blog.csdn.net/zhoudapeng01/article/details/103893014

这种方法有一个问题就是只有绘制的图像，不能针对相关的数据进行后续处理。实际上mne中还有其他一些细节计算的方法，比如：mne.time_frequency.psd_welch()，mne.time_frequency.psd_multitaper() 等，其在官方的实例中均有体现。这里以用psd_multitaper计算不同频率区间的累加和为例：主要的目的是计算不同事件在不同频率区间的能量和，这里只分析一个通道的结果。

import mne
import numpy as np
import matplotlib.pyplot as plt
from mne.time_frequency import psd_multitaper

mne.set_log_level(False)
# 设置分析频率的区间范围，这里有小数主要是为了和后面小波变换方法做对比（暂时未实现）。
iter_freqs = [
    {'name': 'Delta', 'fmin': 0, 'fmax': 3.75},
    {'name': 'Theta', 'fmin': 3.75, 'fmax': 7.5},
    {'name': 'Alpha', 'fmin': 7.5, 'fmax': 12.5},
    {'name': 'Beta',  'fmin': 12.5, 'fmax': 35},
]
# 设置不同事件对应的颜色
color_events = {'cueLeft': (1, 0, 0), 'cueRight': (0, 1, 0), 'cueFoot': (0, 0, 1), 'cueTongue': (1, 1, 0)}

#################psd计算频率区间的能量分布(单通道)############
def EpochPSDEnergy(epochs):
    # 遍历epochs中的不同事件
    for event in epochs.event_id:
        # 计算某个事件对应所有epochs的功率谱密度（返回为ndarray）
        psds, freqs = psd_multitaper(epochs[event], n_jobs=1)
        # 计算事件对应epochs的均值
        psds = np.squeeze(np.average(psds, axis=0))
        #初始化能量矩阵
        eventEnergy =[]
        # 遍历不同频率区间的能量和
        for iter_freq in iter_freqs:
            eventEnergy.append(np.sum(psds[(iter_freq['fmin'] < freqs) & (freqs < iter_freq['fmax'])]))
        # 绘制不同事件、不同频率区间的能量值
        plt.plot([xLabel['name'] for xLabel in iter_freqs], eventEnergy, color=color_events[event], label=event, marker='o', lw=0, ms=5)
    # 设置标题
    plt.title('PSD_SUM')
    # 设置图例
    plt.legend()
    # 绘图显示
    plt.show()


if __name__ == '__main__':
    # 读取筛选好的epoch数据
    epochs = mne.read_epochs(r'F:\BaiduNetdiskDownload\BCICompetition\BCICIV_2a_gdf\Train\Fif\A02T_epo.fif')
    # 这里只分析一个通道的psd
    epochs.pick(['EEG-Cz']).plot_psd()
    # 绘制不同区间的能量分布
    EpochPSDEnergy(epochs.pick(['EEG-Cz']))

只实现了算法，还没具体分析有什么规律。。。

PS:实际上Python中的matplotlib中是自带psd计算方法的，如果不想用mne也可以直接调用matplotlib中的方法。

https://www.osgeo.cn/matplotlib/gallery/lines_bars_and_markers/psd_demo.html

https://www.osgeo.cn/matplotlib/api/_as_gen/matplotlib.pyplot.psd.html?highlight=psd#matplotlib.pyplot.psd