I have been trying to use the bci2000chain tool in Matlab to perform offline analyses of BCI2000 recordings and I've been running in to some difficulties, specifically inconsistencies between what is being returned from the bci2000chain and what I believe is being passed through the real-time processing pipeline of BCI2000. I have tried to distill the problem to the basic elements, so any errors that I'm making will hopefully be very visible to the trained eye. Before getting in to reproduction instructions, here is some information about my system
I'm running Win 7 pro x64, BCI2000 v3.0.5, and building 32-bit binaries using VS2008 pro. Using Matlab R2013a.
The issue: If I run BCI2000 using SignalGenerator | ARSignalProcessing | CursorTask, and compare the ControlSignal values that are arriving in the FeedbackTask::Processing(...) function with what I get out of the same signal processing pipeline run through the bci2000chain tool, I observe large discrepancies when common average re-referencing. Similarly, when using FilePlayback | ARSignalProcessing | CursorTask to play back previously recorded neural data (ECoG), I observe the same issue.
My suspicion is that I'm running in to a numeric precision issue in one of the two versions of the code path, which is exacerbated when using CAR, because the resultant signal that goes through the rest of the signal processing pipeline is very small (especially in the case of the SignalGenerator where the only difference between channels is the additive noise). It is also a distinct possibility that I've misunderstood the way these tools are meant to be used.
At the end of the day, my goal is simply to get a faithful offline reproduction of exactly the same values of ControlSignal that are being passed in to the cursor task. Any advice is greatly appreciated!
Reproduction instructions:
1) In Visual Studio, modify FeedbackTask.cpp (BCI2000FrameworkAppModule/Source/shared/modules/application/FeedbackTask.cpp), specifically FeedbackTask::Process( ... ), after the first if block, by adding the following line:
Code: Select all
bciout << "process: " << Input(0, 0) << " " << Input(1, 0) << endl;
3) From the BCI2000 Launcher launch the following module set: SignalGenerator | ARSignalProcessing | CursorTask
4) In the operator window, click config, select the filtering tab, and change the values in the Adaptation field to "0 0" (without quotes).
5) Click set config, then start.
6) Allow BCI2000 to run for at least 15 seconds, then click pause.
7) Copy everything in the system Log to a file, note the file path
8 ) Determine the path to the file that was recorded: in my case this was the most recent dat file in
<BCI2000_ROOT>\data\Name001\
9) Run the matlab script, Comparator.m (source below) selecting the appropriate log file and BCI2000 dat file that you noted in steps 7 and 8 (respectively).
10) Verify that the output of the real-time signal processing pathway is equivalent to the output of the bci2000chain tool
11) Now, go back to the BCI20000 and click config, select the filtering tab, and change the SpatialFilterType from "full matrix" to "common average reference CAR". Also, verify that the SpatialFilterCAROutput field is empty.
12) Repeat steps 5-9, copying only the new lines from the BCI2000 log window in to a separate text file. Note now that there is a significant difference between the output from the real-time pathway and the bci2000chain tool
Code: Select all
%% Source code for Comparator.m
% J Wander 2/4/2014
%% path setup
try
bci2kRoot = bci2000path;
catch foo
bci2kRoot = uigetdir(pwd, 'Please select the root directory of your bci2000 installation');
end
addpath(fullfile(bci2kRoot, 'tools', 'mex'));
addpath(fullfile(bci2kRoot, 'tools', 'matlab'));
managepath('-addtosystempath', fullfile(bci2kRoot, 'tools', 'cmdline'));
if (isempty(which('bci2000path')))
error('Please add the bci2000 matlab tools directory to your matlab path.');
end
if (isempty(which('load_bcidat')))
error('Please add the bci2000 matlab mex directory to your matlab path.');
end
%% get the paths to the files of interest
[recFilename, recFilepath] = uigetfile('*.dat', 'Select a BCI2000 Recording file', bci2000path);
recFile = fullfile(recFilepath, recFilename);
[logFilename, logFilepath] = uigetfile('*.*', 'Select a BCI2000 console output log file', pwd);
logFile = fullfile(logFilepath, logFilename);
%% load in the parameters used in the recording
[~,~,par] = load_bcidat(recFile);
parstr = make_bciprm(par);
% and write them out to a file
handle = fopen('replay.prm', 'w');
fwrite(handle, parstr);
fclose(handle);
%% run the bci2000chain, note the omission of the ExpressionFilter, should
% be ok because it is setup as a passthrough by default.
output = bci2000chain(recFile, 'TransmissionFilter|SpatialFilter|ARFilter|LinearClassifier|LPFilter|Normalizer', '-2', 'replay.prm');
synthesized = output.Signal;
%% now parse the log file you just created
handle = fopen(logFile,'r');
line = fgetl(handle);
realtime = [];
while line > 0
temp = regexp(line, '.*?ess: ([0-9\. ]+).', 'tokens');
if (~isempty(temp))
realtime(end+1,:) = sscanf(temp{1}{1}, '%f %f %f %f');
end
line = fgetl(handle);
end
fclose(handle);
%% force the same length, for easy comparison
len = min(size(realtime, 1), size(synthesized, 1));
synthesized = synthesized(1:len,:);
realtime = realtime(1:len,:);
%% visual comparison
figure
for c = 1:2
subplot(2,1,c);
plot(realtime(:,c));
hold on;
plot(synthesized(:,c), 'g:');
legend('realtime', 'synthesized');
title(sprintf('channel %d', c));
end
%% numeric comparison
mse = mean((realtime-synthesized).^2,1);
fprintf('Mean Squared Error by channel: %s\n', num2str(mse));
fprintf('As a percent of dynamic range: %s\n', num2str(mse ./ range(realtime,1) * 100));