Skip to content
Permalink
release
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
function [source] = ft_sourcedescriptives(cfg, source)
% FT_SOURCEDESCRIPTIVES computes descriptive parameters of the source
% analysis results.
%
% Use as
% [source] = ft_sourcedescriptives(cfg, source)
%
% where cfg is a structure with the configuration details and source is the
% result from a beamformer source estimation. The configuration can contain
% cfg.cohmethod = 'regular', 'lambda1', 'canonical'
% cfg.powmethod = 'regular', 'lambda1', 'trace', 'none'
% cfg.supmethod = 'chan_dip', 'chan', 'dip', 'none' (default)
% cfg.projectmom = 'yes' or 'no' (default = 'no')
% cfg.eta = 'yes' or 'no' (default = 'no')
% cfg.kurtosis = 'yes' or 'no' (default = 'no')
% cfg.keeptrials = 'yes' or 'no' (default = 'no')
% cfg.keepcsd = 'yes' or 'no' (default = 'no')
% cfg.keepnoisecsd = 'yes' or 'no' (default = 'no')
% cfg.keepmom = 'yes' or 'no' (default = 'yes')
% cfg.keepnoisemom = 'yes' or 'no' (default = 'yes')
% cfg.resolutionmatrix = 'yes' or 'no' (default = 'no')
% cfg.feedback = 'no', 'text' (default), 'textbar', 'gui'
%
% The following option only applies to timecourses.
% cfg.flipori = 'yes' or 'no' (default = 'no')
%
% The following option only applies to single-trial timecourses.
% cfg.fixedori = 'within_trials' or 'over_trials' (default = 'over_trials')
%
% If repeated trials are present that have undergone some sort of
% resampling (i.e. jackknife, bootstrap, singletrial or rawtrial), the mean,
% variance and standard error of mean will be computed for all source
% parameters. This is done after applying the optional transformation
% on the power and projected noise.
%
% 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_SOURCEANALYSIS, FT_SOURCESTATISTICS, FT_MATH
% Copyright (C) 2004-2015, Robert Oostenveld & Jan-Mathijs Schoffelen
%
% 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$
% 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 source
ft_preamble provenance source
% 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
% source = ft_checkdata(source, 'datatype', 'source', 'feedback', 'yes');
% cfg = ft_checkconfig(cfg, 'forbidden', {'trials'}); % trial selection is not implented here, you may want to consider ft_selectdata
% DEPRECATED by roboos on 13 June 2013
% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2199 for more details
% support for this functionality can be removed at the end of 2013
cfg = ft_checkconfig(cfg, 'deprecated', {'transform'}); % please use ft_math instead
% set the defaults
cfg.transform = ft_getopt(cfg, 'transform', []);
cfg.projectmom = ft_getopt(cfg, 'projectmom', 'no'); % if yes -> svdfft
cfg.numcomp = ft_getopt(cfg, 'numcomp', 1);
cfg.powmethod = ft_getopt(cfg, 'powmethod', []); % see below
cfg.cohmethod = ft_getopt(cfg, 'cohmethod', []); % see below
cfg.feedback = ft_getopt(cfg, 'feedback', 'textbar');
cfg.supmethod = ft_getopt(cfg, 'supmethod', 'none');
cfg.resolutionmatrix = ft_getopt(cfg, 'resolutionmatrix', 'no');
cfg.eta = ft_getopt(cfg, 'eta', 'no');
cfg.fa = ft_getopt(cfg, 'fa', 'no');
cfg.kurtosis = ft_getopt(cfg, 'kurtosis', 'no');
cfg.keeptrials = ft_getopt(cfg, 'keeptrials', 'no');
cfg.trials = ft_getopt(cfg, 'trials', 'all');
cfg.keepcsd = ft_getopt(cfg, 'keepcsd', 'no');
cfg.keepmom = ft_getopt(cfg, 'keepmom', 'yes');
cfg.keepnoisecsd = ft_getopt(cfg, 'keepnoisecsd', 'no');
cfg.keepnoisemom = ft_getopt(cfg, 'keepnoisemom', 'yes');
cfg.fwhm = ft_getopt(cfg, 'fwhm', 'no');
cfg.fwhmremovecenter = ft_getopt(cfg, 'fwhmremovecenter', 0);
cfg.fwhmmethod = ft_getopt(cfg, 'fwhmmethod', 'barnes');
cfg.fwhmmaxdist = ft_getopt(cfg, 'fwhmmaxdist', []);
cfg.fixedori = ft_getopt(cfg, 'fixedori', 'over_trials');
cfg.flipori = ft_getopt(cfg, 'flipori', 'no');
% only works for mne
cfg.demean = ft_getopt(cfg, 'demean', 'yes');
cfg.baselinewindow = ft_getopt(cfg, 'baselinewindow', [-inf 0]);
cfg.zscore = ft_getopt(cfg, 'zscore', 'yes');
zscore = strcmp(cfg.zscore, 'yes');
demean = strcmp(cfg.demean, 'yes');
if ischar(cfg.trials) && strcmp(cfg.trials,'all')
% do nothing
elseif ischar(cfg.trials)
ft_error('only ''all'' is allowed for string input for cfg.trials');
else
% check whether there's a trial field in the source structure, and
% subselect, otherwise error
if isfield(source, 'trial')
source.trial = source.trial(cfg.trials);
if isfield(source, 'cumtapcnt'), source.cumtapcnt = source.cumtapcnt(cfg.trials,:); end
else
ft_error('subselecting trials in ft_sourcedescriptives is currently only possible with a ''trial'' field');
end
end
% get desired method from source structure
source.method = ft_getopt(source,'method',[]);
% this is required for backward compatibility with the old sourceanalysis
if isfield(source, 'method') && strcmp(source.method, 'randomized')
source.method = 'randomization';
elseif isfield(source, 'method') && strcmp(source.method, 'permuted')
source.method = 'permutation';
elseif isfield(source, 'method') && strcmp(source.method, 'jacknife')
source.method = 'jackknife';
end
% determine the type of data, this is only relevant for a few specific types
ispccdata = isfield(source, 'avg') && isfield(source.avg, 'csdlabel');
islcmvavg = isfield(source, 'avg') && isfield(source, 'time') && isfield(source.avg, 'mom') && any(size(source.avg.pow)==1);
islcmvtrl = isfield(source, 'trial') && isfield(source, 'time') && isfield(source.trial, 'mom');
ismneavg = isfield(source, 'avg') && isfield(source, 'time') && isfield(source.avg, 'mom') && size(source.avg.pow, 2)==numel(source.time);
% check the consistency of the defaults
if strcmp(cfg.projectmom, 'yes')
if isempty(cfg.powmethod)
cfg.powmethod = 'regular'; % set the default
elseif ~strcmp(cfg.powmethod, 'regular')
ft_error('unsupported powmethod in combination with projectmom');
end
if isempty(cfg.cohmethod)
cfg.cohmethod = 'regular'; % set the default
elseif ~strcmp(cfg.cohmethod, 'regular')
ft_error('unsupported cohmethod in combination with projectmom');
end
else
if isempty(cfg.powmethod)
cfg.powmethod = 'lambda1'; % set the default
end
if isempty(cfg.cohmethod)
cfg.cohmethod = 'lambda1'; % set the default
end
end
% this is required for backward compatibility with an old version of sourcedescriptives
if isfield(cfg, 'singletrial'), cfg.keeptrials = cfg.singletrial; end
% do a validity check on the input data and specified options
if strcmp(cfg.resolutionmatrix, 'yes')
if ~isfield(source.avg, 'filter')
ft_error('The computation of the resolution matrix requires keepfilter=''yes'' in sourceanalysis.');
elseif ~isfield(source, 'leadfield')
ft_error('The computation of the resolution matrix requires keepleadfield=''yes'' in sourceanalysis.');
end
end
if istrue(cfg.fwhm) && ~isfield(source.avg, 'filter')
ft_error('computation of the fwhm requires keepfilter=''yes'' in sourceanalysis.');
end
if istrue(cfg.eta) && strcmp(cfg.cohmethod, 'svdfft')
ft_error('eta cannot be computed in combination with the application of svdfft');
end
if istrue(cfg.keeptrials) && ~strcmp(cfg.supmethod, 'none')
ft_error('you cannot keep trials when you want to partialize something');
end
if istrue(cfg.flipori) && ~istrue(cfg.projectmom)
ft_error('flipori requires projectmom=''yes''');
end
% set some flags for convenience
isnoise = isfield(source, 'avg') && isfield(source.avg, 'noisecsd');
keeptrials = strcmp(cfg.keeptrials, 'yes');
projectmom = strcmp(cfg.projectmom, 'yes');
flipori = strcmp(cfg.flipori, 'yes');
% determine the subfunction used for computing power
switch cfg.powmethod
case 'regular'
powmethodfun = @powmethod_regular;
case 'lambda1'
powmethodfun = @powmethod_lambda1;
case 'trace'
powmethodfun = @powmethod_trace;
case 'none'
powmethodfun = [];
otherwise
ft_error('unsupported powmethod');
end
% represent the selection of sources in the brain as a row-vector with indices
insideindx = find(source.inside(:)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ispccdata
% the source reconstruction was computed using the pcc beamformer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ndipole = size(source.pos,1);
if ischar(source.avg.csdlabel{1}), source.avg.csdlabel = {source.avg.csdlabel}; end
if numel(source.avg.csdlabel)==1
source.avg.csdlabel = repmat(source.avg.csdlabel, [Ndipole 1]);
end
dipsel = find(strcmp(source.avg.csdlabel{1}, 'scandip'));
refchansel = find(strcmp(source.avg.csdlabel{1}, 'refchan'));
refdipsel = find(strcmp(source.avg.csdlabel{1}, 'refdip'));
supchansel = find(strcmp(source.avg.csdlabel{1}, 'supchan'));
supdipsel = find(strcmp(source.avg.csdlabel{1}, 'supdip'));
% cannot handle reference channels and reference dipoles simultaneously
if numel(refchansel)>0 && numel(refdipsel)>0
ft_error('cannot simultaneously handle reference channels and reference dipole');
end
% these are only used to count the number of reference/suppression dipoles and channels
refsel = [refdipsel refchansel];
supsel = [supdipsel supchansel];
% first do the projection of the moment, if requested
if projectmom
source.avg.ori = cell(1, Ndipole);
ft_progress('init', cfg.feedback, 'projecting dipole moment');
for i=insideindx
ft_progress(i/length(insideindx), 'projecting dipole moment %d/%d\n', i, length(insideindx));
if numel(source.avg.csdlabel)>1
dipsel = find(strcmp(source.avg.csdlabel{i}, 'scandip'));
refchansel = find(strcmp(source.avg.csdlabel{i}, 'refchan'));
refdipsel = find(strcmp(source.avg.csdlabel{i}, 'refdip'));
supchansel = find(strcmp(source.avg.csdlabel{i}, 'supchan'));
supdipsel = find(strcmp(source.avg.csdlabel{i}, 'supdip'));
% these are only used to count the number of reference/suppression dipoles and channels
refsel = [refdipsel refchansel];
supsel = [supdipsel supchansel];
end
mom = source.avg.mom{i}(dipsel, :);
ref = source.avg.mom{i}(refdipsel, :);
sup = source.avg.mom{i}(supdipsel, :);
refchan = source.avg.mom{i}(refchansel, :);
supchan = source.avg.mom{i}(supchansel, :);
% compute the projection of the scanning dipole along the direction of the dominant amplitude
if length(dipsel)>1, [mom, rmom] = svdfft(mom, cfg.numcomp, source.cumtapcnt); else rmom = []; end
source.avg.ori{i} = rmom;
% compute the projection of the reference dipole along the direction of the dominant amplitude
if length(refdipsel)>1, [ref, rref] = svdfft(ref, 1, source.cumtapcnt); else rref = []; end
% compute the projection of the supression dipole along the direction of the dominant amplitude
if length(supdipsel)>1, [sup, rsup] = svdfft(sup, 1, source.cumtapcnt); else rsup = []; end
% compute voxel-level fourier-matrix
source.avg.mom{i} = cat(1, mom, ref, sup, refchan, supchan);
% create rotation-matrix
rotmat = zeros(0, length(source.avg.csdlabel{i}));
if ~isempty(rmom)
rotmat = [rotmat; rmom zeros(1,numel(refsel)+numel(supsel))];
end
if ~isempty(rref)
rotmat = [rotmat; zeros(1, numel(dipsel)), rref, zeros(1,numel(refchansel)+numel(supsel))];
end
if ~isempty(rsup)
rotmat = [rotmat; zeros(1, numel(dipsel)+numel(refdipsel)), rsup, zeros(1,numel(refchansel)+numel(supchansel))];
end
for j=1:length(supchansel)
rotmat(end+1,:) = 0;
rotmat(end,numel(dipsel)+numel(refdipsel)+numel(supdipsel)+j) = 1;
end
for j=1:length(refchansel)
rotmat(end+1,:) = 0;
rotmat(end,numel(dipsel)+numel(refdipsel)+numel(supdipsel)+numel(supchansel)+j) = 1;
end
% compute voxel-level csd-matrix
if isfield(source.avg, 'csd'), source.avg.csd{i} = rotmat * source.avg.csd{i} * rotmat'; end
% compute voxel-level noisecsd-matrix
if isfield(source.avg, 'noisecsd'), source.avg.noisecsd{i} = rotmat * source.avg.noisecsd{i} * rotmat'; end
% compute rotated filter
if isfield(source.avg, 'filter'), source.avg.filter{i} = rotmat * source.avg.filter{i}; end
if isfield(source.avg, 'csdlabel')
% remember what the interpretation is of all CSD output components
scandiplabel = repmat({'scandip'}, 1, cfg.numcomp); % only one dipole orientation remains
refdiplabel = repmat({'refdip'}, 1, length(refdipsel)>0); % for svdfft at max. only one dipole orientation remains
supdiplabel = repmat({'supdip'}, 1, length(supdipsel)>0); % for svdfft at max. only one dipole orientation remains
refchanlabel = repmat({'refchan'}, 1, length(refchansel));
supchanlabel = repmat({'supchan'}, 1, length(supchansel));
% concatenate all the labels
source.avg.csdlabel{i} = cat(2, scandiplabel, refdiplabel, supdiplabel, refchanlabel, supchanlabel);
end
% compute rotated leadfield
% FIXME in the presence of a refdip and/or supdip, this does not work; leadfield is Nx3
if isfield(source, 'leadfield')
% FIXME this is a proposed dirty fix
n1 = size(source.leadfield{i},2);
%n2 = size(rotmat,2) - n1;
n2 = size(rotmat,2) - n1 +1; %added 1 JM
source.leadfield{i} = source.leadfield{i} * rotmat(1:n2, 1:n1)';
end
end % for i=insideindx
ft_progress('close');
% update the indices
dipsel = find(strcmp(source.avg.csdlabel, 'scandip'));
refchansel = find(strcmp(source.avg.csdlabel, 'refchan'));
refdipsel = find(strcmp(source.avg.csdlabel, 'refdip'));
supchansel = find(strcmp(source.avg.csdlabel, 'supchan'));
supdipsel = find(strcmp(source.avg.csdlabel, 'supdip'));
refsel = [refdipsel refchansel];
supsel = [supdipsel supchansel];
end % if projectmom
if keeptrials
cumtapcnt = source.cumtapcnt(:);
sumtapcnt = cumsum([0;cumtapcnt]);
Ntrial = length(cumtapcnt);
ft_progress('init', cfg.feedback, 'computing singletrial voxel-level cross-spectral densities');
for triallop = 1:Ntrial
source.trial(triallop).csd = cell(Ndipole, 1); % allocate memory for this trial
source.trial(triallop).mom = cell(Ndipole, 1); % allocate memory for this trial
ft_progress(triallop/Ntrial, 'computing singletrial voxel-level cross-spectral densities %d%d\n', triallop, Ntrial);
for i=insideindx
dat = source.avg.mom{i};
tmpmom = dat(:, sumtapcnt(triallop)+1:sumtapcnt(triallop+1));
tmpcsd = (tmpmom * tmpmom') ./cumtapcnt(triallop);
source.trial(triallop).mom{i} = tmpmom;
source.trial(triallop).csd{i} = tmpcsd;
end % for i=insideindx
end % for triallop
ft_progress('close');
% remove the average, continue with separate trials, but keep track of
% the csdlabel
csdlabel = source.avg.csdlabel;
source = rmfield(source, 'avg');
else
fprintf('using average voxel-level cross-spectral densities\n');
csdlabel = source.avg.csdlabel;
end % if keeptrials
% process the csdlabel for each of the dipoles
hasrefdip = true;
hasrefchan = true;
hassupdip = true;
hassupchan = true;
dipselcell = cell(Ndipole,1);
refdipselcell = cell(Ndipole,1);
refchanselcell = cell(Ndipole,1);
supdipselcell = cell(Ndipole,1);
supchanselcell = cell(Ndipole,1);
for i = insideindx
dipsel = find(strcmp(csdlabel{i}, 'scandip'));
refchansel = find(strcmp(csdlabel{i}, 'refchan'));
refdipsel = find(strcmp(csdlabel{i}, 'refdip'));
supchansel = find(strcmp(csdlabel{i}, 'supchan'));
supdipsel = find(strcmp(csdlabel{i}, 'supdip'));
hasrefdip = ~isempty(refdipsel) && hasrefdip; %NOTE: it has to be true for all dipoles!
hasrefchan = ~isempty(refchansel) && hasrefchan;
hassupdip = ~isempty(supdipsel) && hassupdip;
hassupchan = ~isempty(supchansel) && hassupchan;
dipselcell{i} = dipsel;
refdipselcell{i} = refdipsel;
refchanselcell{i} = refchansel;
supdipselcell{i} = supdipsel;
supchanselcell{i} = supchansel;
end
if keeptrials
% do the processing of the CSD matrices for each trial
if ~strcmp(cfg.supmethod, 'none')
ft_error('suppression is only supported for average CSD');
end
%dipselcell = mat2cell(repmat(dipsel(:)', [Ndipole 1]), ones(Ndipole,1), length(dipsel));
%if hasrefdip, refdipselcell = mat2cell(repmat(refdipsel(:)', [Ndipole 1]), ones(Ndipole,1), length(refdipsel)); end
%if hasrefchan, refchanselcell = mat2cell(repmat(refchansel(:)', [Ndipole 1]), ones(Ndipole,1), length(refchansel)); end
%if hassupdip, supdipselcell = mat2cell(repmat(supdipsel(:)', [Ndipole 1]), ones(Ndipole,1), length(supdipsel)); end
%if hassupchan, supchanselcell = mat2cell(repmat(supchansel(:)', [Ndipole 1]), ones(Ndipole,1), length(supchansel)); end
ft_progress('init', cfg.feedback, 'computing singletrial voxel-level power');
for triallop = 1:Ntrial
%initialize the variables
source.trial(triallop).pow = zeros(Ndipole, 1);
if hasrefdip, source.trial(triallop).refdippow = zeros(Ndipole, 1); end
if hasrefchan, source.trial(triallop).refchanpow = zeros(Ndipole, 1); end
if hassupdip, source.trial(triallop).supdippow = zeros(Ndipole, 1); end
if hassupchan, source.trial(triallop).supchanpow = zeros(Ndipole, 1); end
ft_progress(triallop/Ntrial, 'computing singletrial voxel-level power %d%d\n', triallop, Ntrial);
source.trial(triallop).pow(source.inside) = cellfun(powmethodfun, source.trial(triallop).csd(source.inside), dipselcell(source.inside));
if hasrefdip, source.trial(triallop).refdippow(source.inside) = cellfun(powmethodfun,source.trial(triallop).csd(source.inside), refdipselcell(source.inside)); end
if hassupdip, source.trial(triallop).supdippow(source.inside) = cellfun(powmethodfun,source.trial(triallop).csd(source.inside), supdipselcell(source.inside)); end
if hasrefchan, source.trial(triallop).refchanpow(source.inside) = cellfun(powmethodfun,source.trial(triallop).csd(source.inside), refchanselcell(source.inside)); end
if hassupchan, source.trial(triallop).supchanpow(source.inside) = cellfun(powmethodfun,source.trial(triallop).csd(source.inside), supchanselcell(source.inside)); end
%FIXME kan volgens mij niet
if isnoise && isfield(source.trial(triallop), 'noisecsd')
% compute the power of the noise projected on each source component
source.trial(triallop).noise = cellfun(powmethodfun,source.trial(triallop).csd, dipselcell);
if hasrefdip, source.trial(triallop).refdipnoise = cellfun(powmethodfun,source.trial(triallop).noisecsd, refdipselcell); end
if hassupdip, source.trial(triallop).supdipnoise = cellfun(powmethodfun,source.trial(triallop).noisecsd, supdipselcell); end
if hasrefchan, source.trial(triallop).refchannoise = cellfun(powmethodfun,source.trial(triallop).noisecsd, refchanselcell); end
if hassupchan, source.trial(triallop).supchannoise = cellfun(powmethodfun,source.trial(triallop).noisecsd, supchanselcell); end
end % if isnoise
end % for triallop
ft_progress('close');
if strcmp(cfg.keepcsd, 'no')
source.trial = rmfield(source.trial, 'csd');
end
else
% do the processing of the average CSD matrix
for i=insideindx
switch cfg.supmethod
case 'chan_dip'
supindx = [supdipsel supchansel];
if i==insideindx(1), refsel = refsel - length(supdipsel); end % adjust index only once
refchanselcell{i} = refchanselcell{i} - length(supdipsel);
refdipselcell{i} = refdipselcell{i} - length(supdipsel);
case 'chan'
supindx = supchansel;
case 'dip'
supindx = supdipsel;
if i==insideindx(1), refsel = refsel - length(supdipsel); end
refchanselcell{i} = refchanselcell{i} - length(supdipsel);
refdipselcell{i} = refdipselcell{i} - length(supdipsel);
case 'none'
% do nothing
supindx = [];
end
tmpcsd = source.avg.csd{i};
scnindx = setdiff(1:size(tmpcsd,1), supindx);
tmpcsd = tmpcsd(scnindx, scnindx) - tmpcsd(scnindx, supindx)*pinv(tmpcsd(supindx, supindx))*tmpcsd(supindx, scnindx);
source.avg.csd{i} = tmpcsd;
end % for insideindx
% source.avg.csdlabel = source.avg.csdlabel(scnindx);
if isnoise && ~strcmp(cfg.supmethod, 'none')
source.avg = rmfield(source.avg, 'noisecsd');
end
% initialize the variables
source.avg.pow = nan(Ndipole, 1);
if hasrefdip, source.avg.refdippow = nan(Ndipole, 1); end
if hasrefchan, source.avg.refchanpow = nan(Ndipole, 1); end
if hassupdip, source.avg.supdippow = nan(Ndipole, 1); end
if hassupchan, source.avg.supchanpow = nan(Ndipole, 1); end
if isnoise
source.avg.noise = nan(Ndipole, 1);
if hasrefdip, source.avg.refdipnoise = nan(Ndipole, 1); end
if hasrefchan, source.avg.refchannoise = nan(Ndipole, 1); end
if hassupdip, source.avg.supdipnoise = nan(Ndipole, 1); end
if hassupchan, source.avg.supchannoise = nan(Ndipole, 1); end
end % if isnoise
if hasrefdip||hasrefchan, source.avg.coh = nan(Ndipole, 1); end
if strcmp(cfg.eta, 'yes')
source.avg.eta = nan(Ndipole, 1);
source.avg.ori = cell(1, Ndipole);
end
if strcmp(cfg.eta, 'yes') && ~isempty(refsel)
source.avg.etacsd = nan(Ndipole, 1);
source.avg.ucsd = cell(1, Ndipole);
end
if strcmp(cfg.fa, 'yes')
source.avg.fa = nan(Ndipole, 1);
end
for i=insideindx
dipsel = dipselcell{i};
refchansel = refchanselcell{i};
refdipsel = refdipselcell{i};
refsel = [refchansel refdipsel];
supchansel = supchanselcell{i};
% compute the power of each source component
if strcmp(cfg.projectmom, 'yes') && cfg.numcomp>1
source.avg.pow(i) = powmethodfun(source.avg.csd{i}(dipselcell{i},dipselcell{i}), 1);
else
source.avg.pow(i) = powmethodfun(source.avg.csd{i}(dipselcell{i},dipselcell{i}));
end
if hasrefdip, source.avg.refdippow(i) = powmethodfun(source.avg.csd{i}(refdipsel,refdipsel)); end
%if hassupdip, source.avg.supdippow(i) = powmethodfun(source.avg.csd{i}(supdipsel,supdipsel)); end
if hasrefchan, source.avg.refchanpow(i) = powmethodfun(source.avg.csd{i}(refchansel,refchansel)); end
%if hassupchan, source.avg.supchanpow(i) = powmethodfun(source.avg.csd{i}(supchansel,supchansel)); end
if isnoise
% compute the power of the noise projected on each source component
if strcmp(cfg.projectmom, 'yes') && cfg.numcomp>1
source.avg.noise(i) = powmethodfun(source.avg.noisecsd{i}(dipselcell{i},dipselcell{i}), 1);
else
source.avg.noise(i) = powmethodfun(source.avg.noisecsd{i}(dipselcell{i},dipselcell{i}));
end
if hasrefdip, source.avg.refdipnoise(i) = powmethodfun(source.avg.noisecsd{i}(refdipsel,refdipsel)); end
if hassupdip, source.avg.supdipnoise(i) = powmethodfun(source.avg.noisecsd{i}(supdipsel,supdipsel)); end
if hasrefchan, source.avg.refchannoise(i) = powmethodfun(source.avg.noisecsd{i}(refchansel,refchansel)); end
if hassupchan, source.avg.supchannoise(i) = powmethodfun(source.avg.noisecsd{i}(supchansel,supchansel)); end
end % if isnoise
if ~isempty(refsel)
% compute coherence
csd = source.avg.csd{i};
switch cfg.cohmethod
case 'regular'
% assume that all dipoles have been projected along the direction of maximum power
Pd = abs(csd(dipsel, dipsel));
Pr = abs(csd(refsel, refsel));
Cdr = csd(dipsel, refsel);
source.avg.coh(i) = (Cdr.^2) ./ (Pd*Pr);
case 'lambda1'
% compute coherence the Joachim Gross' way
Pd = lambda1(csd(dipsel, dipsel));
Pr = lambda1(csd(refsel, refsel));
Cdr = lambda1(csd(dipsel, refsel));
source.avg.coh(i) = abs(Cdr).^2 ./ (Pd*Pr);
case 'canonical'
% compute canonical coherence
ccoh = ft_connectivity_cancorr(csd([dipsel refsel],[dipsel refsel]), 'indices', [ones(1,numel(dipsel)) ones(1,numel(refsel))*2]);
source.avg.coh(i) = ccoh(1,2);
otherwise
ft_error('unsupported cohmethod');
end % cohmethod
end
% compute eta
if strcmp(cfg.eta, 'yes')
[source.avg.eta(i), source.avg.ori{i}] = csd2eta(source.avg.csd{i}(dipselcell{i},dipselcell{i}));
if ~isempty(refsel)
%FIXME this only makes sense when only a reference signal OR a dipole is selected
[source.avg.etacsd(i), source.avg.ucsd{i}] = csd2eta(source.avg.csd{i}(dipsel,refsel));
end
end
%compute fa
if strcmp(cfg.fa, 'yes')
source.avg.fa(i) = csd2fa(source.avg.csd{i}(dipsel,dipsel));
end
end % for diplop
if strcmp(cfg.keepcsd, 'no')
source.avg = rmfield(source.avg, 'csd');
end
if strcmp(cfg.keepnoisecsd, 'no') && isnoise
source.avg = rmfield(source.avg, 'noisecsd');
end
end % if keeptrials
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ismneavg
% the source reconstruction was computed using mne and contains an average timecourse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if demean
begsmp = nearest(source.time, cfg.baselinewindow(1));
endsmp = nearest(source.time, cfg.baselinewindow(2));
ft_progress('init', cfg.feedback, 'baseline correcting dipole moments');
for diplop=1:length(insideindx)
ft_progress(diplop/length(insideindx), 'baseline correcting dipole moments %d/%d\n', diplop, length(insideindx));
mom = source.avg.mom{insideindx(diplop)};
mom = ft_preproc_baselinecorrect(mom, begsmp, endsmp);
source.avg.mom{insideindx(diplop)} = mom;
end
ft_progress('close');
end
if projectmom
if isfield(source, 'tri')
nrm = normals(source.pos, source.tri, 'vertex');
source.avg.phi = zeros(size(source.pos,1),1);
end
ft_progress('init', cfg.feedback, 'projecting dipole moment');
for diplop=1:length(insideindx)
ft_progress(diplop/length(insideindx), 'projecting dipole moment %d/%d\n', diplop, length(insideindx));
mom = source.avg.mom{insideindx(diplop)};
[mom, rmom] = svdfft(mom, 1);
source.avg.mom{insideindx(diplop)} = mom;
source.avg.ori{insideindx(diplop)} = rmom;
end
if isfield(source, 'tri')
for diplop=insideindx
source.avg.phi(diplop) = source.avg.ori{diplop}*nrm(diplop,:)';
end
end
if isfield(source.avg, 'noisecov')
source.avg.noise = nan+zeros(size(source.pos,1),1);
for diplop=insideindx
rmom = source.avg.ori{diplop};
source.avg.noise(diplop) = rmom*source.avg.noisecov{diplop}*rmom';
end
end
ft_progress('close');
end % if projectmom
if flipori
tmpmom = cat(1, source.avg.mom{source.inside});
[u, s, v] = svd(tmpmom, 'econ');
flip( source.inside) = sign(u(:,1));
flip(~source.inside) = nan;
for i=1:numel(source.inside)
if source.inside(i)
source.avg.mom{i} = flip(i) * source.avg.mom{i};
source.avg.ori{i} = flip(i) * source.avg.ori{i};
end
end
end % if flipori
if zscore
begsmp = nearest(source.time, cfg.baselinewindow(1));
endsmp = nearest(source.time, cfg.baselinewindow(2));
% zscore using baselinewindow for power
ft_progress('init', cfg.feedback, 'computing power');
%source.avg.absmom = source.avg.pow;
for diplop=1:length(insideindx)
ft_progress(diplop/length(insideindx), 'computing power %d/%d\n', diplop, length(insideindx));
mom = source.avg.mom{insideindx(diplop)};
mmom = mean(mom(:,begsmp:endsmp),2);
smom = std(mom(:,begsmp:endsmp),[],2);
pow = sum(((mom-mmom(:,ones(size(mom,2),1)))./smom(:,ones(size(mom,2),1))).^2,1);
source.avg.pow(insideindx(diplop),:) = pow;
source.avg.mom{insideindx(diplop)} = diag(1./smom)*mom;
end
ft_progress('close');
else
% just square for power
ft_progress('init', cfg.feedback, 'computing power');
%source.avg.absmom = source.avg.pow;
for diplop=1:length(insideindx)
ft_progress(diplop/length(insideindx), 'computing power %d/%d\n', diplop, length(insideindx));
mom = source.avg.mom{insideindx(diplop)};
pow = sum(mom.^2,1);
source.avg.pow(insideindx(diplop),:) = pow;
%source.avg.absmom(insideindx(diplop),:) = sum(mom,1);
end
ft_progress('close');
end % if zscore
if strcmp(cfg.kurtosis, 'yes')
fprintf('computing kurtosis based on dipole timecourse\n');
source.avg.k2 = nan(size(source.pos,1),1);
for diplop=1:length(insideindx)
mom = source.avg.mom{insideindx(diplop)};
if length(mom)~=prod(size(mom))
ft_error('kurtosis can only be computed for projected dipole moment');
end
source.avg.k2(insideindx(diplop)) = kurtosis(mom);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif islcmvavg
% the source reconstruction was computed using the lcmv beamformer and contains an average timecourse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if projectmom
ft_progress('init', cfg.feedback, 'projecting dipole moment');
for diplop=1:length(insideindx)
ft_progress(diplop/length(insideindx), 'projecting dipole moment %d/%d\n', diplop, length(insideindx));
mom = source.avg.mom{insideindx(diplop)};
[mom, rmom] = svdfft(mom, 1);
source.avg.mom{insideindx(diplop)} = mom;
source.avg.ori{insideindx(diplop)} = rmom;
end
ft_progress('close');
end
if flipori
tmpmom = cat(1, source.avg.mom{source.inside});
[u, s, v] = svd(tmpmom, 'econ');
flip( source.inside) = sign(u(:,1));
flip(~source.inside) = nan;
for i=1:numel(source.inside)
if source.inside(i)
source.avg.mom{i} = flip(i) * source.avg.mom{i};
source.avg.ori{i} = flip(i) * source.avg.ori{i};
end
end
end
if ~strcmp(cfg.powmethod, 'none')
fprintf('recomputing power based on dipole timecourse\n')
source.avg.pow = nan(size(source.pos,1),1);
for diplop=1:length(insideindx)
mom = source.avg.mom{insideindx(diplop)};
cov = mom * mom';
source.avg.pow(insideindx(diplop)) = powmethodfun(cov);
end
end
if strcmp(cfg.kurtosis, 'yes')
fprintf('computing kurtosis based on dipole timecourse\n');
source.avg.k2 = nan(size(source.pos,1),1);
for diplop=1:length(insideindx)
mom = source.avg.mom{insideindx(diplop)};
if length(mom)~=prod(size(mom))
ft_error('kurtosis can only be computed for projected dipole moment');
end
source.avg.k2(insideindx(diplop)) = kurtosis(mom);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif islcmvtrl
% the source reconstruction was computed using the lcmv beamformer and contains a single-trial timecourse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ntrial = length(source.trial);
if projectmom && strcmp(cfg.fixedori, 'within_trials')
% the dipole orientation is re-determined for each trial
ft_progress('init', cfg.feedback, 'projecting dipole moment');
for trllop=1:ntrial
ft_progress(trllop/ntrial, 'projecting dipole moment %d/%d\n', trllop, ntrial);
for diplop=1:length(insideindx)
mom = source.trial(trllop).mom{insideindx(diplop)};
[mom, rmom] = svdfft(mom, 1);
source.trial(trllop).mom{insideindx(diplop)} = mom;
source.trial(trllop).ori{insideindx(diplop)} = rmom; % remember the orientation
end
end
ft_progress('close');
elseif projectmom && strcmp(cfg.fixedori, 'over_trials')
ft_progress('init', cfg.feedback, 'projecting dipole moment');
% compute average covariance over all trials
for trllop=1:ntrial
for diplop=1:length(insideindx)
mom = source.trial(trllop).mom{insideindx(diplop)};
if trllop==1
cov{diplop} = mom*mom'./size(mom,2);
else
cov{diplop} = mom*mom'./size(mom,2) + cov{diplop};
end
end
end
% compute source orientation over all trials
for diplop=1:length(insideindx)
[dum, ori{diplop}] = svdfft(cov{diplop}, 1);
end
% project the data in each trial
for trllop=1:ntrial
ft_progress(trllop/ntrial, 'projecting dipole moment %d/%d\n', trllop, ntrial);
for diplop=1:length(insideindx)
mom = source.trial(trllop).mom{insideindx(diplop)};
mom = ori{diplop}*mom;
source.trial(trllop).mom{insideindx(diplop)} = mom;
source.trial(trllop).ori{insideindx(diplop)} = ori{diplop};
end
end
ft_progress('close');
end
if ~strcmp(cfg.powmethod, 'none')
fprintf('recomputing power based on dipole timecourse\n')
for trllop=1:ntrial
for diplop=1:length(insideindx)
mom = source.trial(trllop).mom{insideindx(diplop)};
cov = mom * mom';
source.trial(trllop).pow(insideindx(diplop)) = powmethodfun(cov);
end
end
end
if strcmp(cfg.kurtosis, 'yes')
fprintf('computing kurtosis based on dipole timecourse\n');
for trllop=1:ntrial
source.trial(trllop).k2 = nan(size(source.pos,1),1);
for diplop=1:length(insideindx)
mom = source.trial(trllop).mom{insideindx(diplop)};
if length(mom)~=numel(mom)
ft_error('kurtosis can only be computed for projected dipole moment');
end
source.trial(trllop).k2(insideindx(diplop)) = kurtosis(mom);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % dealing with pcc, lcmv, dics or mne input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfield(source, 'avg') && isfield(source.avg, 'pow') && isfield(source.avg, 'noise') && ~ismneavg
% compute the neural activity index for the average
source.avg.nai = source.avg.pow(:) ./ source.avg.noise(:);
end
if isfield(source, 'trial') && isfield(source.trial, 'pow') && isfield(source.trial, 'noise')
% compute the neural activity index for the trials
ntrials = length(source.trial);
for trlop=1:ntrials
source.trial(trlop).nai = source.trial(trlop).pow ./ source.trial(trlop).noise;
end
end
if strcmp(source.method, 'randomization') || strcmp(source.method, 'permutation')
% compute the neural activity index for the two randomized conditions
source.avgA.nai = source.avgA.pow ./ source.avgA.noise;
source.avgB.nai = source.avgB.pow ./ source.avgB.noise;
for trlop=1:length(source.trialA)
source.trialA(trlop).nai = source.trialA(trlop).pow ./ source.trialA(trlop).noise;
end
for trlop=1:length(source.trialB)
source.trialB(trlop).nai = source.trialB(trlop).pow ./ source.trialB(trlop).noise;
end
end
if ~isempty(cfg.transform)
fprintf('applying %s transformation on the power and projected noise\n', cfg.transform);
% apply the specified transformation on the power
if isfield(source, 'avg' ) && isfield(source.avg , 'pow'), source.avg .pow = feval(cfg.transform, source.avg .pow); end
if isfield(source, 'avgA' ) && isfield(source.avgA , 'pow'), source.avgA.pow = feval(cfg.transform, source.avgA.pow); end
if isfield(source, 'avgB' ) && isfield(source.avgB , 'pow'), source.avgB.pow = feval(cfg.transform, source.avgB.pow); end
if isfield(source, 'trial' ) && isfield(source.trial , 'pow'), for i=1:length(source.trial ), source.trial (i).pow = feval(cfg.transform, source.trial (i).pow); end; end
if isfield(source, 'trialA') && isfield(source.trialA, 'pow'), for i=1:length(source.trialA), source.trialA(i).pow = feval(cfg.transform, source.trialA(i).pow); end; end
if isfield(source, 'trialB') && isfield(source.trialB, 'pow'), for i=1:length(source.trialB), source.trialB(i).pow = feval(cfg.transform, source.trialB(i).pow); end; end
% apply the specified transformation on the projected noise
if isfield(source, 'avg' ) && isfield(source.avg , 'noise'), source.avg .noise = feval(cfg.transform, source.avg .noise); end
if isfield(source, 'avgA' ) && isfield(source.avgA , 'noise'), source.avgA.noise = feval(cfg.transform, source.avgA.noise); end
if isfield(source, 'avgB' ) && isfield(source.avgB , 'noise'), source.avgB.noise = feval(cfg.transform, source.avgB.noise); end
if isfield(source, 'trial' ) && isfield(source.trial , 'noise'), for i=1:length(source.trial ), source.trial (i).noise = feval(cfg.transform, source.trial (i).noise); end; end
if isfield(source, 'trialA') && isfield(source.trialA, 'noise'), for i=1:length(source.trialA), source.trialA(i).noise = feval(cfg.transform, source.trialA(i).noise); end; end
if isfield(source, 'trialB') && isfield(source.trialB, 'noise'), for i=1:length(source.trialB), source.trialB(i).noise = feval(cfg.transform, source.trialB(i).noise); end; end
end
if strcmp(source.method, 'pseudovalue')
% compute the pseudovalues for the beamformer output
avg = source.trial(1); % the first is the complete average
Ntrials = length(source.trial)-1; % the remaining are the leave-one-out averages
pseudoval = [];
if isfield(source.trial, 'pow')
allavg = getfield(avg, 'pow');
for i=1:Ntrials
thisavg = getfield(source.trial(i+1), 'pow');
thisval = Ntrials*allavg - (Ntrials-1)*thisavg;
pseudoval(i).pow = thisval;
end
end
if isfield(source.trial, 'coh')
allavg = getfield(avg, 'coh');
for i=1:Ntrials
thisavg = getfield(source.trial(i+1), 'coh');
thisval = Ntrials*allavg - (Ntrials-1)*thisavg;
pseudoval(i).coh = thisval;
end
end
if isfield(source.trial, 'nai')
allavg = getfield(avg, 'nai');
for i=1:Ntrials
thisavg = getfield(source.trial(i+1), 'nai');
thisval = Ntrials*allavg - (Ntrials-1)*thisavg;
pseudoval(i).nai = thisval;
end
end
if isfield(source.trial, 'noise')
allavg = getfield(avg, 'noise');
for i=1:Ntrials
thisavg = getfield(source.trial(i+1), 'noise');
thisval = Ntrials*allavg - (Ntrials-1)*thisavg;
pseudoval(i).noise = thisval;
end
end
% store the pseudovalues instead of the original values
source.trial = pseudoval;
end
if strcmp(source.method, 'jackknife') || strcmp(source.method, 'bootstrap') || strcmp(source.method, 'pseudovalue') || strcmp(source.method, 'singletrial') || strcmp(source.method, 'rawtrial')
% compute descriptive statistics (mean, var, sem) for multiple trial data
% compute these for as many source parameters as possible
% for convenience copy the trials out of the source structure
dip = source.trial;
% determine the (original) number of trials in the data
if strcmp(source.method, 'bootstrap') %VERANDERD ER ZAT GEEN .RESAMPLE IN SOURCE
Ntrials = size(source.trial,2); % WAS size(source.resample, 2);
else
Ntrials = length(source.trial);
end
fprintf('original data contained %d trials\n', Ntrials);
% allocate memory for all elements in the dipole structure
sumdip = [];
if isfield(dip(1), 'var'), sumdip.var = zeros(size(dip(1).var )); sumdip.var(~source.inside) = nan; end
if isfield(dip(1), 'pow'), sumdip.pow = zeros(size(dip(1).pow )); sumdip.pow(~source.inside) = nan; end
if isfield(dip(1), 'coh'), sumdip.coh = zeros(size(dip(1).coh )); sumdip.coh(~source.inside) = nan; end
if isfield(dip(1), 'rv'), sumdip.rv = zeros(size(dip(1).rv )); sumdip.rv(~source.inside) = nan; end
if isfield(dip(1), 'noise'), sumdip.noise = zeros(size(dip(1).noise)); sumdip.noise(~source.inside) = nan; end
if isfield(dip(1), 'nai'), sumdip.nai = zeros(size(dip(1).nai )); sumdip.nai(~source.inside) = nan; end
sqrdip = [];
if isfield(dip(1), 'var'), sqrdip.var = zeros(size(dip(1).var )); sqrdip.var(~source.inside) = nan; end
if isfield(dip(1), 'pow'), sqrdip.pow = zeros(size(dip(1).pow )); sqrdip.pow(~source.inside) = nan; end
if isfield(dip(1), 'coh'), sqrdip.coh = zeros(size(dip(1).coh )); sqrdip.coh(~source.inside) = nan; end
if isfield(dip(1), 'rv'), sqrdip.rv = zeros(size(dip(1).rv )); sqrdip.rv(~source.inside) = nan; end
if isfield(dip(1), 'noise'), sqrdip.noise = zeros(size(dip(1).noise)); sqrdip.noise(~source.inside) = nan; end
if isfield(dip(1), 'nai'), sqrdip.nai = zeros(size(dip(1).nai )); sqrdip.nai(~source.inside) = nan; end
if isfield(dip(1), 'mom')
sumdip.mom = cell(size(dip(1).mom));
sqrdip.mom = cell(size(dip(1).mom));
for i=1:length(dip(1).mom)
sumdip.mom{i} = zeros(size(dip(1).mom{i}));
sqrdip.mom{i} = zeros(size(dip(1).mom{i}));
end
end
if isfield(dip(1), 'csd')
sumdip.csd = cell(size(dip(1).csd));
sqrdip.csd = cell(size(dip(1).csd));
for i=1:length(dip(1).csd)
sumdip.csd{i} = zeros(size(dip(1).csd{i}));
sqrdip.csd{i} = zeros(size(dip(1).csd{i}));
end
end
for trial=1:length(dip)
% compute the sum of all values
if isfield(dip(trial), 'var'), sumdip.var = sumdip.var + dip(trial).var; end
if isfield(dip(trial), 'pow'), sumdip.pow = sumdip.pow + dip(trial).pow; end
if isfield(dip(trial), 'coh'), sumdip.coh = sumdip.coh + dip(trial).coh; end
if isfield(dip(trial), 'rv'), sumdip.rv = sumdip.rv + dip(trial).rv; end
if isfield(dip(trial), 'noise'), sumdip.noise = sumdip.noise + dip(trial).noise; end
if isfield(dip(trial), 'nai'), sumdip.nai = sumdip.nai + dip(trial).nai; end
% compute the sum of squared values
if isfield(dip(trial), 'var'), sqrdip.var = sqrdip.var + (dip(trial).var ).^2; end
if isfield(dip(trial), 'pow'), sqrdip.pow = sqrdip.pow + (dip(trial).pow ).^2; end
if isfield(dip(trial), 'coh'), sqrdip.coh = sqrdip.coh + (dip(trial).coh ).^2; end
if isfield(dip(trial), 'rv'), sqrdip.rv = sqrdip.rv + (dip(trial).rv ).^2; end
if isfield(dip(trial), 'noise'), sqrdip.noise = sqrdip.noise + (dip(trial).noise).^2; end
if isfield(dip(trial), 'nai'), sqrdip.nai = sqrdip.nai + (dip(trial).nai ).^2; end
% do the same for the cell-array with mom
if isfield(dip(trial), 'mom')
for i=1:length(dip(1).mom)
sumdip.mom{i} = sumdip.mom{i} + dip(trial).mom{i};
sqrdip.mom{i} = sqrdip.mom{i} + (dip(trial).mom{i}).^2;
end
end
% do the same for the cell-array with csd
if isfield(dip(trial), 'csd')
for i=1:length(dip(1).csd)
sumdip.csd{i} = sumdip.csd{i} + dip(trial).csd{i};
sqrdip.csd{i} = sqrdip.csd{i} + (dip(trial).csd{i}).^2;
end
end
end
% compute the mean over all repetitions
if isfield(sumdip, 'var'), dipmean.var = sumdip.var / length(dip); end
if isfield(sumdip, 'pow'), dipmean.pow = sumdip.pow / length(dip); end
if isfield(sumdip, 'coh'), dipmean.coh = sumdip.coh / length(dip); end
if isfield(sumdip, 'rv'), dipmean.rv = sumdip.rv / length(dip); end
if isfield(sumdip, 'noise'), dipmean.noise = sumdip.noise / length(dip); end
if isfield(sumdip, 'nai'), dipmean.nai = sumdip.nai / length(dip); end
% for the cell-array with mom, this is done further below
% for the cell-array with csd, this is done further below
% the estimates for variance and SEM are biased if we are working with the jackknife/bootstrap
% determine the proper variance scaling that corrects for this bias
% note that Ntrials is not always the same as the length of dip, especially in case of the bootstrap
if strcmp(source.method, 'singletrial')
bias = 1;
elseif strcmp(source.method, 'rawtrial')
bias = 1;
elseif strcmp(source.method, 'jackknife')
% Effron gives SEM estimate for the jackknife method in equation 11.5 (paragraph 11.2)
% to get the variance instead of SEM, we also have to multiply with the number of trials
bias = (Ntrials - 1)^2;
elseif strcmp(source.method, 'bootstrap')
% Effron gives SEM estimate for the bootstrap method in algorithm 6.1 (equation 6.6)
% to get the variance instead of SEM, we also have to multiply with the number of trials
bias = Ntrials;
elseif strcmp(source.method, 'pseudovalue')
% note that I have not put any thought in this aspect yet
ft_warning('don''t know how to compute bias for pseudovalue resampling');
bias = 1;
end
% compute the variance over all repetitions
if isfield(sumdip, 'var'), dipvar.var = bias*(sqrdip.var - (sumdip.var .^2)/length(dip))/(length(dip)-1); end
if isfield(sumdip, 'pow'), dipvar.pow = bias*(sqrdip.pow - (sumdip.pow .^2)/length(dip))/(length(dip)-1); end
if isfield(sumdip, 'coh'), dipvar.coh = bias*(sqrdip.coh - (sumdip.coh .^2)/length(dip))/(length(dip)-1); end
if isfield(sumdip, 'rv' ), dipvar.rv = bias*(sqrdip.rv - (sumdip.rv .^2)/length(dip))/(length(dip)-1); end
if isfield(sumdip, 'noise' ), dipvar.noise = bias*(sqrdip.noise - (sumdip.noise .^2)/length(dip))/(length(dip)-1); end
if isfield(sumdip, 'nai' ), dipvar.nai = bias*(sqrdip.nai - (sumdip.nai .^2)/length(dip))/(length(dip)-1); end
% compute the SEM over all repetitions
if isfield(sumdip, 'var'), dipsem.var = (dipvar.var /Ntrials).^0.5; end
if isfield(sumdip, 'pow'), dipsem.pow = (dipvar.pow /Ntrials).^0.5; end
if isfield(sumdip, 'coh'), dipsem.coh = (dipvar.coh /Ntrials).^0.5; end
if isfield(sumdip, 'rv' ), dipsem.rv = (dipvar.rv /Ntrials).^0.5; end
if isfield(sumdip, 'noise' ), dipsem.noise = (dipvar.noise /Ntrials).^0.5; end
if isfield(sumdip, 'nai' ), dipsem.nai = (dipvar.nai /Ntrials).^0.5; end
% compute the mean and SEM over all repetitions for the cell-array with mom
if isfield(dip(trial), 'mom')
for i=1:length(dip(1).mom)
dipmean.mom{i} = sumdip.mom{i}/length(dip);
dipvar.mom{i} = bias*(sqrdip.mom{i} - (sumdip.mom{i}.^2)/length(dip))/(length(dip)-1);
dipsem.mom{i} = (dipvar.mom{i}/Ntrials).^0.5;
end
end
% compute the mean and SEM over all repetitions for the cell-array with csd
if isfield(dip(trial), 'csd')
for i=1:length(dip(1).csd)
dipmean.csd{i} = sumdip.csd{i}/length(dip);
dipvar.csd{i} = bias*(sqrdip.csd{i} - (sumdip.csd{i}.^2)/length(dip))/(length(dip)-1);
dipsem.csd{i} = (dipvar.csd{i}/Ntrials).^0.5;
end
end
if strcmp(source.method, 'pseudovalue')
% keep the trials, since they have been converted to pseudovalues
% and hence the trials contain the interesting data
elseif keeptrials
% keep the trials upon request
else
% remove the original trials
source = rmfield(source, 'trial');
% assign the descriptive statistics to the output source structure
source.avg = dipmean;
source.var = dipvar;
source.sem = dipsem;
end
end
if strcmp(cfg.resolutionmatrix, 'yes')
% this is only implemented for pcc and no refdips/chans at the moment
Nchan = size(source.leadfield{insideindx(1)}, 1);
Ninside = length(insideindx);
allfilter = zeros(Ninside,Nchan);
allleadfield = zeros(Nchan,Ninside);
dipsel = match_str(source.avg.csdlabel, 'scandip');
ft_progress('init', cfg.feedback, 'computing resolution matrix');
for diplop=1:length(insideindx)
ft_progress(diplop/length(insideindx), 'computing resolution matrix %d/%d\n', diplop, length(insideindx));
% concatenate all filters
allfilter(diplop,:) = source.avg.filter{insideindx(diplop)}(dipsel,:);
% concatenate all leadfields
allleadfield(:,diplop) = source.leadfield{insideindx(diplop)};
end
ft_progress('close');
% multiply the filters and leadfields to obtain the resolution matrix
% see equation 1 and 2 in De Peralta-Menendez RG, Gonzalez-Andino SL: A critical analysis of linear inverse solutions to the neuroelectromagnetic inverse problem. IEEE Transactions on Biomedical Engineering 45: 440-448, 1998.
source.resolution = nan(Ndipole, Ndipole);
source.resolution(insideindx, insideindx) = allfilter*allleadfield;
end
% compute fwhm
if strcmp(cfg.fwhm, 'yes')
switch cfg.fwhmmethod
case 'barnes'
if ~isfield(source, 'dim')
ft_error('computation of fwhm is not possible with method ''barnes'' is not possible when the dipoles are not defined on a regular 3D grid');
end
fprintf('computing fwhm of spatial filters using method ''barnes''\n');
source = estimate_fwhm1(source, cfg.fwhmremovecenter);
case 'gaussfit'
fprintf('computing fwhm of spatial filters using method ''gaussfit''\n');
source = estimate_fwhm2(source, cfg.fwhmmaxdist);
otherwise
ft_error('unknown method for fwhm estimation');
end
end
% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous source
ft_postamble provenance source
ft_postamble history source
ft_postamble savevar source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute eta from a csd-matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [eta, u] = csd2eta(csd)
[u,s,v] = svd(real(csd));
eta = s(2,2)./s(1,1);
u = u'; %orientation is defined in the rows
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute fa from a csd-matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fa] = csd2fa(csd)
s = svd(real(csd));
ns = rank(real(csd));
s = s(1:ns);
ms = mean(s);
fa = sqrt( (ns./(ns-1)) .* (sum((s-ms).^2))./(sum(s.^2)) );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute power
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = powmethod_lambda1(x, ind)
if nargin==1
ind = 1:size(x,1);
end
s = svd(x(ind,ind));
p = s(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute power
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = powmethod_trace(x, ind)
if nargin==1
ind = 1:size(x,1);
end
p = trace(x(ind,ind));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute power
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = powmethod_regular(x, ind)
if nargin==1
ind = 1:size(x,1);
end
p = abs(x(ind,ind));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to obtain the largest singular value or trace of the
% source CSD matrices resulting from DICS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = lambda1(x)
s = svd(x);
s = s(1);