function [grandavg] = freqgrandaverage(cfg, varargin);

% FREQGRANDAVERAGE computes the average powerspectrum or time-frequency spectrum
% over multiple subjects
%
% Use as
%   [grandavg] = freqgrandaverage(cfg, freq1, freq2, freq3...)
%
% The input data freq1..N are obtained from either FREQANALYSIS with
% keeptrials=no or from FREQDESCRIPTIVES. The configuration structure
% can contain
%   cfg.keepindividual = 'yes' or 'no' (default)
%
% See also TIMELOCKGRANDAVERAGE, FREQANALYSIS, FREQDESCRIPTIVES

% Copyright (C) 2005, Robert Oostenveld
%
% $Log: freqgrandaverage.m,v $
% Revision 1.4  2006/02/23 10:28:16  roboos
% changed dimord strings for consistency, changed toi and foi into time and freq, added fixdimord where neccessary
%
% Revision 1.3  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.2  2005/05/17 17:50:37  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
%
% Revision 1.1  2005/01/03 12:23:17  roboos
% initial implementation
%

% set the defaults
if ~isfield(cfg, 'keepindividual'), cfg.keepindividual = 'no'; end

% for backward compatibility with old data structures
for i=1:length(varargin)
  varargin{i} = fixdimord(varargin{i});
end

% collect information about the input data, assume that they behave nicely
haspow   = isfield(varargin{1}, 'powspctrm');  % this is always true
hascoh   = isfield(varargin{1}, 'cohspctrm');
Nsubj    = length(varargin);
dimord   = varargin{1}.dimord;
dim      = size(varargin{1}.powspctrm);
if hascoh, 
  cohdim = size(varargin{1}.cohspctrm);
end
  
% for backward compatibility with old data structures
if isfield(varargin{1}, 'sgn') && ~isfield(varargin{1}, 'label')
  warning('renaming "sng" field into label for backward compatibility');
  for s=1:Nsubj
    varargin{s}.label = varargin{s}.sgn;
  end
end

% for backward compatibility with old data structures
if isfield(varargin{1}, 'sgncmb') && ~isfield(varargin{1}, 'labelcmb')
  warning('renaming "sngcmb" field into labelcmb for backward compatibility');
  for s=1:Nsubj
    varargin{s}.labelcmb = varargin{s}.sgncmb;
  end
end

% check whether the input data is suitable
if ~isempty(findstr(dimord, 'rpt'))
  error('the input data of each subject should be an average, use FREQDESCRIPTIVES first');
end
if ~isempty(findstr(dimord, 'tap'))
  error('multiple tapers in the input are not supported');
end

% give some feedback on the screen
if strcmp(cfg.keepindividual, 'no')
  if haspow, fprintf('computing average power over %d subjects\n', Nsubj); end
  if hascoh, fprintf('computing average coherence over %d subjects\n', Nsubj); end
else
  if haspow, fprintf('not computing grand average, but keeping individual power for %d subjects\n', Nsubj); end
  if hascoh, fprintf('not computing grand average, but keeping individual coherence for %d subjects\n', Nsubj); end
end

% allocate memory to hold the data
if strcmp(cfg.keepindividual, 'no')
  if haspow, s_pow = zeros(   dim); end
  if hascoh, s_coh = zeros(cohdim); end
else
  if haspow, s_pow = zeros([Nsubj    dim]); end
  if hascoh, s_coh = zeros([Nsubj cohdim]); end
end

for s=1:Nsubj
  if strcmp(cfg.keepindividual, 'no')
    % add this subject to the total sum
    if haspow, s_pow = s_pow + varargin{s}.powspctrm; end
    if hascoh, s_coh = s_coh + varargin{s}.cohspctrm; end
  else
    % concatenate this subject to the rest
    if haspow, s_pow(s,:) = varargin{s}.powspctrm(:); end
    if hascoh, s_coh(s,:) = varargin{s}.cohspctrm(:); end
  end
end

% collect the output data
grandavg.label = varargin{1}.label;
grandavg.freq  = varargin{1}.freq;
if isfield(varargin{1}, 'time')
  % remember the time axis
  grandavg.time = varargin{1}.time;
end
if isfield(varargin{1}, 'grad')
  % this is not averaged
  warning('discarding gradiometer information');
end
if isfield(varargin{1}, 'elec')
  % this is not averaged
  warning('discarding electrode information');
end

if strcmp(cfg.keepindividual, 'no')
  grandavg.dimord = dimord;
  grandavg.powspctrm = s_pow/Nsubj;
  if hascoh,
    grandavg.labelcmb = varargin{1}.labelcmb;
    grandavg.cohspctrm = s_coh/Nsubj;
  end
else
  grandavg.dimord = ['subj_' dimord];
  grandavg.powspctrm = s_pow;
  if hascoh,
    grandavg.labelcmb  = varargin{1}.labelcmb;
    grandavg.cohspctrm = s_coh;
  end
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: freqgrandaverage.m,v 1.4 2006/02/23 10:28:16 roboos Exp $';
% remember the configuration details of the input data
cfg.previous = [];
for i=1:length(varargin)
  try, cfg.previous{i} = varargin{i}.cfg; end
end
% remember the exact configuration details in the output 
grandavg.cfg = cfg;

