Introduction: Classifying Raw EEG

by Chuck Anderson
Colorado State University
November, 2005

Create the Data

We will first build EEG sample matrices and correct class labels to be used in the following examples. We select data from Subject 1, all five tasks, and trials 1 through 5, all from data available at http://www.cs.colostate.edu/eeg.


load_eegdata;
subject = 1
tasks = [1 2 3 4 5]
%tasks = [2 3]
trials = 1:5
trials = squeeze(eegdata(subject,tasks,trials));
% Remove EOG channel
[nr,nc] = size(trials);
for r = 1:nr
  for c = 1:nc
    trials{r,c} = trials{r,c}(1:6,:);
  end
end

classes = repmat((1:length(tasks))',1,length(trials));

Create Function for Representing Data as Features

Now let's say we want to try to classify the raw EEG data, possibly lagged. Write a function that transforms samples into features. Same function is used to calculate parameters needed for transformation by passing in third argument of 'make'. For this case, 'make' does nothing. For features that are just the raw EEG data with lags, this function will suffice.


function featuresOrCV = makeRawEEGFeatures(trials,CV,command)
%%%  featuresOrCV = makeRawEEGFeatures(trials,CV,command)
%%%   trials: ntasks x ntrials cell array of data
%%%   CV.parms contains
%%%           .nlags
%%%   featuresOrCV: ntasks x ntrials cell array of featurized data
%%%         or  set to CV if command='make'

if nargin > 2 && strcmp(command,'make')
  featuresOrCV = CV;
  return;
end

nlags = CV.parms.nlags;
featuresOrCV = lagize(trials,nlags);

This function requires


function lagged = lagize(data,nlags,delta)
% lagged = lagize(datacellofwindows,nlags,delta)
if nargin < 3
  delta = 1;
end

nlags = nlags + 1;

[r,c] = size(data);
lagged = cell(r,c);
for ri = 1:r
  for ci = 1:c
    datacols = data{ri,ci};
    ncols = size(datacols,2);
    starts = [0:nlags-1]*delta+1;
    laggedWindow = [];
    for i = starts
      laggedWindow = [laggedWindow; datacols(:,i:(ncols-starts(end)+i))];
    end
    lagged{ri,ci} = laggedWindow;
  end
end

Test it by


>> CV.parms.nlags = 1;
>> features = makeRawEEGFeatures(trials,CV);
>> size(features)
ans =
     5     5
>> size(features{1,1})
ans =
          12        2499

We have a 5 x 5 cell array of matrices of features, where each matrix is 12 x 2499, because we have voltages from 6 channels augmented by the same voltages delayed by one sample. Starting with 2500 samples (10 seconds of EEG sampled 250 times per second), we get 2499 total samples.

Create Function for Selecting Feature Samples

The cross-validation code handles the random selection of a fraction of samples from each trial. However, it does need a function that combines the requested samples into a matrix of input data, X, and target data, Y. For the raw representation, the following function works fine.


function [X,Y] = selectSamples(oneTrial,indices,class,CV)

X = oneTrial(:,indices);
Y = repmat(class,1,length(indices));

Create Function for Classifying Features

Now we need a classifier. Like the representation function, if the last argument is 'make', the classifier's parameters are chosen by performing cross-validation to evaluate every set of specified parameter values.

First, let's discuss the parameters. For the raw EEG representation, the only paramter is the number of lags. The parameter sets are specified as a vector of structs, like this.


lagsraw = 0:16;
clear parameterSets;
k = 0;
for nl = lagsraw
  k = k + 1;
  parameterSets(k).nlags = nl;
end

Cross-validation is performed by choosing one trial to be the test trial. The remaining trials form the training data. We must also specify the fraction of samples from each trial that are to be randomly selected from the training trials for the validation set, and the number of times we repeat this random selection. Let's call the training data remaining after removing the validation data the tuning data, used to tune the classifier. For each repetition, a classifier is formed on the tuning data, applied to the validation data, and the resulting RMS error is averaged over these representations to evaluate each parameter set. The parameter set that is determined to give the lowest average validation data error is saved as the best parameters in CV.parms.

Then, a classifier is constructed on all the training data using these best parameters and applied to the test trial. This error is saved. Then the whole procedure is repeated after each trial is selected as the test trial. With five trials, this procedure is repeated five times.

Anyway, here is an example classifier function that implements linear discriminant analysis (LDA). It should work for any representation you come up with, so for now you can choose to not read through this code and skip to the next section.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Make a Linear Discriminant Classifier 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fractionCorrect,classifierOrClasses] = LDA(data,CV,selectInOneTrialF,command)
%%% [fractionCorrect,classifierOrClasses] = LDA(train,CV)
%%%   train.features: ntasks x ntrials cell array of featurized data
%%%        .classes: ntasks x ntrials matrix of class labels
%%%   CV contains for collectTrainValidation()
%%%     .parms.firstMode
%%%     .parms.nModes
%%%     .validateIndices: ntasks x ntrials cell array of indices
%%%     .classes:  ntasks x ntrials cell array of class labels (ints)
%%%   command:  Optional. If present, must be 'make' and classifier is created.
%%% Returns
%%%   classifier.weights  
%%%             .bias
%%%   validateCorrect: fraction of validation samples correct
%%%   X: training data

if nargin > 3 && strcmp(command,'make')

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% make the classifier
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if ~isfield(CV.parms,'firstMode')
    CV.parms.firstMode = 1;  %so this works if data.features are not cell array of
    CV.parms.nModes = 1;    %windows, but just nchannels x nsamples matrix
  end

  [X,Y,Xval,Yval] = collectPartitions(data,CV,selectInOneTrialF,'train');

  [p,N] = size(X);
  K = size(data.features,1);

  %% Calculate class priors, means and sum of covariance matrices.
  priors = zeros(1,K);
  means = zeros(K,p);
  covsum = zeros(p,p);
  classNames = unique(Y);
  for i = 1:length(classNames)
    mask = Y == classNames(i);
    Nthisclass = sum(mask);
    priors(i) = Nthisclass / N;
    Xmat = X(:,mask);
    means(i,:) = mean(Xmat');
    XmeanZero = Xmat - repmat(means(i,:)',1,Nthisclass);
    covsum = covsum + XmeanZero * XmeanZero';
  end
  covariance = covsum / (N-K);
  invcov = inv(covariance);

  %% Calculate weights and biases for each class
  weights = zeros(p,K);
  bias = zeros(K,1);
  for k = 1:K
    weights(:,k) = invcov * means(k,:)';
    bias(k) = - 0.5 * means(k,:) * weights(:,k) + log(priors(k));
  end

  classifier.weights = weights;
  classifier.bias = bias;

  if ~isempty(CV.validateIndices)
    disc = classifier.weights' * Xval + repmat(classifier.bias,1,size(Xval,2));
    [junk,whichmax] = max(disc);
    classes = classNames(whichmax);
    validateCorrect = sum(classes == Yval) / length(Yval);
  else
    validateCorrect = [];
  end

  classifierOrClasses = classifier;
  fractionCorrect = validateCorrect;

elseif nargin > 3 && ~strcmp(command,'make')

  error('Command to LDA must be ''make'' or absent.');

else
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% use the classifier
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  [X,Y] = collectPartitions(data,CV,selectInOneTrialF);

  disc = CV.classifier.weights' * X + repmat(CV.classifier.bias,1,size(X,2));
  [junk,whichmax] = max(disc);
  classNames = unique(Y);
  classes = classNames(whichmax);
  testCorrect = sum(classes == Y) / length(Y);

  classifierOrClasses = classes;
  fractionCorrect = testCorrect;

end

Put It All Together by Calling `runCV`

Now the cross-validation procedure is performed by calling the main function runCV. Here is a complete example, including some of the code shown above and using the functions defined above.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Making data
clear trials classes;

load_eegdata;

subject = 1
tasks = [1 2 3 4 5]
%tasks = [2 3]
trials = 1:5

trials = squeeze(eegdata(subject,tasks,trials));
% Remove EOG channel
[nr,nc] = size(trials);
for r = 1:nr
  for c = 1:nc
    trials{r,c} = trials{r,c}(1:6,:);
  end
end

classes = repmat((1:length(tasks))',1,length(trials));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Specify common paramters
validationRepetitions = 5;
validationFraction = 0.2;
lags = 0:16;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Raw lagged representation

clear parameterSets;
k = 0;
for nl = lags
  k = k + 1;
  parameterSets(k).nlags = nl;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Run all experiments
CV.raw = runCV(trials,classes,parameterSets,...
    validationFraction,validationRepetitions,...
    @makeRawEEGFeatures,@LDA,@selectSamples);

save allCV CV; % Save result to file named 'allcv.mat'

Classifying MSF Bases

Representation

This representation take a little more work, because the samples must be collected into windows, the MSF decomposition must be calculated, then a subset of the basis vectors must be selected and concatenated for the feature representation of a window. The following function performs all these steps.


function features = makeSVDModes(trials,CV,command)
%%%  features = makeSVDModes(trials,CV)
%%%   trials: ntasks x ntrials cell array of data
%%%   CV.parms contains
%%%           .nlags
%%%           .winSize
%%%           .winShift
%%%   features: ntasks x ntrials cell array of featurized data

if nargin > 2 && strcmp(command,'make')
  features = CV;
  return;
end

nlags = CV.parms.nlags;
winSize = CV.parms.winSize;
winShift = CV.parms.winShift;


features = msfize(windowize(lagize(trials,nlags),winSize,winShift));

It is based on lagize given above plus the following function in windowize.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function allwindows = windowize(data,wsize,winshift)

if nargin == 2
  winshift = 1;
end

[r,c] = size(data);
allwindows = cell(r,c);
for ri = 1:r
  for ci = 1:c
    oneTrial = data{ri,ci};
    nwin = 0;
    for wi = 1:winshift:size(oneTrial,2)-wsize+1
      nwin = nwin + 1;
    end
    for w = 1:nwin
      start = (w-1)*winshift+1;
      windows{w} = oneTrial(:,start:start+wsize-1);
    end
    allwindows{ri,ci} = windows;
  end
end

and the following function in msfize


function msfwindows = msfize(data)

[r,c] = size(data);
msfwindows = cell(r,c);
for ri = 1:r
  for ci = 1:c
    windows = data{ri,ci};
    nwin = length(windows);
    msfwindowsTrial = cell(1,nwin);
    for w = 1:nwin
      windows{w} = windows{w} - repmat(mean(windows{w},2),1,size(windows{w},2));
      [u,s] = mnf2(windows{w});
      u = s';
      for col = 1:size(u,2)
	if u(1,col) < 0
	  u(:,col) = u(:,col) * -1;
	end
      end
      msfwindowsTrial{w} = u;
    end
    msfwindows{ri,ci} = msfwindowsTrial;
  end
end

which depends on mnf2.m:


function [phi,psi] = mnf2(X)
% function [psi,phi] = mnf(X)
%    X - an p x N data matrix  (N is number of samples)
%
% Out:
%    psi - the mixing coefficients (pxp)
%   phi - basis  p x N
%
% Example:
%    s = sin(0:99);
%    n = [zeros(1,50) rand(1,50)];
%    M = [1 0.5; 0.2 0.8];
%    X = M * [s;n];
%    subplot(2,1,1);
%    plot(X');
%    subplot(2,1,2);
%    [psi,phi] = mnf2(X);
%    plot(phi');


[p,N] = size(X);

dX = X(:,1:end-1) - X(:,2:end);

[psi,D] = eig(X*X', dX*dX');

phi = psi' * X;
psi = psi'; % to agree with output of mnf.m

Selecting Samples

We must also write a new selection function to deal with windows of basis vectors.


function [X,Y] = selectModes(oneTrial,indices,class,CV)

if ~iscell(oneTrial)
  error('selectModes must be applied to cella array, but it is not.');
else
  X = [];
  Y = [];
  
  modespan = CV.parms.firstMode : CV.parms.firstMode + ...
      CV.parms.nModes - 1;
  for w = indices
    modes = oneTrial{w}(:,modespan);
    modes = reshape(modes,size(modes,1)*size(modes,2),1);
    X = [X modes];
    Y = [Y class];
  end
end

Now Do It

Now, put it all together as before, still using LDA to classify, and we can run this experiment by doing the following.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% MSF basis vector representation

clear parameterSets;
k = 0;
for winSize = winSizes
  for nl = lags
    for nModes = 1:nchannels*(nl+1)
      for firstMode = 1:(nchannels*(nl+1))-nModes+1
	k = k + 1;
	parameterSets(k).nlags = nl;
	parameterSets(k).winSize = winSize;
	parameterSets(k).winShift = winShift;
	parameterSets(k).firstMode = firstMode;
	parameterSets(k).nModes = nModes;
      end
    end
  end
end

%%% Run all experiments
CV.msf = runCV(trials,classes,parameterSets,...
    validationFraction,validationRepetitions,...
    @makeMSFModes,@LDA,@selectModes);

save allCV CV