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