FiledTrip is a matlab toolbox for MEG and EEG analysis , of Course, it mainly focus on the MEG analysis if you are a EEG user, that will give you some confusion, however , you’d better to go though slowly with the tutorials , herein I will present my way to keep  away obstacles  from me .

用中文标注将EEG信号溯源的流程,

如有问题请随时电邮至jiahua.xu@med.ovgu.de

这个教程比较早的时候做的  后来逐步发现需要几点需要改善  :

1. 没有把EGI channel position maping到 template MRI 上

2. 没有注意template mri 和 channel 以及 headmodel 的 unit 不同

3. 溯源时候 没有归一化 ,有可能导致噪音集中在脑袋中心

下面代码已经做改善,请留意其中部分说明, 采用了新的surface画图 效果很明显

 

 

数据简介:

脑部中风病人的10分钟静息态EEG数据

使用软件:

matlab+filedtrip

处理流程:

1.数据读取和剪切
cfg=[];
cfg.dataset='01_AD 20140317 0753001.RAW';%%读取原数据
cfg.trialdef.triallength=60;%%每个epoch 60s
cfg.trialdef.ntrials=Inf;%%有多少数据剪切多少trials
cfg=ft_definetrial(cfg);  %%定义数据结构


2.预处理和操作
cfg.continuous='yes';%%将数据转换成连续的
cfg.bsfilter ='yes';
cfg.bsfreq = [49 51; 99 100; 149 151];%%过滤外部电流影响
cfg.detrend = 'yes'; %%remove linear trend from the data
cfg.demean = 'yes';%% apply for the baseline correcttion
cfg.bpfilter = 'yes';
cfg.bpfreq = [1 48];%%带通过滤器,只需要1-48hz数据
cfg.reref = 'yes';
cfg.refchannel = 'all';
cfg.refmethod = 'avg';  %% common aveger method 进行数据参考
data=ft_preprocessing(cfg);%% 操作数据

3. 准备EEG channle的数据
cfg=[];
cfg.layout ='GSN128.sfp';
lay = ft_prepare_layout(cfg); 

load('F:\Daten-jxu\fieldtrip\fieldtrip-20160814\template\headmodel\standard_mri.mat');%standard MRI colin27

%% load elec and realign to the template 
elec = ft_read_sens('GSN-HydroCel-128.sfp'); 

elec=ft_convert_units(elec,'mm'); %consist with the mri mm

% get these positions in the ctf coordinate system using the transformation matrix of the mri 
%and the ft_warp_apply function.
nas=mri.hdr.fiducial.mri.nas;
lpa=mri.hdr.fiducial.mri.lpa;
rpa=mri.hdr.fiducial.mri.rpa;
 
transm=mri.transform;
 
nas=ft_warp_apply(transm,nas, 'homogenous');
lpa=ft_warp_apply(transm,lpa, 'homogenous');
rpa=ft_warp_apply(transm,rpa, 'homogenous');

%align the position of the fiducials in the electrode structure (defined with labels 'FidNz';'FidT9';'FidT10') 
%to their ctf-coordinates that we acquired from the anatomical mri (nas, lpa, rpa).
fid.elecpos = [nas; lpa; rpa]; % ctf-coordinates of fiducials
fid.label = {'FidNz';'FidT9';'FidT10'}; % same labels as in elec 
fid.unit = 'mm'; 


cfg = [];
cfg.method = 'fiducial'; 
cfg.target = fid; % see above
cfg.elec = elec;
cfg.fiducial = {'FidNz';'FidT9';'FidT10'}; % labels of fiducials in fid and in elec
elec = ft_electroderealign(cfg); %%这里就是将EGI channel position 放到template Mri的最后点电极分布

因为后续的画图都需要电极的坐标,所以根据不同系统有不同坐标文件,具体在fieldtrip/templete/

4. 重新抽样和选择感兴趣时间点
cfg=[];
cfg.resamplefs=200;
data_ressa=ft_resampledata(cfg,data);

cfg=[];
cfg.toilim=[15 60];
datafinal=ft_redefinetrial(cfg,data_ressa);

5.去除伪迹
cfg = [];
cfg.method = 'summary';
cfg.layout = lay;
data_clean = ft_rejectvisual(cfg, data_ressa);
这步骤是肉眼审查trials的数据分散程度,做ICA前一定要先做,否则严重影响ICA的准确度
% 
% %% run ica to remove the EOG EMG ECG
cfg=[];
cfg.method='fastica';%% runica
datacpm=ft_componentanalysis(cfg,data_clean);

这三行就是ICA的执行,可以选择多种方法,如 fastica  runica(from EEGlab)

6.肉眼审核ICA的component 去除EOG ECG EMG 
cfg=[];
cfg.viewmode='component';
cfg.layout=lay;
ft_databrowser(cfg,datacpm);

7.记下comp号 然后去除坏的component
cfg=[];
cfg.component=[1 2 3];
dataica=ft_rejectcomponent(cfg,datacpm);

这样我们就得到了干净的EEG data

下面进入重点 也是就是fieldTrip EEG信号溯源的过程
至于原理和方法,不在此叙述,当然简单的来说,在头皮记录的信号并不能反映内脑信号的分布情况,
我们需要构建也是常用BEM方式的volume conduct model 主要分为 brain skull scalp三层,根据其不同的导电率整合
多种算法来预估最终脑部活动区域,鉴于大部分EEG信号记录时都没有对应的mri解剖模型,fieldtrip提供了标准的
mri模型以及其他各种sourcemode skinmodel,都位于filetrip/templete/文件下,你可以直接采用也可以rebuild
FT提供的headmodel模型是使用dipole方法跑出来的BEM(bundary element method),可惜的是在windows系统上无法运行
于是我自己动手跑了一个 BEMCP模型,具体流程如下:
load('standard_mri.mat');%standard MRI 下载标准mri数据

%% volmesegment or load('standard_seg.mat') here with BEM method 

cfg = [];
cfg.output = {'brain','skull','scalp'};
segmentedmri = ft_volumesegment(cfg, mri);

cfg=[];
cfg.tissue={'brain','skull','scalp'};
cfg.numvertices = [3000 2000 1000];
bnd=ft_prepare_mesh(cfg,segmentedmri);


%% headmodel volume conduct or load('standard_bem.mat')

cfg = [];
cfg.method = 'bemcp';
headmodel = ft_prepare_headmodel(cfg, bnd);

% keep the same units with the leadfield matrix 
headmodel=ft_convert_units(headmodel,'cm');  统一单位

%% caculate the cross spectual density计算交叉频谱密度
cfg = [];
cfg.method = 'mtmfft';
cfg.output = 'powandcsd';
cfg.tapsmofrq = 1;
cfg.foilim = [10 10];
freqfinal = ft_freqanalysis(cfg, datafinal);
这里我们对apla频率进行溯源,以便观察10hz发生的脑部区域

%% Compute lead field
cfg = [];
cfg.normalize       ='yes';% 将mridata进行normallize so that the data could be well behavored
cfg.elec = elec;
cfg.headmodel = headmodel;
cfg.reducerank = 2;
cfg.grid.resolution = 1; % use a 3-D grid with a 1 cm resolution
cfg.grid.unit = 'cm';
grid= ft_prepare_leadfield(cfg);
这部使用建好的headmodel和channel location 进行定位溯源

%% Source Analysis: without contrasting condition
cfg = []; 
cfg.method = 'dics';
cfg.elec = elec;
cfg.frequency = [10 10]; 
cfg.grid = grid; 
cfg.headmodel = headmodel;
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda = 0;
sourcePost_nocon = ft_sourceanalysis(cfg, freqfinal);

上述代码使用dics的方法进行脑电信号溯源,得出的输出就是指定频率脑部的发生区域

如下图:
viresad110hzsr
上图其实是有噪音的,定位偏移到了中心,cfg.normalize       ='yes'我们计算leadfield matrix的时候将其归一化
另外 需要注意的是 如果有对比,可以不用neuro activity index ,如有没有最好用NAI来显示最后参考图
untitled
备注:
里面有几个函数需要注意:
1.

ft_volumesegment
2.
ft_prepare_mesh
3.
ft_prepare_headmodel
4.
ft_prepare_leadfield
5.
ft_sourceanalysis