csp

[catha]

Hey,

I noticed that the CSP box in Openvibe does not give exactly the same coefficients as I get when I calculate them in Matlab. I used 2 different Matlab codes to calculate CSP coefficients, and they correlate with 1, therefore I assume the Matlab codes working proper. But the correlation of Openvibe coefficients and Matlab coefficients are around 0.94. Does anybody have an explanation why they are not exactly the same?

Thanks
Catha

[fabien.lotte]

Hi Catha,

There are actually different ways to compute the CSP, so I guess it depends which method was used in your matlab codes. For instance, some CSP implementations normalise the covariance matrices, some don't. Some implementations use generalized eigen value decomposition and some other use matrix inversion and (standard) eigen value decomposition. Some implementations use the ratio of the variance of one class to the variance of the other class, while some other use the ratio of the difference of the variance of the two classes to the sum of the variance of the two classes. I didn't designed the CSP box in OpenViBE but I guess these different ways of designing the CSP may explain the differences you observed. Do you know which approach was used in your matlab codes?

Best,
Fabien

[catha]

Dear Fabian

Thank you for your answer. We used different implementations.

We normalise the covariance matrices. We used both implementations based on generalized eigen value decomposition and matrix inversion + (standard) eigen value decomposition, but this should give in theory, and gives in practise the same result. Concerning the ratio, we checked both the eigen value decomposition from the difference of one class to the other, as the difference of one class to the sum of both classes.

All our Matlab results are correlated over 0.99 %. Because openvibe results correlate to these results "only" with 0.94 %, we would be interested to see which approach is used there. Anyone can give more details on this?

Thanks a lot,
Catha


In case it helps, we attach the implementations:

function [unmixing] = csp_new(dat1, dat2)

% CSP Common spatial pattern decomposition
%
% Use as
% [unmixing] = csp(dat1, dat2)
%
% This implements Ramoser, H., Gerking, M., and Pfurtscheller, G. "Optimal
% spatial filtering of single trial EEG during imagined hand movement."% IEEE Trans. Rehab. Eng 8 (2000), 446, 441.

% The initial version was coded by James Ethridge and William Weaver.
% See http://www.mathworks.com/matlabcentral/ ... l-patterns
% Some cleanups by Robert Oostenveld, 2012

R1 = dat1*dat1';
R1 = R1/trace(R1);
R2 = dat2*dat2';
R2 = R2/trace(R2);

% Ramoser equation (2)
Rsum = R1+R2;

% Find Eigenvalues and Eigenvectors of RC
% Sort eigenvalues in descending order
[EVecsum,EValsum] = eig(Rsum);
[EValsum,ind] = sort(diag(EValsum),'descend');
EVecsum = EVecsum(:,ind);

% Find Whitening Transformation Matrix - Ramoser Equation (3)
W = sqrt(pinv(diag(EValsum))) * EVecsum';

% Whiten Data Using Whiting Transform - Ramoser Equation (4)
S1 = W * R1 * W';
S2 = W * R2 * W';

% Ramoser equation (5)
% [U1,Psi1] = eig(S1);
% [U2,Psi2] = eig(S2);

%generalized eigenvectors/values
[B,D] = eig(S1,S2);

% Simultanous diagonalization
% Should be equivalent to [B,D]=eig(S1);

% verify algorithim
%disp('test1:Psi1+Psi2=I')
%Psi1+Psi2

% sort ascending by default
%[Psi1,ind] = sort(diag(Psi1)); U1 = U1(:,ind);
%[Psi2,ind] = sort(diag(Psi2)); U2 = U2(:,ind);
[D,ind]=sort(diag(D));
B=B(:,ind);


%keyboard

% Resulting Projection Matrix-these are the spatial filter coefficients
unmixing = B'*W;

=============================================================


function [V,d] = csp(ECM,arg2,arg3)
% CSP computes common spatial patterns
% this version supports multiple classes using a One-vs-Rest scheme
%
% [V] = csp(ECM)
% [V] = csp(X,Y)
% [V] = csp(...,Mode)
%
% ECM(k,:,:) is the extended covariance matrices (see COVM) for class k
% X,Y are matrices of the two classes (one channel per column)
% the number of columns must be the same for X and Y
% Mode = 'CSP0' uses common diagonalization
% = 'CSP3' solves generalized eigenvalue problem
% V each column represents one CSP component.
%
% REFERENCES:
% [1] Koles ZJ, Soong AC.
% EEG source localization: implementing the spatio-temporal decomposition approach.
% Electroencephalogr Clin Neurophysiol. 1998 Nov;107(5):343-52
% [2] Ramoser, H.; Muller-Gerking, J.; Pfurtscheller, G.;
% Optimal spatial filtering of single trial EEG during imagined hand movement
% Rehabilitation Engineering, IEEE Transactions on [see also IEEE Trans. on Neural Systems and Rehabilitation]
% Volume 8, Issue 4, Dec. 2000 Page(s):441 - 446
% [3] Dornhege, G.; Blankertz, B.; Curio, G.; Muller, K.-R.;
% Boosting bit rates in noninvasive EEG single-trial classifications by feature combination and multiclass paradigms
% Biomedical Engineering, IEEE Transactions on
% Volume 51, Issue 6, June 2004 Page(s):993 - 1002
% Digital Object Identifier 10.1109/TBME.2004.827088
% [4] Lemm, S.; Blankertz, B.; Curio, G.; Muller, K.-R.;
% Spatio-Spectral Filters for Improving the Classification of Single Trial EEG
% Biomedical Engineering, IEEE Transactions on
% Volume 52, Issue 9, Sept. 2005 Page(s):1541 - 1548
% Digital Object Identifier 10.1109/TBME.2005.851521

% $Id: csp.m 2205 2009-10-27 12:18:15Z schloegl $
% Copyright (C) 2007,2008,2009 by Alois Schloegl <a.schloegl@ieee.org>
% This is part of the BIOSIG-toolbox http://biosig.sf.net/

p = 2;

Mode =[];
sz = size(ECM);
if (length(sz)<3) && isnumeric(arg2),
X = ECM;
ECM = permute(cat(3,covm(X,'E'),covm(arg2,'E')),[3,1,2]);
end;
sz = size(ECM);
if (nargin>1) && ischar(arg2)
Mode = arg2;
end;
if (nargin>1) && ischar(arg2)
Mode = arg3;
end;
if isempty(Mode),
Mode = 'CSP3';
end;
%keyboard;




%COV = zeros(size(ECM)-[1,1,0]);
for k=1:sz(1),
[mu,sd,COV(:,:,k),xc,N,R2]=decovm(squeeze(ECM(k,:,:)));
end;


V = repmat(NaN,sz(2)-1,2*sz(1));
if 0,

elseif strcmpi(Mode,'CSP0');
% common diagonalization
[P,D] = eig(squeeze(sum(COV,3)));
P = diag(sqrt(1./diag(D)))*P';

% class specific components
V = repmat(NaN,sz(2)-1,2*sz(1));
d = V(1,:);
for k = 1:sz(1),
C = P * squeeze(COV(:,:,k)) * P';
[R,d1] = eig((C+C')/2);
[d1,ix] = sort(diag(d1));

V(:,2*k+[-1:0]) = P'*R(:,ix([1,end]));
d(1,2*k+[-1:0]) = d1(ix([1,end]))';
end;

elseif strcmpi(Mode,'CSP3');
%% do actual CSP calculation as generalized eigenvalues
%% R = permute(COV,[2,3,1]);
for k = 1:sz(1),
[W,D] = eig(COV(:,:,k),sum(COV,3));
V(:,2*k+[1-p:0]) = W(:,[1,end]);
end;

end;

[ddvlamin]

Probably this comes a bit late.

Anyway, I wrote that code a while ago. There are indeed some things that have been done differently and could be improved/changed. Each epoch of each class is first centred, then the covariance is computed and then this is trace normalized. However in OpenViBE I thought it was not the best idea to keep all epochs in memory and then do a batch covariance computation on all the data. Instead this should be done incrementally as the data comes in. Because of this I just averaged over all separate covariance matrices which is not really a valid way to do it. When both classes contain the same amount of samples this will be very similar, definitely for a big amount of samples (I think). Then just the standard eigenvalue decomposition is applied on inv(simga_2)*sigma_1. So no it probably does not exactly follow the implementation that you gave. But it should be easy to change. There are better/valid ways of doing the incremental computation of the covariance matrix

[nsprt]

hi;

please help me. I really don't know whether the input dat of csp function 'csp_new(dat1, dat2)' is [channels×samples] or [ samples ×channels].