Spectrally weighted CSP OV algorithm box: help to validate ?

Making & changing box plugins and external apps
Post Reply
Jeff_B
Posts: 37
Joined: Wed Apr 11, 2012 1:31 pm
Location: Nice - Alpes Maritimes

Spectrally weighted CSP OV algorithm box: help to validate ?

Post by Jeff_B »

Hi Openvibers,
Here is the pre-beta version (well... a code which still needs to be validated functionnally) of the spectrally weighted CSP along the lines of:

[1] Tomioka, R., Dornhege, G., Aihara, K., and Mueller, K.-R. "An iterative algorithm for spatio-temporal filter optimization."
In Proceedings of the 3rd International Brain-Computer Interface Workshop and Training Course 2006.

[2] Ryota Tomioka, Guido Dornhege, Guido Nolte, Benjamin Blankertz, Kazuyuki Aihara, and Klaus-Robert Mueller
"Spectrally Weighted Common Spatial Pattern Algorithm for Single Trial EEG Classification",
Mathematical Engineering Technical Reports (METR-2006-40), July 2006.

This OV implementation closely follows the (bcilab) Matlab code by Christian Kothe, UCSD (29-04-2010) available at:
https://github.com/sccn/BCILAB/blob/mas ... mSpecCSP.m

Code: Select all

function model = specCSP(X, Fs, pp, qp, prior, n_of, steps, ratio_mode, csp_mode);

% INPUT PARAMETERS :
%
% X{2}          matrix of epoched signals for the two classes, expect the classic dims: [C,S,T] = size(X{1}) = size(X{2})
% Fs            sampling frequency
% pp and qp     regularization parameters p and q
% prior         spectral prior: '@(f) f>=prior(1) & f<=prior(2)'
% n_of          half the number of CSP patterns
% steps         number of iterations for the spectral weighting estimation
% ratio_mode    (bool)use the spectral prior built from CSP on the individual cross-spectrum matrices and ietrative ratio optimization
% MI_mode       (bool)use the spectral prior built from MI labels vs features for each channel and frequency bin and iterative MI optimization 
%
% OUTPUT PARAMETERS
%
% model         a structure (can be used with my_topoplot) with the following fields: 
%               'W'     (spatial weights),
%               'P'     (inverse of the spatial 'projection'),
%               'alpha' (spectral weights),
%               'freqs' (fft frequency bins)        and 'bands'         (retained frequency bin indexes),
%               'lambda'(composite CSP eigenvalues) and 'lambda_cross'  (subband cross-spectrum CSP eigenvalues)

% re-parameterize the hyper-parameters p' and q' into p and q (pp=0 & qp=1 => p=1 & q=1)
p  = pp+qp;
q  = qp;
% spectral prior
if isnumeric(prior) && length(prior) == 2
    prior = @(f) f >= prior(1) & f <= prior(2);
end
% number of C=Channels, S=Samples and T=Trials
[C,S,T] = size(X{1});
% build a frequency table (one per DFT bin)
freqs = (0:S-1)*Fs/S;
% evaluate the prior I
I = prior(freqs);
% and find table indices that are supported by the prior
bands = find(I);
% preprocessing

for c=1:2
    % compute the per-class epoched data X and its Fourier transform (along time), Xfft
    [C,S,T] = size(X{c});
    Xfft{c} = fft(X{c},[],2);
    % the full spectrum F of covariance matrices per every DFT bin and trial of the data
    F{c} = single(zeros(C,C,max(bands),T));
    for k=bands
        for t=1:T
            F{c}(:,:,k,t) = 2*real(Xfft{c}(:,k,t)*Xfft{c}(:,k,t)');
        end
    end
    % compute the cross-spectrum V as an average over trials
    V{c} = mean(F{c},4);
end

% 1. initialize the filter set alpha and the number of filters J
 J = 1; alpha{J}(bands) = 1;                     % original specCSP

% 2. For each step
for step=1:steps
    % 3. For each set of spectral coefficients alpha{j} (j=1,...,J)
    for j=1:J
        % 4. Calculate sensor covariance matrices for each class weighted by alpha{j}
        for c = 1:2
            Sigma{c} = zeros(C);
            for b=bands
                Sigma{c} = Sigma{c} + alpha{j}(b)*V{c}(:,:,b);
            end
        end
        if csp_mode == 'CSP0'
            % 5. Common diagonalization version of the CSP problem Eq. (2)
            [Ph,D] = eig(Sigma{1}+Sigma{2});
            Ph = diag(sqrt(1./diag(D)))*Ph';
            C1 = Ph * Sigma{1} * Ph';     C2 = Ph * Sigma{2} * Ph';
            C1 = (C1+C1')/2 ;             C2 = (C2+C2')/2 ;
            [R,DD]  = eig(C1,C2);
            VV      = Ph'*R;
        elseif csp_mode == 'CSP3'
            % 5. Solve the generalized eigenvalue problem Eq. (2)
            [VV,DD] = eig(Sigma{1},Sigma{1}+Sigma{2});
        end
        % and retain n_of top eigenvectors at both ends of the eigenvalue spectrum...
        W{j} = { VV(:,1:n_of),  VV(:,end-n_of+1:end)};
        iVV  = inv(VV)';
        P{j} = {iVV(:,1:n_of), iVV(:,end-n_of+1:end)};
        % as well as the top eigenvalue for each class
        lambda(j,:) = [DD(1), DD(end)];
    end
    % 7. Set W{c} from all W{j}{c} such that lambda(j,c) is minimal/maximal over j
    W = {W{argmin(lambda(:,1))}{1}, W{argmax(lambda(:,2))}{2}};
    P = {P{argmin(lambda(:,1))}{1}, P{argmax(lambda(:,2))}{2}};
    % 8. for each projection w in the concatenated [W{1},W{2}]...
    Wcat = [W{1} W{2}]; J = 2*n_of;
    Pcat = [P{1} P{2}];
    for j=1:J
        w = Wcat(:,j);
        % 9. Calculate (across trials within each class) mean and variance of the w-projected cross-spectrum components
        for c=1:2
            % Part of Eq. (3)
            s{c} = zeros(size(F{c},4),max(bands));
            for k=bands
                for t = 1:size(s{c},1)
                    s{c}(t,k) = w'*F{c}(:,:,k,t)*w;
                end
            end
            mu_s{c} = mean(s{c});
            var_s{c} = var(s{c});
        end
        % 10. Update alpha{j} according to Eqs. (4) and (5)
        for c=1:2            
            for k=bands                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % [JeffAdd] instead of Eq. (4), trying to maximize a Rayleigh ratio <s{+}> / <s{-}> not a (normalized) difference <s{+}> - <s{-}> %          
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                if ratio_mode,
                % Eq. (4) modified for a Rayleigh ratio optim.
                    alpha_opt{c}(k) = max( 0 ,   mu_s{c}(k) ./ mu_s{3-c}(k) );   % should be positive ! (no need for max(0,...)) 
                else,
                % Eq. (4) (original specCSP):  
                    alpha_opt{c}(k) = max( 0 , ( mu_s{c}(k) -  mu_s{3-c}(k) ) / ( var_s{1}(k) + var_s{2}(k) ) );
                end
                % Eq. (5), with prior from Eq. (6)
                alpha_tmp{c}(k) = (I(k)*alpha_opt{c}(k)).^q * ((mu_s{1}(k) + mu_s{2}(k))/2).^p;
            end
        end
        % ... as the maximum for both classes
        alpha{j} = max( alpha_tmp{1} , alpha_tmp{2} );
        % and normalize alpha{j} so that it sums to unity
        alpha{j} = alpha{j} / sum(alpha{j});
    end
    alpha = [vertcat(alpha{:})'; zeros(S-length(alpha{1}),length(alpha))];
    model = struct('W',{Wcat},'P',{Pcat},'alpha',{alpha_tot},'freqs',{freqs},'bands',{bands},'lambda',{lambda});
    waitbar(step/steps,h);
end
Implementation hints:

Compare to the original CSP, it also implements a (Rayleigh) "ratio mode" for spectral optimization.

This code reuses the existing OV implementation for CSP (the ovpCBoxAlgorithmRegularizedCSPTrainer) for the linear algebra needed for the CSP problem resolution, but it also requires to epoch signals as 3d real and complex matrix (FFTs of the epoched time domain signals) and to store instantaneous cross spectrum matrix as 4d real arrays. In most of the implementation, the class dimension could have been added to avoid the double declaration of arrays for the two classes (but adding one more dimension to all the arrays).

There is a slight difference of implementation with the OV and bcilab algorithms though: in OV, the eigenvalues an vectors are computed for the two classes and the best (highest eigenvalues) are selected for the two classes, whereas in bcilab, they are computed for only one class, the best eigenvectors are associated to the current class and the lowest vectors are projectors for the other class. This is why the step 7 is slightly different (two argmax for OV, one argmax an one argmin for bcilab). It should be completely equivalent though ? I'll check this on Matlab.

This code should compile, run and generate the output spatiospectral filters config file, but it still needs to be validated functionnally (verify step by step that it behaves accordingly to the Matlab/bcilab model in both modes, ratio mode or original CSP mode). Convergence indicators of the iterative algorithm that are not implemented in the bcilab/matlab model should also be added to help the user to select the proper number of steps for the spatiospectral filters optimisation.

Since I won't have much time in the coming days for this validation (I will continue though), if anyone want to help, welcome ! :) Beyond functional validation, if there are better implementations for the data structures, please let me know !

Here is the code (with many unecessary file writings or variable storages for debugging purposes).
Jeff
Attachments
ovpCBoxAlgorithmSpecCSPTrainer.cpp
(34.75 KiB) Downloaded 242 times
ovpCBoxAlgorithmSpecCSPTrainer.h
(9.11 KiB) Downloaded 247 times

Jeff_B
Posts: 37
Joined: Wed Apr 11, 2012 1:31 pm
Location: Nice - Alpes Maritimes

Re: Spectrally weighted CSP OV algorithm box: help to valida

Post by Jeff_B »

I forgot:

the connectivity of the SpecCSPtrainer box is the same than the RegularizedCSPTrainer box: the stimulation signal is the 1st input, the stimulation based signals for class #1 and #2 are the inputs 2 an 3... and it generates a train_completed stimulation signal on the only output (stimulation signal).

Hence the connections & setting parameters are as such:
Image

Jeff_B
Posts: 37
Joined: Wed Apr 11, 2012 1:31 pm
Location: Nice - Alpes Maritimes

Re: Spectrally weighted CSP OV algorithm box: help to valida

Post by Jeff_B »

Hi OVibers,
The use of the boost multidimensional array library for all the 2d, 3d, 4d (real or complex) data containers
http://www.boost.org/doc/libs/1_60_0/li ... /user.html
... could be a better choice for the algorithm data structures ! (instead of my nested vectors). I am going to rewrite that.
Jeff

Jeff_B
Posts: 37
Joined: Wed Apr 11, 2012 1:31 pm
Location: Nice - Alpes Maritimes

Re: Spectrally weighted CSP OV algorithm box: help to valida

Post by Jeff_B »

Hi OVibers,
There is an error for the eq. 4 of the implementation, it was:

Code: Select all

alpha_opt{c}(k) = max( 0 , ( mu_s{c}(k) -  mu_s{3-c}(k) ) / ( var_s{1}(k) + var_s{2}(k) ) );
... and should be:

Code: Select all

alpha_opt{c}(k) = max( 0 , ( mu_s{c}(k) -  mu_s{3-c}(k) ) / sqrt( var_s{1}(k) + var_s{2}(k) ) );
I corrected the OV code and let the EEGLAB team know about this bug.

Given this, on the default EEG signals budlded with the OV motor imagery scenario, the first outputs of the algorithm seems relevant for the "ratio mode" (see below) but unexpected for the default mode.

the spectral weighting seems to peak in the beta band for the best (the #1 and #6 for 10 spatiospectral filters) spectiospectral projectors:
Image

the spatial weighting seems relevant to MI (right left hand)
Image

and, like in the original article, the convergence of the (max) spectral weighting is very fast (3 iterations should do it):
Image

... which is also visible on the convergence of the max eigenvalues
Image

but for the default mode, the spectral weighting is inviting us to look deeper :shock:
Image

I will thus focus on the original mode and begin the step by step comparison between OV and Matlab.
Jeff

Jeff_B
Posts: 37
Joined: Wed Apr 11, 2012 1:31 pm
Location: Nice - Alpes Maritimes

Re: Spectrally weighted CSP OV algorithm box: help to valida

Post by Jeff_B »

Hi there,

Here is my last version (version 0.6) of the SpecCSP algorithm.

- Data structures (2d, 3d and 4d real or complex matrix) migrated toward boost::multi_array

- Dual CSP mode added: CSP1 is based on the eigendecomposition eig(inv(R2)*R1) while CSP2 is eig(White*R1*White') where White is the whitenning matrix based on the eigendecomposition of the composite covariance R1+R2. CSP1 is the mode used in the 2 other OV CSP filter trainers (ovpCBoxAlgorithmCSPTrainer and ovpCBoxAlgorithmRegularizedCSPTrainer) while CSP2 is the default here.

- Two simple convergence indicators (l_pMaxLambda and l_pMaxAlpha) added to the output config. file


I realized a bit too late that boost::multi_array is incompatible with MSVC2010 in debug mode :oops:
https://svn.boost.org/trac/boost/ticket/4874
I'll ask Jussi whether it is redhibitory (for an OV inclusion) and if I reverse to the nested standard vectors to build the 2d, 3d and 4d matrix.

The BCILAB CSP implementation is based on the generalized eigendecomposition eig( R1 , R1+R2 ) but if the EIGEN library (used in OV) has a generalized eigensolver, it does not produces any eigenvectors :shock: (only the eigenvalues are accessible, the eigenvectors are in the EIGEN todo list).

With this CSP2 mode, the original SpecCSP algorithm seems to produce spectral weightings closer to our expectations for a right or left hand imagery (note: I am testing it on the standard OV motor imagery signals: ${Path_Samples}/signals/bci-motor-imagery.ov)

Image

Again with the CSP2 mode, the (Rayleigh) Ratio mode produces similar spatial weights with narrower spectral weightings:

Image

With the previous CSP1 mode (and still in Ratio mode), while the spectral weights are nearly the same (the two classes are inverted between CSP1 and CSP2) the spatial weights are very close but slightly differs for left hand movement (closer to the head top or closer to the electrode C1 with the CSP1 mode):

Image

And here is the last configuration (original spectral mode, CSP1 mode):

Image
Attachments
ovpCBoxAlgorithmSpecCSPTrainer.cpp
(34.24 KiB) Downloaded 246 times
ovpCBoxAlgorithmSpecCSPTrainer.h
(8.95 KiB) Downloaded 241 times

jtlindgren
Posts: 778
Joined: Tue Dec 04, 2012 3:53 pm
Location: INRIA Rennes, FRANCE

Re: Spectrally weighted CSP OV algorithm box: help to valida

Post by jtlindgren »

Hello Jeff, absolutely super work and analysis, and thanks for the code!

I had cranked a few of the OV motor imagery files through Matlab analysis earlier (they were recorded long before my time...), and I could notice a discriminating spectral characteristics between the conditions and laterality around 12hz, matching your results -- though I didn't make a spatial heatmap.

About boost, our preferred approach is to use Eigen where applicable, but its in the works to support VS2013 in the next release (1.2.0), so if the boost problem doesn't appear with that VS, its just fine. In that case people will just need to compile your plugin with a newer VS. And btw, for OV 1.2.0, the boost used and provided will be 1.55.

Would you like the plugin to be integrated in openvibe? BCILab license seems to be GPL, so there's no problem from that direction.


Cheers & thanks,
Jussi

Jeff_B
Posts: 37
Joined: Wed Apr 11, 2012 1:31 pm
Location: Nice - Alpes Maritimes

Re: Spectrally weighted CSP OV algorithm box: help to valida

Post by Jeff_B »

Thank you Jussi !
The generated spatial heatmaps are pretty much consistent between each others and with the spatial filters generated by the other OV spatial filter trainers. But I still want to validate it step by step (the ones numbered in the code) along with Matlab.
The boost library is a pain in the ass for me (using MSVC2010) ...as I am coding/debugging the companion box of this SpecCSP algorithm box, a box to apply the spatiospectral filters it estimates (Spatio Spectral filter box) to the signals. So I would like to get rid of boost for Eigen, but is it still possible to create multidimensional arrays (4d real matrix and 3d complex matrix) with Eigen without complicated pointers arithmetics ?
Jeff
Image

Post Reply