英文版路径 http://www.fieldtriptoolbox.org/tutorial/preprocessing_erp

中文版是笔者在做完整个流程后的心得和翻译综合而成,FT没有中文教程,太可惜了,这个代码库强大到辣眼睛,尤其是在脑信号分析领域,也支持脑机接口和实时数据传输,最期望的是FiledTrip 能出Python版本,届时真的是雄霸天下,和EEGLAB完全不在一个层次,当然eeglab也很牛。

这篇教程是针对ERP脑信号的预处理,数据有FT 提供,主要识别人或动物的正面性和负面性。

下面就根据我自己写的代码做注释,请以官方为权威。

数据内容截图

微信截图_20160825131217

 

数据读取:

FT数据读取主要有两种方式,一种是全部读取到内存,然后处理,另外一种只读取所需的,FT当然选择后面一种。


 

%% reference the channel with RM REF
cfg=[];
% cfg.datafile='s04.eeg';
% cfg.headerfile='s04.vhdr'; 注释掉的部分和用dataset效果一样,
cfg.dataset='s04.vhdr';

这一步将eeg文件读取到matlab

 

%%
cfg.demean='yes';是否需要基线校正
cfg.baselinewindow=[-0.2 0];如果要,设置基线校正的长度

这一步是对ERP信号类型做基线校正,一般选择事件发生前的200ms,


%% lowpass filter
cfg.lpfilter='yes';
cfg.lpfreq=100;

这一步就是传说中的低通过滤器 lowpass 因为EEG的有效成分频率都在100以下,所以做了100的低通过率


%% refer the channle 

cfg.refer='yes';
cfg.implicitref='REF';
cfg.refchannel={'RM','REF'};

cfg.trialfun='trialfun_affcog';
data_erp=ft_definetrial(cfg);

这一步是对全部信号进行参考电位处理,根据数据的描述
In this raw BrainVision dataset, 
the signal from all electrodes is monopolar and referenced to the left mastoid. 
We want the signal to be referenced to linked (left and right) mastoids. 
During the acquisition the 'RM' electrode (number 32) had been placed on the right mastoid. 
In order to re-reference the data (e.g. including also the right mastoid in the reference)
we added implicit channel 'REF' to the channels (which represents the left mastoid), 
and assigned two reference channels ('REF' and 'RM', channels of the left and right mastoids).
在录制数据的时候,所有电极都是单级的,并且用了左乳突为参考,
其中31号电极被放置到右乳突,系统里面记做RM,
在这里为了重新做参考点位,包含右乳突
我们加上了隐藏信道‘REF’,就是左乳突。
并把这两个信道作为参考电位。
这里的时候  笔者是不明白怎么回事的,
但是FT的代码里面基本都是这么干的,
一般把一个乳突作为参考电极(无数据或者皆为0),另外一个正常记录,
在做重新参考是将这两个电位平均即可,不同系统是不同的设置
比如我们所用EGI的,基本上都是用Cz作为参考电位,记录后是没有数据的,
用全部平均avg后,cz才会有数据。



%% preprocessing 
data_raw=ft_preprocessing(data_erp);将上面配置的参数进行处理

%%这里处理的是眼电,一般情况像我们所里,只是用来观察眼点的存在,实际上没有用到这两个数据
但是有些需要水平和垂直眼电,下面就是对眼电的处理
cfg=[];
cfg.refer='yes';是否需要参考点位
cfg.channel={'53','LEOG'};应用的channel,注意只选这两个
cfg.refchannel={'53'}; 53作为参考电极,53自己参考自己,其实就是0电位  LEOG减去了53的值
cfg.implicitref=[];这里是空的 
eogv=ft_preprocessing(cfg,data_raw);

%% renmae eogv
cfg=[];
cfg.channel={'LEOG'};
eogv=ft_selectdata(cfg,eogv);
eogv.label={'eogv'};


%% eogv
cfg=[];水平眼电
cfg.refer='yes';
cfg.channel={'57','25'};57作为参考点位   25的值被减去了57的值
cfg.refchannel={57};
cfg.implicitref=[];
eogh=ft_preprocessing(cfg,data_raw);

%% renmae eogv 
cfg=[];
cfg.channel={'25'};
eogh=ft_selectdata(cfg,eogh);
eogh.label={'eogh'};

这样就是完成了 LEOG 和 25channel的眼电计算,被用过的57 和53从数据里面剔除。


cfg=[];
cfg.channel=setdiff(1:60,[53 57 25]);
data=ft_selectdata(cfg,data_raw);选择数据
data=ft_appenddata(cfg,data,eogv,eogh);重新把数据放在一起

去除杂讯
cfg=[];
cfg.layout='mpi_customized_acticap64.mat';电极的分布数据
cfg.channel=[1:60];选择有用的  
cfg.method='summary'; 
data_clean=ft_rejectvisual(cfg,data);右眼手动去除有杂音信号 如下:
微信截图_20160825143114

cfg=[];
cfg.method='viewmode';
Data_mode=ft_databrowser(cfg,data_clean);画出所有信号


cfg=[];
cfg.trials=find(data_clean.trialinfo==1);选出event
task1=ft_timelockanalysis(cfg,data_clean);计算ERP


cfg=[];
cfg.trials=find(data_clean.trialinfo==2);选出event
task2=ft_timelockanalysis(cfg,data_clean);计算ERP

cfg=[];
cfg.layout='mpi_customized_acticap64.mat';
cfg.interactive='yes'; 可以放大  zoom in 
cfg.showoutline='yes';
ft_multiplotER(cfg,task1,task2);画出ERP trial时序图
微信截图_20160825144157
% 
cfg=[];计算波形差
cfg.operation='subtract';减去
cfg.parameter='avg';求平均
difference=ft_math(cfg,task1,task2); 波形差值

cfg=[];
cfg.layout='mpi_customized_acticap64.mat';
cfg.interactive='yes';
cfg.showoutline='yes';
ft_multiplotER(cfg,difference);画出波形差 
微信截图_20160825144234


总结一下:

有几个函数要注意
cfg.trialfun     = 'trialfun_affcog';
这个参数配置,就已经将所有trial截取,执行ft_preprocessing的时候,注意过程

ft_definetrial的时候  取出了所有trials的开始时间 结束时间  和offset 
执行preprocessing  系统按照起始时间 将信号自动处理。


data        = ft_selectdata(cfg, data); 

这个命令 其实是根据最新的参数 来过滤数据 得到最新的数据。

cfg.channel = setdiff(1:60, [53, 57, 25]); 

这个函数1:60 之间的数  去除 53  57  25 

cfg.method   = 'summary';

data_clean   = ft_rejectvisual(cfg, data);

手动去除伪迹信号  其实就根据计算方差 然后删除outlier  效果杠杠的  


cfg = [];
cfg.operation = 'subtract';
cfg.parameter = 'avg';
difference = ft_math(cfg, task1, task2);

ft_math函数 根据 operation 和avg来计算task1和 task2

其实就是 60个channel 每个channel有184个trials   
每个trial为1的做平均  为2的做平均 然后向减得到的值

ft_multiplotER()画出ER时序图