function [source] = sourcedescriptives(cfg, source)

% SOURCEDESCRIPTIVES computes descriptive parameters of the beamformer source
% analysis results.
%
% Use as: 
%   [source] = 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.transform  = string describing the transformation (default is [])
%
% The nai, i.e. neural activity index (power divided by projected noise),
% is computed prior to applying the optional transformation.
%
% Subsequently, the transformation is applied on the power and on the projected
% noise using "feval". A usefull transformation is for example 'log' or 'log10'.
%
% If repeated trials are present that have undergone some sort of
% resampling (i.e. jacknife, 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.
%
% See also SOURCEANALYSIS, SOURCESTATISTICS

% undocumented options
%   cfg.cohmethod = 'dics', 'svdfft' or 'canonical'

% Copyright (C) 2004, Robert Oostenveld
%
% $Log: sourcedescriptives.m,v $
% Revision 1.21  2006/01/30 17:10:33  roboos
% replaced fprintf with progress() indicator
%
% Revision 1.20  2006/01/24 14:31:05  roboos
% added cfg.cohmethod as undocumented option (not in main help, but just below it)
%
% Revision 1.19  2005/10/18 12:21:27  roboos
% fixed bug in power estimation from noise csd
% removed some obsolete code (was already commented out)
%
% Revision 1.18  2005/09/16 07:37:20  jansch
% in case of svd, only do svd on the real part
%
% Revision 1.17  2005/09/05 06:39:20  jansch
% experimental inclusion of partialisation in svdfft-option. no guarantees
% for success!
%
% Revision 1.16  2005/09/01 12:52:07  jansch
% included possibility to do svd on the single-trial ffts.
%
% Revision 1.15  2005/06/29 12:46:29  roboos
% the following changes were done on a large number of files at the same time, not all of them may apply to this specific file
% - added try-catch around the inclusion of input data configuration
% - moved cfg.version, cfg.previous and the assignment of the output cfg to the end
% - changed help comments around the configuration handling
% - some changes in whitespace
%
% Revision 1.14  2005/06/15 07:10:40  roboos
% fixed obvious typing mistake that caused error for mom field
%
% Revision 1.13  2005/06/13 08:04:20  roboos
% fixed a important bug in the computation of jacknife estimate of variance: the estimate was sqrt(N) times too small (i.e. it was the SEM instead of VAR)
% made the computation of the variance more readible by separating the bias term from the 1/(N-1) division
%
% Revision 1.12  2005/05/20 16:50:19  roboos
% added normalisation of the CSD matrix computed from the fourier spectra
%
% Revision 1.11  2005/05/19 16:15:19  roboos
% source does not always contain an avg, do conditional isfield for source.avg.csdlabel
%
% Revision 1.10  2005/05/17 17:58:28  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
% added preliminary implementation for postprocessing of pcc source csd matrices
%
% Revision 1.9  2005/03/03 13:57:29  roboos
% implemented pseudovalue resampling, i.e. a combination of jacknife and the complete average
%
% Revision 1.8  2004/10/27 13:13:02  roboos
% included source.trial.noise and source.trial.nai in the computation of avg, var and sem
%
% Revision 1.7  2004/10/13 15:54:05  roboos
% implemented trandformation on pow and noise (cfg.transform) prior to computing descriptive stats
%
% Revision 1.6  2004/09/22 14:09:52  roboos
% renamed method=randomized into method=randomization, with backward compatibility
%
% Revision 1.5  2004/09/01 17:59:28  roboos
% added copyright statements to all filed
% added cfg.version to all functions that give configuration in their output
% added cfg.previous to all functions with input data containing configuration details
%
% Revision 1.4  2004/08/04 16:10:41  roboos
% fixed bug in randomized data when different number of trials, added support for permuted resamplings
%
% Revision 1.3  2004/05/24 14:53:56  roberto
% fixed a typing mistake in the help
%
% Revision 1.2  2004/02/10 13:02:26  roberto
% fixed (obvious) bug due to not looping over trials
%
% Revision 1.1  2004/01/21 13:18:56  roberto
% new function that computes nai for all sourcereconstructions and
% var/avg/sem for sourcereconstructions that have been jacknifed or
% otherwise resampled over trials
%

% set the defaults
if ~isfield(cfg, 'transform'), cfg.transform = [];        end
if ~isfield(cfg, 'cohmethod'), cfg.cohmethod = 'dics';    end
if ~isfield(cfg, 'feedback'),  cfg.feedback  = 'textbar'; end

% 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';
end

if isfield(source, 'avg') && isfield(source.avg, 'csdlabel') 
  % it looks like a pcc beamformer was used
  Ndipole = length(source.inside) + length(source.outside);
  dipsel = match_str(source.avg.csdlabel, 'scandip');
  refchansel = match_str(source.avg.csdlabel, 'refchan');
  refdipsel  = match_str(source.avg.csdlabel, 'refdip');
  supchansel = match_str(source.avg.csdlabel, 'supchan');
  supdipsel  = match_str(source.avg.csdlabel, 'supdip');
  if length(refchansel)>0 && length(refdipsel)>0
    error('cannot simultaneously handle reference channels and reference dipole');
  elseif length(refchansel)>0
    refsel = refchansel;
  elseif length(refdipsel)>0
    refsel = refdipsel;
  else
    refsel = [];
  end
  if length(supchansel)>0 && length(supdipsel)>0
    error('cannot simultaneously handle suppressing channels and suppressing dipole');
  elseif length(supchansel)>0
    supsel = supchansel;
  elseif length(supdipsel)>0
    supsel = supdipsel;
  else
    supsel = [];
  end
  if ~isfield(source.avg, 'coh') && ~isempty(refsel) && strcmp(cfg.cohmethod, 'svdfft'),
    % FIXME this is not keeping track of the tapers!!!!!
    source.avg.pow = zeros(Ndipole, 1);
    source.avg.coh = zeros(Ndipole, 1);
    progress('init', cfg.feedback, 'performing svd on dipole moment');
    for diplop=1:length(source.inside)
      progress(diplop/length(source.inside), 'performing svd on dipole moment %d/%d\n', diplop, length(source.inside));
      i   = source.inside(diplop);
      mom = source.avg.mom{i}(dipsel,:);
      if length(dipsel>1)
        % compute the projection of the scanning dipole along the direction of the dominant amplitude
        mom = svdfft(mom, 1);
      end
      ref = source.avg.mom{i}(refsel,:);
      if length(refsel>1)
        % compute the projection of the reference dipole along the direction of the dominant amplitude
        ref = svdfft(ref, 1);
      end
      sup = source.avg.mom{i}(supsel,:);
      if length(supsel>1)
        % compute the projection of the superence dipole along the direction of the dominant amplitude
        sup = svdfft(sup, 1);
      end
      dat = [mom; ref; sup];
      csd = dat*ctranspose(dat);
      source.avg.csd{i} = csd; %this does not take the tapers into account
      source.avg.mom{i} = dat;
      if size(dat,1) == 2,
        source.avg.coh(i) = (csd(1,2).^2)./(abs(csd(1,1)).*abs(csd(2,2)));
      elseif size(dat,1) == 3,
        source.avg.coh(i)  = (csd(1,2).^2)./(abs(csd(1,1)).*abs(csd(2,2)));
        pcsd               = csd([1 2],[1 2]) - (csd([1 2],3)*inv(csd(3,3))*csd(3,[1 2]));
        source.avg.pcsd{i} = pcsd;
        source.avg.pcoh(i) = (pcsd(1,2).^2)./(abs(pcsd(1,1)).*abs(pcsd(2,2)));
      end
      source.avg.pow(i) = abs(csd(1,1));
    end
    progress('close');
    if length(refchansel)>0 && length(supchansel)>0,
      source.avg.csdlabel = {'scandip' 'refchan' 'supchan'};
    elseif length(refchansel)>0 && length(supdipsel)>0,
      source.avg.csdlabel = {'scandip' 'refchan' 'supdip'};
    elseif length(refchansel)>0,
      source.avg.csdlabel = {'scandip' 'refchan'};
    elseif length(refdipsel)>0 && length(supchansel)>0,
      source.avg.csdlabel = {'scandip' 'refdip' 'supchan'};
    elseif length(refdipsel)>0 && length(supdipsel)>0,
      source.avg.csdlabel = {'scandip' 'refdip' 'supdip'};
    elseif length(refdipsel)>0,
      source.avg.csdlabel = {'scandip' 'refdip'};
    end
    if isfield(source.avg, 'noisecsd'), source.avg = rmfield(source.avg, 'noisecsd'); end; 
  elseif ~isfield(source.avg, 'coh') && ~isempty(refsel) && strcmp(cfg.cohmethod, 'canonical'),
    fprintf('computing canonical coherence for pcc beamformer\n');
    % compute power and canonical coherence
    source.avg.pow = zeros(Ndipole, 1);
    source.avg.coh = zeros(Ndipole, 1);
    for i=source.inside(:)'
      source.avg.coh(i) = max(diag(cancorr(source.avg.csd{i}, dipsel, refsel)));
    end
    source.avg.pow(source.outside) = nan;
    source.avg.coh(source.outside) = nan;
  elseif ~isfield(source.avg, 'coh') && ~isempty(refsel) && strcmp(cfg.cohmethod, 'dics'),
    fprintf('computing power and coherence for pcc beamformer\n');
    % compute power and coherence on Joachim Gross' way
    source.avg.pow = zeros(Ndipole, 1);
    source.avg.coh = zeros(Ndipole, 1);
    for i=source.inside(:)'
      Pd  = lambda1(source.avg.csd{i}(dipsel,dipsel));
      Pr  = lambda1(source.avg.csd{i}(refsel,refsel));
      Cdr = lambda1(source.avg.csd{i}(dipsel,refsel));
      source.avg.pow(i) = Pd;
      source.avg.coh(i) = abs(Cdr).^2 ./ (Pd*Pr);
    end
    source.avg.pow(source.outside) = nan;
    source.avg.coh(source.outside) = nan;
  elseif ~isfield(source.avg, 'pow')
    fprintf('computing power for pcc beamformer\n');
    % compute power on Joachim Gross' way
    source.avg.pow = zeros(Ndipole, 1);
    for i=source.inside(:)'
      source.avg.pow(i) = lambda1(source.avg.csd{i});
    end
    source.avg.pow(source.outside) = nan;
  end
  if ~isfield(source.avg, 'noise') && isfield(source.avg, 'noisecsd')
    % compute power of projected noise
    fprintf('computing power of projected noise for pcc beamformer\n');
    source.avg.noise = zeros(Ndipole, 1);
    % compute projected noise on Joachim Gross' way
    for i=source.inside(:)'
      source.avg.noise(i) = lambda1(source.avg.noisecsd{i});
    end
    source.avg.noise(source.outside) = nan;
  end
end

if isfield(source, 'avg') && isfield(source.avg, 'pow') && isfield(source.avg, 'noise')
  % 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, 'jacknife') || 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')
    Ntrials = 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  )); end
  if isfield(dip(1), 'pow'),    sumdip.pow    = zeros(size(dip(1).pow  )); end
  if isfield(dip(1), 'coh'),    sumdip.coh    = zeros(size(dip(1).coh  )); end
  if isfield(dip(1), 'rv'),     sumdip.rv     = zeros(size(dip(1).rv   )); end
  if isfield(dip(1), 'noise'),  sumdip.noise  = zeros(size(dip(1).noise)); end
  if isfield(dip(1), 'nai'),    sumdip.nai    = zeros(size(dip(1).nai  )); end
  sqrdip = [];
  if isfield(dip(1), 'var'),    sqrdip.var    = zeros(size(dip(1).var  )); end
  if isfield(dip(1), 'pow'),    sqrdip.pow    = zeros(size(dip(1).pow  )); end
  if isfield(dip(1), 'coh'),    sqrdip.coh    = zeros(size(dip(1).coh  )); end
  if isfield(dip(1), 'rv'),     sqrdip.rv     = zeros(size(dip(1).rv   )); end
  if isfield(dip(1), 'noise'),  sqrdip.noise  = zeros(size(dip(1).noise)); end
  if isfield(dip(1), 'nai'),    sqrdip.nai    = zeros(size(dip(1).nai  )); 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 jacknife/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, 'jacknife')
    % Effron gives SEM estimate for the jacknife 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
    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
  else
    % remove the original trials
    source = rmfield(source, 'trial');
  end
  % assign the output source structure
  source.avg = dipmean;
  source.var = dipvar;
  source.sem = dipsem;
end

% add version information to the configuration
try
  % get the full name of the function
  cfg.version.name = mfilename('fullpath');
catch
  % required for compatibility with Matlab versions prior to release 13 (6.5)
  [st, i] = dbstack;
  cfg.version.name = st(i);
end
cfg.version.id = '$Id: sourcedescriptives.m,v 1.21 2006/01/30 17:10:33 roboos Exp $';
% remember the configuration details of the input data
try, cfg.previous = source.cfg; end
% remember the exact configuration details in the output 
source.cfg = cfg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);

