function [comp] = ft_componentanalysis(cfg, data)

% FT_COMPONENTANALYSIS performs independent component analysis or other

% spatio-temporal decompositions of EEG or MEG data. This function computes

% the topography and timecourses of the components. The output of this

% function can be further analyzed with FT_TIMELOCKANALYSIS or

% FT_FREQANALYSIS.

%

% Use as

% [comp] = ft_componentanalysis(cfg, data)

% where cfg is a configuration structure and the input data is obtained from

% FT_PREPROCESSING or from FT_TIMELOCKANALYSIS.

%

% The configuration should contain

% cfg.method = 'runica', 'fastica', 'binica', 'pca', 'svd', 'jader',

% 'varimax', 'dss', 'cca', 'sobi', 'white' or 'csp'

% (default = 'runica')

% cfg.channel = cell-array with channel selection (default = 'all'),

% see FT_CHANNELSELECTION for details

% cfg.split = cell-array of channel types between which covariance

% is split, it can also be 'all' or 'no' (default = 'no')

% cfg.trials = 'all' or a selection given as a 1xN vector (default = 'all')

% cfg.numcomponent = 'all' or number (default = 'all')

% cfg.demean = 'no' or 'yes', whether to demean the input data (default = 'yes')

% cfg.updatesens = 'no' or 'yes' (default = 'yes')

% cfg.feedback = 'no', 'text', 'textbar', 'gui' (default = 'text')

%

% The runica method supports the following method-specific options. The

% values that these options can take can be found with HELP RUNICA.

% cfg.runica.extended

% cfg.runica.pca

% cfg.runica.sphering

% cfg.runica.weights

% cfg.runica.lrate

% cfg.runica.block

% cfg.runica.anneal

% cfg.runica.annealdeg

% cfg.runica.stop

% cfg.runica.maxsteps

% cfg.runica.bias

% cfg.runica.momentum

% cfg.runica.specgram

% cfg.runica.posact

% cfg.runica.verbose

% cfg.runica.logfile

% cfg.runica.interput

%

% The fastica method supports the following method-specific options. The

% values that these options can take can be found with HELP FASTICA.

% cfg.fastica.approach

% cfg.fastica.numOfIC

% cfg.fastica.g

% cfg.fastica.finetune

% cfg.fastica.a1

% cfg.fastica.a2

% cfg.fastica.mu

% cfg.fastica.stabilization

% cfg.fastica.epsilon

% cfg.fastica.maxNumIterations

% cfg.fastica.maxFinetune

% cfg.fastica.sampleSize

% cfg.fastica.initGuess

% cfg.fastica.verbose

% cfg.fastica.displayMode

% cfg.fastica.displayInterval

% cfg.fastica.firstEig

% cfg.fastica.lastEig

% cfg.fastica.interactivePCA

% cfg.fastica.pcaE

% cfg.fastica.pcaD

% cfg.fastica.whiteSig

% cfg.fastica.whiteMat

% cfg.fastica.dewhiteMat

% cfg.fastica.only

%

% The binica method supports the following method-specific options. The

% values that these options can take can be found with HELP BINICA.

% cfg.binica.extended

% cfg.binica.pca

% cfg.binica.sphering

% cfg.binica.lrate

% cfg.binica.blocksize

% cfg.binica.maxsteps

% cfg.binica.stop

% cfg.binica.weightsin

% cfg.binica.verbose

% cfg.binica.filenum

% cfg.binica.posact

% cfg.binica.annealstep

% cfg.binica.annealdeg

% cfg.binica.bias

% cfg.binica.momentum

%

% The dss method requires the following method-specific option and supports

% a whole lot of other options. The values that these options can take can

% be found with HELP DSS_CREATE_STATE.

% cfg.dss.denf.function

% cfg.dss.denf.params

%

% The sobi method supports the following method-specific options. The

% values that these options can take can be found with HELP SOBI.

% cfg.sobi.n_sources

% cfg.sobi.p_correlations

%

% The csp method implements the common-spatial patterns method. For CSP, the

% following specific options can be defined:

% cfg.csp.classlabels = vector that assigns a trial to class 1 or 2.

% cfg.csp.numfilters = the number of spatial filters to use (default: 6).

%

% The icasso method implements icasso. It runs fastica a specified number of

% times, and provides information about the stability of the components found

% The following specific options can be defined, see ICASSOEST:

% cfg.icasso.mode

% cfg.icasso.Niter

%

% Instead of specifying a component analysis method, you can also specify

% a previously computed unmixing matrix, which will be used to estimate the

% component timecourses in this data. This requires

% cfg.unmixing = NxN unmixing matrix

% cfg.topolabel = Nx1 cell-array with the channel labels

%

% You may specify a particular seed for random numbers called by

% rand/randn/randi, or the random state used by a previous call to this

% function to replicate results. For example:

% cfg.randomseed = integer seed value of user's choice

% cfg.randomseed = comp.cfg.callinfo.randomseed (from previous call)

%

% To facilitate data-handling and distributed computing you can use

% cfg.inputfile = ...

% cfg.outputfile = ...

% If you specify one of these (or both) the input data will be read from a *.mat

% file on disk and/or the output data will be written to a *.mat file. These mat

% files should contain only a single variable, corresponding with the

% input/output structure.

%

% See also FT_TOPOPLOTIC, FT_REJECTCOMPONENT, FASTICA, RUNICA, BINICA, SVD,

% JADER, VARIMAX, DSS, CCA, SOBI, ICASSO

% Copyright (C) 2003-2012, Robert Oostenveld

%

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org

% for the documentation and details.

%

% FieldTrip is free software: you can redistribute it and/or modify

% it under the terms of the GNU General Public License as published by

% the Free Software Foundation, either version 3 of the License, or

% (at your option) any later version.

%

% FieldTrip is distributed in the hope that it will be useful,

% but WITHOUT ANY WARRANTY; without even the implied warranty of

% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

% GNU General Public License for more details.

%

% You should have received a copy of the GNU General Public License

% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.

%

% $Id$

% undocumented cfg options:

% cfg.cellmode = string, 'no' or 'yes', allows to run in cell-mode, i.e.

% no concatenation across trials is needed. This is based on experimental

% code and only supported for 'dss', 'fastica' and 'bsscca' as methods.

% these are used by the ft_preamble/ft_postamble function and scripts

ft_revision = '$Id$';

ft_nargin = nargin;

ft_nargout = nargout;

% do the general setup of the function

ft_defaults

ft_preamble init

ft_preamble debug

ft_preamble loadvar data

ft_preamble provenance data

ft_preamble randomseed

% the ft_abort variable is set to true or false in ft_preamble_init

if ft_abort

return

end

% check if the input data is valid for this function

istimelock = ft_datatype(data, 'timelock');

data = ft_checkdata(data, 'datatype', 'raw', 'feedback', 'yes');

% check if the input cfg is valid for this function

cfg = ft_checkconfig(cfg, 'forbidden', {'channels', 'trial'}); % prevent accidental typos, see issue 1729

cfg = ft_checkconfig(cfg, 'forbidden', {'detrend'});

cfg = ft_checkconfig(cfg, 'renamed', {'blc', 'demean'});

cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'predetermined mixing matrix', 'predetermined unmixing matrix'});

cfg = ft_checkconfig(cfg, 'deprecated', {'topo'});

% set the defaults

cfg.method = ft_getopt(cfg, 'method', 'runica');

cfg.demean = ft_getopt(cfg, 'demean', 'yes');

cfg.trials = ft_getopt(cfg, 'trials', 'all', 1);

cfg.channel = ft_getopt(cfg, 'channel', 'all');

cfg.split = ft_getopt(cfg, 'split', 'no');

cfg.numcomponent = ft_getopt(cfg, 'numcomponent', 'all');

cfg.normalisesphere = ft_getopt(cfg, 'normalisesphere', 'yes');

cfg.cellmode = ft_getopt(cfg, 'cellmode', 'no');

cfg.doscale = ft_getopt(cfg, 'doscale', 'yes');

cfg.updatesens = ft_getopt(cfg, 'updatesens', 'yes');

cfg.feedback = ft_getopt(cfg, 'feedback', 'text');

% select channels, has to be done prior to handling of previous (un)mixing matrix

cfg.channel = ft_channelselection(cfg.channel, data.label);

if istrue(cfg.cellmode)

ft_hastoolbox('cellfunction', 1);

end

if isfield(cfg, 'topo') && isfield(cfg, 'topolabel')

ft_warning(['Specifying cfg.topo (= mixing matrix) to determine component '...

'timecourses in specified data is deprecated; please specify an '...

'unmixing matrix instead with cfg.unmixing. '...

'Using cfg.unmixing=pinv(cfg.topo) for now to reproduce old behavior.']);

cfg.unmixing = pinv(cfg.topo);

cfg = rmfield(cfg, 'topo');

end

if isfield(cfg, 'unmixing') && isfield(cfg, 'topolabel')

% use the previously determined unmixing matrix on this dataset

% test whether all required channels are present in the data

[datsel, toposel] = match_str(cfg.channel, cfg.topolabel);

if length(toposel)~=length(cfg.topolabel)

ft_error('not all channels that are required for the unmixing are present in the data');

end

% ensure that all data channels not used in the unmixing should be removed from the channel selection

tmpchan = match_str(cfg.channel, cfg.topolabel);

cfg.channel = cfg.channel(tmpchan);

% update some settings where there is no further choice to be made by the user

cfg.numcomponent = 'all';

cfg.method = 'predetermined unmixing matrix';

end

% add the options for the specified methods to the configuration, only if needed

switch cfg.method

case 'icasso'

cfg.icasso = ft_getopt(cfg, 'icasso', []);

cfg.icasso.mode = ft_getopt(cfg.icasso, 'mode', 'both');

cfg.icasso.Niter = ft_getopt(cfg.icasso, 'Niter', 15);

cfg.icasso.method = ft_getopt(cfg.icasso, 'method', 'fastica');

cfg.fastica = ft_getopt(cfg, 'fastica', []);

case 'fastica'

% additional options, see FASTICA for details

cfg.fastica = ft_getopt(cfg, 'fastica', []);

case 'runica'

% additional options, see RUNICA for details

cfg.runica = ft_getopt(cfg, 'runica', []);

cfg.runica.lrate = ft_getopt(cfg.runica, 'lrate', 0.001);

case 'binica'

% additional options, see BINICA for details

cfg.binica = ft_getopt(cfg, 'binica', []);

cfg.binica.lrate = ft_getopt(cfg.binica, 'lrate', 0.001);

case 'dss'

% additional options, see DSS for details

cfg.dss = ft_getopt(cfg, 'dss', []);

cfg.dss.denf = ft_getopt(cfg.dss, 'denf', []);

cfg.dss.denf.function = ft_getopt(cfg.dss.denf, 'function', 'denoise_fica_tanh');

cfg.dss.denf.params = ft_getopt(cfg.dss.denf, 'params', []);

cfg.dss.preprocf = ft_getopt(cfg.dss, 'preprocf', []);

cfg.dss.preprocf.function = ft_getopt(cfg.dss.preprocf, 'function', 'pre_sphere');

cfg.dss.preprocf.params = ft_getopt(cfg.dss.preprocf, 'params', []);

case 'csp'

% additional options, see CSP for details

cfg.csp = ft_getopt(cfg, 'csp', []);

cfg.csp.numfilters = ft_getopt(cfg.csp, 'numfilters', 6);

cfg.csp.classlabels = ft_getopt(cfg.csp, 'classlabels');

case 'bsscca'

% additional options, see BSSCCA for details

cfg.bsscca = ft_getopt(cfg, 'bsscca', []);

cfg.bsscca.refdelay = ft_getopt(cfg.bsscca, 'refdelay', 1);

cfg.bsscca.chandelay = ft_getopt(cfg.bsscca, 'chandelay', 0);

if strcmp(cfg.cellmode, 'no')

ft_error('cfg.mehod = ''bsscca'' requires cfg.cellmode = ''yes''');

end

otherwise

% do nothing

end

% select trials of interest

tmpcfg = keepfields(cfg, {'trials', 'channel', 'tolerance', 'showcallinfo', 'trackcallinfo', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo'});

data = ft_selectdata(tmpcfg, data);

% restore the provenance information

[cfg, data] = rollback_provenance(cfg, data);

% deal with different chantypes if requested

if isequal(cfg.split, 'no')

chantype = {};

elseif isequal(cfg.split, 'all')

chantype = unique(ft_chantype(data.label));

else

chantype = cfg.split;

end

if numel(chantype)>0

% recurse per specified chantype

tmpdata = cell(1, numel(chantype));

for k = 1:numel(chantype)

tmpcfg = cfg;

tmpcfg.channel = data.label(ft_chantype(data.label, lower(chantype{k})));

tmpcfg.split = 'no';

tmpcfg.chantype = lower(chantype{k}); % makes the output labels unique, to allow appending later on

tmpdata{1,k} = ft_componentanalysis(tmpcfg, data);

end

comp = ft_appenddata([], tmpdata{:});

return;

else

%

end

Ntrials = length(data.trial);

Nchans = length(data.label);

if Nchans==0

ft_error('no channels were selected');

end

% default is to compute just as many components as there are channels in the data

if strcmp(cfg.numcomponent, 'all')

defaultNumCompsUsed = true(1);

cfg.numcomponent = length(data.label);

else

defaultNumCompsUsed = false(1);

end

% determine the size of each trial, they can be variable length

Nsamples = zeros(1,Ntrials);

for trial=1:Ntrials

Nsamples(trial) = size(data.trial{trial},2);

end

if strcmp(cfg.demean, 'yes')

% optionally perform baseline correction on each trial

ft_info('baseline correcting data \n');

for trial=1:Ntrials

data.trial{trial} = ft_preproc_baselinecorrect(data.trial{trial});

end

end

if strcmp(cfg.doscale, 'yes')

% determine the scaling of the data, scale it to approximately unity

% this will improve the performance of some methods, esp. fastica

trlidx = 1;

tmp = data.trial{trlidx};

while all(isnan(tmp(:))) % if all data in this trial is NaN

trlidx = trlidx + 1; % try next trial

tmp = data.trial{trlidx}; % overwrite tmp with next trial

end

tmp(~isfinite(tmp)) = 0; % ensure that the scaling is a finite value

scale = norm((tmp*tmp')./size(tmp,2)); clear tmp;

scale = sqrt(scale);

if scale ~= 0

ft_info('scaling data with 1 over %f\n', scale);

for trial=1:Ntrials

data.trial{trial} = data.trial{trial} ./ scale;

end

else

ft_info('no scaling applied, since factor is 0\n');

end

else

ft_info('no scaling applied to the data\n');

end

if strcmp(cfg.method, 'sobi')

% concatenate all the data into a 3D matrix respectively 2D (sobi)

ft_info('concatenating data');

Nsamples = Nsamples(1);

dat = zeros(Ntrials, Nchans, Nsamples);

% all trials should have an equal number of samples

% and it is assumed that the time axes of all trials are aligned

for trial=1:Ntrials

ft_info('.');

dat(trial,:,:) = data.trial{trial};

end

ft_info('\n');

ft_info('concatenated data matrix size %dx%dx%d\n', size(dat,1), size(dat,2), size(dat,3));

if Ntrials == 1

dummy = 0;

[dat, dummy] = shiftdim(dat);

else

dat = shiftdim(dat,1);

end

elseif strcmp(cfg.method, 'csp')

% concatenate the trials into two data matrices, one for each class

sel1 = find(cfg.csp.classlabels==1);

sel2 = find(cfg.csp.classlabels==2);

if min(length(sel1), length(sel2)) == 0

ft_error('CSP requires class labels!');

end

if length(sel1)+length(sel2)~=length(cfg.csp.classlabels)

ft_warning('not all trials belong to class 1 or 2');

end

dat1 = cat(2, data.trial{sel1});

dat2 = cat(2, data.trial{sel2});

ft_info('concatenated data matrix size for class 1 is %dx%d\n', size(dat1,1), size(dat1,2));

ft_info('concatenated data matrix size for class 2 is %dx%d\n', size(dat2,1), size(dat2,2));

elseif ~strcmp(cfg.method, 'predetermined unmixing matrix') && strcmp(cfg.cellmode, 'no')

% concatenate all the data into a 2D matrix unless we already have an

% unmixing matrix or unless the user request it otherwise

ft_info('concatenating data');

dat = zeros(Nchans, sum(Nsamples));

ft_progress('init', cfg.feedback, 'concatenating trials...');

for trial=1:Ntrials

ft_progress(trial/Ntrials, 'Concatenating trial %d from %d', trial, Ntrials);

begsample = sum(Nsamples(1:(trial-1))) + 1;

endsample = sum(Nsamples(1:trial));

dat(:,begsample:endsample) = data.trial{trial};

end

ft_progress('close')

ft_info('concatenated data matrix size %dx%d\n', size(dat,1), size(dat,2));

hasdatanans = any(~isfinite(dat(:)));

if hasdatanans && strcmp(cfg.method, 'dss')

ft_error('DSS does not work with nans or inf in the data');

elseif hasdatanans

ft_info('data contains nan or inf, only using the samples without nan or inf\n');

finitevals = sum(~isfinite(dat))==0;

if ~any(finitevals)

ft_error('no samples remaining');

else

dat = dat(:,finitevals);

end

end

else

ft_info('not concatenating data\n');

dat = data.trial;

% FIXME cellmode processing is not nan-transparent yet

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% perform the component analysis

ft_info('starting decomposition using %s\n', cfg.method);

switch cfg.method

case 'icasso'

% check whether the required low-level toolboxes are installed

ft_hastoolbox('icasso', 1);

if strcmp(cfg.icasso.method, 'fastica')

ft_hastoolbox('fastica', 1);

cfg.fastica.numOfIC = cfg.numcomponent;

optarg = ft_cfg2keyval(cfg.(cfg.icasso.method));

sR = icassoEst(cfg.icasso.mode, dat, cfg.icasso.Niter, optarg{:});

elseif strcmp(cfg.icasso.method, 'dss')

% recurse into ft_componentanalysis and do some post processing

tmpcfg = rmfield(cfg, 'icasso');

tmpcfg.method = cfg.icasso.method;

tmpdata = data;

% initialize the variables to hold the output

sR.W = cell(cfg.icasso.Niter, 1);

sR.A = cell(cfg.icasso.Niter, 1);

sR.index = zeros(0,2);

for k = 1:cfg.icasso.Niter

tmp = ft_componentanalysis(tmpcfg, tmpdata);

sR.W{k} = tmp.unmixing;

sR.A{k} = tmp.topo;

sR.index = cat(1, sR.index, [k*ones(size(tmp.topo,2),1) (1:size(tmp.topo,2))']);

sR.whiteningMatrix = tmp.cfg.dss.V;

sR.dewhiteningMatrix = tmp.cfg.dss.dV;

end

sR.signal = dat;

sR.mode = cfg.icasso.mode;

sR.rdim = size(tmp.topo,2);

else

ft_error('only ''fastica'' or ''dss'' is supported as method for icasso');

end

% do the rest of the icasso related processing

sR = icassoCluster(sR, 'strategy', 'AL', 'simfcn', 'abscorr', 's2d', 'sim2dis', 'L',cfg.numcomponent);

sR = icassoProjection(sR, 'cca', 's2d', 'sqrtsim2dis', 'epochs', 75);

[Iq, mixing, unmixing, dum, index2centrotypes] = icassoResult(sR,cfg.numcomponent);

% this step is done, because in icassoResult mixing is determined to be

% pinv(unmixing), which yields strange results. Better take it from the

% individual iterations. NOTE: as a consequence unmixing*mixing is not

% necessarily identity anymore !!!

for k = 1:size(mixing,2)

ix = sR.index(index2centrotypes(k),:);

mixing(:,k) = sR.A{ix(1)}(:,ix(2));

end

%[Iq, mixing, unmixing, dat] = icassoShow(sR, 'estimate', 'off', 'L', cfg.numcomponent);

% sort the output according to Iq

[srt, ix] = sort(-Iq); % account for NaNs

mixing = mixing(:, ix);

unmixing = unmixing(ix, :);

cfg.icasso.Iq = Iq(ix);

cfg.icasso.sR = rmfield(sR, 'signal'); % keep the rest of the information

case 'fastica'

% check whether the required low-level toolboxes are installed

ft_hastoolbox('fastica', 1); % see http://www.cis.hut.fi/projects/ica/fastica

if ~defaultNumCompsUsed &&...

(~isfield(cfg, 'fastica') || ~isfield(cfg.fastica, 'numOfIC'))

% user has specified cfg.numcomponent and not specified

% cfg.fastica.numOfIC, so copy cfg.numcomponent over

cfg.fastica.numOfIC = cfg.numcomponent;

elseif ~defaultNumCompsUsed &&...

isfield(cfg, 'fastica') && isfield(cfg.fastica, 'numOfIC')

% user specified both cfg.numcomponent and cfg.fastica.numOfIC,

% unsure which one to use

ft_error('you can specify either cfg.fastica.numOfIC or cfg.numcomponent (they will have the same effect), but not both');

end

try

% construct key-value pairs for the optional arguments

optarg = ft_cfg2keyval(cfg.fastica);

[mixing, unmixing] = fastica(dat, optarg{:});

catch

% the "catch me" syntax is broken on MATLAB74, this fixes it

me = lasterror;

% give a hopefully instructive error message

ft_info(['If you get an out-of-memory in fastica here, and you use fastica 2.5, change fastica.m, line 482: \n' ...

'from\n' ...

' if ~isempty(W) %% ORIGINAL VERSION\n' ...

'to\n' ...

' if ~isempty(W) && nargout ~= 2 %% if nargout == 2, we return [A, W], and NOT ICASIG\n']);

% forward original error

rethrow(me);

end

case 'runica'

% check whether the required low-level toolboxes are installed

% see http://www.sccn.ucsd.edu/eeglab

ft_hastoolbox('eeglab', 1);

if ~defaultNumCompsUsed &&...

(~isfield(cfg, 'runica') || ~isfield(cfg.runica, 'pca'))

% user has specified cfg.numcomponent and not specified

% cfg.runica.pca, so copy cfg.numcomponent over

cfg.runica.pca = cfg.numcomponent;

elseif ~defaultNumCompsUsed &&...

isfield(cfg, 'runica') && isfield(cfg.runica, 'pca')

% user specified both cfg.numcomponent and cfg.runica.pca,

% unsure which one to use

ft_error('you can specify either cfg.runica.pca or cfg.numcomponent (they will have the same effect), but not both');

end

% construct key-value pairs for the optional arguments

optarg = [ft_cfg2keyval(cfg.runica) {'reset_randomseed' 0}]; % let FieldTrip deal with the random seed handling

[weights, sphere] = runica(dat, optarg{:});

% scale the sphering matrix to unit norm

if strcmp(cfg.normalisesphere, 'yes')

sphere = sphere./norm(sphere);

end

unmixing = weights*sphere;

mixing = [];

case 'binica'

% check whether the required low-level toolboxes are installed

% see http://www.sccn.ucsd.edu/eeglab

ft_hastoolbox('eeglab', 1);

if ~defaultNumCompsUsed &&...

(~isfield(cfg, 'binica') || ~isfield(cfg.binica, 'pca'))

% user has specified cfg.numcomponent and not specified

% cfg.binica.pca, so copy cfg.numcomponent over

cfg.binica.pca = cfg.numcomponent;

elseif ~defaultNumCompsUsed &&...

isfield(cfg, 'binica') && isfield(cfg.binica, 'pca')

% user specified both cfg.numcomponent and cfg.binica.pca,

% unsure which one to use

ft_error('you can specify either cfg.binica.pca or cfg.numcomponent (they will have the same effect), but not both');

end

% construct key-value pairs for the optional arguments

optarg = ft_cfg2keyval(cfg.binica);

[weights, sphere] = binica(dat, optarg{:});

% scale the sphering matrix to unit norm

if strcmp(cfg.normalisesphere, 'yes')

sphere = sphere./norm(sphere);

end

unmixing = weights*sphere;

mixing = [];

case 'jader'

% check whether the required low-level toolboxes are installed

% see http://www.sccn.ucsd.edu/eeglab

ft_hastoolbox('eeglab', 1);

unmixing = jader(dat, cfg.numcomponent);

mixing = [];

case 'varimax'

% check whether the required low-level toolboxes are installed

% see http://www.sccn.ucsd.edu/eeglab

ft_hastoolbox('eeglab', 1);

unmixing = varimax(dat);

mixing = [];

case 'cca'

% check whether the required low-level toolboxes are installed

% see http://www.sccn.ucsd.edu/eeglab

ft_hastoolbox('cca', 1);

[y, w] = ccabss(dat);

unmixing = w';

mixing = [];

case 'pca'

% compute data cross-covariance matrix

if iscell(dat)

C = zeros(size(dat{1},1));

nC = 0;

for k = 1:numel(dat)

C = C + (dat{k}*dat{k}');

nC = nC + size(dat{k},2);

end

C = C./(nC-1);

else

C = (dat*dat')./(size(dat,2)-1);

end

% eigenvalue decomposition (EVD)

[E,D] = eig(C);

% sort eigenvectors in descending order of eigenvalues

d = cat(2,(1:1:Nchans)',diag(D));

d = sortrows(d, -2);

% return the desired number of principal components

unmixing = E(:,d(1:cfg.numcomponent,1))';

mixing = [];

clear C D E d

case 'kpca'

% linear kernel (same as normal covariance)

%kern = @(X,y) (sum(bsxfun(@times, X, y),2));

% polynomial kernel degree 2

%kern = @(X,y) (sum(bsxfun(@times, X, y),2).^2);

% RBF kernel

kern = @(X,y) (exp(-0.5* sqrt(sum(bsxfun(@minus, X, y).^2, 2))));

% compute kernel matrix

C = zeros(Nchans,Nchans);

ft_progress('init', cfg.feedback, 'computing kernel matrix...');

for k = 1:Nchans

ft_progress(k/Nchans, 'computing kernel matrix %d from %d', k, Nchans);

C(k,:) = kern(dat, dat(k,:));

end

ft_progress('close');

% eigenvalue decomposition (EVD)

[E,D] = eig(C);

% sort eigenvectors in descending order of eigenvalues

d = cat(2,(1:1:Nchans)',diag(D));

d = sortrows(d, -2);

% return the desired number of principal components

unmixing = E(:,d(1:cfg.numcomponent,1))';

mixing = [];

clear C D E d

case 'svd'

% it is more memory efficient to use the (non-scaled) covariance

if cfg.numcomponent<Nchans

% compute only the first components

[u, s, v] = svds(dat*dat', cfg.numcomponent);

else

% compute all components

[u, s, v] = svd(dat*dat', 0);

end

clear s v % not needed

unmixing = u';

mixing = [];

case 'dss'

% check whether the required low-level toolboxes are installed

% see http://www.cis.hut.fi/projects/dss

ft_hastoolbox('dss', 1);

params = removefields(struct(cfg.dss), {'V' 'dV' 'W' 'indx'});

params.denf.h = str2func(cfg.dss.denf.function);

params.preprocf.h = str2func(cfg.dss.preprocf.function);

if isfield(cfg.dss, 'wdim') && ~isempty(cfg.dss.wdim)

params.wdim = cfg.dss.wdim;

end

if ~ischar(cfg.numcomponent)

params.sdim = cfg.numcomponent;

if isfield(params, 'wdim')

params.sdim = min(params.sdim, params.wdim);

end

end

if isfield(params.denf, 'params') && isfield(params.denf.params, 'artifact')

% this may require the sampleinfo in the params structure, to keep the sampling bookkeeping correct

params.denf.params.sampleinfo = data.sampleinfo;

end

% create the state

state = dss_create_state(dat, params);

if isfield(cfg.dss, 'V') && ~isempty(cfg.dss.V)

state.V = cfg.dss.V;

state.Y = cfg.dss.V*dat;

end

if isfield(cfg.dss, 'dV') && ~isempty(cfg.dss.dV)

state.dV = cfg.dss.dV;

end

if isfield(cfg.dss, 'W') && ~isempty(cfg.dss.W)

state.W = cfg.dss.W;

end

if isfield(cfg.dss, 'indx') && ~isempty(cfg.dss.indx)

state.indx = cfg.dss.indx; %may be needed for dss_core_mim

end

% increase the amount of information that is displayed on screen

% state.verbose = 3;

% start the decomposition

state = denss(state); % this is for the DSS toolbox version 1.0

mixing = state.A;

unmixing = state.B;

% remember the updated configuration details

cfg.dss.denf = state.denf;

cfg.dss.orthof = state.orthof;

cfg.dss.preprocf = state.preprocf;

cfg.dss.stopf = state.stopf;

cfg.dss.W = state.W;

cfg.dss.V = state.V;

cfg.dss.dV = state.dV;

if isfield(state, 'D'), cfg.dss.D = state.D(1:min([state.sdim size(state.dV)])); end

cfg.numcomponent = min([state.sdim size(state.dV)]);

case 'sobi'

% check whether the required low-level toolboxes are installed

% see http://www.sccn.ucsd.edu/eeglab

ft_hastoolbox('eeglab', 1);

% check for additional options, see SOBI for details

if ~isfield(cfg, 'sobi')

mixing = sobi(dat, cfg.numcomponent);

elseif isfield(cfg.sobi, 'n_sources') && isfield(cfg.sobi, 'p_correlations')

mixing = sobi(dat, cfg.sobi.n_sources, cfg.sobi.p_correlations);

elseif isfield(cfg.sobi, 'n_sources')

mixing = sobi(dat,cfg.sobi.n_sources);

else

ft_error('unknown options for SOBI component analysis');

end

unmixing = [];

case 'predetermined unmixing matrix'

% check which labels from the cfg are identical to those of the data

% this gives us the rows of cfg.topo (the channels) and of

% data.trial (also channels) that we are going to use later

[datsel, chansel] = match_str(data.label, cfg.topolabel);

% ensure 1:1 corresponcence between cfg.topolabel & data.label

% otherwise we cannot compute the components (if source channels are

% missing) or will have a problem when projecting it back (because we

% dont have a marker to say that there are channels in data.label

% which we did not use and thus can't recover from source-space)

if length(cfg.topolabel)<length(chansel)

ft_error('cfg.topolabels do not uniquely correspond to data.label, please check')

end

if length(data.label)<length(datsel)

ft_error('cfg.topolabels do not uniquely correspond to data.label, please check')

end

% reorder the mixing matrix so that the channel order matches the order in the data

cfg.unmixing = cfg.unmixing(:,chansel);

cfg.topolabel = cfg.topolabel(chansel);

unmixing = cfg.unmixing;

mixing = [];

case 'white'

% compute the covariance matrix and an unmixing matrix that makes the data white

c = dat*dat';

c = c./(size(dat,2)-1);

[u, s] = svd(c);

% split the singular values into half

for i=1:size(s)

if (s(i,i)/s(1,1))>(100*eps)

s(i,i) = 1./sqrt(s(i,i));

else

s(i,i) = 0;

end

end

unmixing = s * u';

mixing = [];

case 'csp'

C1 = cov(dat1');

C2 = cov(dat2');

unmixing = csp(C1, C2, cfg.csp.numfilters);

mixing = []; % will be computed below

case 'bsscca'

% this method relies on time shifting of the original data, in much the

% same way as ft_denoise_tsr. as such it is more natural to represent

% the data in the cell-array, because the trial-boundaries are clear.

% if represented in a concatenated array one has to keep track of the

% trial boundaries

optarg = ft_cfg2keyval(cfg.bsscca);

optarg = cat(2,optarg, {'time', data.time});

[unmixing, mixing, rho, compdata, time] = bsscca(dat, optarg{:});

data.trial = mixing*compdata;

data.time = time;

data = removefields(data, 'sampleinfo');

if size(mixing,1)>numel(data.label)

for m = 1:(size(mixing,1)-numel(data.label))

data.label{end+1} = sprintf('refchan%03d',m);

end

end

% remember the canonical correlations

cfg.bsscca.rho = rho;

case 'parafac'

ft_error('parafac is not supported anymore in ft_componentanalysis');

otherwise

ft_error('unknown method for component analysis');

end % switch method

% make sure we have both mixing and unmixing matrices

% if not, compute (pseudo-)inverse to go from one to the other

if isempty(unmixing) && ~isempty(mixing)

if (size(mixing,1)==size(mixing,2))

unmixing = inv(mixing);

else

unmixing = pinv(mixing);

end

elseif isempty(mixing) && ~isempty(unmixing)

if (size(unmixing,1)==size(unmixing,2)) && rank(unmixing)==size(unmixing,1)

mixing = inv(unmixing);

else

mixing = pinv(unmixing);

end

elseif isempty(mixing) && isempty(unmixing)

% this sanity check is needed to catch convergence problems in fastica

% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1519

ft_error('the component unmixing failed');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% collect the results and construct data structure

comp = keepfields(data, {'time' 'fsample'});

% make sure we don't return more components than were requested

% (some methods respect the maxcomponent parameters, others just always

% return a fixed (i.e., numchans) number of components)

if size(unmixing,1) > cfg.numcomponent

unmixing(cfg.numcomponent+1:end,:) = [];

end

if size(mixing,2) > cfg.numcomponent

mixing(:,cfg.numcomponent+1:end) = [];

end

% compute the activations in each trial

if strcmp(cfg.doscale, 'yes')

for trial=1:Ntrials

comp.trial{trial} = scale * unmixing * data.trial{trial};

end

else

for trial=1:Ntrials

comp.trial{trial} = unmixing * data.trial{trial};

end

end

% store mixing/unmixing matrices in structure

comp.topo = mixing;

comp.unmixing = unmixing;

% get the labels

if strcmp(cfg.method, 'predetermined unmixing matrix')

prefix = 'component';

else

prefix = cfg.method;

end

st = dbstack;

if numel(st)>1 && isequal(st(2).name, 'ft_componentanalysis')

% this is a recursive call, as per the cfg.split option, add something

% extra to the prefix

chantype = ft_getopt(cfg, 'chantype', '');

prefix = [prefix chantype];

end

for k = 1:size(comp.topo,2)

comp.label{k,1} = sprintf('%s%03d', prefix, k);

end

comp.topolabel = data.label(:);

sensfield = cell(0,1);

if isfield(data, 'grad')

sensfield{end+1} = 'grad';

end

if isfield(data, 'elec')

sensfield{end+1} = 'elec';

end

if isfield(data, 'opto')

sensfield{end+1} = 'opto';

end

% apply the linear projection also to the sensor description

if ~isempty(sensfield)

if strcmp(cfg.updatesens, 'yes')

% construct a montage and apply it to the sensor description

montage = [];

montage.labelold = data.label;

montage.labelnew = comp.label;

montage.tra = unmixing;

for m = 1:numel(sensfield)

ft_info('also applying the unmixing matrix to the %s structure\n', sensfield{m});

comp.(sensfield{m}) = ft_apply_montage(data.(sensfield{m}), montage, 'balancename', 'comp', 'keepunused', 'yes');

% The output sensor array cannot simply be interpreted as the input

% sensor array, hence the type should be removed to allow autodetection

% See also http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1806

if isfield(comp.(sensfield{m}), 'type')

comp.(sensfield{m}) = rmfield(comp.(sensfield{m}), 'type');

end

end

else

for m = 1:numel(sensfield)

ft_info('not applying the unmixing matrix to the %s structure\n', sensfield{m});

% simply copy it over

comp.(sensfield{m}) = data.(sensfield{m});

end

end

end % if sensfield

% copy the sampleinfo into the output

if isfield(data, 'sampleinfo')

comp.sampleinfo = data.sampleinfo;

end

% copy the trialinfo into the output

if isfield(data, 'trialinfo')

comp.trialinfo = data.trialinfo;

end

% convert back to input type if necessary

if istimelock

comp = ft_checkdata(comp, 'datatype', 'timelock+comp');

end

% do the general cleanup and bookkeeping at the end of the function

ft_postamble debug

ft_postamble randomseed

ft_postamble previous data

ft_postamble provenance comp

ft_postamble history comp

ft_postamble savevar comp