Automatic artifact rejection

Before further analysis in any of the other tutorials, it is best to have artifact free data. Within FieldTrip you can choose to do visual/manual or automatic artifact detection and rejection.

MEG data from the VSM/CTF system can be recorded either continuously or as trial based. The dataset used in the other tutorials is collected by trial based acquisition. The default settings in the automatic artifact detection routine are not optimal for trial based recordings. Therefore, to demonstrate automatic artifact detection and rejection we will use a continuously recorded dataset. How to do the trial selection is described here.

For a successful analysis of EEG or MEG signals, “clean” data is required. That means that you should try to reduce the amount of variance in the data due to factors that you cannot influence. One of the factors that is difficult to control are the presence of artifacts in the data. These artifact are physiological or can result from the acquisition electronics. The strongest physiological artifacts stem from eye blinks, eye movements and head movements. Muscle artifact from swallowing and neck contraction can be a problem as well. Artifacts related to the electronics are 'SQUID jumps' or spikes seen in several channels.

To start with, it is best to avoid those artifacts during the recording. You can instruct the subject not to blink during the trial, but instead give him some well-defined time between the trials in which he is allowed to blink. But of course there will always be some artifacts in the raw data.

The following steps are taken to detect artifacts and reject them from the data:

• Define segments of data of interest (the trial definition) using ft_definetrial
• Detect EOG, muscle and jump artifacts using the respective artifact detection functions: ft_artifact_eog, ft_artifact_muscle and ft_artifact_jump
• Reject the artifacts from the trial definition using ft_rejectartifact
• Read the data into Matlab using ft_preprocessing

The Matlab code will have the following outline:

% define the trials
cfg = [];
cfg.trialdef = ...
cfg = ft_definetrial(cfg)

% detect the artifacts
cfg = ft_artifact_eog(cfg)
cfg = ft_artifact_muscle(cfg)
cfg = ft_artifact_jump(cfg)

% reject the artifacts
cfg = ft_rejectartifact(cfg)

% read and preprocess the data
data = ft_preprocessing(cfg)

Note that for each subsequent step you continue with the output cfg of the previous step. Both the trial and artifact definitions are kept in the cfg.

The data used in this example is the MEG dataset ArtifactMEG.ds (available from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/ArtifactMEG.zip). This dataset was acquired continuously with trials of 10 seconds, and during the recording the experimenter instructed the subject to make artifacts and wrote down the trial number in which the artifacts were made. We will use these manual trial numbers later on.

First we will extract all trials from the dataset:

cfg                    = [];
cfg.dataset            = 'ArtifactMEG.ds';
cfg.headerformat       = 'ctf_ds';
ctf.dataformat         = 'ctf_ds'; 
cfg.trialdef.eventtype = 'trial';
cfg                    = ft_definetrial(cfg);

% remember the complete trial definition, we need it again later
trl = cfg.trl;
trl([1 end],:) = []; % remove first and last trial - otherwise padding is not possible in this dataset

The trial definition “trl” is an Nx3 matrix, N is the number of trials. The first column contains the sample-indices of the begin of each trial relative to the beginning of the raw data, the second column contains the sample-indices of the end of each trial, and the third column contains the offset of the trigger with respect to the trial. An offset of 0 means that the first sample of the trial corresponds to the trigger. A positive offset indicates that the first sample is later than the trigger, a negative offset indicates that the trial begins before the trigger.

In order to detect artifacts in the data, the artifact detection functions ft_artifact_eog, ft_artifact_muscle and ft_artifact_jump will be used. These functions process the data in a way optimized to detect the respective artifacts.

For EOG, jump and muscle artifact detection the same z-score-based algorithm is used, which is based on the idea that artifacts are anomalies with more prominent amplitudes than general brain activity. The processing steps are:

• Channel selection: limits the data to the channels closest to the artifact sources.
• Bandpass filtering: limits the data to the frequencies in which the artifacts are most dominant.
• Hilbert analytic amplitude: yields the envelope of the signal (for each selected channel).
• Normalization by calculating the z-scores for each channel
• Summation for obtaining one summed z-value for each moment in time by adding the z-scores of all selected channels, and normalizing the sum by dividing it by the root of the number of channels summed.
• Threshold comparison is done to determine if an artifact is detected.

The formulas for calculating the z-scores are:

z_ch,t = (x_ch,t - mu_ch)/(sigma_ch)

where: mu_ch = (1/N) sum_t=1^N (x_ch,t)

with: N = the total number of time samples and

sigma_ch = sqrt( (1/N) sum_t=1^N(x_ch,t - \mu_ch)² )

In the code this formula is formed such as to optimize the calculation of the channel means and standard deviations.

The summation is performed like:

zsum_t = (sum_ch^C(z_ch,t) / (sqrt(C))

with: C = the number of channels.

Furthermore, data padding is used. Here we distinguish between three different types of padding:

• Trial padding

You can include segments of the data preceding or following the trial (e.g. if you want to avoid trials with an eye blink right before the trial onset). For this purpose trial padding (indicated in cfg.artfctdef.xxx.trlpadding) is used. Trial padding pads data on both sides of the trial with the specified length, such that artifact detection and rejection is also performed for those padded segments.

• Filter padding

Each artifact type can best be detected with filters of a certain pass band (e.g. pass band of 1–15 Hz for eye blinks, and 110–140 Hz for muscle artifacts). However, filters are known to produce edge effects which can also be detected by the artifact-detection routine and mistaken for real artifacts. To avoid that, in addition to trial padding, filter padding (specified in cfg.artfctdef.xxx.fltpadding) is used. Filter padding of a specified amount is appended to the already existing trial padding (if cfg.artfctdef.xxx.trlpadding had been used) on both sides. After filtering has been performed, the filter padding is removed.

Figure 1; Trial and filter padding. Filter padding is only used during filtering and removed afterwards. Artifacts contained in the trial padding will also be detected.

• Artifact padding

When the artifact has been detected, artifact padding (cfg.artfctdef.xxx.artpadding) is used to pad the artifact on both sides (Fig. 2). This has two purposes. The first purpose is to link the artifact to the trial with which it should be associated to ensure that the trial (or its part) will be subsequently rejected. The second purpose is to define the extent of a long artifact (e.g. of a blink).

Figure 2; Artifact padding extends the segment of data that will be marked as artifact.

EOG artifacts can be detected by the function ft_artifact_eog. By indicating feedback='yes', you enter an interactive mode where you can browse through the data and adjust the cut-off value accordingly (i.e., the z-value used for thresholding).

cfg                        = [];
cfg.trl                    = trl;                  % use the previously determined trial definition
cfg.dataset                = 'ArtifactMEG.ds';
cfg.continuous             = 'yes';
cfg.artfctdef.eog.feedback = 'yes';
cfg.artfctdef.eog.channel  = 'EOG';                % only use the EOG channel
                                                   % note nomenclature of channel (not necessarily 'EOG')
cfg = ft_artifact_eog(cfg);

Figure 3; Summary of artifact detection by ft_artifact_eog, showing the cumulative z-value for each time sample. The current cut-off value is indicated by the dashed red line, all samples exceeding this will be marked as artifacts.

In the command window you see the following question:

would you like to page through the data [y/n]? 

If you type “y” you can browse through all the trials to see whether the default threshold is suitable to detect all the eye artifacts in this dataset.

Figure 4; You can browse through all the trials using the feedback='yes' mode of ft_artifact_eog. The current cut-off value is indicated by the dashed red line, all samples exceeding this will be marked as artifacts. Notice the artifact padding! The upper panel shows the raw data of the indicated channel and trial, the lower panel shows the z-scored data.

If you press the stop button the figure disappears and you see the following question in the command window:

give new cutoff value, or press enter to accept current value [4]: 

You can change the threshold and browse through the data again to see whether this threshold is better to detect all the EOG artifacts. You can repeat this procedure of adjusting the threshold until you are satisfied that most EOG artifact are detected.

If you type enter, all the pieces of data that were colored red in the figure are marked as artifacts and the output of ft_artifact_eog is generated. This output is a cfg which includes a field cfg.artfctdef.eog.artifact that contains the begin and end samples of the detected artifacts.

The same way as ft_artifact_eog is used to detect eye artifacts, ft_artifact_muscle and ft_artifact_jump can be used to detect muscle and jump artifacts. Note that here all MEG channels are scanned for artifacts (instead of only the EOG channel in the eye artifacts case), which will take more time. Make sure that for each step you continue with the output cfg of the previous step, because this contains both the trial and artifact definitions.

% muscle
cfg.artfctdef.muscle.feedback = 'yes';
cfg.artfctdef.muscle.channel  = 'MEG';
cfg = ft_artifact_muscle(cfg);

For detection of jump artifacts it is important to specify the padding (cfg.padding) that will be used for subsequent preprocessing! The data that will be used for padding during preprocessing - especially relevant when you plan to filter your data - is than also checked for artifacts. This is important, as you do not want to pad your artifact free trials with trials containing a jump (filtering this would be detrimental).

% jump
cfg.padding                   = 14; % length to which the trials are padded for filtering
cfg.artfctdef.jump.feedback   = 'yes';
cfg.artfctdef.jump.channel    = 'MEG';
cfg = ft_artifact_jump(cfg);

Figure 5; Using ft_artifact_jump, a SQUID jump can be observed on channel MLT42, trial 17.

Reject those data segments which are contaminated by an artifact. One can either reject the complete trial or reject only the part with the artifact (this leads to partial trials with variable trial lengths). To remove the detected artifacts from the trial definition (trl) use the following code:

cfg.artfctdef.reject = 'complete';   % this rejects complete trials, use 'partial' if you want to do partial artifact rejection
cfg = ft_rejectartifact(cfg);

The output contains cfg.trl which is the cleaned trial definition, and cfg.trlold which contains the old trl (the output of ft_definetrial).

Now we are ready to read the data into Matlab using ft_preprocessing. Only the artifact-free trials will be read.

cfg.channel = 'MEG';
data = ft_preprocessing(cfg);

This tutorial was last tested with version 20090330 of FieldTrip.