Programming Tutorial:Implementing a simple Matlab-based Filter

From BCI2000 Wiki
Jump to navigation Jump to search

Online Algorithm Verification

In the field of BCI signal processing research, novel methods are often proposed and tested on the basis of existing data. While this is adequate as a research strategy, it is important to keep in mind that timely feedback of brain signal classification is an essential element of a BCI. From this, it is clear that, ideally, a proposed BCI signal processing or classification method should be verified with respect to its viability and usefulness in a true online setting, providing feedback to the participant causally and in real-time.

BCI2000 facilitates the transformation of an existing offline data analysis method into a functional online system.

  • For signal processing components, it provides a convenient, straightforward programming interface.
  • A signal processing component may be implemented as a set of Matlab scripts.

Still, it is easy to underestimate the effort required to transform an existing offline implementation of a signal processing algorithm into a functional online implementation. While BCI2000 tries to make the transformation as simple as possible, it cannot remove the effort required to deal with chunks of data, which implies the need of

  • buffering -- rather than having immediate access to a continuous data set, it may be necessary to maintain an additional data buffer;
  • with the Matlab interface, maintaining a consistent state between subsequent calls to the processing script.

An Example Algorithm

In this scenario, we use a simple, straightforward BCI signal processing algorithm, to show the steps necessary to modify the algorithm such that it may be used to build an online system in BCI2000.

In the example, signal processing consists of IIR bandpass filtering, followed with RMS envelope computation, and linear classification. A typical Matlab implementation of that algorithm might consist of about ten lines:

function class_result = classify( data, band, classifier );

% Use as
%   class_result = classify( data, band, classifier )
%
% This function takes raw data as a [channels x samples] vector in the 'data'
% input variable. 
%
% Then, it computes bandpower features for the band specified in the 'band'
% input variable, which is a number between 0 and 0.5, specifying
% center frequency in terms of the sampling rate.
%
% As a last step, it applies the 'classifier' matrix to the features in order
% to obtain a single classification result for each sample. The 'classifier' 
% vector specifies a classification weight for each processed channel.
% 
% The result is a single classification result for each sample.
%
% This requires the Matlab signal processing toolbox.

% Design bandpass filters and apply them to the data matrix.
% The filtered data will contain bandpass filtered data as channels.
[n,Wn]=buttord(band*[0.9 1.1]/2,band*[0.7 1.4]/2,1,60);
[b,a]=butter(n,Wn);
processed_data = filter(b,a,data);

% For demodulation, rectify and apply a low pass.
[b,a]=butter(1,band/4);
processed_data = filter(b,a,abs(processed_data));

% Finally, apply the linear classifier.
class_result = processed_data * classifier;

Note that, to be viable in an online environment, an algorithm must operate on its signal in a causal way, i.e. it may not use future samples in order to process present samples. (Still, a certain amount of non-causality may be possible by windowed operation on buffered data, although this will increase the effective delay between input and output data.)

Also note that the classify function omits spatial filtering -- in Matlab, this may be done easily by pre-multiplying the data with a spatial filter matrix if desired.

Translating the existing MATLAB code into BCI2000 functions

For use with BCI2000, the signal processing code needs to be cast into a form that suits the BCI2000 filter interface. In this event-based model, portions of code are called at certain times to configure a filter component's internal state, and to act upon the signal in chunks on its way through the BCI2000 chain of filters.

Process

The most central event in the filter interface is the Process event. The Process event handler receives a signal input, will process this input according to the filter component's role in the signal processing chain, and return the result of processing in a signal output variable. It is important to understand that the Process handler is called separately for each chunk of data, and thus does not see the signal in its entirety. Also, the size of data blocks (chunks) is freely configurable by the user, although restricted to limits stemming from recording hardware properties. This implies that Process scripts may not depend on a certain data block size, and will sometimes need to maintain their own data buffers when the algorithm in question operates on windows of data rather than continuously. In the current example, this is not the case, so we need not maintain an internal data buffer. Still, we need to maintain an internal state between calls to the Process handler in order to preserve the state of IIR filter delays. This will allow continuous operation on the signal with the same processing result as in the offline version of the algorithm.

The algorithm's online version in the bci_Process script will thus be:

function out_signal = bci_Process( in_signal )

% Process a single block of data by applying a filter to in_signal, 
% and return the result in out_signal.
% Signal dimensions are ( channels x samples ).

% We use global variables to store classifier, 
% filter coefficients and filter state.
global a b z lpa lpb lpz classifier;

[out_signal, z] = filter( b, a, in_signal, z );
out_signal = abs( out_signal );
[out_signal, lpz] = filter( lpb, lpa, out_signal, lpz );
out_signal = out_signal * classifier;

Initialize

Determination of filter coefficients is part of per-run initialization, and occurs in the Initialize event handler:

function bci_Initialize( in_signal_dims, out_signal_dims )

% Perform configuration for the bci_Process script.

% Parameters and states are global variables.
global bci_Parameters bci_States;

% We use global variables to store classifier vector, 
% filter coefficients and filter state.
global a b z lpa lpb lpz classifier;

% Configure the Bandpass filter
band = str2double( bci_Parameters.Passband ) / str2double( bci_Parameters.SamplingRate );
[n,Wn]=buttord(band*[0.9 1.1]/2,band*[0.7 1.4]/2,1,60);
[b,a]=butter(n,Wn);
z=zeros(max(length(a),length(b))-1,in_signal_dims(1));

% Configure the Lowpass filter
[lpb,lpa]=butter(1,band/4);
lpz=zeros(max(length(lpa),length(lpb))-1,in_signal_dims(1));

% Configure the Classifier vector
classifier = str2double( bci_Parameters.ClassVector );

StartRun

In addition, we need to reset filter state at the beginning of each run, using the StartRun event handler:

function bci_StartRun

% Reset filter state at the beginning of a run.

global z lpz;
z   = zeros(size(z));
lpz = zeros(size(lpz));

Constructor

To complete the Matlab filter code, we need to declare the Band and Classifier parameters in a Constructor event handler:

function [ parameters, states ] = bci_Construct

% Request BCI2000 parameters by returning parameter definition 
% lines as demonstrated below.

parameters = { ...
  [ 'BPClassifier float Passband= 10 10 0 % % // Bandpass frequency in Hz' ] ...
  [ 'BPClassifier matrix ClassVector= 1 1 1 0 % % // Linear classifier vector' ] ...
};

Preflight

Finally, we need to check parameters for consistency in a Preflight script, and declare the size of our filter's output signal:

function [ out_signal_dim ] = bci_Preflight( in_signal_dim )

% Check whether parameters are accessible, and whether
% parameters have values that allow for safe processing by the 
% bci_Process function.
% Also, report output signal dimensions in the 'out_signal_dim' argument.

% Parameters and states are global variables.
global bci_Parameters bci_States;

band = str2double( bci_Parameters.Passband ) / str2double( bci_Parameters.SamplingRate );
if( band <= 0 )
  error( 'The passband parameter must be greater zero' );
elseif( band > 0.5 / 1.4 )
  error( 'The Passband parameter conflicts with the sampling rate' ); 
end

out_signal_dim = [ 1, size( in_signal_dim, 2 ) ];
if( in_signal_dim( 1 ) ~= size( bci_Parameters.ClassVector, 2 ) )
  error( 'ClassVector length must match the input signal''s number of channels' );
end

See also

Programming Reference:MatlabFilter, User Reference:MatlabFilter, Contributions:FieldTripBuffer