Category Archives: Matlab

First steps with a €20 single-channel EEG system

My friend Vladimir recently demonstrated a single-channel EEG system that he got at a hackathon in London. When he mentioned that it only costs €20 (or actually 20 GBP to be more precisely) I immediately decided to order one myself. The ICI-BCI system is a low cost open source brain computer interface.

The bag clearly and rightfully indicates that it is a totally experimental system, and that it should be used with caution.

The basic idea of the amplifier is that it takes an 1000 Hz analog audio signal from the computer or mobile phone, which is amplitude modulated by the ExG signal and subsequently fed back as microphone signal. So the system is fully analog and requires the DAC to be done by the audio input of the computer or phone.

Continue reading

Peer-to-peer distributed Matlab computing – update

After discussing in detail with colleagues at the Donders and at the FIL, I have implemented the peer-to-peer distributed computing toolbox for MATLAB. Most of the desired functionality is now in place, and it seems to work robustly and efficiently.

The peer toolbox allows you to do something like this in MATLAB

a = randn(400,400);
tic; cellfun('pinv', {a, a, a, a, a}, 'UniformOutput', false); toc
tic; peercellfun('pinv', {a, a, a, a, a}, 'UniformOutput', false); toc

Continue reading

Peer-to-peer distributed Matlab computing

In a recent meeting with the SPM developers, we discussed parallel computing using the Matlab distributed computing toolbox, Star-P, Sun Grid Engine, and other batch systems that can be linked to Matlab. These are all limited in their usefulness for the typical neuroimaging research setting in that they are based on a centralized job distribution system. That may work fine on a large cluster with a centralized configuration and system administration, but even then the usefullness is limited because all input and output data (which are typically large) have to be send over the network twice: first to the job manager, then to the compute node (and vice versa for the results).

To resolve some of these problems, I came up with the idea of peer-to-peer distributed computing in Matlab. The full description can be found on http://fieldtrip.fcdonders.nl/development/peer

Automatically compile a missing mex file on the fly

I am maintaining a few open source Matlab toolboxes. Although most functions are plain Matlab, some functions are implemented as mex file. The mex files have to be compiled on each platform, i.e. 32 and 64 Linux, 32 and 64 bit Windows, ppc and intel Mac OS X, etc. Since I don’t have access to all those platforms, I have devised a trick to automatically compile the mex files on the end-user’s platform.

Let’s say that you have a function mymexfile, then the file mymexfile.c would contain the c-code implementation of the Matlab mex file. You should save the following code as mymexfile.m. If the compiled mex file (e.g. mymexfile.mexw32) does not exist upon the first call to the function, the m-file will be executed. The m-file will automatically compile the c-file and will subsequently call the compiled mex file. On all subsequent calls, the compiled mex file will be executed directly.

function [varargout] = autocompile(varargin)

% AUTOCOMPILE compile the missing mex file on the fly

% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

try
% try to compile the mex file on the fly
warning('trying to compile MEX file from %s', mexsrc);
cd(mexdir);
mex(mexsrc);
cd(pwdir);
success = true;

catch
% compilation failed
disp(lasterr);
error('could not locate MEX file for %s', mexname);
cd(pwdir);
success = false;
end

if success
% execute the mex file that was juist created
funname = mfilename;
funhandle = str2func(funname);
[varargout{1:nargout}] = funhandle(varargin{:});
end

Multithreading in a MATLAB mex file

Together with colleagues from the F.C. Donders Centre and the NICI, I am currently working on a project on real-time analysis of EEG and MEG data for Brain-Computer Interfacing (BCI). Since we already have a lot of relevant code in MATLAB and since a lot of the target end-users are familiar with it, we would like to do the real-time processing in MATLAB.

The idea is to read the data that is streaming from the EEG amplifier, do the computation on the data (e.g. spatial filtering, spectral decomposition and classification) and then continue with the next section of the data. The challenge is not to lose any of the streaming EEG data while MATLAB is busy with the computation. A potential solution would be to implement the reading/acquisition of the data into matlab in a mex file which is multi-threaded.

For this purpose I wrote some code for testing the multithreading in MATLAB. The test code below demonstrates that it works, and shows how easy it actually is.

I am working on an Apple PowerBook G4 with OS X 10.4.11 and with MATLAB Version 7.4 (R2007a). The example code uses the POSIX Pthread libraries and I would expect the same code to work on most Unix and Linux systems. It will not work on a windows system, but might work with pthreads-win32.

Put the following code in a file with the name threads.c

#include < pthread.h > /* for threading */
#include < unistd.h > /* for sleep */
#include "mex.h"
#include "matrix.h"

void *do_thread(void *threadid)
{
int tid;
tid = (int)threadid;
mexPrintf("In thread: just started thread #%dn", tid);
sleep(4); /* pretend that the code is busy for a few seconds */
mexPrintf("In thread: about to exit from thread #%dn", tid);
pthread_exit(NULL);
}

void mexFunction (int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
int num;

if (nrhs<1)
mexErrMsgTxt("not enough input arguments");
else
num = mxGetScalar(prhs[0]);

pthread_t threads[num];
int rc, t;
for(t=0; t < num ; t++){ mexPrintf("In main: creating thread %dn", t); rc = pthread_create(&threads[t], NULL, do_thread, (void *)t); if (rc){ mexErrMsgTxt("problem with return code from pthread_create()"); } sleep(1); /* wait some time before making the next thread */ } return; }

Subsequently, in MATLAB you should compile the mex file using the following command:

>> mex threads.c -lpthread

After that, you can run the mex function. It takes a single argument, corresponding to the number of threads to create. Each thread will live for three seconds (pretending that it is busy for some time), and the time between threads being created is one second. Some example output is below.

>> threads(6)
In main: creating thread 0
In thread: just started thread #0
In main: creating thread 1
In thread: just started thread #1
In main: creating thread 2
In thread: just started thread #2
In main: creating thread 3
In thread: just started thread #3
In thread: about to exit from thread #0
In main: creating thread 4
In thread: just started thread #4
In thread: about to exit from thread #1
In main: creating thread 5
In thread: just started thread #5
In thread: about to exit from thread #2
>> In thread: about to exit from thread #3
In thread: about to exit from thread #4
In thread: about to exit from thread #5
>>

You can see each thread starting and stopping. The last threads stop after the control has been returned to the MATLAB prompt.

Here you can find some comments from Mathworks, including some windows specific example code. Important to note is that MATLAB itself is single-threaded, and not thread safe. Hence you should take care in interfacing your threaded c-code in the mex file with the MATLAB core.

I would appreciate it to get some feedback on this. Please leave a comment below if you tested it successfully or if it does not work for you, with some details on your configuration (like operating system and MATLAB version).

EEProbe

EEProbe is a Linux software package for EEG acquisition and analysis that was originally develloped at the Max Planck Institute in Leipzig, Germany (there known as EEP), and that is now further developed as a commercial product by ANT Software in Enschede, The Netherlands.

EEProbe stores its continuous data in a file with the extension *.cnt (not to be confused with the Neuroscan *.cnt format). The *.cnt file is a binary file and the data in it is compressed, which makes reading the data in your own software quite difficult. The file format for averaged ERPs (*.avr) is also binary, but not compressed. Furthermore there are files that contain event information such as triggers (*.trg), rejection marks (*.rej) and some others.

ANT has released a library with c-code (cntopenlib) to its customers, containing functions to read the header, the binary data and to decompress the data. I have used that library to make two mex files to import the binary data into Matlab, you can download those functions here. Furthermore, at the bottom of this page you can also download the source code of the mex files.

Note: since the release of EEProbe 3.2.4 in March 2004 the binary format of the *.avg files has been slightly changed (they now also contain history information). No new version of cntopenlib has yet been released, and therefore the mex function cannot always read the averages properly. You can use the “avrstrip” program from EEProbe to remove the history information, the resulting file can be read into Matlab.

Plain Matlab functions:

Mex files for Mac OS X (created using Matlab 6.5):

Mex files for Solaris (created using Matlab 6.1):

Mex files for Linux (created using Matlab 6.1):

Mex files for Windows (created using Matlab 6.1):

Source code of the mex files, released under the BSD license:

Or you can download them all together in a single package here.

Matlab compiler #2

I discovered another problem with the matlab compiler (version 3, which is included with Matlab 6.5), namely that it will not translate an m-file to c-code when that m-file tries to read a variable from a mat-file.

Mathworks does specify that one of the things that is not supported by the compiler (see here) are m-files that dynamically name variables to be loaded or saved. They give this example which is disallowed by the compiler:

x= 'f';
load('foo.mat',x);

However, this function also will not compile:

function testload(cfg);
b = load(cfg.filename);
c = getfield(b, 'a');
disp(c);

giving the following error

>> mcc -m testload
testload.c: In function `Mtestload':
testload.c:103: error: `mlxLoadStruct' undeclared (first use in this function)
testload.c:103: error: (Each undeclared identifier is reported only once
testload.c:103: error: for each function it appears in.)
mbuild: compile of 'testload.c' failed.

If I slightly rewrite the function it does work. The alternative function that has the desired functionality and that does compile correctly is

function testload(cfg);
filename = cfg.filename;
b = load(filename);
c = getfield(b, 'a');
disp(c);

At least this gives me a handle on how to modify my code with little effort to make it compatible with the Matlab compiler.

Matlab compiler

Currently I am testing out the use of my mentat toolbox for parallel computing within Matlab in combination with GridEngine on our linux cluster. Initially it worked fine, but that was because I was testing it with simple level Matlab functions like “rand” and “fft”. However, I am now trying to get it to work with more complex real-world scripts for EEG and MEG analysis that use the FieldTrip toolbox and I am running into a few problems:

I had to edit the file /opt/matlab-6.5.1/toolbox/matlab/iofun/private/readavi.m and rename the output argument “varargin” into “varargout”. I think that this is a plain bug in that function, which only becomes apparent after translating the code to c. The original file resulted in a compilation error, changing the output variable name fixed that.

Somehow, the code wants to include the Matlab function wk1read in the stand-alone executable. I cannot determine where that happens, because all functions are very deeply nested: mcc generates ~350 c-files on which my analysis function depends. The problem is that, after compilation, wk1read results in a unresolved symbol error. Since I am not interested in reading Lotus database stuff (which is what that particular function does), I have simply replaced it on my own Matlab search path with an alternative (empty) implementation for wk1read. That solved the unresolved symbol error.

Furthermore, I encountered a problems with the exist function. In the stand-alone executable, it does not seem to detect the existence of a variable. Since I am using that function on multiple locations throughout the code, I have to think of a nice workaround. Related to this problem I found this page which gives some relevant hints on using the Matlab compiler:

  1. The following things cause a warning in compilation and may cause a runtime error if the offending statement is executed: objects, most toolboxes, inline, whos, some graphics functions that are passed functions by string.
  2. Program compiled with graphics will not run without X windows. This includes running on beowulf nodes from PBS scripts.
  3. Program containing graphics commands and compiled without graphics will produce compile time warning but will run fine as long as the graphics commands are not executed.
  4. The following commands work differently: eval, feval, exist, global.
  5. Use handles to pass functions rather than strings.
  6. Do not use exist to check for the presence of arguments in a function.