BCI2000 + MATLAB = Valence and Arousal

Forum for discussion on different signal processing algorithms
Post Reply
lavincent
Posts: 11
Joined: 10 Mar 2016, 06:25

BCI2000 + MATLAB = Valence and Arousal

Post by lavincent » 16 May 2016, 15:01

Hello all,

I'm a software engineering student and this is the first time I use BCI2000. I've often read usefull tips on this forum so i would ask help to you:
as in the topic, I'm trying to create a BCI that recognizes brain signal and obtains valence and arousal.

I'm using EMOTIV Epoc headset as signal source and I found convenient to use batch file. To start I've downloaded Emotiv.exe from contrib that I've used to replace signal source in the cursorTask file batch. I've fixed some issues (channel number, and other) of headset and currently i can successfully use a batch file with my params.

It contains this lines:

...

Change directory $BCI2000LAUNCHDIR
Show window; Set title ${Extract file base $0}
Reset system
Startup system localhost
Start executable Emotiv --local
Start executable SpectralSignalProcessing --local
Start executable CursorTask --local
Wait for Connected
Load parameterfile "../parms/LaV/EmotivEpoc.prm"


Everything works.
My question, now is, how can I obtain valence and arousal from that? I've read on "A Pratical Guide to Brain-Computer Interfacing with BCI2000" some info about classifier but I'm not sure how to use it. For this step, I would use MATLAB but, reading the book, I've made confusion between MEX file and MatlabSignalProcessing.

Can you help me how to procede now?

Then I'll ask you also some info about AppConnector (exactly paragr. 6.4.5)

Thanks

pbrunner
Posts: 344
Joined: 17 Sep 2010, 12:43

Re: BCI2000 + MATLAB = Valence and Arousal

Post by pbrunner » 17 May 2016, 10:16

Lavincen,

the following steps are necessary to achieve your goals:

(1) Literature research on Valence and Arousal to figure out what features in the EEG represent these states.

(2) Record screening data in which you pace a few subjects through the different states (e.g., Valence and Arousal). Make sure they remain still during the task, as otherwise the data will be poisoned with EMG data, which can create a confound.

(3) Load the data into MATLAB for post-hoc analysis using the load_bcidat function. Extract your EEG features and train a linear classifier/regression between the two states. I would recommend using either a regularized least square (lasso) or elastic net regularization (lassoglm) in MATLAB.

(4) Write the feature extraction parameters and linear classifier weights into a BCI2000 parameter structure. Write out the parameters into a BCI2000 parameter fragment using the convert_bciprm mex function.

(5) Setup a BCI2000 signal processing filter chain that implements the same steps as you have done in your post-hoc MATLAB analysis.

(6) Load these parameters into your online system and use the CursorTask as the application.

Regards, Peter

lavincent
Posts: 11
Joined: 10 Mar 2016, 06:25

Re: BCI2000 + MATLAB = Valence and Arousal

Post by lavincent » 18 May 2016, 16:16

Thank you very much for your reply,
but i'm so inexpert to do what you said :(

I'll ask to my thesis teacher some info about valence and arousal to figure them out from the EGG.
In addition, I need to do an online analysis. For this, I'm trying to use MatlabSignalProcessing. This module can correctly open the MATLAB console but I don't know hot to use it.
I would use a SVM classifier, can I do something like in page 114?

Then i need to export the results in a C# application.
Have I use AppConnector (TCP or UDP) to do this?


Thanks and sorry for my inexperience.. I'm really interested in this world but currently Iìm groping in the dark.

pbrunner
Posts: 344
Joined: 17 Sep 2010, 12:43

Re: BCI2000 + MATLAB = Valence and Arousal

Post by pbrunner » 18 May 2016, 17:04

Lavincent,

I would suggest that you first study the literature on what is known about features in the EEG that represent these emotional states. Once you know the features and their scalp location, you will know whether you can even capture them with your EMOTIV Epoc headset. If the EMOTIV headset is suitable, you then will try to reproduce the effect in these features using your own experiment. Only then you will start thinking about how to translate this into an online experiment.

Please follow this link for your literature research:
http://www.ncbi.nlm.nih.gov/pubmed/?ter ... rousal+eeg

Regards, Peter

lavincent
Posts: 11
Joined: 10 Mar 2016, 06:25

Re: BCI2000 + MATLAB = Valence and Arousal

Post by lavincent » 21 May 2016, 10:28

Thanks Peter,
I can't believe to find someone so patient! :)

I've done some research finding that the values of valence and arousal are (meanly) in the alpha band peaks, so:

1) I need to "clean" my EMOTIV Epoc headset signal with
-a pass low filter at 0.1Hz
-and a pass high filter at 30Hz. (can be this right? tell me if is this wrong in your experience)

2) I have to apply a CAR (Indipendent Component analysis)

3) the channels I'm interested in are:
AF3, AF4, F3, F4, FC5, FC6

Supposing that my results are correct, how can I continue, now?


Thanks a lot.

pbrunner
Posts: 344
Joined: 17 Sep 2010, 12:43

Re: BCI2000 + MATLAB = Valence and Arousal

Post by pbrunner » 22 May 2016, 17:32

Lavincent,

the next step is to design an experiment that paces the subject through different emotions. For this you can use BCI2000 with your source module, DummySignalProcessing and StimulusPresentation. Implement an interleaved block design with at least 50 trials for each condition. Each trial should be at least 10 seconds long and there should be a 10 second long inter-stimulus interval. Load the recorded data into MATLAB and perform the signal processing and classification as you have decided on from your research. If you can find a reliable effect, i.e., one that survives randomization tests and cross validation, you can proceed to translating the task into an online BC2000 experiment.

Please note that you would preferable use a recording system that is capable to record more channels and most importantly EMG (to rule out a confound) before you decide to use the Emotiv headset.

Regards, Peter

lavincent
Posts: 11
Joined: 10 Mar 2016, 06:25

Re: BCI2000 + MATLAB = Valence and Arousal

Post by lavincent » 25 May 2016, 14:41

Thanks again Peter.

I've load my experiments (done using DummySignalProcessing and StimulusPresentation) into Matlab using load_bcidat function, getting a lot of data.
But I'm so inexperienced with Matlab and now I don't know how to execute yours second tip: "...and perform the signal processing and classification as you have decided on from your research".

I've read about fitcsvm for Matlab, is this functions relevant?

If useful, as I think you already know, exists the IAPS (International Affective Picture System) database for emotive image.
Anyway, in my thesis work, is not strictly requested to train a classifier but, if needed, I would be grateful if you could help me.

Exactly, how can I now procede?

Thank you so much.

pbrunner
Posts: 344
Joined: 17 Sep 2010, 12:43

Re: BCI2000 + MATLAB = Valence and Arousal

Post by pbrunner » 25 May 2016, 16:38

Lavincent,

I assume you loaded the data in the following way:

Code: Select all

[ signal, states, parameters ] = load_bcidat( settings.filename);
What is relevant for your analysis are signal and states.StimulusCode. What you probably want is the band-power for different frequency bands in your signal to calculate a linear model that predicts the emotional state coded in the states.StimulusCode. I will paste a few code snippets that should help you.

Regards, Peter

Code: Select all

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1, '> CALCULATE FILTER COEFFICIENTS FOR IIR BANDPASS FILTER BANK \n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f_modulated = zeros(size(settings.frequency_terms));

for idx = 1:length(settings.frequency_terms),
    eval(['f_modulated(idx) = ',settings.frequency_terms{idx},';']);
end

for idx = 1:length(f_modulated),

    % define passband, stopband and ripples
    bandpass{idx}.Wp = [f_modulated(idx)-settings.bw_pass f_modulated(idx)+settings.bw_pass]/(parameters.SamplingRate.NumericValue/2); 
    bandpass{idx}.Ws = [f_modulated(idx)-settings.bw_stop f_modulated(idx)+settings.bw_stop]/(parameters.SamplingRate.NumericValue/2);
    bandpass{idx}.Rp = 3; 
    bandpass{idx}.Rs = 30;

    % calculate the minimum filter order
    [bandpass{idx}.n,bandpass{idx}.Wn] = buttord(bandpass{idx}.Wp,bandpass{idx}.Ws,bandpass{idx}.Rp,bandpass{idx}.Rs);
    bandpass{idx}.n = bandpass{idx}.n + rem(bandpass{idx}.n,2);

    % caclulate the filter coefficients in Zero-Pole-Gain design
    [bandpass{idx}.z, bandpass{idx}.p, bandpass{idx}.k] = butter(bandpass{idx}.n,bandpass{idx}.Wn,'bandpass');
    [bandpass{idx}.sos,bandpass{idx}.g]=zp2sos(bandpass{idx}.z,bandpass{idx}.p,bandpass{idx}.k);
    bandpass{idx}.h=dfilt.df2sos(bandpass{idx}.sos,bandpass{idx}.g);
    
end

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1, '> VISUALIZE IIR FILTER BANK \n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for idx=1:length(bandpass),
    vis.bandpass.h(idx) = bandpass{idx}.h;
    vis.bandpass.legend{idx} = sprintf('%2.2f-%2.2f Hz',f_modulated(idx)-settings.bw_pass,f_modulated(idx)+settings.bw_pass);
end

h=fvtool(vis.bandpass.h,'FrequencyScale','linear','Fs',parameters.SamplingRate.NumericValue);
set(h,'NumberofPoints',2^16);
xlim([0,0.1]);
legend(vis.bandpass.legend,'FontSize',22);
title('Band-pass filters','FontSize',18,'FontWeight','bold');

Code: Select all

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1, '> EXTRACT FEATURES \n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% decimate data by a factor
decimation_stepping = round(parameters.SampleBlockSize.NumericValue/settings.decimation_factor);

% fix issue at the beginning of the StimulusCode state
states.StimulusCode(1:parameters.SampleBlockSize.NumericValue) = states.StimulusCode(parameters.SampleBlockSize.NumericValue+1);
StimulusCode = states.StimulusCode(1:decimation_stepping:end);

% create correlation variable
index_task   = find(StimulusCode == 2 | StimulusCode == 4);
corr_var(StimulusCode == 2) = -1;
corr_var(StimulusCode == 4) = 1;

% limit correlation variable to only task period
corr_var_task = corr_var(index_task);

% determine start and stop of each block
index_start = find(diff(double(StimulusCode == 2 | StimulusCode == 4)) >= 1)+1;
index_stop  = find(diff(double(StimulusCode == 2 | StimulusCode == 4)) <= -1);

% create correlation variable for task
corr_var_trial = double(StimulusCode(index_stop) - 3);

% determine size of feature space
num_frequencies = length(f_modulated);
num_features    = num_channels * num_frequencies;

% allocate feature variable
signal_power   = zeros(length(index_task),num_features);

% allocate markers for mapping between feature index and corresponding channel and frequency
list_channel      = zeros(num_features,1);
list_frequency    = zeros(num_features,1);
list_channel_used = zeros(num_features,1);

% calculate envelope for each filter 
for idx_channel = 1:num_channels,
    
    signal_channel = signal(:,idx_channel);
    
    for idx_frequency = 1:length(f_modulated), 
        
        % create decimated envelope of band-pass filtered signal for this channel
        signal_var = decimate(abs(hilbert(filtfilt(bandpass{idx_frequency}.sos,bandpass{idx_frequency}.g,double(signal_channel-signal_channel(1))))),decimation_stepping);
 
        % determine index of the feature (i.e., channels by frequency)
        idx_feature                 = (idx_frequency-1)*num_channels+idx_channel;
        
        % mark which frequency and channel this feature corresponds to
        list_channel(idx_feature)      = idx_channel;
        list_frequency(idx_feature)    = idx_frequency;
        
        % mark if this feature is in the set of channels that should be use for the model 
        if (~isempty(intersect(idx_channel,settings.list_channel_index))),
            list_channel_used(idx_feature) = 1;
        end
        
        % store features for this feature/channel combination
        signal_power(:,idx_feature) = signal_var(index_task);
     

    end 
end

Code: Select all

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1, '> CALCULATE LINEAR MODEL AND PLOT MODEL PREDICTION \n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

index = find(list_channel_used);

switch settings.method
    
    case 'least squares fit'
        b = pinv(signal_power(:,index))*corr_var_task;
        stats.intercept = -mean(signal_power(:,index)*b);   
        inmodel = ones(size(b));
        
    case 'stepwise regression'       
        [b,se,pval,inmodel,stats,nextstep,history] = stepwisefit(signal_power(:,index),corr_var_task,'penter',1e-8,'premove',0.01);
        b(~inmodel) = 0;        

    case 'regularized least square'
        [b,FitInfo] = lasso(signal_power(:,index),corr_var_task,'CV',10,'NumLambda',20);
        lassoPlot(b,FitInfo,'plottype','CV');
        idx = FitInfo.IndexMinMSE;
        b = b(:,idx);
        stats.intercept = FitInfo.Intercept(idx);
        inmodel = b~=0;
        
    case 'elastic net regularization'
        [b,FitInfo] = lassoglm(signal_power(:,index),corr_var_task,'normal','CV',10,'NumLambda',20);
        lassoPlot(b,FitInfo,'plottype','CV');
        idx = FitInfo.IndexMinDeviance;
        b = b(:,idx);
        stats.intercept = FitInfo.Intercept(idx);
        inmodel = b~=0;

    otherwise

        error('settings.method = %s is unkown',settings.method);
        
end

response = signal_power(:,index)*b+stats.intercept;
plot(response), hold on
plot(corr_var_task,'r-');
legend({'response','task'},'FontSize',18);

Code: Select all

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1, '> WRITE BCI2000 PARAMETER FILE FOR REAL-TIME EXPERIMENT \n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Classifier
prm.Classifier.Section      = 'Filtering';
prm.Classifier.Type         = 'matrix';
prm.Classifier.DefaultValue = '';
prm.Classifier.LowRange     = '';
prm.Classifier.HighRange    = '';
prm.Classifier.Comment      = 'classifier settings';
prm.Classifier.Value        = cell(sum(inmodel)*(parameters.SampleBlockSize.NumericValue/decimation_stepping),4);
prm.Classifier.RowLabels    = cell(sum(inmodel),1);
prm.Classifier.ColumnLabels = cell(1,4);

prm.Classifier.ColumnLabels{1}  = 'input channel';
prm.Classifier.ColumnLabels{2}  = 'input element (bin)';
prm.Classifier.ColumnLabels{3}  = 'output channel';
prm.Classifier.ColumnLabels{4}  = 'weight';

index = find(inmodel);

counter = 1;
for idx=1:length(index),
   
    for idx_sample = decimation_stepping:decimation_stepping:parameters.SampleBlockSize.NumericValue,
    
        prm.Classifier.RowLabels{counter} = sprintf('%d',idx);
        prm.Classifier.Value{counter,1}      = sprintf('%d',index(idx));
        prm.Classifier.Value{counter,2}      = sprintf('%d',idx_sample);
        prm.Classifier.Value{counter,3}      = '1';
        prm.Classifier.Value{counter,4}      = sprintf('%1.8d',b(index(idx))/settings.decimation_factor);

        counter = counter + 1;

    end
    
end

parameter_lines = convert_bciprm( prm );

fid = fopen('filter_bank.prm','w');

for idx=1:length(parameter_lines),
    fprintf(fid,'%s',parameter_lines{idx});
    fprintf(fid,'\r\n');
end

fclose(fid);

lavincent
Posts: 11
Joined: 10 Mar 2016, 06:25

Re: BCI2000 + MATLAB = Valence and Arousal

Post by lavincent » 25 May 2016, 19:15

Thank you Peter,
yes, I've used load_bcidat but what "settings" is?

You have been infinitely kind. I don't want disturb you over.
I will try to use your code and solve all future issues.

Best regards,
LaV

pbrunner
Posts: 344
Joined: 17 Sep 2010, 12:43

Re: BCI2000 + MATLAB = Valence and Arousal

Post by pbrunner » 26 May 2016, 10:03

Lavincent,

the settings is just a structure where you keep your global variables.

Code: Select all

settings.filename            = 'bci2000_data_file.dat';
settings.method              = 'least squares fit';                 %'stepwise regression' | 'least squares fit' | 'regularized least square' | 'elastic net regularization' 
settings.list_reference      = [6];
settings.frequency_terms     = {'f1','f2','2*f1','2*f2'};           %[f1,f2,f1+f2,abs(f1-f2),2*f1-f2,2*f2-f1,2*f1,3*f1,2*f2,3*f2]
settings.bw_pass             = 0.1; % Hz 
settings.bw_stop             = 0.5; % Hz
settings.decimation_factor   = 10;
settings.list_channel_names = {'P3','Pz','P4','PO7','PO8','Oz'};    %  {'F3','Fz','F4','T7','C3','Cz','C4','T8','CP3','CP4','P3','Pz','P4','PO7','PO8','Oz'};
Regards, Peter

lavincent
Posts: 11
Joined: 10 Mar 2016, 06:25

Re: BCI2000 + MATLAB = Valence and Arousal

Post by lavincent » 27 May 2016, 19:42

Some news,
my teacher says that mine, is a three-year thesis so I've not implement the classifier.

Thanks really a lot, Peter, for your help, but I need an already done classifier. Suggestion from where can i take it? (If exists, 'f course).
In addition, I've found that someone have used WEKA. Do you know it? Do you think this tool would be easier for me?

Then I don't know how to train a classifier, neither with your tips.


Thanks again.

pbrunner
Posts: 344
Joined: 17 Sep 2010, 12:43

Re: BCI2000 + MATLAB = Valence and Arousal

Post by pbrunner » 01 Jun 2016, 16:48

Lavincent,

if you look in code fragment above for the section labeled "CALCULATE LINEAR MODEL AND PLOT MODEL PREDICTION" you will find examples for how to use an 'existing' classifier.

Regards, Peter

Post Reply

Who is online

Users browsing this forum: No registered users and 1 guest