FieldTrip toolbox for EEG/MEG (Matlab)

Tutorial for FieldTrip
(If you have a problem to add path of fieldtrip toolbox, command >> rehash toolbox)

1) Introduction

The knowledge or assumptions incorporated by the researcher in the protocol might be:

• At a certain time after the presentation of a stimulus, the brain will process the stimulus.
• The brain will process repeated stimuli similarly.
• Electrical activity generated by the neurons in the brain has certain properties that distinguish it from artifacts, such as amplitude (eye movements result in much larger scalp potentials) or distribution (signals originating from the brain are picked up simultaneously by multiple neighboring electrodes).
• The scalp distribution of the EEG/MEG depends on the location of the underlying activity.
• Oscillatory activity is produced by certain brain regions.

 These tools can be combined in an analysis protocol that for example looks like Figure 1.

Figure 1; Analysis protocol for Event-Related Potentials (ERPs)

• Functions for preprocessing, reading and coverting data (e.g. ft_preprocessing)
• Functions for analyzing event-related fields or potentials (ERF/ERP) (e.g. ft_timelockanalysis)
• Functions for frequency and time-frequency analysis (e.g. ft_freqanalysis)
• Functions for source analysis (e.g. ft_sourceanalysis)
• Functions for statistical analysis (e.g. ft_timelockstatistics)
• Functions for plotting and displaying the data (e.g. ft_multiplotER)
• Function for real-time acquisition and analysis (e.g. ft_realtime_average)

The full list of functions can be found here. 

2) Preprocessing - Reading continuous data (EEG)

•  read the data for the EEG channels using ft_preprocessing, apply a filter and refererence to linked mastoids
• read the data for the horizontal and vertical EOG channels using ft_preprocessing, and compute the horizontal and vertical bipolar EOG derivations
• combine the EEG and EOG into a single data representation using ft_appenddata
• determine interesting pieces of data based on the trigger events using ft_definetrial
• segment the continuous data into trials using ft_redefinetrial
• segment the continuous data into one-second pieces using ft_redefinetrial

%%%%%%%%% define  and preprocessing data %%%%%%%%%%%%%
cfg = [];

cfg.dataset = '*.eeg';

data_org = ft_preprocessing(cfg); % raw data

cfg.trialdef.eventtype = 'Stimulus';

% cfg.trialdef.eventtype  = '?'; % get the eventvalues

cfg.trialdef.eventvalue = {'S 16', 'S 32'};
% 2 stimulus

cfg.trialdef.prestim = 0.1;
cfg.trialdef.poststim = 0.8; % -100 ms ~ 800 ms

cfg_stimulus = ft_definetrial(cfg);

data_stimulus = ft_redefinetrial(cfg_stimulus,data_org);

% new data by trials

Subsequently we could do artifact detection with ft_rejectvisual to remove trials with artifacts and average the trials using ft_timelockanalysis to get the ERP, or use ft_freqanalysis to obtain an averaged time-frequency representation of the data in both conditions. If you use ft_timelockanalysis or ft_frequencyanalysis with the option cfg.keeptrials='yes', you subsequently could use ft_timelockstatistics or ft_freqstatistics for statistical comparison of the animal-tool contrast in these stimuli. 

3) Introduction: dealing with artifacts

3.1) rejecting the piece of data containing the artifact (e.g. for a short-lived artifact)   

Manual/visual detection 

• ft_rejectvisual
• ft_databrowser (prefer to use this method)

%%%%%%%%% artifact rejection %%%%%%%%%%%%%%%%%%
cfg_stimulus.method = 'trial';
cfg_stimulus.alim = 200;
% set scaling to 200 micro Vlot
data_stimulus_rejectArtifact = ft_rejectvisual (cfg_stimulus,data_stimulus); 

%%%%%%%%% artifact rejection %%%%%%%%%%%%%%%%%%

cfg.blocksize = 100; %100s per segment
cfg_reject = ft_databrowser(cfg,data_org); 
data_stimulus_rejectArtifact = ft_rejectartifact(cfg_reject, data_stimulus); 

Automatic detection

• ft_artifact_clip
• ft_artifact_ecg
• ft_artifact_threshold
• ft_artifact_eog
• ft_artifact_jump
• ft_artifact_muscle
• ft_artifact_zvalue

3.2) subtracting the spatio-temporal contribution of the artifact from the data (e.g. for line noise)

Using ICA to identify artifacts

4) Re-referencing

This can be done at first step.

%%%%%%%%%%%%%%%%%%%%% re-referencing %%%%%%%%%%%%%%%%

cfg = [];

cfg.reref = 'yes';

cfg.refchannel = {'2'}; %the channels '2'is used as the new reference

data_stimulus_rejectArtifact_reref = ft_preprocessing(cfg,data_stimulus_rejectArtifact);

5) Filtering

This can be done at first step.

%%%%%%%%%%%%%%%%% Filtering %%%%%%%%%%%%%%%%%%

cfg_filt = [];
lowPass = 250;
highPass = 1;
cfg_filt.lpfilter = 'yes';
cfg_filt.lpfreq = [lowPass];
cfg_filt.hpfilter = 'yes';
cfg_filt.hpfreq = [highPass];
data_filt = ft_preprocessing(cfg_filt, data_reref);


if (exist('data_stimulus_rejectArtifact','var')||exist('data_stimulus_rejectArtifact_reref','var'))
    if reref==1
        data_filt =    ft_preprocessing                       (cfg_filt,data_stimulus_rejectArtifact_reref);
    else
        data_filt = ft_preprocessing (cfg_filt,data_stimulus_rejectArtifact);
    end
elseif exist('data_stimulus','var')
    data_filt = ft_preprocessing  (cfg_filt,data_stimulus);    
end

6)Time-Frequency analysis

    
%%%%%%%%% Time-Frequency analysis %%%%%%%%%%%%%%%%%%
cfg = [];
lowFreq = 1;
highFreq = 30;
stepFreq = 1;
widthFactor = 1;
    
cfg.method = 'wavelet';
cfg.width = widthFactor;
cfg.output = 'pow';
cfg.foi = lowFreq : stepFreq : highFreq;
cfg.toi = -0.1 : 0.05 : 0.8;
    
cfg.trials = find(data_filt.trialinfo(:,1) == 'triggerCode');    
TFRwave_singletrigger{triggerCode} = ft_freqanalysis(cfg, data_filt); 

%%%%%%%%% plot Time-Frequency analysis %%%%%%%%%%%%% 
temp = TFRwave_singletrigger{triggerCode}.powspctrm (channel,:,:);
temp = reshape(temp,size(temp,2),size(temp,3));
                    pcolor(TFRwave_singletrigger{triggerCode}.time',TFRwave_singletrigger{triggerCode}.freq',temp);  

7) Timelock analysis

The function ft_timelockanalysis makes averages of all the trials in a data structure. It requires preprocessed data.

reference: http://fieldtrip.fcdonders.nl/tutorial/eventrelatedaveraging

%%%%%%%%%%%%% timelock analysis %%%%%%%%%%%%%%%%
cfg = [];
cfg.trials=find(data_filt.trialinfo(:,1)==triggerCode);
data_temp = ft_timelockanalysis (cfg, data_filt);

%%%%%%%%%%%%%%%% set baseline %%%%%%%%%%%%%%%%%%%
cfg = [];
cfg.baseline = [-0.1,0]; 
data_singletrigger{i}=ft_timelockbaseline(cfg,data_temp);

%%%%%%%%% plot Timelock analysis %%%%%%%%%%%%% 
plot(data_singletrigger{triggerCode}.time, data_singletrigger{triggerCode}.avg(channel,:));