Discrepancies between bci2000chain and real-time BCI2000

Forum for discussion on different user applications
Locked
jdwander
Posts: 6
Joined: 04 Apr 2012, 15:55

Discrepancies between bci2000chain and real-time BCI2000

Post by jdwander » 04 Feb 2014, 19:30

All,

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;
2) Build CursorTask

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));

boulay
Posts: 382
Joined: 25 Dec 2011, 21:14

Re: Discrepancies between bci2000chain and real-time BCI2000

Post by boulay » 08 Feb 2014, 11:07

I'm sorry I don't have time to test everything myself. I think the first step is to try to isolate the problem.

Is it a discrepancy on data load?
Write a simple logging filter.
Write, make & build a command line filter tool based on this logging filter.
Write, make & build a signal processing module using only this filter.
Test your method using the above. If the data logged by the module are different to the data loaded by the bci2000 filter chain then you can stop here and please report that. If they are the same, then continue.

Is it a problem with the filtering pipeline?
Repeat the above process with increasingly complex Signal Processing modules, starting with the above then add filters one at a time.
Use BCI2000\src\core\SignalProcessing\Spectral\PipeDefinition.cpp as a guideline.
Make sure your filter chain includes only the filters used in your module.

Did I understand correctly that the output data are identical when you don't use CAR? If that's true then it is likely a problem with the spatial filter, but please do the above process to isolate it down to that filter.

By the way, ARSignalProcessing is deprecated since its functionality is entirely contained within SpectralSignalProcessing. Simply use the SpectralSignalProcessing module and choose to use the auto regressive spectral estimation instead of FFT. (AR might even be default). This should not affect your problem.

Locked

Who is online

Users browsing this forum: No registered users and 28 guests