This example script has been promoted to a tutorial under the name Preprocessing of EEG data and compute ERPs. All information below might be outdated. Please check the tutorial instead.

Preprocessing of EEG data and computing ERPs

Introduction

This tutorial shows how to read and pre-process EEG data in Fieldtrip. In this example we'll be using a dataset recorded with the BrainProducts EEG system. The relevant differences with datasets from some other EEG systems are briefly discussed.

If you are interested in how to do preprocessing on your data prior to segmenting it into trials, you can check the Preprocessing - Reading continuous data tutorial. If you are interested in how to preprocess MEG data, see the Preprocessing - Trigger-based trial selection tutorial.

Background

In FieldTrip the preprocessing of data refers to the reading of the data, segmenting the data around interesting events such as triggers, temporal filtering and optionally rereferencing. The ft_preprocessing function takes care of all these steps, i.e., it reads the data and applies the preprocessing options.

There are largely two alternative approaches for preprocessing, which especially differ in the amount of memory required. The first approach is to read all data from the file into memory, apply filters, and subsequently cut the data into interesting segments. The second approach is to first identify the interesting segments, read those segments from the data file and apply the filters to those segments only. The remainder of this tutorial explains the second approach, as that is the most appropriate for large data sets such as the MEG data used in this tutorial. The approach for reading and filtering continuous data and segmenting afterwards is explained in another tutorial.

Preprocessing involves several steps including identifying individual trials from the dataset, filtering and artifact rejections. This tutorial covers how to identify trials using the trigger signal. Defining data segments of interest can be done

  • according to a specified trigger channel
  • according to your own criteria when you write your own trial function

Examples for both ways are described in this tutorial, and both ways depend on ft_definetrial.

The output of ft_definetrial is a configuration structure containing the field cfg.trl. This is a matrix representing the relevant parts of the raw datafile which are to be selected for further processing. Each row in the trl-matrix represents a single epoch-of-interest, and the trl-matrix has at least 3 columns. The first column defines (in samples) the beginpoint of each epoch with respect to how the data are stored in the raw datafile. The second column defines (in samples) the endpoint of each epoch, and the third column specifies the offset (in samples) of the first sample within each epoch with respect to timepoint 0 within that epoch.

Dataset

The EEG dataset used in this script is available here. In the experiment, subjects made positive/negative or animal/human judgments on nouns. The nouns were either positive animals (puppy), negative animals (maggot), positive humans (princess), or negative humans (murderer). The nouns were presented visually (written words). The task cue (which judgement to make) was given with each word.

Procedure

Defining trials

Make sure that all files that you have downloaded from the ftp link are unzipped and are located in the present working directory in MATLAB. In the command window, you can type pwd to see what the present directory is, and you can type dir to see the content of the working directory.

For memory efficiency (especially relevant for large MEG datasets), with FieldTrip we commonly use the strategy to only read in those segments of data that are of interest. This requires first to define the segments of interest (the trials) and subsequently to read then in and preprocess them. It is also possible to read in the whole continuous data, and segment the data in memory (see here).

Instead of using the default 'trialfun_general' function with ft_definetrial, we will use a custom 'trialfun_affcog' that has been written specifically for this experiment. This custom function reads markers from the EEG record and identifies trials that belong to condition 1 (positive-negative judgement) or 2 (animal-human judgement). The function is available along with the data.

The custom trial function is available from here or can be found at the end in the appendix of this example script. Please save it to a local file with the name trialfun_affcog.m.

cfg = [];
cfg.trialfun     = 'trialfun_affcog';
cfg.headerfile   = 's04.vhdr';
cfg = ft_definetrial(cfg);

After the call to ft_definetrial, the cfg now not only stores the dataset name, but also the definition of the segments of data that will be used for further processing and analysis. The first column is the begin sample, the second the end sample, the third the offset and the fourth contains the condition for each trial (1=affective, 2=ontological).

>> disp(cfg.trl)
ans =
       52441       53041        -100           2
       56740       57340        -100           1
       61845       62445        -100           1
       66383       66983        -100           2
       70402       71002        -100           1
       74747       75347        -100           1
       ...

Pre-processing and re-referencing

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).

Now call pre-processing:

% Baseline-correction options
cfg.demean          = 'yes';
cfg.baselinewindow  = [-0.2 0];

% Fitering options
cfg.lpfilter        = 'yes';
cfg.lpfreq          = 100;

% Re-referencing options - see explanation below
cfg.reref         = 'yes';
cfg.implicitref   = 'REF';
cfg.refchannel    = {'RM' 'REF'};

data = ft_preprocessing(cfg);

Try ft_databrowser now to visualize the data segments that were read into memory.

 
cfg = [];  % use only default options                 
ft_databrowser(cfg, data);

You can also use ft_databrowser to visualize the continuous data that is stored on disk.

 
cfg         = [];
cfg.dataset = 's04.vhdr'; 
ft_databrowser(cfg);

Extracting the EOG signals

We now continue with re-referencing to extract the bipolar EOG signal from the data. In the BrainAmp acquisition system, all channels are measured relative to a common reference. For the horizonal EOG we will compute the potential difference between channels 57 and 25 (see the plot of the layout and the figure below). For the vertical EOG we will use channel 53 and channel “LEOG” which was placed below the subjects' left eye (not pictured on the layout).

Some acquisition systems, such as Biosemi, allow for direct bipolar recording of EOG. The re-referencing step to obtain the EOG is therefore not required when working with Biosemi data.

% EOGV channel
cfg              = [];
cfg.channel      = {'53' 'LEOG'};
cfg.reref        = 'yes';
cfg.implicitref  = [];
cfg.refchannel   = {'53'};
eogv             = ft_preprocessing(cfg, data);

% only keep one channel, and rename to eogv
cfg              = [];
cfg.channel      = 'LEOG';
eogv             = ft_selectdata(cfg, eogv); 
eogv.label       = {'eogv'};

% EOGH channel
cfg              = [];
cfg.channel      = {'57' '25'};
cfg.reref        = 'yes';
cfg.implicitref  = [];
cfg.refchannel   = {'57'};
eogh             = ft_preprocessing(cfg, data);

% only keep one channel, and rename to eogh
cfg              = [];
cfg.channel      = '25';
eogh             = ft_selectdata(cfg, eogh); 
eogh.label       = {'eogh'};

We now discard extra channels from the data and add the bipolar-referenced EOGv and EOGh channels that we have just created:

% only keep all non-EOG channels
cfg         = [];
cfg.channel = setdiff(1:60, [53, 57, 25]);              % you can use either strings or numbers as selection
data        = ft_selectdata(cfg, data); 

% append the EOGH and EOGV channel to the 60 selected EEG channels 
cfg = [];
data = ft_appenddata(cfg, data, eogv, eogh);

You can check the channel labels that are now present in the data and use ft_databrowser to look at all data combined.

disp(data.label')
  Columns 1 through 12

    '1'    '2'    '3'    '4'    '5'    '6'    '7'    '8'    '9'    '10'    '11'    '12'

  Columns 13 through 23

    '13'    '14'    '15'    '16'    '17'    '18'    '19'    '20'    '21'    '22'    '23'

  Columns 24 through 34

    '24'    '26'    '27'    '28'    '29'    '30'    '31'    'RM'    '33'    '34'    '35'

  Columns 35 through 45

    '36'    '37'    '38'    '39'    '40'    '41'    '42'    '43'    '44'    '45'    '46'

  Columns 46 through 56

    '47'    '48'    '49'    '50'    '51'    '52'    '54'    '55'    '56'    '58'    '59'

  Columns 57 through 59

    '60'    'eogv'    'eogh'

Channel layout

For topoplotting and sometimes for analysis it is necessary to know how the electrodes were positioned on the scalp. In contrast to the sensor arrangement from a given MEG manufacturer, the topographical arrangement of the channels in EEG is not fixed. Different acquisition systems are designed for different electrode montages, and the number and position of electrodes can be adjusted depending on the experimental goal. In the current experiment, so-called 64-electrodes equidistant montage (ActiCap, BrainVision) was used:

The channel positions are not stored in the EEG dataset. You have to use a layout file; this is a *.mat file that contains the 2-D positions of the channels. FieldTrip provides a number of default layouts for BrainVision EEG caps in the fieldtrip/template/layout directory. It is also possible to create custom layouts (see ft_prepare_layout and the layout tutorial). In this example we will use an existing layout file that is included with the example data.

cfg = [];
cfg.channel  = data.label;
cfg.layout   = 'mpi_customized_acticap64';
cfg.feedback = 'yes';
lay = ft_prepare_layout(cfg);
disp(lay)

lay = 
        pos: [59x2 double]
      label: {59x1 cell}
      width: [59x1 double]
     height: [59x1 double]
       mask: {[72x2 double]}
    outline: {[123x2 double]}

Note that the layout should contain correct channel labels that match the channel labels in the data.

Artifacts

An next important step of EEG preprocessing is detection (and rejection) of artifacts. Different approaches of dealing with artifacts are presented in details in the introductory tutorial on artifacts, the visual artifact removal tutorial and the automatic artifact rejection removal tutorial. In this example script, we will use ft_rejectvisual function to visually inspect the data and reject the trials or channels that contain artifacts. We first will try the “channel” mode. In this mode all channels are displayed at once allowing paging through the trials. Then we will try the “summary” mode.

Channel mode

cfg        = [];
cfg.method = 'channel';
ft_rejectvisual(cfg, data)

You can scroll to the vertical EOG channel ('veog', number 61) and confirm to yourself that trials 22, 42, 126, 136 and 150 contain blinks. You can exclude a trial from the data by clicking on it. Note, however, that in this example we do not assign any output to the function. Matlab will create the default output “ans” variable. All the changes (rejections) that you make will be applied to the “ans”. The “data” will remain the same, no trials will be removed!

In ft_rejectvisual with cfg.method='channel' you can go to channel '43' (note that the channel name is '43' and its number is also 43). There you will see that in trials 138 to 149 this channel is a bit more noisy, suggesting that the electrode contact on this side of the cap was temporarily bad. Neighboring channels also suggest that at trial 138 something happened, perhaps a movement of the electrode cap. We are not going to deal with this now, but it is something that you might want to keep in mind for optional cleaning of the data with ft_componentanalysis and ft_rejectcomponent

Summary mode

The data can be also displayed in a “summary” mode, in which case the variance (or another metric) in each channel and each trial is computed. Close the “channel” mode figure and try the “summary” mode. Note, that a new variable “data_clean” will be created now.

cfg = [];
cfg.method   = 'summary';
cfg.layout   = lay;       % this allows for plotting
cfg.channels = [1:60];    % do not show EOG channels
data_clean   = ft_rejectvisual(cfg, data);

The left lower box of Figure 4 shows the variance of the signal in each trial. By dragging the mouse over the trials in this box you can remove them from the plot and reject them from the data. You will see the numbers of the rejected trials in the box on the right. You can undo the rejection by typing the trial's number in “Toggle trial” box. You can also plot the signal in a specific trial with “Plot trial” box. Here, we have plotted the trial 90 - the one with the highest variance. On the topoplot you can see a drift in the channel 48. You can zoom in to this channel by dragging the mouse over it.

Rejection of trials based on visual inspection is somewhat arbitrary. Sometimes it is not easy to decide if a trial has to be rejected or not. In this exercise we suggest that you remove 8 trials with the highest variance (trial numbers 22, 42, 89, 90, 92, 126, 136 and 150). As you see, the trials with blinks that we saw in the “Channel” mode are among them. To complete the rejection press “Quit” button. You get the data_clean variable that will be used for subsequent analyses.

After removing data segments that contain artifacts, you might want to do a last visual inspection of the EEG traces.

cfg          = [];
cfg.viewmode = 'vertical';
ft_databrowser(cfg, data_clean);

Computing and plotting the ERP's

We now would like to compute the ERP's for two conditions: positive-negative judgement and human-animal judgement. For each trial, the condition is assigned by the trialfun that we used in the beginning when defined the trials, this information is kept with the data in data.trialinfo.

disp(data.trialinfo)

 Columns 1 through 19
   2 1 1 2 1 1 2 1 1 2 1 1 1 2 1 1 2 2 2   

 Columns 20 through 38
   1 2 2 2 2 2 1 2 1 2 1 2 2 1 2 1 2 1 2   

 ...

 Columns 172 through 184
   2 1 1 2 2 2 1 2 1 1 1 1 2

FieldTrip automatically kept track of the trials with artifacts that were rejected: the data_clean.trialinfo field contains the condition code for the 184 clean trials, whereas the data.trialinfo field contained the information for the original 192 trials.

We now select the trials with conditions 1 and 2 and compute ERP's.

% use ft_timelockanalysis to compute the ERPs 
cfg = [];
cfg.trials = find(data_clean.trialinfo==1);
task1 = ft_timelockanalysis(cfg, data_clean);

cfg = [];
cfg.trials = find(data_clean.trialinfo==2);
task2 = ft_timelockanalysis(cfg, data_clean);

cfg = [];
cfg.layout = lay;
cfg.interactive = 'yes';
ft_multiplotER(cfg, task1, task2)

Note, that we use the layout file for plotting the results. With the cfg.interactive = 'yes' option you can select channels and zoom in.

The following code allows you to look at the ERP difference waves.

 
difference     = task1;                   % copy one of the structures
difference.avg = task1.avg - task2.avg;   % compute the difference ERP

cfg = [];
cfg.layout      = 'mpi_customized_acticap64.mat';
cfg.interactive = 'yes';
ft_multiplotER(cfg, difference)

Appendix: the trialfun used in this example

function [trl, event] = trialfun_affcog(cfg)

%% the first part is common to all trial functions
% read the header (needed for the samping rate) and the events
hdr        = ft_read_header(cfg.headerfile);
event      = ft_read_event(cfg.headerfile);

%% from here on it becomes specific to the experiment and the data format
% for the events of interest, find the sample numbers (these are integers)
% for the events of interest, find the trigger values (these are strings in the case of BrainVision)
EVsample   = [event.sample]';
EVvalue    = {event.value}';

% select the target stimuli
Word = find(strcmp('S141', EVvalue)==1);

% for each word find the condition
for w = 1:length(Word)
  % code for the judgement task: 1 => Affective; 2 => Ontological;
  if strcmp('S131', EVvalue{Word(w)+1}) == 1
    task(w,1) = 1;
  elseif strcmp('S132', EVvalue{Word(w)+1}) == 1
    task(w,1) = 2;
  end
end

PreTrig   = round(0.2 * hdr.Fs);
PostTrig  = round(1 * hdr.Fs);

begsample = EVsample(Word) - PreTrig;
endsample = EVsample(Word) + PostTrig;

offset = -PreTrig*ones(size(endsample));

%% the last part is again common to all trial functions
% return the trl matrix (required) and the event structure (optional)
trl = [begsample endsample offset task];

end % function

Suggested further reading

After having finished this tutorial on EEG data, you can look at the event related averaging tutorial for MEG data or continue with the time-frequency analysis tutorial.

If you have more questions about preprocessing or timelocked-analysis, you can also read the following faq-s:

Or you can also read the example scripts:


This tutorial was last tested with version 20030101 of FieldTrip using Matlab 1.0a on a 16-bit OS/2 platform.

example/preprocessing_eeg.txt · Last modified: 2014/01/09 08:59 by jmhorschig

You are here: startexamplepreprocessing_eeg
This DokuWiki features an Anymorphic Webdesign theme, customised by Eelke Spaak and Stephen Whitmarsh.
www.chimeric.de Valid CSS Driven by DokuWiki do yourself a favour and use a real browser - get firefox!! Recent changes RSS feed Valid XHTML 1.0