function [output] = freqdescriptives(cfg, freq)

% FREQDESCRIPTIVES computes descriptive statistics of the frequency
% or time-frequency decomposition of the EEG/MEG signal, such as
% average power and coherence.
%
% Use as
%   [freq] = frqanastatistic(cfg, freq)
% where cfg can be empty.
%
% The input freq should be organised in a structure as obtained from
% from the FREQANALYSIS function.
%
% The output depends on the input:
%
% If the input contains only power-spectra, then only the statistics of the power is computed.
% If the input contains power- and cross-spectral densities, then power and coherence are computed.
%
% If results from trials and/or tapers have not been kept, only coherence will be calculated.
% If results from trials and/or tapers have been kept, the jackknife estimate of the SEM of
% power and eventually coherence are calculated. In that case, also the degrees of freedom (dof)
% are computed.
%
% See also FREQANALYSIS, FREQSTATISTICS

% Copyright (C) 2004, F.C. Donders Centre
%
% $Log: freqdescriptives.m,v $
% Revision 1.27  2006/02/27 13:39:03  roboos
% fixed some dimord problems
%
% Revision 1.26  2006/02/27 11:22:45  roboos
% fixed incorrect variable name, data should be freq
%
% Revision 1.25  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.24  2006/02/09 10:06:27  roboos
% changed a help comment
%
% Revision 1.23  2006/02/01 12:26:00  roboos
% made all uses of dimord consistent with the common definition of data dimensions, see the fixdimord() function
%
% Revision 1.22  2005/09/05 06:37:51  jansch
% small change in feedback
%
% Revision 1.21  2005/08/16 11:19:04  jansch
% moved some stuff in 'ro's implementation to a separate function
% fourier2crsspctrm to facilitate re-use
%
% Revision 1.20  2005/08/09 13:43:16  jansch
% made some changes in the handling of fourier-input in the ro-implementation.
%
% Revision 1.19  2005/08/08 10:36:03  jansch
% included the beta-constant of 23/20 in the z-transformation of the coherence
%
% Revision 1.18  2005/08/08 09:44:58  jansch
% fixed some small bugs in ro-implementation, added hidden option of
% cfg.labelcmb in the case of fft-input.
%
% Revision 1.17  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.16  2005/05/27 09:40:40  jansch
% fixed bug in combinechan
%
% Revision 1.15  2005/05/18 12:46:02  jansch
% fixed bug in correct conjugation of cross-spectral densities
%
% Revision 1.14  2005/05/17 17:53:18  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.13  2005/05/17 09:31:10  jansch
% fixed bug in plv
%
% Revision 1.12  2005/05/13 08:15:07  jansch
% fixed some small bugs
%
% Revision 1.11  2005/05/13 08:00:21  jansch
% included jacknife, pseudovalue, and z-transformation. moved partialisation
% part to a subfunction
%
% Revision 1.10  2005/05/09 11:38:43  jansch
% improved ro-implementation
%
% Revision 1.9  2005/04/13 08:52:25  roboos
% completely new implementation, separate from old code by Pascal
% the two implementations co-exists as subfunctions  with cfg.implementation
% the new implementation is not completely ready yet, but is sufficient for testing by some users
% new functionality is canonical coherence, phase locking value, redesign of jacknifing, complex or absolute output, coherence between planar gradiometers
%
% Revision 1.8  2005/04/12 09:43:12  roboos
% added many comments, improved the help
% no functionality changes, although some lines of code have been rearranged
%
% Revision 1.7  2005/01/27 08:47:21  roboos
% fixed bug in finding channel indices  (added 'exact' to strmatch)
% added version details to the configuration
%
% Revision 1.6  2005/01/20 11:03:34  jansch
% changed sgncmb into labelcmb, allowing for backwardcompatibility
%
% Revision 1.5  2004/10/08 10:17:16  jansch
% *** empty log message ***
%
% Revision 1.4  2004/10/08 06:50:18  jansch
% fixed small bug which caused a crash when cfg.keeptrials has been 'no'
%

% set the defaults
if ~isfield(cfg, 'implementation'), cfg.implementation = 'pf'; end  % pf, ro

% for backward compatibility with old data structures
freq = fixdimord(freq);

% call the appropriate subfunction
if strcmp(cfg.implementation, 'pf')
  [output] = freqdescriptives_pf(cfg, freq);
elseif strcmp(cfg.implementation, 'ro')
  [output] = freqdescriptives_ro(cfg, freq);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION with Pascal Fries' original implementation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [output] = freqdescriptives_pf(cfg, freq);

csdflg = isfield(freq, 'labelcmb');       % boolean, indicates presence of cross-spectra
rptflg = strcmp(freq.dimord(1:3),'rpt');  % boolean, indicates presence of single trials

if csdflg
  % find the indices into the power spectra for each cross-spectra, i.e. 2 elements
  for labelcmblop = 1:size(freq.labelcmb,1)
    labelcmbind(labelcmblop,1) = strmatch(freq.labelcmb(labelcmblop,1), freq.label, 'exact');
    labelcmbind(labelcmblop,2) = strmatch(freq.labelcmb(labelcmblop,2), freq.label, 'exact');
  end
end

if ~rptflg
  % there are no individual trials in the input, only compute average
  newdimord = freq.dimord;
  if ~csdflg
    disp('Input was averaged PSD only, allowing no statistic.')
  elseif csdflg
    if strcmp(freq.dimord(end-3:end),'freq')
      output.coh = abs(freq.crsspctrm ./ ...
        sqrt(freq.powspctrm(labelcmbind(:,1),:) .* ...
        freq.powspctrm(labelcmbind(:,2),:)));
    elseif strcmp(freq.dimord(end-3:end),'time')
      output.coh = abs(freq.crsspctrm ./ ...
        sqrt(freq.powspctrm(labelcmbind(:,1),:,:) .* ...
        freq.powspctrm(labelcmbind(:,2),:,:)));
    end
  end
  output.pow = freq.powspctrm;
elseif rptflg
  % there are individual trials in the input, compute jacknife variance estimates
  newdimord = freq.dimord(5:end);
  if strcmp(freq.dimord(end-3:end),'freq')
    siz = size(freq.powspctrm);
    numrpt       = siz(1);
    output.dof   = numrpt;
    numlabel     = siz(2);
    numboi       = siz(3);
    jckavgpowarr = zeros(numrpt,numlabel,numboi);
    sumpowspctrm = sum(freq.powspctrm,1);
    for rptlop = 1:numrpt
      jckavgpowarr(rptlop,:,:) = (sumpowspctrm - freq.powspctrm(rptlop,:,:)) ./ (numrpt -1);
    end
    avgjckavgpow = mean(jckavgpowarr,1);
    jckavgpowdum = zeros(1,numlabel,numboi);
    for rptlop = 1:numrpt
      jckavgpowdum = jckavgpowdum + ((jckavgpowarr(rptlop,:,:) - avgjckavgpow) .^ 2);
    end
    jckavgpowsem = sqrt((numrpt-1) .* jckavgpowdum ./ numrpt);
    avgpow = sumpowspctrm ./ numrpt;
    output.pow = reshape(avgpow,[numlabel,numboi]);
    output.jckpowsem = reshape(jckavgpowsem,[numlabel,numboi]);
    if csdflg
      numlabelcmb   = size(freq.labelcmb,1);
      jckcoharr = zeros(numrpt,numlabelcmb,numboi);
      sumcrsspctrm = sum(freq.crsspctrm,1);
      for rptlop = 1:numrpt
        jckcoharr(rptlop,:,:) = abs((sumcrsspctrm - freq.crsspctrm(rptlop,:,:)) ./ ...
          sqrt((sumpowspctrm(1,labelcmbind(:,1),:) - freq.powspctrm(rptlop,labelcmbind(:,1),:)) .* ...
          (sumpowspctrm(1,labelcmbind(:,2),:) - freq.powspctrm(rptlop,labelcmbind(:,2),:))));
      end
      avgjckcoh = mean(jckcoharr,1);
      jckcohdum = zeros(1,numlabelcmb,numboi);
      for rptlop = 1:numrpt
        jckcohdum = jckcohdum + ((jckcoharr(rptlop,:,:) - avgjckcoh) .^ 2);
      end
      jckcohsem = sqrt((numrpt-1) .* jckcohdum ./ numrpt);
      coh = abs(sumcrsspctrm ./ ...
        sqrt(sumpowspctrm(:,labelcmbind(:,1),:) .* ...
        sumpowspctrm(:,labelcmbind(:,2),:)));
      output.coh = reshape(coh,[numlabelcmb,numboi]);
      output.jckcohsem = reshape(jckcohsem,[numlabelcmb,numboi]);
    end
  elseif strcmp(freq.dimord(end-3:end),'time')
    siz = size(freq.powspctrm);
    numrpt       = siz(1);
    numlabel     = siz(2);
    numfoi       = siz(3);
    numtoi       = siz(4);
    output.dof   = zeros(numfoi,numtoi);
    jckavgpowarr = zeros(numrpt,numlabel,numfoi,numtoi);
    for toilop = 1:numtoi
      for foilop = 1:numfoi
        finiteind = find(isfinite(freq.powspctrm(:,1,foilop,toilop)));
        numfiniteind = length(finiteind);
        output.dof(foilop,toilop) = numfiniteind;
        if numfiniteind < 3
          jckavgpowarr(:,:,foilop,toilop) = nan;
        else
          sumpowspctrm = sum(freq.powspctrm(finiteind,:,foilop,toilop),1);
          for rptlop = 1:length(finiteind)
            jckavgpowarr(rptlop,:,foilop,toilop) = (sumpowspctrm - ...
              freq.powspctrm(finiteind(rptlop),:,foilop,toilop)) ./ (numfiniteind -1);
          end
          jckavgpowarr(numfiniteind+1:end,:,foilop,toilop) = nan;
        end
      end
    end
    jckavgpowdum = zeros(1,numlabel,numfoi,numtoi);
    for toilop = 1:numtoi
      for foilop = 1:numfoi
        finiteind = find(isfinite(jckavgpowarr(:,1,foilop,toilop)));
        numfiniteind = length(finiteind);
        if numfiniteind < 3
          jckavgpowdum(1,:,foilop,toilop) = nan;
        else
          avgjckavgpow = mean(jckavgpowarr(finiteind,:,foilop,toilop),1);
          for rptlop = 1:numfiniteind
            jckavgpowdum(1,:,foilop,toilop) = jckavgpowdum(1,:,foilop,toilop) + ...
              ((jckavgpowarr(finiteind(rptlop),:,foilop,toilop) - avgjckavgpow) .^ 2);
          end
        end
      end
    end
    jckavgpowsem = sqrt((numrpt-1) .* jckavgpowdum ./ numrpt);
    avgpow = zeros(1,numlabel,numfoi,numtoi);
    for toilop = 1:numtoi
      for foilop = 1:numfoi
        finiteind = find(isfinite(freq.powspctrm(:,1,foilop,toilop)));
        numfiniteind = length(finiteind);
        if numfiniteind < 3
          avgpow(:,:,foilop,toilop) = nan;
        else
          avgpow(1,:,foilop,toilop) = mean(freq.powspctrm(finiteind,:,foilop,toilop),1);
        end
      end
    end
    output.pow = reshape(avgpow,[numlabel,numfoi,numtoi]);
    output.jckpowsem = reshape(jckavgpowsem,[numlabel,numfoi,numtoi]);
    if csdflg
      numlabelcmb   = size(freq.labelcmb,1);
      jckcoharr = zeros(numrpt,numlabelcmb,numfoi,numtoi);
      for toilop = 1:numtoi
        for foilop = 1:numfoi
          finiteind = find(isfinite(freq.crsspctrm(:,1,foilop,toilop)));
          if length(finiteind) < 3
            jckcoharr(:,:,foilop,toilop) = nan;
          else
            sumcrsspctrm = sum(freq.crsspctrm(finiteind,:,foilop,toilop),1);
            sumpowspctrm = sum(freq.powspctrm(finiteind,:,foilop,toilop),1);
            for rptlop = 1:length(finiteind)
              jckcoharr(rptlop,:,foilop,toilop) = abs(...
                (sumcrsspctrm - freq.crsspctrm(finiteind(rptlop),:,foilop,toilop)) ./ ...
                sqrt((sumpowspctrm(1,labelcmbind(:,1)) - ...
                freq.powspctrm(finiteind(rptlop),labelcmbind(:,1),foilop,toilop)) .* ...
                (sumpowspctrm(1,labelcmbind(:,2)) - ...
                freq.powspctrm(finiteind(rptlop),labelcmbind(:,2),foilop,toilop))));
            end
            jckcoharr(length(finiteind)+1:end,:,foilop,toilop) = nan;
          end
        end
      end
      jckcohdum = zeros(1,numlabelcmb,numfoi,numtoi);
      for toilop = 1:numtoi
        for foilop = 1:numfoi
          finiteind = find(isfinite(jckcoharr(:,1,foilop,toilop)));
          if length(finiteind) < 3
            jckcohdum(1,:,foilop,toilop) = nan;
          else
            avgjckcoh = mean(jckcoharr(finiteind,:,foilop,toilop),1);
            for rptlop = 1:length(finiteind)
              jckcohdum(1,:,foilop,toilop) = jckcohdum(1,:,foilop,toilop) + ...
                ((jckcoharr(finiteind(rptlop),:,foilop,toilop) - avgjckcoh) .^ 2);
            end
          end
        end
      end
      jckcohsem = sqrt((numrpt-1) .* jckcohdum ./ numrpt);
      coh = zeros(1,numlabelcmb,numfoi,numtoi);
      for toilop = 1:numtoi
        for foilop = 1:numfoi
          finiteind = find(isfinite(freq.crsspctrm(:,1,foilop,toilop)));
          if length(finiteind) < 3
            coh(:,:,foilop,toilop) = nan;
          else
            sumcrsspctrm = sum(freq.crsspctrm(finiteind,:,foilop,toilop),1);
            sumpowspctrm = sum(freq.powspctrm(finiteind,:,foilop,toilop),1);
            coh(1,:,foilop,toilop) = abs(sumcrsspctrm ./ ...
              sqrt((sumpowspctrm(1,labelcmbind(:,1))) .* ...
              (sumpowspctrm(1,labelcmbind(:,2)))));
          end
        end
      end
      output.coh = reshape(coh,[numlabelcmb,numfoi,numtoi]);
      output.jckcohsem = reshape(jckcohsem,[numlabelcmb,numfoi,numtoi]);
    end
  end
end

output.powspctrm = output.pow;
output = rmfield(output,'pow');

if rptflg
  % rename the output from jckpowsem to powspctrmsem
  output.powspctrmsem = output.jckpowsem;
  output = rmfield(output,'jckpowsem');
end

if csdflg
  % rename the output from coh to cohspctrm
  output.cohspctrm = output.coh;
  output = rmfield(output,'coh');
  output.labelcmb = freq.labelcmb;
  if rptflg
    % rename the output from jckcohsem to cohspctrmsem
    output.cohspctrmsem = output.jckcohsem;
    output = rmfield(output,'jckcohsem');
  end
end

% collect everything into a self-descriptive output structure
output.label    = freq.label;
output.dimord   = newdimord;
output.freq     = freq.freq;
if strcmp(freq.dimord(end-3:end),'time')
  output.time    = freq.time;
end
if isfield(freq, 'grad')
  output.grad      = freq.grad;
end
if ~rptflg && isfield(freq,'dof'),
  % copy the degrees of freedom from the input
  output.dof = freq.dof;
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: freqdescriptives.m,v 1.27 2006/02/27 13:39:03 roboos Exp $';
% remember the configuration details of the input data
try, cfg.previous = freq.cfg; end
% remember the exact configuration details in the output 
output.cfg = cfg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION with Robert Oostenveld's implementation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [output] = freqdescriptives_ro(cfg, freq);

% TODO:
%  jacknife in output
%  output for absangle and realimag
%  pseudovalue
%  coh and plv using svd
%  partialisation
%  replace squared coherence by non-squared coherence

% set the defaults
if ~isfield(cfg, 'feedback'),          cfg.feedback    = 'textbar';    end
if ~isfield(cfg, 'cohmethod'),         cfg.cohmethod   = 'canonical';  end    % svd, canonical, plain
if ~isfield(cfg, 'coh_ztransform'),    cfg.coh_ztransform = 'no';      end
if ~isfield(cfg, 'plv'),               cfg.plv         = 'no';         end
if ~isfield(cfg, 'plvmethod'),         cfg.plvmethod   = 'plain';      end    % svd, canonical, plain
if ~isfield(cfg, 'complex'),           cfg.complex     = 'no';         end    % no, yes, absangle, realimag
if ~isfield(cfg, 'combinechan'),       cfg.combinechan = 'no';         end    % no, planar, ...
if ~isfield(cfg, 'jacknife'),          cfg.jacknife    = 'no';         end
if ~isfield(cfg, 'pseudovalue'),       cfg.pseudovalue = 'no';         end    % yes, no
if ~isfield(cfg, 'variance'),          cfg.variance    = 'no';         end %FIXME, include sanity check on config.
if ~isfield(cfg, 'partchan'),          cfg.partchan    = 'no';         end

output   = [];
hascsd   = isfield(freq, 'crsspctrm');
hasdft   = isfield(freq, 'fourierspctrm');
haspow   = isfield(freq, 'powspctrm');
hasrpt   = ~isempty(strfind(freq.dimord, 'rpt'));
hastim   = ~isempty(strfind(freq.dimord, 'time'));
jckflg   = strcmp(cfg.jacknife,    'yes');
psdflg   = strcmp(cfg.pseudovalue, 'yes');
varflg   = strcmp(cfg.variance,    'yes');
plvflg   = strcmp(cfg.plv,         'yes');
 
if ~hasrpt && (strcmp(cfg.jacknife, 'yes') || strcmp(cfg.pseudovalue, 'yes')),
  error('a jacknife-estimate or pseudovalue-estimate without repeated observations in the input is not possible');
end

if ~(hascsd || hasdft) && ~strcmp(cfg.partchan, 'no'),
  error('cannot do partialisation without cross-spectra in the input-data');
end

if ~isfield(freq, 'cumtapcnt'),
  freq = freq2cumtapcnt(freq, cfg.fsample);
end

% determine the labels of the input data
for i=1:length(freq.label)
  inputindex{i} = i;
end

%put the powerspectra and cross-spectra into one matrix,
if ~hascsd, 
  freq.crsspctrm = []; 
  freq.labelcmb  = {};
end 

if hasdft,
  if ~hasrpt,
    error('not yet implemented without having kept trials');
  else
    if isfield(cfg, 'channel'),    config.channel    = cfg.channel;    end;
    if isfield(cfg, 'channelcmb'), config.channelcmb = cfg.channelcmb; end;
    if isfield(cfg, 'foilim'),     config.foilim     = cfg.foilim;     end;
    if isfield(cfg, 'keeptrials'), 
      config.keeptrials = cfg.keeptrials; 
    elseif ~(jckflg || psdflg || varflg || plvflg), 
      config.keeptrials = 'no'; 
    else 
      config.keeptrials = 'yes';
    end
    if isfield(cfg, 'crsfrq') && strcmp(cfg.crsfrq, 'yes'),
      error('not yet implemented');
    else
      fprintf('computing cross-spectra from fourier-input\n');
      freq   = fourier2crsspctrm(config, freq);
      hasrpt = strcmp(config.keeptrials, 'yes');   
      hascsd = 1;
    end
  end
end

if haspow, %concatenate powerspectrum and cross-spectrum
           %eventually this will be done at a different place
   if ~hasrpt, 
      freq.crsspctrm = cat(1,freq.powspctrm,freq.crsspctrm);
      freq.labelcmb  = [freq.label freq.label; freq.labelcmb];
   else
      freq.crsspctrm = cat(2,freq.powspctrm,freq.crsspctrm);
      freq.labelcmb  = [freq.label freq.label; freq.labelcmb];
   end
end

if ~hasrpt,
    dum(1,:,:,:)   = freq.crsspctrm;
    freq.crsspctrm = dum;
    Nrpt           = 1;
    Nsgn           = length(freq.label);
    Nfrq           = size(freq.crsspctrm,3);
    Ntim           = size(freq.crsspctrm,4);
else
    Nrpt           = size(freq.crsspctrm,1);
    Nsgn           = length(freq.label);
    Nfrq           = size(freq.crsspctrm,3);
    Ntim           = size(freq.crsspctrm,4);
end

if plvflg,
  freq.ncrsspctrm = freq.crsspctrm./abs(freq.crsspctrm);
end

if jckflg || psdflg, %create the cross-spectral density matrix containing the leave-one-out
                    %averages, with the average of all trials appended to it.
  sumcrsspctrm = sum(freq.crsspctrm,1);
  sumtap       = sum(freq.cumtapcnt,1);
  if plvflg, sumncrsspctrm = sum(freq.ncrsspctrm,1); end;
  progress('init', cfg.feedback, 'computing the leave-one-out averages');
  for j = 1:Nrpt %do it the inefficient way, save memory, to be sure
     progress(j/Nrpt, 'computing coherence %d/%d\n', j, Nrpt);
     freq.crsspctrm(j,:,:,:) = (sumcrsspctrm - freq.crsspctrm(j,:,:,:))./(Nrpt - 1);  
     freq.cumtapcnt(j,:)     = sumtap - freq.cumtapcnt(j,:);
     if plvflg, freq.ncrsspctrm(j,:,:,:) = (sumncrsspctrm - freq.ncrsspctrm(j,:,:,:))./(Nrpt - 1); end; 
  end
  freq.crsspctrm(end+1,:,:,:) = sumcrsspctrm./Nrpt;
  freq.cumtapcnt(end+1,:)     = sumtap;
  if plvflg, freq.ncrsspctrm(end+1,:,:,:) = sumncrsspctrm./Nrpt; end;
  progress('close');
elseif varflg,
  %do nothing
  warning('variance computation only interpretable for power-estimate');
else
  %freq.crsspctrm = mean(freq.crsspctrm,1);
  %freq.cumtapcnt = sum(freq.cumtapcnt,1);
  %FIXME this leads to crashes if there is only one freq-bin
  %quick and dirty fix
  if Nfrq == 1,
    dum = mean(freq.crsspctrm, 1);
    freq.crsspctrm = cat(3, dum, dum);
    freq.cumtapcnt = sum(freq.cumtapcnt,1);
    Nfrq = 2;
    freq.freq = [freq.freq freq.freq];
  else
    freq.crsspctrm = mean(freq.crsspctrm,1);
    freq.cumtapcnt = sum(freq.cumtapcnt,1);
  end
 
  if plvflg, freq.ncrsspctrm = mean(freq.ncrsspctrm,1); end
end


%if hasrpt && ~(jckflg || psdflg),
%  freq.crsspctrm = mean(freq.crsspctrm(1:end-1,:,:,:),1);
%  if plvflg, freq.ncrsspctrm = mean(freq.ncrsspctrm(1:end-1,:,:,:),1); end
if ~hasrpt,
  freq.crsspctrm(1,:,:,:) = freq.crsspctrm;
  if plvflg, freq.ncrsspctrm(1,:,:,:) = freq.ncrsspctrm; end
end

%---------------------------- 
%partialise the cross-spectra

if ~strcmp(cfg.partchan, 'no'),
  [freq] = partialisation(cfg, freq, 'crsspctrm');
  if plvflg,
    [freq] = partialisation(cfg, freq, 'ncrsspctrm');
  end 
end

%end of partialisation procedure
%-------------------------------

%-----------------------------
%recombination of the channels

[outputlabel, outputindex] = inputlabel2outputlabel(cfg, freq);

%end of recombination of the channels
%------------------------------------

%create signalcombination-matrix, each non-nan entry corresponding
%with the index of the signalcombination in the cross-spectrum
%negative entries should be conjugated before further processing.
fprintf('indexing the channels and channel combinations\n');
sgnindx    = [1:length(freq.label)]';
if hascsd,
  sgncmbindx = zeros(size(freq.labelcmb));
  for i=1:size(freq.labelcmb,1)
    sgncmbindx(i,1) = strmatch(freq.labelcmb{i,1}, freq.label, 'exact');
    sgncmbindx(i,2) = strmatch(freq.labelcmb{i,2}, freq.label, 'exact');
  end
  sgncmbmat  = nan*zeros(length(sgnindx),length(sgnindx));
  for i=1:size(sgncmbindx,1)
    sgncmbmat(sgncmbindx(i,2),sgncmbindx(i,1)) = -i;
    sgncmbmat(sgncmbindx(i,1),sgncmbindx(i,2)) = i;
  end
end

%--------------------------------------------
%compute the recombined and partialised power

pow    = zeros(size(freq.crsspctrm,1),length(outputlabel), Nfrq, Ntim);

progress('init', cfg.feedback, 'computing the combined power in the required output-channels');
for i = 1:length(outputlabel)
  progress(i/length(outputlabel), 'computing power %d/%d\n', i, length(outputlabel));
  for j = 1:length(outputindex{i})
    indx{i}(j) = find(sgncmbindx(:,1)==outputindex{i}(j) & sgncmbindx(:,2)==outputindex{i}(j));
    if isempty(indx{i}(j)),
      indx{i}(j) = find(sgncmbindx(:,2)==outputindex{i}(j) & sgncmbindx(:,1)==outputindex{i}(j));
    end; %find the auto-spectra
  end
  pow(:,i,:,:) = abs(sum(freq.crsspctrm(:,indx{i},:,:),2));
end
progress('close');

freq.powspctrm = pow;  %this now contains the combined power-spectra

%end of power combination
%------------------------

%---------------------
%compute the coherence

if hascsd,
  progress('init', cfg.feedback, 'computing the rearrangement of the cross-spectra');
  m   = 0;
  k   = 0;
  cmb = {};
  for i=1:length(outputlabel)
    for j=i:length(outputlabel)
      m = m + 1;
      cmbtmp1 = [];
      cmbtmp2 = [];
      progress(m/((length(outputlabel).^2)./2));
      if i==j
        continue
      end
      % construct the complete CSD matrix for this combination of channels
      sel1 = outputindex{i};
      sel2 = outputindex{j};
      if isempty(sel1) || isempty(sel2)
        continue
      end
      if ~isempty(intersect(sel1,sel2)),%computing multiple coherence between sets with overlapping sensors will give 1.
        continue
      end
      if strcmp(cfg.combinechan(1), 'pseudomeg') && length(sel1)==1 && length(sel2)==1,
        continue
      end
      sel    = [sel1 sel2];
      dum    = sgncmbmat(sel,sel);
      if any(isnan(dum))
        continue
      else
         k = k + 1;
         cmb{k} = dum(:);
         cmbind{k,1} = sel1;
         cmbind{k,2} = sel2;
         outputlabelcmb{k,1} = outputlabel{i};
         outputlabelcmb{k,2} = outputlabel{j};
      end
    end
  end
  progress('close');
 
  Nrpt           = size(freq.crsspctrm,1); 
  Ncmb           = length(cmb);
  coh            = complex(zeros(Nrpt,Ncmb,Nfrq,Ntim),zeros(Nrpt,Ncmb,Nfrq,Ntim));
  if plvflg, plv = complex(zeros(Nrpt,Ncmb,Nfrq,Ntim),zeros(Nrpt,Ncmb,Nfrq,Ntim)); end

  cbin = 1; 
  progress('init', cfg.feedback, 'computing the coherence');
  for rbin=1:Nrpt 
    for cbin=1:Ncmb
      progress((Ncmb*(rbin-1)+cbin)./(Ncmb*Nrpt));
      nsgn = sqrt(length(cmb{cbin}));
      csd  = squeeze(freq.crsspctrm(rbin,abs(cmb{cbin}),:,:));
      csd(find(cmb{cbin}<0),:,:) = conj(csd(find(cmb{cbin}<0),:,:));
      csd    = reshape(csd, [nsgn,nsgn,Nfrq,Ntim]);
      if plvflg,
        ncsd = squeeze(freq.ncrsspctrm(rbin,abs(cmb{cbin}),:,:));
        ncsd(find(cmb{cbin}<0),:,:) = conj(ncsd(find(cmb{cbin}<0),:,:));
        ncsd = reshape(ncsd, [nsgn,nsgn,Nfrq,Ntim]);
      end
      for fbin=1:Nfrq
        for tbin=1:Ntim
          n1 = length(cmbind{cbin,1});
          n2 = length(cmbind{cbin,2});
          coh(rbin,cbin,fbin,tbin)       = coherence(squeeze(csd(:,:,fbin,tbin)), n1, n2, cfg.cohmethod);
          %if plvflg, plv(rbin,cbin,fbin,tbin) = phaselocking(squeeze(ncsd(:,:,fbin,tbin)), n1, n2, cfg.plvmethod); end
          if plvflg, plv(rbin,cbin,fbin,tbin) = coherence(squeeze(ncsd(:,:,fbin,tbin)), n1, n2, cfg.plvmethod); end
        end
      end
    end 
  end
  progress('close');
  
  if ~strcmp(cfg.coh_ztransform, 'no'),
    zcoh = zeros(size(coh));
    if plvflg, zplv = zeros(size(plv)); end 
    for rbin = 1:size(coh,1)
       for cbin = 1:size(coh,2)
                    zcoh(rbin,cbin,:,:) = zcoherence(squeeze(abs(coh(rbin,cbin,:,:))), freq.cumtapcnt(rbin,:).*2, cfg.coh_ztransform);
         if plvflg, zplv(rbin,cbin,:,:) = zcoherence(squeeze(abs(plv(rbin,cbin,:,:))), freq.cumtapcnt(rbin,:).*2, cfg.coh_ztransform); end
       end
    end
  end
end

%end of coherence computation
%----------------------------

%--------------------------------------
%put output coherence in correct format

if hascsd,
   if exist('zcoh', 'var'),
      if jckflg || psdflg, jckzcoh = zcoh(1:end-1,:,:,:); end
      zcoh = zcoh(end,:,:,:);
   end

  switch cfg.complex
    case 'yes'
      % leave all values complex
      if jckflg || psdflg, jckcoh = coh(1:end-1,:,:,:); end   
      coh = coh(end,:,:,:);
      if plvflg,
        if jckflg || psdflg,
          jckplv = plv(1:end-1,:,:,:);
        end
        plv = plv(end,:,:,:); 
      end
    case 'no'
      % convert to absolute values
      if jckflg || psdflg, jckcoh = abs(coh(1:end-1,:,:,:)); end
      coh = abs(coh(end,:,:,:));
      if plvflg,      
        if jckflg || psdflg,
          jckplv = abs(plv(1:end-1,:,:,:));
        end
        plv = abs(plv(end,:,:,:));
      end
    case 'realimag'
      % convert to real values
      if jckflg || psdflg, jckcoh_real = real(coh(1:end-1,:,:,:));
                          jckcoh_imag = imag(coh(1:end-1,:,:,:)); end
      coh_real = real(coh(end,:,:,:));
      coh_imag = imag(coh(end,:,:,:));
      
      if plvflg,
        if jckflg || psdflg,
          jckplv_real = real(plv(1:end-1,:,:,:));
          jckplv_imag = imag(plv(1:end-1,:,:,:));
        end
        plv_real = real(plv(end,:,:,:));
        plv_imag = imag(plv(end,:,:,:));
      end
    case 'absangle'
      % convert to absolute values
      if jckflg || psdflg, jckcoh_abs = abs(coh(1:end-1,:,:,:)); 
                          jckcoh_angle = angle(coh(1:end-1,:,:,:)); end
      coh_abs   = abs(coh(end,:,:,:));
      coh_angle = angle(coh(end,:,:,:));
      if plvflg,
        if jckflg || psdflg,
          jckplv_abs   = abs(plv(1:end-1,:,:,:));
          jckplv_angle = angle(plv(1:end-1,:,:,:));
        end
        plv_abs   = abs(plv(end,:,:,:));
        plv_angle = angle(plv(end,:,:,:));
      end 
    otherwise
      error(sprintf('method ''%s'' is not implemented for cfg.complex', cfg.complex));
  end
end

%---------------------------
%create the output-structure

output.avg.powspctrm   = squeeze(mean(freq.powspctrm,1));
output.freq            = freq.freq;
if hastim, output.time = freq.time; end
output.label           = outputlabel(:);

if varflg,
  output.var.powspctrm = squeeze( std(freq.powspctrm,[],1)).^2;
end

if jckflg,
  output.jckavg.powspctrm = squeeze(mean(freq.powspctrm));
  output.jckvar.powspctrm = squeeze( std(freq.powspctrm,[],1)).^2;
end

if psdflg,
  Nrpt = size(pow,1) - 1;
  output.pseudo.powspctrm = repmat(Nrpt.*pow(end,:,:,:), [Nrpt 1 1 1]) - (Nrpt-1).*pow(1:end-1,:,:,:);
end

if exist('coh','var'),
  output.avg.cohspctrm = squeeze(coh);
  if jckflg,
    output.jckavg.cohspctrm = squeeze(mean(jckcoh));
    output.jckvar.cohspctrm = squeeze( std(jckcoh,[],1).^2);
  end

  if psdflg,
    output.pseudo.cohspctrm = repmat(Nrpt.*coh, [Nrpt 1 1 1]) - (Nrpt-1).*jckcoh;
  end
end

if exist('zcoh','var'),
  output.avg.zcohspctrm = squeeze(zcoh);
  if jckflg,
    output.jckavg.zcohspctrm = squeeze(mean(jckzcoh));
    output.jckvar.zcohspctrm = squeeze( std(jckzcoh,[],1).^2);
  end
  
  if psdflg,
    output.pseudo.zcohspctrm = repmat(Nrpt.*zcoh, [Nrpt 1 1 1]) - (Nrpt-1).*jckzcoh;
  end
end

if exist('plv','var'),
  output.avg.plvspctrm = squeeze(plv);
  if jckflg,
    output.jckavg.plvspctrm = squeeze(mean(jckplv));
    output.jckvar.plvspctrm = squeeze( std(jckplv,[],1).^2);
  end

  if psdflg,
    output.pseudo.plvspctrm = repmat(Nrpt.*plv, [Nrpt 1 1 1]) - (Nrpt-1).*jckplv;
  end
end        % FIXME, bias correction needed

if exist('outputlabelcmb','var'), output.labelcmb = outputlabelcmb; 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: freqdescriptives.m,v 1.27 2006/02/27 13:39:03 roboos Exp $';
try, cfg.previous = freq.cfg; end

% remember the configuration details
output.cfg = cfg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that partialises cross-spectra
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [freq] = partialisation(cfg, freq, field)

crsspctrm = getfield(freq, field);
Nrpt      = size(crsspctrm,1);
Nsgn      = size(crsspctrm,2);
Nfrq      = size(crsspctrm,3);
Ntim      = size(crsspctrm,4);
%create signalcombination-matrix, each non-nan entry corresponding
%with the index of the signalcombination in the cross-spectrum
%negative entries should be conjugated before further processing.
fprintf('indexing the channels and channel combinations\n');
sgnindx    = [1:length(freq.label)]';
%if hascsd,
  sgncmbindx = zeros(size(freq.labelcmb));
  for i=1:size(freq.labelcmb,1)
    sgncmbindx(i,1) = strmatch(freq.labelcmb{i,1}, freq.label, 'exact');
    sgncmbindx(i,2) = strmatch(freq.labelcmb{i,2}, freq.label, 'exact');
  end
  %sgncmbindx = [sgnindx sgnindx; sgncmbindx];
  sgncmbmat  = nan*zeros(length(sgnindx),length(sgnindx));
  for i=1:size(sgncmbindx,1)
    sgncmbmat(sgncmbindx(i,2),sgncmbindx(i,1)) = -i;
    sgncmbmat(sgncmbindx(i,1),sgncmbindx(i,2)) = i;
  end
%end

fprintf('checking whether partialisation of the requested channel(s) is possible\n');
prtindx = [];
for j = 1:length(cfg.partchan)
  prtindx = [prtindx; strmatch(cfg.partchan{j}, freq.label, 'exact')];
end
if isempty(prtindx), warning('no partialisation will be performed, requested channel(s) not present in data'); end;

prtindx     = sort(prtindx);
sgnindx     = setdiff(sgnindx, prtindx);
prtsgnmat   = sgncmbmat(prtindx,prtindx);
rstsgnmat   = sgncmbmat(sgnindx,sgnindx);
Nprtsgn     = length(prtindx);
Nrstsgn     = length(sgnindx);  
prtsgnlabel = freq.label(prtindx); 
label       = freq.label(sgnindx);

%the reformatting of the cross-spectral densities leads to a re-ordering of the signal-combinations
[ind1, ind2] = find(~isnan(rstsgnmat));
labelcmb     = [label(ind1) label(ind2)];
sgncmbindx   = [ind1 ind2];

%these cross-spectral density matrices should be complete
if any(isnan(prtsgnmat)),    error('partialisation of the requested channel(s) is not possible'); end;
prtrstsgnmat = sgncmbmat(prtindx,sgnindx);
if any(isnan(prtrstsgnmat)), error('partialisation of the requested channel(s) is not possible'); end;
rstprtsgnmat = sgncmbmat(sgnindx,prtindx);
if any(isnan(rstprtsgnmat)), error('partialisation of the requested channel(s) is not possible'); end;

%allocate memory
csdrr = nan+complex(zeros(Nrstsgn.^2,Nfrq,Ntim)      , zeros(Nrstsgn.^2,Nfrq,Ntim));
csdpp = nan+complex(zeros(Nprtsgn.^2,Nfrq,Ntim)      , zeros(Nprtsgn.^2,Nfrq,Ntim));
csdrp = nan+complex(zeros(Nrstsgn.*Nprtsgn,Nfrq,Ntim), zeros(Nrstsgn.*Nprtsgn,Nfrq,Ntim));
csdpr = nan+complex(zeros(Nprtsgn.*Nrstsgn,Nfrq,Ntim), zeros(Nprtsgn.*Nrstsgn,Nfrq,Ntim));

progress('init', cfg.feedback, 'partialising out the requested channel(s)')
for m = 1:size(crsspctrm,1)
  progress(m/size(crsspctrm,1), 'computing coherence %d/%d\n', m, size(crsspctrm,1));
  %a bit of hocus-pocus with the cross-spectral densities
  csdrr(find(~isnan(rstsgnmat)),:,:)  = squeeze(crsspctrm(m,abs(rstsgnmat(find(~isnan(rstsgnmat)))),:,:));
  csdpp(:) = squeeze(crsspctrm(m,abs(prtsgnmat),:,:));  
  csdrp(:) = squeeze(crsspctrm(m,abs(rstprtsgnmat),:,:));
  csdpr(:) = squeeze(crsspctrm(m,abs(prtrstsgnmat),:,:));
  
  %take care of the conjugates
  csdrr(find(rstsgnmat<0),:,:)     = conj(csdrr(find(rstsgnmat<0),:,:));
  csdpp(find(prtsgnmat<0),:,:)     = conj(csdpp(find(prtsgnmat<0),:,:));
  csdrp(find(rstprtsgnmat<0),:,:)  = conj(csdrp(find(rstprtsgnmat<0),:,:));
  csdpr(find(prtrstsgnmat<0),:,:)  = conj(csdpr(find(prtrstsgnmat<0),:,:));
  
  %perform the partialisation
  for j = 1:Nfrq
    for k = 1:Ntim
      rr  = reshape(squeeze(csdrr(:,j,k)),[Nrstsgn Nrstsgn]);
      rp  = reshape(squeeze(csdrp(:,j,k)),[Nrstsgn Nprtsgn]);
      pp  = reshape(squeeze(csdpp(:,j,k)),[Nprtsgn Nprtsgn]);
      pr  = reshape(squeeze(csdpr(:,j,k)),[Nprtsgn Nrstsgn]);
      rrp = rr - rp*pinv(pp)*pr;    
      csdp(m,:,j,k) = rrp(find(~isnan(rrp)));
    end
  end
end
progress('close');

freq.label     = label;
freq.labelcmb  = labelcmb;
%freq.crsspctrm = csdp;
freq           = setfield(freq, field, csdp);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that computes multichannel coherence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function coh = coherence(csd, n1, n2, method)

switch(method)
  case 'plain'
    %coh = (csd(1,2) ./sqrt(csd(1,1)*csd(2,2))); %this outputs the "normal" coherence
    coh = (csd(1,2).^2) ./(csd(1,1)*csd(2,2));
  case 'canonical' 
    coh = max(diag(cancorr(csd, 1:n1, n1+(1:n2)))); %this outputs the "squared" coherence
  case 'canonical_complex'
    dum      = diag(cancorr(csd, 1:n1, n1+(1:n2)));
    [m, ind] = max(dum);
    coh      = dum(ind);
  otherwise
    error(sprintf('method ''%s'' is not implemented for coherence', method));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that applies bias-correction to the coherence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function zcoh = zcoherence(coh, dof, method)

switch(method)
  case 'mitra'  
    if size(coh)~=size(dof),
      coh = coh';
    end 
    zcoh = sqrt( -(dof - 2) .* log(1 - coh) );
    beta = 23/20;
    zcoh = beta .* (zcoh - beta); 
  otherwise
    error(sprintf('method ''%s'' is not yet implemented for coherence', method));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that computes multichannel phase locking value
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plv = phaselocking(csd, n1, n2, method)
% this requires multiple repetitions
switch(method)
  case 'plain'
    plv = mean(csd(1,2) ./ abs(csd(1,2)));
  case 'canonical'
    csd = csd ./ abs(csd);
    plv = max(diag(cancorr(csd, 1:n1, n1+(1:n2))));
  otherwise
    error(sprintf('method ''%s'' is not implemented for phaselocking', method));
end
