MNE预处理脑电数据

1.mne安装和导入数据

import os
import numpy as np
import pandas as pd
import mne
import matplotlib.pyplot as plt
#可交互
%matplotlib qt 

#将数据读进来
raw = mne.io.read_raw_eeglab('data/eeg/f006.set')
#raw1 = mne.io.read_raw_edf('CHB-MIT\chb-mit-scalp-eeg-database-1.0.0\chb01\chb01_01.edf')

 

#查看raw类型
type(raw)
#type(raw1)

#画一下数据图从第20秒开始，持续20秒
raw.plot(start=20,duration=20)
plt.show()

 

 2.剔除坏导

#导入数据并剔除无用channels
raw = mne.io.read_raw_eeglab('data/eeg/f006.set')
#告诉它哪个是eeg,哪个是eog
raw.set_channel_types({'VEOG':'eog'})
raw.pick(picks='all',exclude=['HEOG','EKG','EMG','Trigger']) #选出所有的导，不包括这几个，没有记录的无用数据
raw.plot(start=20,duration=1)
plt.show()

 3.裁剪片段

#找到并裁剪出我们需要的片段，比如说标记了的事件相关信号
events=mne.events_from_annotations(raw)  #将raw里面的annotation读取出来
# print(events) #
# print(events[0][:,2]) #  0索引中所有行第三列的数值
# print(events[1]['14']) # 1索引中关键字为14的值
# print( np.where(events[0][:,2] == 3)[0])  #输出值为3对应的关键字
ind = np.where(events[0][:,2] == events[1]['14'])[0].min()  
# print(ind)
start_t = events[0][ind,0] - 3000
ind = np.where(events[0][:,2] == events[1]['18'])[0].max()
end_t = events[0][ind,0] + 6500
#print(events[0][ind,0])
print(start_t,end_t)  #从start裁到end时间

#因为crop裁剪操作不可逆，先进行copy一份进行裁剪处理。
raw_cropped = raw.copy()
#把数据从min裁剪到max
raw_cropped.crop(tmin = start_t/1000,tmax = end_t/1000)  
raw_cropped.plot(start=20,duration=1)
plt.show()


#因为数据集和视频老师的不一样，所以出现错误，自己改了长短

 

#定义坏导，在图像中点击fz
raw_cropped.info['bads']

 4.定位电极

#定位电极，用坏导周围的电极来平均坏导的值
raw_cropped.set_montage('standard_1020',on_missing='warn')  #定义eeg 10-20模板   周围的四个电极在模板里找不到，位置太下了，所以就警告

5.插值坏导

#插值坏导
raw_cropped.load_data()
raw_cropped.interpolate_bads(exclude=['F11','F12','FT11','FT12'])  #对raw对象进行插值坏导
raw_cropped.plot(start=20,duration=1,n_channels=33,block=True,title='坏导插值完成，如无误请关闭窗口')
plt.show()

 6.低通滤波

#低通滤波
raw_filter = raw_ref.copy()
raw_filter.filter(l_freq=1,h_freq=None)  #只做了低通滤波
#raw_filter.notch_filter(freq=50)  #一般做一个50hz凹陷滤波，事件干扰工频干扰
raw_filter.plot_psd(fmax=60)
plt.show(block=False)
raw_filter.plot(start=20,duratioon=1,block = True,title='低通滤波完成，准备ICA，无误请关闭窗口')

7.去伪迹，跑ICA

#去伪迹，选用ICA
#先定义ica方法
ica = mne.preprocessing.ICA(n_components=10,method='picard',max_iter=800) #10个一组，训练800次
ica.fit(raw_filter) #训练
raw_filter.load_data()
ica.plot_sources(raw_filter,show_scrollbars=True,block=True,title='请选择需要除去的成分')   #主成分分析跑出来，
plt.show()  #在图中选中需要去除的成分
# print(ica)
# raw_recons = raw_filter.copy()
# raw_recons = ica.apply(raw_recons)
# raw_filter.plot(start=20,duration=1,n_channels=33,title='ICA处理前，确认请关闭')
# raw_recons.plot(start=20,duration=1,n_channels=33,title='ICA处理前，确认请关闭')
# plt.show(block=True)

#跑完ICA,画前后对比图
print(ica)
raw_recons = raw_filter.copy()
raw_recons = ica.apply(raw_recons)
raw_filter.plot(start=20,duration=1,n_channels=33,title='ICA处理前，确认请关闭')
raw_recons.plot(start=20,duration=1,n_channels=33,title='ICA处理前，确认请关闭')
plt.show(block=True)  #关闭窗口程序可以接着往下走

8.高通滤波凹陷滤波

#高通滤波与凹陷滤波
raw_recons.filter(l_freq=None,h_freq=30)
raw_recons.notch_filter(freq=50)

9.提取epochs

#提取epochs   把感兴趣的epochs提取出来
#先将标记的注提取出来，放到events里面
events = mne.events_from_annotations(raw_recons)
events  #可以看一下，events有5种注释    tupe类型 [1]有三列，时间，不用看，对应的值。
event_dic = ['pos':events[1]['20'],'neg':events[1]['22']]  #感兴趣的是正值和负值，对应的是20和22标签的值
reject_criteria = dict(eeg=100e-6) #拒绝振幅大于100uV的eeg,振幅过大。可选可不选


#raw数据是一整段，epoch就是提取的一小段一小段数据，与你感兴趣的事件相关的那一段eeg
#                    数据读进来   取array,  感兴趣的事件，  加载数据      提取事件前200毫秒到事件后1秒     振幅过大不要
epochs = mne.Epochs(raw_recons,events[0],event_id=event_dic,preload=True,tmax=1,tmin=-0.2,reject=reject_criteria)
epochs.plot(events=events[0],block=True,title='请目视挑选出来坏EPOCHES')
plt.show()  #点击坏的片段

10.平均epochs

#对每个片段平均  对pos和neg分别做平均
#先复制
pos_epochs = epochs['pos']
neg_epochs = epoch['neg']

pos_evoked = pos_epochs.average()
neg_evoked = neg_epochs.average()
neg_evoked.plot()   #最后画出 负200毫秒到1秒之间32个导的平均值，异常的为坏导

#将32个导进行平均
# 两种epoch之间的对比，取全部32个导的平均                                                                       上正下负，
mne.viz.plot_compare_evoked(dict(positive=pos_evoked,negative=neg_evooked),combine='mean',picks=np.arange(32)，invert_y=True,legend=True,ci=True)