Applying independent component analysis (ICA) and time/frequency analysis to collections of single-trial EEG (or MEG) data, or to collections of averaged event-related potential (ERP) or field (ERF) epochs.
This tutorial will go through the process of decomposing
collections of averaged ERP and unaveraged single-trial EEG data
using the ICA EEG toolbox of MATLAB functions for ICA, single-trial visualization
and time/frequency analysis, with some explanatory
comments and considerations.
We assume you know the basics of using MATLAB.
Currently, inter-window communication works better using Internet Explorer (5.0)
than Netscape (4.74).
There is also
a FAQ page
on ICA decomposition of EEG/ERP data,
and a set of pages on
other ICA research
at CNL/Salk Institute.
To navigate, you may use the popup
Tutorial Outline window.
For help in installing the toolbox, click here.
Open the function help window.
Loading and preprocessing data
Before you start to use the toolbox, you must edit one file,
icadefs.m,
which stores global toolbox parameters including
the location of the toolbox functions.
See its
help file for details.
Next, we define some constants used in the function call examples below.
Note our terminology: a vector of all channels at one time point is termed
a "frame" of data, and a single contiguous data time series, a data
"epoch." A "data matrix" may be composed of
one or more concatenated data epochs.
Note: The MATLAB code and plots below are mostly taken from
icademo().
Next, we load an ascii data matrix consisting of
two concatenated averaged event-related potential (ERP) epochs from
different experimental conditions (in this case, responses to auditory noise
stimuli, first those responded to (Hits), and next those not responded to (Lapses)
by the drowsy subject - see Makeig et al, PNAS 1997).
Note: Another supported input format is native machine floats,
multiplexed by channel, using toolbox function
floatread().
Now for convenience filter the data to remove 60-Hz line noise, other high-frequency noise, and 40-Hz steady-state responses from the data as not of interest for this tutorial.
eegfilt())
is one of only two toolbox functions that use the MATLAB signal processing toolbox.
Note: Click on any of the underlined
function calls to view its help message.
The assumed format for toolbox data is two-dimensional matrices of size
(channels,frames*epochs).
The toolbox includes a function,
matsel(),
that facilitates extracting channels,frames, and/or epochs from such matrices.
For example, to plot channels 1-5 of the first data epoch on the same timebase:
Above, the prime ()' is needed to transpose the submatrix returned by
matsel(),
since plot() plots columns instead of rows.
Note: A 3-D matrix M, as defined in MATLAB5,
of size (channels, frames, epochs) can be transformed to the two-dimensional
toolbox format by M = M(:,:);
Another simple preprocessing function is
blockave(),
which averages specified epochs in a data matrix:
Additional preprocessing tools for extracting and averaging single trials
from continuous or multi-epoch data are described in
the section on single-trial data decomposition.
The scrolling plotting tool
eegplot()
is useful for examining large data sets of any sort.
Here we plot our data as two successive 1-second data epochs.
When you use
eegplot(),
try using the on-screen buttons and pull-down menu (not shown below)
to adjust the display.
Note: Click on any of the plots for
a close-up view in a new window.
Function
timtopo()
displays a data epoch by plotting
all the single channel traces on a single axis
and showing interpolated scalp maps of the field at specified time points.
Notice how the field pattern changes across the response epoch.
Function
timtopo()
can be thought of as showing a few "still frames" from a movie of
the changing electric field across the data epoch. (To make and view an actual
field movie, use
eegmovie()
and seemovie()).
Note:
Given a default '0' instead of the four specified plot times,
timtopo()
would plot the map at the
power maximum.
timtopo()
plots each scalp map using
topoplot().
Below, we use
eegplot()
again twice to plot both data epochs in separate figure axes.
With the 'noui' argument,
eegplot()
can be used to embed plots (without gui buttons) in larger figures.
(Toolbox function
sbplot(),
used below,
is more flexible than standard Matlab "subplot()").
Next, we overplot the two 312-point data epochs at each channel using
plotdata().
The epoch end times are 0 ms and 995 ms.
Read four-character channel identifiers from the file 'chan.nam.' Use the default (0) color order and
default (no arg) ydir (positive-up).
Another way to plot the same data -- in two-dimensional topographic order
of the locations of the recording electrodes.
Below, the channels are arrayed from the back of the
head (lowest traces) to the front (highest traces).
plottopo()
reads the channel locations from file "eloc.loc".
To see the format of this file, say
topoplot example
In addition to the two-dimensional scalp map plotting function
topoplot(),
the toolbox contains a 3-D scalp map plotting function,
headplot().
Use of
headplot()
requires a wire frame head model. Though the current toolbox release does not contain the
model we have available (taken from MRI data), as we are not completely satisfied with it,
this model may be downloaded from
our ftp site (98k).
The
headplot()
script requires building an interpolation model for each electrode montage,
a lengthy computation.
Once this is constructed, however, making individual data plots is quick
(see small example below using another head model),
although the resulting plot files may be large. Within MATLAB, the 3-D scalp
map images may be rotated using the mouse.
The post-motor P3 subcomponent Pmp (discussed in Makeig et al., JNS 1999)
Next, decompose the data using ICA ...
Display the tutorial outline
Open the function help window
Download the ICA EEG toolbox
View log of changes to this tutorial.
chans = 14; % data channels (rows in data matrix)
frames = 312; % frames per data epoch (columns in data submatrix)
epochs = 2; % epochs in simulated data
srate = 312.5; % data sampling rate in Hz
limits = [0 995 0 0]; % epoch time limits 0-995 msec, amp. limits unspecified
lowpass = 30; % filter band lowpass freq in Hz
load pnas.adt -ascii
% Load an ascii ERP data matrix and reshape as two
% successive (14,312) 14-channel responses.
% This data is included in the toolbox.
data = [pnas(15:28,:) pnas(1:14,:)];
% Shape the data matrix. Make each column a data frame.
data = data/24.2;
% Scale to uV
% data = eegfilt(data,srate,lowpassHz,highpassHz,frames);
[data,filtwts] = eegfilt(data,srate,0,lowpass,frames);
% Lowpass filter below 30 Hz; filter each epoch similarly
% [dataout] = matsel(data,frames,framelist,chanlist,epochlist);
figure; plot(matsel(data,312,0,1:5,1)');
% dataout = blockave(data,frames,epochs);
figure; plot(blockave(data,312)'); % Plot the mean of the two data epochs
Visualizing the data
% eegplot(data,srate,spacing,eloc_file,windowlength,title)
eegplot(data,srate,0,'eloc.loc',1,'Two data epochs');
% eegplot() is the only toolbox function that creates its own figure.
% timtopo(data,'chan_locs',[limits],[plottimes]','title',[plotchans],[voffsets]);
figure; timtopo(matsel(data,frames,0,0,1),'eloc.loc',[0 995 -12 12],...
[250 320 390 500],'Target Hits'); % Image scalp maps at 250, .. 500 ms
% eegplot('noui',data,srate,spacing,eloc_file,startsec,color)
figure
sbplot(1,2,1);
eegplot('noui',data(:,1:312),srate,0,'eloc.loc',0,'r'); title('Lapses')
sbplot(1,2,2);
eegplot('noui',data(:,313:624),srate,0,'eloc.loc',0,'b'); title('Hits')
% plotdata(data,frames,limits,title,channames,colors,rtitle,ydir)
figure; plotdata(data,frames,[0 995 -10 10],'ERP Data',...
'chan.nam',0,'(2 conditions)');
% plottopo(data,'chan_locs',frames,limits,title,...
% channels,axsize,colors,ydir)
figure; plottopo(data,'eloc.loc',frames,[0 995 -10 10],'ERP Data');
Comments and suggestions on this tutorial
are welcome. Email scott@salk.edu