Minimalist BCI

This section is associated with the book chapter

Delorme, A., Kothe, C., Bigdely, N., Vankov, A., Oostenveld, R., Makeig, S. Matlab Tools for BCI Research? In "human-computer interaction and brain-computer interfaces". Editors : Tan, D. and Nijholt, A. To appear in 2010. Springer Publishing. Download PDF here.

The minimalistic BCI script was programmed by C. Kothe and A. Delorme. The top function (train_bci) performs temporal filtering and training of a CSP filter. The second function (test_bci) applies the model to incoming blocks of raw data, and can be used for online processing. The scripts load data (here from the BCI Competition III), perform training and testing, and display results. Assuming a function “get_rawdata” allows asynchronous data collection, the fourth script performs real-time classification and displays an evolving time course of classification over the past 2 seconds with a refresh rate of 100 ms.

Training BCI function

The first function extract data epochs from continuous data and trains the BCI using class-grouping, CSP

function [S,T,w,b] = train_bci(EEG,Fs,mrk,wnd,f,nof,n)
% [S,T,w,b] = train_bci(Raw-Signal, Sample-Rate, Markers, 
%     Epoch-Wnd, Spectral-Flt, Flt-Number, Flt-Length)
% In:  
%       Raw-Signal   :   raw EEG/MEG calibration session data [Samples x Channels]
%       Sample-Rate  :   sampling frequency of the signal in Hz, e.g. 100
%       Markers      :   marker channel: entries in {1,2} indicate different mental conditions,
%                        others are ignored, e.g. sparse(1,[1132 6462 8234 2342],[1 2 2 1])
%       Epoch-Wnd    :   window, in seconds, relative to each condition marker to 
%                        indicate the duration of the conditions, e.g. [0.5 3.5]
%       Spectral-Flt :   spectral filter to be used; function of frequency in Hz, e.g. @(f) f>7 & f<30
%       Flt-Number   :   number of spatial filters to compute per condition via CSP, e.g. 3
%       Flt-Length   :   length, in samples, of the spectral filter, e.g. 200
% Out:
%       S,T,w,b     :   parameters for test_bci
% Example:
%       a = <load via EEGLAB, with markers for condition A as '1' and markers for condition B as '2'>
%       [S,T,w,b] = train_bci(a.data, a.srate, sparse(1,[a.event.latency], strcmp({a.event.type},'1') + 2*strcmp({a.event.type},'2'), ...
%                             [0.5 3.5], @(f) f>7 & f<30, 3, 200);
 
% frequency filtering and temporal filter estimation
[t,c] = size(EEG); idx = reshape(1:t*c-mod(t*c,n),n,[]);
FLT = real(ifft(fft(EEG).*repmat(f(Fs*(0:t-1)/t)',1,c)));
T = FLT(idx)/EEG(idx);
 
% data epoching, class-grouping and CSP
wnd = round(Fs*wnd(1)):round(Fs*wnd(2));
for k = 1:2
    EPO{k} = FLT(repmat(find(mrk==k),length(wnd),1) + repmat(wnd',1,nnz(mrk==k)),:);
end
[V,D] = eig(cov(EPO{2}),cov(EPO{1})+cov(EPO{2}));
S = V(:,[1:nof end-nof+1:end]);
 
% log-variance feature extraction and LDA
for k = 1:2
    X{k} = squeeze(log(var(reshape(EPO{k}*S, length(wnd),[],2*nof))));
end
w = ((mean(X{2})-mean(X{1}))/(cov(X{1})+cov(X{2})))';
b = (mean(X{1})+mean(X{2}))*w/2;

Testing BCI function

The second function test the BCI. It applies the temporal and spatial filter given by the training function above

function y = test_bci(X,S,T,w,b)
% Prediction = test_bci(Raw-Block, Spatial-Flt, Temporal-Flt, Weights, Bias)
% A linear online detector for oscillatory processes.
%
% In:
%   X : incoming raw sample data [Samples x Channels]
%   S,T: spatio-temporal filter [Channels x Filters], [Samples x Samples]
%   w,b: linear classifier [Filters x 1], [1 x 1]
%
% Out:
%   y : the processed output
%	
% Notes:	
%	y can be post-processed by sign(y) for pure classification, or 1./(1+exp(-y)) for logistic regression
%
 
global B; % B is the buffer
if any(size(B) ~= [length(T),length(S)]) 
    B = zeros(length(T),length(S)); 
end
B = [B;X]; B = B(end-length(T)+1:end,:);
y = log(var(T*(B*S)))*w - b;

Script running the training and testing BCI function

The script below reads training data, train the BCI. It then reads some test data to test the BCI. The data is available here (50 Mb).

load data_set_IVb_al_train
flt = @(f)(f>7&f<30).*(1-cos((f-(7+30)/2)/(7-30)*pi*4));
[S,T,w,b] = train_bci(single(cnt), nfo.fs, ...
    sparse(1,mrk.pos,(mrk.y+3)/2),[0.5 3.5],flt,3,200);
 
load data_set_IVb_al_test
for x=1:length(cnt)
    y(x) = test_bci(single(cnt(x,:)),S,T,w,b); 
end
 
load true_labels
plot((1:length(cnt))/nfo.fs,[y/sqrt(mean(y.*y)); true_y']);
xlabel('time (seconds)'); ylabel('class');

The figure plotted by the script above is shown below

Script for online BCI testing and result plotting

Finally, assuming offline training training, the small script below performs online classification. It assumes the existence of a function "get_rawdata" that delivers streaming data in real time.

y = []; 
start(timer('R = get_rawdata; y(end+1) = test_bci(R,S,T,w,b); plot(y(max(1,length(y)-20):end);', 'InstantPeriod', 0.1));

Obtaining the data

The data can be obtained from the BCI competition III webpage, dataset IVb.