MNE_EEG

2024-05-15 EEG►ERSP ERSP Comments

目的：介绍用python的MNE预处理EEG数据，并绘制ERSP

EEG analysis with MNE-Python
This tutorial covers the basic EEG pipeline for event-related analysis: loading data, preprocessing, epoching, averaging, plotting, time-frequency analysis, and estimating cortical activity from sensor data.

一、实验范式

二、预处理

1. 加载package

import os
import mne
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mne.preprocessing import EOGRegression
from scipy.ndimage import gaussian_filter1d
from mne.time_frequency import tfr_morlet

import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter('ignore')

%matplotlib qt

2. 加载数据

path = '.\\EEG\\set\\'
fileName = "clench_right"
raw = mne.io.read_raw_eeglab(path + fileName + '.set') 
raw.info

3. 设置眼电、心电、去除不要的channel

这里根据你自己的电极来设定

1
2
3

pd.DataFrame(raw.ch_names)                               # 查看通道名称
raw.set_channel_types({'HEOG':'eog', 'VEOG':'eog'})      # 32导联，VEO垂直眼电，HEO水平眼电，ECG心电 # 
raw.pick(picks='all', exclude=['TRIGGER', 'ECG', 'EMG']) # 去除无用电极

4. 查看raw数据波形，关注信噪比，该数据是否可用【坏导插值】

4.0 查看raw波形，绘制psd

raw.plot(start=10, duration=1)  # 从10s开始，窗口显示时长1s
plt.show()
raw.plot_psd(fmin=0, fmax=250)
plt.show()

4.1 挑选坏导

把认为是坏导的，鼠标左键点击一下变灰即可

1 2	raw.plot(start=10, duration=1, title='请选择需要插值的坏导') # 从20s开始，窗口显示时长1s plt.show()

4.2 插值坏导

1
2

raw.load_data()  # 默认情况下，MNE不会将数据加载到主内存以节约资源。插值需要加载原始数据。在构造函数或raw.load_data()中使用preload=True（或字符串）。
raw.interpolate_bads()

5. 带通滤波：2 ~ 40 *

1
2
3

raw_filter = raw.copy()
raw_filter.load_data()
raw_filter.filter(l_freq=2, h_freq=40, n_jobs=1)

6. 重参考 *

1	raw_ref = raw_filter.copy()

6.1 双侧乳突参考

因为我在采集数据时，用的参考电极是M2，这里再次选择双侧乳突参考的原理参考这篇文章

''' 1. 添加 M2，值为0 '''
# add new reference channel (all zero) 添加一个新的参考电极，值为0，这样选择M1，M2重参考，就相当于 1/2 M1啦
raw_new_ref = mne.add_reference_channels(raw_ref, ref_channels=["M2"])
# raw_new_ref.plot()
''' 2. 双侧乳突重参考''' 
raw_new_ref.set_eeg_reference(ref_channels=['M1', 'M2'])

6.2 平均参考

平均参考的话，一般放到ICA之后

1	raw_ref.set_eeg_reference(ref_channels='average')

7. ICA去眼电 *

7.1 选取需要去除的ICA成分

''' 1. ICA，独立成分分析去除眼电和心电伪迹'''
ica = mne.preprocessing.ICA(n_components=20, max_iter=800)  # 定义一个ICA方法
ica.fit(raw_new_ref)                                         # 训练ICA模型
raw_new_ref.load_data()
ica.plot_sources(raw_new_ref, show_scrollbars=True, title='请选择需要去除的成分')  
ica.plot_components()
''' 2. 对ICA之后的成分进行绘制，跟前面选择坏导方法一致 ''' 
plt.show()

7.2 ICA成分剔除

''' 3. 查看选择去除的成分个数 ''' 
print(ica)
raw_ICA = raw_new_ref.copy()
''' 4. 去除刚刚选择的伪迹成分 '''
raw_ICA = ica.apply(raw_ICA)  
''' 5. ICA去除前后波形比较''' 
raw_new_ref.plot(start=1, duration=5, title='ICA处理前，如无误请关闭窗口')
raw_ICA.plot(start=1, duration=5, title='ICA处理后，如无误请关闭窗口')
plt.show()

8. 通道定位 *

这一步一般不需要，只不过是因为我的电极重新布置了，所以得重新定位
自定义电极位置，需要自定义电极布局文件对数据进行通道定位可以查看这篇文章。

# 查看通道名称
raw_ICA_rename = raw_ICA.copy()
raw_ICA_rename.pick(picks='all', exclude=['M1', 'M2','HEOG','VEOG'])
# pd.DataFrame(raw_ICA_rename.ch_names)
raw_chName = raw_ICA_rename.ch_names

rename_chName = {}
for ch in raw_chName:
    if ch == 'T7':
        rename_chName[ch] = 'C5'
    elif ch == 'TP7':
        rename_chName[ch] = 'CP5'
    elif ch == 'T8':
        rename_chName[ch] = 'C6'
    elif ch == 'TP8':
        rename_chName[ch] = 'CP6'
    else:
        rename_chName[ch] = ch

''' 1. 通道重命名 '''
raw_ICA_rename.rename_channels(rename_chName)

''' 2. 加载定位文件'''
myStardart1020 = pd.read_excel('sensor_dataFrame_32.xlsx', index_col=0) 
 
ch_names = np.array(myStardart1020.index)       # 电极名称
position = np.array(myStardart1020)             # 电极坐标位置
sensorPosition = dict(zip(ch_names, position))  # 制定为字典的形式
montage = mne.channels.make_dig_montage(ch_pos=sensorPosition)
''' 3. 通道定位 '''
raw_ICA_rename.set_montage(montage)

9. 保存预处理后的数据

1 2	sub = 'sub1_' raw_ICA_rename.save(".//1_preProcessData//"+ sub +fileName+"_M1reRef_preP-eeg.fif", overwrite = False)

三、提取epoch，evoked *

1. 绘制events图

events, event_dict = mne.events_from_annotations(raw_ICA_rename)
event_id = {'Fix':event_dict['1'], 'task':event_dict['3'],'Rest':event_dict['5']}  # marker类型由自己的范式任务决定
plt.rcParams['figure.figsize']=(7, 5)
mne.viz.plot_events(events, event_id=event_id, sfreq=raw.info["sfreq"])

2. 根据marker提取epoch

events = mne.events_from_annotations(raw_ICA_rename)
event_dic = {'Fix':events[1]['1'], 
             'task':events[1]['3'],
             'Rest':events[1]['5']} 
# reject_criteria = dict(eeg=100e-6)  # 幅度大于100微幅的不要，这个可要可不要，即振幅过大就去除
epochs = mne.Epochs(raw_ICA_rename, 
                    events[0], 
                    event_id=event_dic, 
                    preload=True,
                    tmax=20, tmin=-10,          # 刺激前10s，刺激后20s
                    # reject=reject_criteria,
                    baseline = None,            # 要做ERSP，提取epoch时不需要基线矫正  baseline= (-2,0)
                   ) 
epochs

3. 目视并手动去除伪迹依然严重的epoch

1
2
3

task_epochs = epochs['task']
task_epochs.plot(title="请目视挑出坏epochs")
plt.show()

4. 绘制evoked时域波形

1
2
3

task_evoked = task_epochs.average()
task_evoked.plot(picks=['C3','C4'], xlim=(-2,2))
plt.show()

5. 绘制evoked的psd

task_evoked.crop(0,12).plot_psd(fmin=2, fmax=30,
                     picks=['C3', 'C4'],
                     # method = 'welch',
                     xscale='linear',
                     dB=False,
                     estimate='amplitude',

                    )
plt.title(sub + fileName)

6. 绘制evoked某一时刻的脑地形图

只能选择某一时刻
参考：https://mne.tools/dev/auto_examples/visualization/evoked_topomap.html#sphx-glr-auto-examples-visualization-evoked-topomap-py

1 2	times = np.arange(-2, 14, 2) task_evoked.plot_topomap(times, ch_type="eeg", ncols=4, nrows="auto")

四. 绘制ERSP *

要计算ERSP，提取epoch时不需要基线矫正。
ERSP：反映了信号中不同频率的功率相对于特定时间点（例如信号出现）的变化程度。
事件相关同步化/去同步化：event-related synchronization / desynchronization （ERS/ERD）

ERD 对应于特定频带内的功率相对于基线的降低。类似地，ERS 对应于功率的增加。
ERDS 图是 ERD/ERS 在一定频率范围内的时间/频率表示形式。
ERDS 地图也称为 ERSP（事件相关频谱扰动：event-related spectral perturbation， ERSP）。
参考：https://zh-1-peng.gitbook.io/eeg-analysis-note/ersp

1. 定义小波函数计算时频谱

参考：https://mne.tools/stable/auto_tutorials/time-freq/20_sensors_time_frequency.html#sphx-glr-auto-tutorials-time-freq-20-sensors-time-frequency-py

from mne.time_frequency import tfr_morlet
freqs = np.logspace(*np.log10([2, 40]), num=80)
n_cycles = freqs / 2.0  # different number of cycle per frequency
# 1. 定义函数
def getPSD(dataEpochs):
    power1, itc = tfr_morlet(
            dataEpochs,
            freqs=freqs,
            n_cycles=n_cycles,
            use_fft=True,
            return_itc=True,
            decim=3,
            n_jobs=1,)
    return power1, itc

# 2. 计算
task_power, itc = getPSD(task_epochs)

2. 绘制单通道ERSP

subtracting the mean of baseline values (‘mean’)
dividing by the mean of baseline values (‘ratio’)
dividing by the mean of baseline values and taking the log (‘logratio’)
subtracting the mean of baseline values followed by dividing by the mean of baseline values (‘percent’)
subtracting the mean of baseline values and dividing by the standard deviation of baseline values (‘zscore’)
dividing by the mean of baseline values, taking the log, and dividing by the standard deviation of log baseline values (‘zlogratio’)

def plotTF(taskPower, tmin1, tmax1, ch, fn):
    # ch：需要绘制的导联名称
    # 绘制时频图
    taskPower.plot(
        picks = task_power.ch_names.index(ch),
        baseline=(-7, -3),                     # 这里我选择基线是（-7~-3），后面的power都减去基线的均值。表示：所有特定频带内的功率相对于基线（-7，-3）的变化程度
        mode="percent",   # 具体含义参考上面的解释
        tmin=tmin1, tmax=tmax1, 
        vmin=-1, vmax=1,
        show=False,
        colorbar=True,
        # cmap='jet',
    )
    plt.ylim([0, 40])
    plt.suptitle(ch)
    plt.savefig(".//1_preProcessData//TFA//"+fn+ch, dpi=600)   # 保存绘制的TF

# 绘制
fn = sub+ fileName
plotTF(task_power, -10, 20, 'C3', fn)

可以发现在任务时期，0~8s内，alpha频段的能量相较于基线明显下降，这是明显的ERD现象。

3. 绘制所有通道的ERSP

1	task_power.plot_topo(baseline=(-7, -3), mode='percent', vmin=-1, vmax=1, title="Average power") # mode="logratio"

4. 保存计算的时频谱

import pickle

with open('task_power.pkl', 'wb') as f:
    pickle.dump(task_power, f)

5. 加载保存好的时频谱

1 2	with open('task_power.pkl', 'rb') as f: task_power = pickle.load(f)

6. 时间进程上特定频段power的变化

def getPowerDataFrameTime(power_task, ch_pair):
    freqs = np.arange(2, 40)  # frequencies from 2-35Hz
    # vmin, vmax = -1, 1.5  # set min and max ERDS values in plot
    baseline = (-7, -3)  # baseline interval (in s)
    tmin, tmax = -7, 20
    ''' ERSD： percent： 减去基线值的平均值，再除以基线值的平均值（"百分比） ''' 
    power_task.crop(tmin, tmax).apply_baseline(baseline, mode="percent")
    
    ''' 0. 字典形式 '''
    df = power_task.to_data_frame(time_format=None, long_format=True)
    
    ''' 1. 选择感兴趣的通道 '''
    new_categories = ch_pair
    df = df[df.channel.isin(new_categories)]
    df["channel"] = df["channel"].cat.remove_unused_categories()
    df["channel"] = df["channel"].cat.reorder_categories(new_categories, ordered=True)

    ''' 2. 选择感兴趣的频段 '''
    freq_bounds = {"_": 0, "delta": 4, "theta": 8, "alpha": 13, "beta": 30, "gamma":40}
    df["band"] = pd.cut(
        df["freq"], list(freq_bounds.values()), labels=list(freq_bounds)[1:]
    )
    freq_bands_of_interest = ["alpha","beta"]  # "delta", "theta",, "beta", "gamma"
    df = df[df.band.isin(freq_bands_of_interest)]
    df["band"] = df["band"].cat.remove_unused_categories()
    return df

# 1. 选择感兴趣的通道
ch_pairs = ['C3', 'CZ', 'C4']

# 2. 获取统计的DataFrame
task_power_time = getPowerDataFrameTime(task_power.copy(), ch_pairs)

# 3. 绘图
g = sns.FacetGrid(task_power_time, row="band", margin_titles=True)
g.map(sns.lineplot, "time", "value", "channel", n_boot=10, linewidth=2)
axline_kw = dict(color="black", linestyle="dashed", linewidth=0.5, alpha=0.5)

g.map(plt.axhline, y=0, **axline_kw)
g.map(plt.axvline, x=0, **axline_kw)
g.refline(x=0, color='red')
g.refline(y=0, color='red')
g.set(ylim=(-0.9, 0.9))
g.set_axis_labels("Time (s)", "ERDS")
g.set_titles(col_template="{col_name}", row_template="{row_name}")
g.add_legend(ncol=3, loc="lower center")
g.fig.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.08)

另一种画法

g = sns.FacetGrid(task_power_time, row="band", col="channel", margin_titles=True)
g.map(sns.lineplot, "time", "value",  n_boot=10, linewidth=2)
axline_kw = dict(color="black", linestyle="dashed", linewidth=0.5, alpha=0.5)

g.map(plt.axhline, y=0, **axline_kw)
g.map(plt.axvline, x=0, **axline_kw)
g.refline(x=0, color='red')
g.refline(y=0, color='red')
g.set(ylim=(-0.9, 0.9))
g.set_axis_labels("Time (s)", "ERDS")
g.set_titles(col_template="{col_name}", row_template="{row_name}")
g.add_legend(ncol=2, loc="lower center")
g.fig.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.08)

7. 频段进程上特定时间段power的变化

def getPowerDataFrameBand(power_task, ch_pair):
    freqs = np.arange(2, 40)  # frequencies from 2-35Hz
    # vmin, vmax = -1, 1.5  # set min and max ERDS values in plot
    baseline = (-7, -3)  # baseline interval (in s)
    tmin, tmax = -7, 20
    
    ''' 基线  ''' 
    power_task.crop(tmin, tmax).apply_baseline(baseline, mode="percent")
    
    ''' 0. 字典形式 '''
    df = power_task.to_data_frame(time_format=None, long_format=True)
    
    ''' 1. 选择感兴趣的通道 '''
    new_categories = ch_pair
    df = df[df.channel.isin(new_categories)]
    df["channel"] = df["channel"].cat.remove_unused_categories()
    df["channel"] = df["channel"].cat.reorder_categories(new_categories, ordered=True)

    ''' 2. 选择感兴趣的时间段 '''
    df_mean = (df.query("time > 1 & time < 7")
                .groupby(["channel","freq"], observed=False)[["value"]]
                .mean()
                .reset_index()
                )
    return df_mean

# 计算
ch_pairs = ['C3', 'C4']
task_power_band = getPowerDataFrameBand(task_power.copy(), ch_pairs)
task_power_band

# 绘制
g = sns.FacetGrid(task_power_band, margin_titles=True)
g.map(sns.lineplot, "freq", "value", "channel",n_boot=10, linewidth=2)
axline_kw = dict(color="black", linestyle="dashed", linewidth=0.5, alpha=0.5)

g.map(plt.axhline, y=0, **axline_kw)
g.map(plt.axvline, x=0, **axline_kw)
g.refline(x=0, color='red')
g.refline(y=0, color='red')
# g.set(ylim=(None, 0.5))
g.set_axis_labels("Frequency (Hz)", "ERDS")
g.set_titles(col_template="{col_name}", row_template="{row_name}")

g.add_legend(ncol=3, loc="upper center")
g.fig.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.08)

8. 绘制ERSD Map

可以选择频段和时间区间

def plotTopoMap(taskPower, tmin1, tmax1, fn):
    # topomap_args = dict(extrapolate="local")
    ''' 1. 选择感兴趣的频段 '''
    plot_dict = dict(Alpha=dict(fmin=9, fmax=12))  
    ''' 2. 选择感兴趣的时间区间：tmin1, tmax1 '''
    topomap_kw = dict(tmin=tmin1, tmax=tmax1, baseline=(-7, -3),  mode="percent",  show=False)
    
    for (title, fmin_fmax) in plot_dict.items():
        taskPower.plot_topomap(**fmin_fmax, 
                               # cmap='Pastel1',
                               # vlim=(-0.7, 0.7),                                        
                               **topomap_kw,
                               # **topomap_args
                              )
          
    
    plt.suptitle(str(tmin1) +'_'+ str(tmax1)+'s_8-11Hz')
    plt.savefig(".//1_preProcessData//TopoMap//"+fn, dpi=600)
    
fn1 = sub + fileName
plotTopoMap(taskPower=task_power, tmin1=0, tmax1=7, fn=fn1)

以上就是从EEG预处理到绘制ERSP的全部过程，这里考虑了很多情况，步骤有点多，如果想要修改但不知道怎么改的可以给我留言哦。

本文链接： https://weisihong9.github.io/2024/05/15/MNE_EEG/

版权声明： 本博客所有文章除特别声明外，均采用 CC BY 4.0 CN协议许可协议。转载请注明出处！

WeisihongPhD in progress & BCI

Happy Researching!