std_ersp.m

% STD_ERSP - Compute ERSP and/or ITC transforms for ICA components 
%              or data channels of a dataset. Save results into Matlab 
%              float files. 
%
% Function description:
%              The function computes the mean ERSP or ITC for the selected 
%              dataset ICA components or data channels in the requested 
%              frequency range and time window (the two are dependent). 
%              Frequencies are equally log spaced. Options specify component 
%              numbers, desired frequency range,  time window length, 
%              frequency resolution, significance level, and wavelet
%              cycles. See >> help newtimef and >> timef details 
%
%              Two Matlab files are saved (for ERSP and ITC). These contain 
%              the ERSP|ITC image, plus the transform parameters 
%              used to compute them. Saves the computed dataset mean images 
%              in dataset-name files with extensions '.icaersp' and '.icaitc'
%              for ICA components or '.datersp', '.datitc' for data channels.
% Usage:  
%              >> [X times logfreqs ] = std_ersp(EEG, 'key', 'val', ...);
% Inputs:
%   EEG          - a loaded epoched EEG dataset structure. May be an array
%                  of such structure containing several datasets.
%
% Other inputs:
%   'trialindices' - [cell array] indices of trials for each dataset.
%                  Default is EMPTY (no trials). NEEDS TO BE SET.
%   'components' - [numeric vector] components of the EEG structure for which 
%                  activation spectrum will be computed. Note that because 
%                  computation of ERP is so fast, all components spectrum are
%                  computed and saved. Only selected component 
%                  are returned by the function to Matlab
%                  {default|[] -> all}
%   'channels'   - [cell array] channels of the EEG structure for which 
%                  activation spectrum will be computed. Note that because 
%                  computation of ERP is so fast, all channels spectrum are
%                  computed and saved. Only selected channels 
%                  are returned by the function to Matlab
%                  {default|[] -> none}
%   'recompute'  - ['on'|'off'] force recomputing data file even if it is 
%                  already on disk.
%   'rmcomps'    - [integer array] remove artifactual components (this entry
%                  is ignored when plotting components). This entry contains 
%                  the indices of the components to be removed. Default is none.
%   'interp'     - [struct] channel location structure containing electrode
%                  to interpolate ((this entry is ignored when plotting 
%                  components). Default is no interpolation.
%   'fileout'    - [string] name of the file to save on disk. The default
%                  is the same name (with a different extension) as the 
%                  dataset given as input.
%  'savetrials'  - ['on'|'off'] save single-trials ERSP. Requires a lot of disk
%                  space (dataset space on disk times 10) but allow for refined
%                  single-trial statistics. This option is obsolete as
%                  trials are now always saved.
%  'savefile'    - ['on'|'off'] save file or simply return measures.
%                  Default is to save files ('on').
%  'getparams'   - ['on'|'off'] return optional parameters for the newtimef 
%                  function (and do not compute anything). This argument is
%                  obsolete (default is 'off').
%
% ERSP optional inputs:
%   'type'       - ['ersp'|'itc'|'ersp&itc'] save ERSP, ITC, or both data 
%                  types to disk {default: 'ersp'}
%   'freqs'      - [minHz maxHz] the ERSP/ITC frequency range to compute 
%                  and return. {default: 3 to EEG sampling rate divided by 3}
%   'timelimits' - [minms maxms] time window (in ms) to compute.
%                  {default: whole input epoch}.
%   'cycles'     - [wavecycles (factor)]. If 0 -> DFT (constant window length 
%                  across frequencies).
%                  If >0 -> the number of cycles in each analysis wavelet. 
%                  If [wavecycles factor], wavelet cycles increase with 
%                  frequency, beginning at wavecyles. (0 < factor < 1) 
%                  factor = 0 -> fixed epoch length (DFT, as in FFT). 
%                  factor = 1 -> no increase (standard wavelets)
%                  {default: [0]}
%   'padratio'   - (power of 2). Multiply the number of output frequencies 
%                  by dividing their frequency spacing through 0-padding.
%                  Output frequency spacing is (low_freq/padratio).
%   'alpha'      - If in (0, 1), compute two-tailed permutation-based 
%                  probability thresholds and use these to mask the output 
%                  ERSP/ITC images {default: NaN}
%   'powbase'    - Deprecated. Note that baseline can be readjusted after 
%                  computation as single trial spectral decompositions are stored.
%
% Other optional inputs:
%   This function will take any of the NEWTIMEF optional inputs (for instance
%   to compute log-space frequencies)...
%
% Outputs:
%   X         - the masked log ERSP/ITC of the requested ICA components/channels 
%               in the selected frequency and time range. Note that for
%               optimization reasons, this parameter is now empty or 0. X
%               thus must be read from the datafile saved on disk.
%   times     - vector of time points for which the ERSPs/ITCs were computed. 
%   logfreqs  - vector of (equally log spaced) frequencies (in Hz) at which the 
%               log ERSP/ITC was evaluated. 
%   parameters - parameters given as input to the newtimef function.
%
% Files written or modified:     
%              [dataset_filename].icaersp   <-- saved component ERSPs
%              [dataset_filename].icaitc    <-- saved component ITCs
%              [dataset_filename].icatimef  <-- saved component single
%                                               trial decompositions.
%  OR for channels
%              [dataset_filename].datersp   <-- saved channel ERSPs
%              [dataset_filename].datitc    <-- saved channel ITCs
%              [dataset_filename].dattimef  <-- saved channel single
%                                               trial decompositions.
% Example: 
%            % Create mean ERSP and ITC images on disk for all comps from 
%            % dataset EEG use three-cycle wavelets (at 3 Hz) to more than 
%            % three-cycle wavelets at 50 Hz. See >> help newtimef 
%            % Return the (equally log-freq spaced, probability-masked) ERSP.
%            >> [Xersp, times, logfreqs] = std_ersp(EEG, ...
%                       'type', 'ersp', 'freqs', [3 50], 'cycles', [3 0.5]);
%
% See also: NEWTIMEF, STD_ERP, STD_SPEC, STD_TOPO, STD_PRECLUST
%
% Authors: Arnaud Delorme, Hilit Serby, SCCN, INC, UCSD, January, 2005-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, October 11, 2004, arno@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [X, times, logfreqs, parameters] = std_ersp(EEG, varargin)

if nargin < 1
    help std_ersp;
    return;
end

X = [];
options = {};
if length(varargin) > 1 
    if ~ischar(varargin{1})
        if length(varargin) > 0, options = { options{:} 'components' varargin{1} }; end
        if length(varargin) > 1, options = { options{:} 'freqs'      varargin{2} }; end
        if length(varargin) > 2, options = { options{:} 'timewindow' varargin{3} }; end
        if length(varargin) > 3, options = { options{:} 'cycles'     varargin{4} }; end
        if length(varargin) > 4, options = { options{:} 'padratio'   varargin{5} }; end
        if length(varargin) > 5, options = { options{:} 'alpha'      varargin{6} }; end
        if length(varargin) > 6, options = { options{:} 'type'       varargin{7} }; end
        if length(varargin) > 7, options = { options{:} 'powbase'    varargin{8} }; end
    else
        options = varargin;
    end
end

[g, timefargs] = finputcheck(options, { ...
                        'components'    'integer'               []          [];
                        'channels'      { 'cell','integer' }    { [] [] }   {};
                        'powbase'       'real'                  []          [];
                        'trialindices' { 'integer','cell' }     []          [];
                        'savetrials'    'string'      { 'on','off' }        'off';
                        'plot'        'string'      { 'on','off' }    'off';
                        'recompute'   'string'      { 'on','off' }    'off';
                        'getparams'     'string'      { 'on','off' }        'off';
                        'savefile'      'string'      { 'on','off' }        'on';
                        'parallel'      'string'      { 'on','off' }        'off';
                        'timewindow',    'real',    [], [];
                        'fileout',       'string',  [], '';
                        'timelimits'    'real'                  []          [EEG(1).xmin EEG(1).xmax]*1000;
                        'cycles'        'real'                  []          [3 .5];
                        'padratio'      'real'                  []          1;
                        'trialinfo'     'struct'                []          struct([]);
                        'freqs'         'real'                  []          [0 EEG(1).srate/2];
                        'rmcomps'       'cell'                  []          cell(1,length(EEG));
                        'interp'        'struct'                { }         struct([]);
                        'freqscale'     'string'                []         'log';
                        'alpha'         'real'                  []          NaN;
                        'baseline'      'real'                  []          0;
                        'type'          'string'      { 'ersp','itc','both','ersp&itc' }  'both'}, 'std_ersp', 'ignore');
if ischar(g), error(g); end
if isempty(g.trialindices), g.trialindices = cell(length(EEG)); end
if ~iscell(g.trialindices), g.trialindices = { g.trialindices }; end
if ~isempty(g.powbase), disp('''powbase'' parameter is no longer supported at computation time'); end

% checking input parameters
% -------------------------
if isempty(g.components) && isempty(g.channels)
    if isempty(EEG(1).icaweights)
        error('EEG.icaweights not found');
    end
    g.components = 1:size(EEG(1).icaweights,1);
    disp('Computing ERSP with default values for all components of the dataset');
end

% select ICA components or data channels
% --------------------------------------
if isempty(g.fileout), g.fileout = fullfile(EEG(1).filepath, EEG(1).filename(1:end-4)); end
if ~isempty(g.components)
    g.indices = g.components;
    prefix = 'comp';
    filenameersp   = [ g.fileout '.icaersp'  ];
    filenameitc    = [ g.fileout '.icaitc'   ];
    filenametrials = [ g.fileout '.icatimef' ];    
    if ~isempty(g.channels)
        error('Cannot compute ERSP/ITC for components and channels at the same time');
    end
elseif ~isempty(g.channels)
    if iscell(g.channels)
        if ~isempty(g.interp)
            g.indices = eeg_chaninds(g.interp, g.channels, 0);
        else
            g.indices = eeg_chaninds(EEG(1), g.channels, 0);
            for ind = 2:length(EEG)
                if ~isequal(eeg_chaninds(EEG(ind), g.channels, 0), g.indices)
                    error([ 'Channel information must be consistent when ' 10 'several datasets are merged for a specific design' ]);
                end
            end
        end
    else
        g.indices = g.channels;
    end
    prefix = 'chan';
    filenameersp   = [ g.fileout '.datersp'  ];
    filenameitc    = [ g.fileout '.datitc'   ];
    filenametrials = [ g.fileout '.dattimef' ];    
end

% Check if ERSP/ITC information found in datasets and if fits requested parameters 
% ----------------------------------------------------------------------------
if exist( filenametrials ) && strcmpi(g.recompute, 'off')
    fprintf('Use existing file for ERSP: %s; check the ''recompute checkbox'' to force recomputing.\n', filenameersp);
    return;
end

% Compute ERSP parameters
% -----------------------
parameters = { 'cycles', g.cycles, 'padratio', g.padratio, ...
               'alpha', g.alpha, 'freqscale', g.freqscale, timefargs{:} };
defaultlowfreq = 3;
[time_range] = compute_ersp_times(g.cycles,  EEG(1).srate, ...
                                 [EEG(1).xmin EEG(1).xmax]*1000 , defaultlowfreq, g.padratio); 
if time_range(1) < time_range(2) && g.freqs(1) == 0
     g.freqs(1) = defaultlowfreq; % for backward compatibility
end
parameters = { parameters{:} 'freqs' g.freqs };
if strcmpi(g.plot, 'off')
    parameters = { parameters{:} 'plotersp', 'off', 'plotitc', 'off', 'plotphase', 'off' };
end
parameters{end+1} = 'baseline';
parameters{end+1} = g.baseline;

% return parameters
% -----------------
if strcmpi(g.getparams, 'on')
    X = []; times = []; logfreqs = [];
    if strcmpi(g.savetrials, 'on')
        parameters = { parameters{:} 'savetrials', g.savetrials };
    end
    return;
end

options = {};
if ~isempty(g.rmcomps), options = { options{:} 'rmcomps' g.rmcomps }; end
if ~isempty(g.interp),  options = { options{:} 'interp' g.interp }; end
if isempty(g.channels)
     X = eeg_getdatact(EEG, 'component', g.indices, 'trialindices', g.trialindices );
else X = eeg_getdatact(EEG, 'channel'  , g.indices, 'trialindices', g.trialindices, 'rmcomps', g.rmcomps, 'interp', g.interp);
end
if size(X, 3) == 1
    error('The data is continuous for one of the dataset. ERSP can only be computed when data trials are present');
end

% frame range
% -----------
pointrange1 = round(max((g.timelimits(1)/1000-EEG(1).xmin)*EEG(1).srate, 1));
pointrange2 = round(min(((g.timelimits(2)+1000/EEG(1).srate)/1000-EEG(1).xmin)*EEG(1).srate, EEG(1).pnts));
pointrange = [pointrange1:pointrange2];

% Compute ERSP && ITC
% ------------------
allTrialsTmp   = cell(1,length(g.indices));
allTrialsTime  = cell(1,length(g.indices));
allTrialsFreqs = cell(1,length(g.indices));
eeglab_options;
usesingle = option_single;

disp('Computing time/frequency decomposition...');
parfor k = 1:length(g.indices)
    tmpparams = parameters;
    if length(g.indices) > 1
        tmpparams{end+1} = 'verbose';
        tmpparams{end+1} = 'off';
    end
    
    % Run timef() to get ERSP
    % ------------------------
    timefdata  = reshape(X(k,pointrange,:), 1, length(pointrange)*size(X,3));
    mytimes = [];
    mylogfreqs = [];
    alltfX   = [];
    if ~isempty(timefdata)
        [logersp,logitc,logbase,mytimes,mylogfreqs,logeboot,logiboot,alltfX] ...
              = newtimef( timefdata, length(pointrange), g.timelimits, EEG(1).srate, tmpparams{2:end});
        %figure; newtimef( TMP.data(32,:), EEG.pnts, [EEG.xmin EEG.xmax]*1000, EEG.srate, cycles, 'freqs', freqs);
        %figure; newtimef( timefdata, length(pointrange), g.timelimits, EEG.srate, cycles, 'freqs', freqs);
    end
    %if strcmpi(g.plot, 'on'), return; end
    if usesingle
        alltfX = single(alltfX);
    end

    allTrialsTmp{k}   = single( alltfX );
    allTrialsTime{k}  = mytimes;
    allTrialsFreqs{k} = mylogfreqs;
end
all_trials = [];
for k = 1:length(g.indices)  % for each (specified) component/channel
    all_trials = setfield( all_trials, [ prefix int2str(g.indices(k)) ], allTrialsTmp{k});
end
X = allTrialsTmp{1};

% Save ERSP into file
% -------------------
logfreqs             = allTrialsFreqs{1};
times                = allTrialsTime{1};
all_trials.freqs     = allTrialsFreqs{1};
all_trials.times     = allTrialsTime{1};
all_trials.parameters = { options{:} parameters{:} };
all_trials.datatype   = 'TIMEF';
all_trials.datafiles  = computeFullFileName( { EEG.filepath }, { EEG.filename });
all_trials.datatrials = g.trialindices;

all_trials.parameters = parameters;
if ~isempty(g.channels)
    if ~isempty(g.interp)
        all_trials.labels = { g.interp(g.indices).labels };
    elseif ~isempty(EEG(1).chanlocs)
        tmpchanlocs = EEG(1).chanlocs;
        all_trials.labels = { tmpchanlocs(g.indices).labels };
    end
end
all_trials.trialinfo = g.trialinfo;

if strcmpi(g.savefile, 'on')
    std_savedat( filenametrials , all_trials );
end

% compute full file names
% -----------------------
function res = computeFullFileName(filePaths, fileNames);
for index = 1:length(fileNames)
    res{index} = fullfile(filePaths{index}, fileNames{index});
end