ft_freqanalysis.m

function [freq] = ft_freqanalysis(cfg, data)

% FT_FREQANALYSIS performs frequency and time-frequency analysis
% on time series data over multiple trials
%
% Use as
%   [freq] = ft_freqanalysis(cfg, data)
%
% The input data should be organised in a structure as obtained from
% the FT_PREPROCESSING or the FT_MVARANALYSIS function. The configuration
% depends on the type of computation that you want to perform.
%
% The configuration should contain:
%   cfg.method      = different methods of calculating the spectra
%                     'mtmfft', analyses an entire spectrum for the entire data
%                       length, implements multitaper frequency transformation.
%                     'mtmconvol', implements multitaper time-frequency
%                       transformation based on multiplication in the
%                       frequency domain.
%                     'wavelet', implements wavelet time frequency
%                       transformation (using Morlet wavelets) based on
%                       multiplication in the frequency domain.
%                     'tfr', implements wavelet time frequency
%                       transformation (using Morlet wavelets) based on
%                       convolution in the time domain.
%                     'mvar', does a fourier transform on the coefficients
%                       of an estimated multivariate autoregressive model,
%                       obtained with FT_MVARANALYSIS. In this case, the
%                       output will contain a spectral transfer matrix,
%                       the cross-spectral density matrix, and the
%                       covariance matrix of the innovatio noise.
%                     'superlet', combines Morlet-wavelet based
%                       decompositions, see below.
%                     'irasa', implements Irregular-Resampling Auto-Spectral 
%                       Analysis (IRASA), to separate the fractal components 
%                       from the periodicities in the signal.
%   cfg.output      = 'pow'       return the power-spectra
%                     'powandcsd' return the power and the cross-spectra
%                     'fourier'   return the complex Fourier-spectra
%                     'fractal'   (when cfg.method = 'irasa'), return the
%                       fractal component of the spectrum (1/f)
%                     'original'  (when cfg.method = 'irasa'), return the
%                       full power spectrum
%                     'fooof'     returns a smooth power-spectrum,
%                       based on a parametrization of a mixture of aperiodic and periodic
%                       components (only works with cfg.method = 'mtmfft')
%                     'fooof_aperiodic' returns a power-spectrum with the
%                       fooof based estimate of the aperiodic part of the signal.
%                     'fooof_peaks' returns a power-spectrum with the fooof
%                       based estimate of the aperiodic signal removed,
%                       it's expressed as
%                       10^(log10(fooof)-log10(fooof_aperiodic))
%   cfg.channel     = Nx1 cell-array with selection of channels (default = 'all'),
%                       see FT_CHANNELSELECTION for details
%   cfg.channelcmb  = Mx2 cell-array with selection of channel pairs (default = {'all' 'all'}),
%                       see FT_CHANNELCOMBINATION for details
%   cfg.trials      = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.keeptrials  = 'yes' or 'no', return individual trials or average (default = 'no')
%   cfg.keeptapers  = 'yes' or 'no', return individual tapers or average (default = 'no')
%   cfg.pad         = number, 'nextpow2', or 'maxperlen' (default), length
%                      in seconds to which the data can be padded out. The
%                      padding will determine your spectral resolution. If you
%                      want to compare spectra from data pieces of different
%                      lengths, you should use the same cfg.pad for both, in
%                      order to spectrally interpolate them to the same
%                      spectral resolution.  The new option 'nextpow2' rounds
%                      the maximum trial length up to the next power of 2.  By
%                      using that amount of padding, the FFT can be computed
%                      more efficiently in case 'maxperlen' has a large prime
%                      factor sum.
%   cfg.padtype     = string, type of padding (default 'zero', see
%                      ft_preproc_padding)
%   cfg.polyremoval = number (default = 0), specifying the order of the
%                      polynome which is fitted and subtracted from the time
%                      domain data prior to the spectral analysis. For
%                      example, a value of 1 corresponds to a linear trend.
%                      The default is a mean subtraction, thus a value of 0.
%                      If no removal is requested, specify -1.
%                      see FT_PREPROC_POLYREMOVAL for details
%
%
% METHOD SPECIFIC OPTIONS AND DESCRIPTIONS
%
% MTMFFT performs frequency analysis on any time series trial data using a
% conventional single taper (e.g. Hanning) or using the multiple tapers based on
% discrete prolate spheroidal sequences (DPSS), also known as the Slepian
% sequence.
%   cfg.taper      = 'dpss', 'hanning' or many others, see WINDOW (default = 'dpss')
%                     For cfg.output='powandcsd', you should specify the channel combinations
%                     between which to compute the cross-spectra as cfg.channelcmb. Otherwise
%                     you should specify only the channels in cfg.channel.
%   cfg.foilim     = [begin end], frequency band of interest
%       OR
%   cfg.foi        = vector 1 x numfoi, frequencies of interest
%   cfg.tapsmofrq  = number, the amount of spectral smoothing through
%                    multi-tapering. Note that 4 Hz smoothing means
%                    plus-minus 4 Hz, i.e. a 8 Hz smoothing box.
%
% MTMCONVOL performs time-frequency analysis on any time series trial data using
% the 'multitaper method' (MTM) based on Slepian sequences as tapers.
% Alternatively, you can use conventional tapers (e.g. Hanning).
%   cfg.tapsmofrq  = vector 1 x numfoi, the amount of spectral smoothing
%                     through multi-tapering. Note that 4 Hz smoothing means
%                     plus-minus 4 Hz, i.e. a 8 Hz smoothing box.
%    cfg.foi        = vector 1 x numfoi, frequencies of interest
%    cfg.taper      = 'dpss', 'hanning' or many others, see WINDOW (default = 'dpss')
%                      For cfg.output='powandcsd', you should specify the channel combinations
%                      between which to compute the cross-spectra as cfg.channelcmb. Otherwise
%                      you should specify only the channels in cfg.channel.
%    cfg.t_ftimwin  = vector 1 x numfoi, length of time window (in seconds)
%    cfg.toi        = vector 1 x numtoi, the times on which the analysis
%                     windows should be centered (in seconds), or a string
%                     such as '50%' or 'all' (default).  Both string options
%                     use all timepoints available in the data, but 'all'
%                     centers a spectral estimate on each sample, whereas
%                     the percentage specifies the degree of overlap between
%                     the shortest time windows from cfg.t_ftimwin.
%
% WAVELET performs time-frequency analysis on any time series trial data using the
% 'wavelet method' based on Morlet wavelets. Using mulitplication in the frequency
% domain instead of convolution in the time domain.
%   cfg.foi        = vector 1 x numfoi, frequencies of interest
%       OR
%   cfg.foilim     = [begin end], frequency band of interest
%   cfg.toi        = vector 1 x numtoi, the times on which the analysis
%                    windows should be centered (in seconds)
%   cfg.width      = 'width', or number of cycles, of the wavelet (default = 7)
%   cfg.gwidth     = determines the length of the used wavelets in standard
%                    deviations of the implicit Gaussian kernel and should
%                    be chosen >= 3; (default = 3)
%
% The standard deviation in the frequency domain (sf) at frequency f0 is
% defined as: sf = f0/width
% The standard deviation in the temporal domain (st) at frequency f0 is
% defined as: st = 1/(2*pi*sf)
%
% SUPERLET performs time-frequency analysis on any time series trial data using the
% 'superlet method' based on a frequency-wise combination of Morlet wavelets of varying cycle
% widths (see Moca et al. 2019, https://doi.org/10.1101/583732).
%   cfg.foi                 = vector 1 x numfoi, frequencies of interest
%       OR
%   cfg.foilim              = [begin end], frequency band of interest
%   cfg.toi                 = vector 1 x numtoi, the times on which the analysis
%                             windows should be centered (in seconds)
%   cfg.width  = 'width', or number of cycles, of the base wavelet (default = 3)
%   cfg.gwidth     = determines the length of the used wavelets in standard
%                             deviations of the implicit Gaussian kernel and should
%                             be chosen >= 3; (default = 3)
%   cfg.combine    = 'additive', 'multiplicative' (default = 'additive')
%                             determines if cycle numbers of wavelets comprising a superlet
%                             are chosen additively or multiplicatively
%   cfg.order      = vector 1 x numfoi, superlet order, i.e. number of combined
%                             wavelets, for individual frequencies of interest.
%
% The standard deviation in the frequency domain (sf) at frequency f0 is
% defined as: sf = f0/width
% The standard deviation in the temporal domain (st) at frequency f0 is
% defined as: st = 1/(2*pi*sf)
%
% HILBERT performs time-frequency analysis on any time series data using a frequency specific
% bandpass filter, followed by the Hilbert transform.
%   cfg.foi              = vector 1 x numfoi, frequencies of interest
%   cfg.toi              = vector 1 x numtoi, the time points for which the estimates will be returned (in seconds)
%   cfg.width            = scalar, or vector (default: 1), specifying the half bandwidth of the filter;
%   cfg.edgartnan        = 'no' (default) or 'yes', replace filter edges with nans, works only for finite impulse response (FIR) filters, and
%                          requires a user specification of the filter order
%   
% For the bandpass filtering the following options can be specified, the default values are as in FT_PREPROC_BANDPASSFILTER, for more 
% information see the help of FT_PREPROCESSING
%   cfg.bpfilttype
%   cfg.bpfiltord        = (optional) scalar, or vector 1 x numfoi;
%   cfg.bpfiltdir        
%   cfg.bpinstabilityfix
%   cfg.bpfiltdf
%   cfg.bpfiltwintype
%   cfg.bpfiltdev
%
% TFR performs time-frequency analysis on any time series trial data using the
% 'wavelet method' based on Morlet wavelets. Using convolution in the time domain
% instead of multiplication in the frequency domain.
%   cfg.foi        = vector 1 x numfoi, frequencies of interest
%       OR
%   cfg.foilim     = [begin end], frequency band of interest
%   cfg.width      = 'width', or number of cycles, of the wavelet (default = 7)
%   cfg.gwidth     = determines the length of the used wavelets in standard
%                    deviations of the implicit Gaussian kernel and should
%                    be choosen >= 3; (default = 3)
%
%
% 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_FREQSTATISTICS, FT_FREQDESCRIPTIVES, FT_CONNECTIVITYANALYSIS

% Guidelines for use in an analysis pipeline:
% after FT_FREQANALYSIS you will have frequency or time-frequency
% representations (TFRs) of the data, represented as power-spectra,
% power and cross-spectra, or complex fourier-spectra, either for individual
% trials or an average over trials.
% This usually serves as input for one of the following functions:
%    * FT_FREQDESCRIPTIVES  to compute descriptive univariate statistics
%    * FT_FREQSTATISTICS    to perform parametric or non-parametric statistical tests
%    * FT_FREQBASELINE      to perform baseline normalization of the spectra
%    * FT_FREQGRANDAVERAGE  to compute the average spectra over multiple subjects or datasets
%    * FT_CONNECTIVITYANALYSIS to compute various measures of connectivity
% Furthermore, the data can be visualised using the various plotting
% functions, including:
%    * FT_SINGLEPLOTTFR     to plot the TFR of a single channel or the average over multiple channels
%    * FT_TOPOPLOTTFR       to plot the topographic distribution over the head
%    * FT_MULTIPLOTTFR      to plot TFRs in a topographical layout

% Undocumented local options:
% cfg.method = 'hilbert'. Keeping this as undocumented as it does not make
%              sense to use in ft_freqanalysis unless the user is doing his
%              own filter-padding to remove edge-artifacts
% cfg.correctt_ftimwin (set to yes to try to determine new t_ftimwins based
%                       on correct cfg.foi)

% Copyright (C) 2003-2006, F.C. Donders Centre, Pascal Fries
% Copyright (C) 2004-2006, F.C. Donders Centre, Markus Siegel
% Copyright (C) 2007-2012, DCCN, The FieldTrip team
%
% 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/>.

% 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 data
ft_preamble provenance data
ft_preamble trackconfig

% 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
data = ft_checkdata(data, 'datatype', {'raw', 'raw+comp', 'mvar'}, 'feedback', 'yes', 'hassampleinfo', 'yes');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels', 'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'renamed',    {'label', 'channel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'sgn',   'channel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'labelcmb', 'channelcmb'});
cfg = ft_checkconfig(cfg, 'renamed',    {'sgncmb',   'channelcmb'});
cfg = ft_checkconfig(cfg, 'required',   {'method'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'fft',    'mtmfft'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'convol', 'mtmconvol'});
cfg = ft_checkconfig(cfg, 'forbidden',  {'latency'}); % see bug 1376 and 1076
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'wltconvol', 'wavelet'});

% set the defaults
cfg.feedback    = ft_getopt(cfg, 'feedback',   'text');
cfg.inputlock   = ft_getopt(cfg, 'inputlock',  []);  % this can be used as mutex when doing distributed computation
cfg.outputlock  = ft_getopt(cfg, 'outputlock', []);  % this can be used as mutex when doing distributed computation
cfg.trials      = ft_getopt(cfg, 'trials',     'all', 1);
cfg.channel     = ft_getopt(cfg, 'channel',    'all');

% select channels and trials of interest, by default this will select all channels and trials
tmpcfg = keepfields(cfg, {'trials', 'channel', 'tolerance', 'showcallinfo', 'trackcallinfo', 'trackconfig', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo'});
data = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

% some proper error handling
if isfield(data, 'trial') && numel(data.trial)==0
  ft_error('no trials were selected'); % this does not apply for MVAR data
end

if numel(data.label)==0
  ft_error('no channels were selected');
end

% switch over method and do some of the method specfic checks and defaulting
switch cfg.method
  
  case 'mtmconvol'
    cfg.taper = ft_getopt(cfg, 'taper', 'dpss');
    if isequal(cfg.taper, 'dpss') && ~isfield(cfg, 'tapsmofrq')
      ft_error('you must specify a smoothing parameter with taper = dpss');
    end
    % check for foi above Nyquist
    if isfield(cfg, 'foi')
      if any(cfg.foi > (data.fsample+100*eps(data.fsample))/2)
        % add a small number to allow for numeric tolerance issues
        ft_error('frequencies in cfg.foi are above Nyquist')
      end
      if isequal(cfg.taper, 'dpss') && not(isfield(cfg, 'tapsmofrq'))
        ft_error('you must specify a smoothing parameter with taper = dpss');
      end
    end
    cfg = ft_checkconfig(cfg, 'required', {'toi', 't_ftimwin'});
    if ischar(cfg.toi)
      begtim  = min(cellfun(@min,data.time));
      endtim  = max(cellfun(@max,data.time));
      if strcmp(cfg.toi, 'all') % each data sample gets a time window
        cfg.toi = linspace(begtim, endtim, round((endtim-begtim) ./ ...
          mean(diff(data.time{1})))+1);
      elseif strcmp(cfg.toi(end), '%') % percent overlap between smallest time windows
        overlap = str2double(cfg.toi(1:(end-1)))/100;
        cfg.toi = linspace(begtim, endtim, round((endtim-begtim) ./ ...
          (overlap * min(cfg.t_ftimwin))) + 1);
      else
        ft_error('cfg.toi should be either a numeric vector or a string: can be ''all'' or a percentage (e.g., ''50%'')');
      end
    end
    
  case 'mtmfft'
    cfg.taper = ft_getopt(cfg, 'taper', 'dpss');
    if isequal(cfg.taper, 'dpss') && not(isfield(cfg, 'tapsmofrq'))
      ft_error('you must specify a smoothing parameter with taper = dpss');
    end
    % check for foi above Nyquist
    if isfield(cfg, 'foi')
      if any(cfg.foi > (data.fsample/2))
        ft_error('frequencies in cfg.foi are above Nyquist')
      end
    end
    if isequal(cfg.taper, 'dpss') && not(isfield(cfg, 'tapsmofrq'))
      ft_error('you must specify a smoothing parameter with taper = dpss');
    end
    
  case 'irasa'
    cfg.taper  = ft_getopt(cfg, 'taper',  'hanning');
    cfg.output = ft_getopt(cfg, 'output', 'fractal');
    cfg.pad    = ft_getopt(cfg, 'pad',    'nextpow2');
    if ~isequal(cfg.taper, 'hanning')
      ft_error('the irasa method supports hanning tapers only');
    end
    if ~isequal(cfg.pad, 'nextpow2')
      ft_warning('consider using cfg.pad=''nextpow2'' for the irasa method');
    end
    % check for foi above Nyquist
    if isfield(cfg, 'foi')
      if any(cfg.foi > (data.fsample/2))
        ft_error('frequencies in cfg.foi are above Nyquist')
      end
    end
    
  case 'wavelet'
    cfg.width  = ft_getopt(cfg, 'width',  7);
    cfg.gwidth = ft_getopt(cfg, 'gwidth', 3);
    
  case 'superlet'
    % reorganize the cfg, a nested cfg is not consistent with the othe methods
    cfg = ft_checkconfig(cfg, 'createtopcfg', 'superlet');
    cfg = removefields(cfg, 'superlet');
    cfg = ft_checkconfig(cfg, 'renamed', {'basewidth', 'width'});
    
    cfg.width   = ft_getopt(cfg, 'width',   3);
    cfg.gwidth  = ft_getopt(cfg, 'gwidth',  3);
    cfg.combine = ft_getopt(cfg, 'combine', 'additive');
    cfg.order   = ft_getopt(cfg, 'order',   ones(1, numel(cfg.foi)));
    if numel(cfg.order) == 1
      cfg.order = cfg.order.*length(cfg.foi);
    elseif numel(cfg.order)~= numel(cfg.foi)
      ft_error('cfg.foi must have the same number of elements as cfg.foi, or must be a scalar');
    end
    
  case 'tfr'
    cfg = ft_checkconfig(cfg, 'renamed', {'waveletwidth', 'width'});
    cfg = ft_checkconfig(cfg, 'unused',  {'downsample'});
    cfg.width  = ft_getopt(cfg, 'width',  7);
    cfg.gwidth = ft_getopt(cfg, 'gwidth', 3);
    
  case 'hilbert'
    ft_warning('method = hilbert may require user action to deal with filtering-artifacts')
    cfg = ft_checkconfig(cfg, 'renamed', {'filttype',  'bpfilttype'});
    cfg = ft_checkconfig(cfg, 'renamed', {'filtorder', 'bpfiltord'});
    cfg = ft_checkconfig(cfg, 'renamed', {'filtdir',   'bpfiltdir'});
    cfg.bpfilttype       = ft_getopt(cfg, 'bpfilttype');
    cfg.bpfiltord        = ft_getopt(cfg, 'bpfiltord');
    cfg.bpfiltdir        = ft_getopt(cfg, 'bpfiltdir');
    cfg.bpinstabilityfix = ft_getopt(cfg, 'bpinstabilityfix');
    cfg.bpfiltdf         = ft_getopt(cfg, 'bpfiltdf');
    cfg.bpfiltwintype    = ft_getopt(cfg, 'bpfiltwintype');
    cfg.bpfiltdev        = ft_getopt(cfg, 'bpfiltdev');
    cfg.width            = ft_getopt(cfg, 'width', 1);
    cfg.edgartnan        = istrue(ft_getopt(cfg, 'edgeartnan', 'no'));
    
    fn = fieldnames(cfg);
    bpfiltoptions = ft_cfg2keyval(keepfields(cfg, fn(startsWith(fn, 'bp'))));
    
  case 'mvar'
    if isfield(cfg, 'inputfile')
      freq = ft_freqanalysis_mvar(cfg);
    else
      freq = ft_freqanalysis_mvar(cfg, data);
    end
    return
    
  case 'neuvar'
    cfg.order  = ft_getopt(cfg, 'order',  1); % order of differentiation
    
  otherwise
    ft_error('specified cfg.method is not supported')
end

% set all the defaults
cfg.pad               = ft_getopt(cfg, 'pad',       []);
if isempty(cfg.pad)
  ft_notice('Default cfg.pad=''maxperlen'' can run slowly. Consider using cfg.pad=''nextpow2'' for more efficient FFT computation.')
  cfg.pad = 'maxperlen';
end
cfg.padtype           = ft_getopt(cfg, 'padtype',   'zero');
cfg.output            = ft_getopt(cfg, 'output',    'pow'); % the default for irasa is set earlier
cfg.calcdof           = ft_getopt(cfg, 'calcdof',   'no');
cfg.channel           = ft_getopt(cfg, 'channel',   'all');
cfg.precision         = ft_getopt(cfg, 'precision', 'double');
cfg.foi               = ft_getopt(cfg, 'foi',       []);
cfg.foilim            = ft_getopt(cfg, 'foilim',    []);
cfg.correctt_ftimwin  = ft_getopt(cfg, 'correctt_ftimwin', 'no');
cfg.polyremoval       = ft_getopt(cfg, 'polyremoval', 0);

% keeptrials and keeptapers should be conditional on cfg.output,
% cfg.output = 'fourier' should always output tapers
if strcmp(cfg.output, 'fourier')
  cfg.keeptrials = ft_getopt(cfg, 'keeptrials', 'yes');
  cfg.keeptapers = ft_getopt(cfg, 'keeptapers', 'yes');
  if strcmp(cfg.keeptrials, 'no') || strcmp(cfg.keeptapers, 'no')
    ft_error('cfg.output = ''fourier'' requires cfg.keeptrials = ''yes'' and cfg.keeptapers = ''yes''');
  end
else
  cfg.keeptrials = ft_getopt(cfg, 'keeptrials', 'no');
  cfg.keeptapers = ft_getopt(cfg, 'keeptapers', 'no');
end

% set flags for keeping trials and/or tapers
if strcmp(cfg.keeptrials, 'no') &&  strcmp(cfg.keeptapers, 'no')
  keeprpt = 1;
elseif strcmp(cfg.keeptrials, 'yes') &&  strcmp(cfg.keeptapers, 'no')
  keeprpt = 2;
elseif strcmp(cfg.keeptrials, 'no') &&  strcmp(cfg.keeptapers, 'yes')
  ft_error('There is currently no support for keeping tapers WITHOUT KEEPING TRIALS.');
elseif strcmp(cfg.keeptrials, 'yes') &&  strcmp(cfg.keeptapers, 'yes')
  keeprpt = 4;
end
if strcmp(cfg.keeptrials, 'yes') && strcmp(cfg.keeptapers, 'yes')
  if ~strcmp(cfg.output, 'fourier')
    ft_error('Keeping trials AND tapers is only possible with fourier as the output.');
  end
end

% Set flags for output
if ismember(cfg.output, {'pow','fractal','original','fooof','fooof_peaks','fooof_aperiodic'})
  powflg = 1;
  csdflg = 0;
  fftflg = 0;
elseif strcmp(cfg.output, 'powandcsd')
  powflg = 1;
  csdflg = 1;
  fftflg = 0;
elseif strcmp(cfg.output, 'fourier')
  powflg = 0;
  csdflg = 0;
  fftflg = 1;
else
  ft_error('Unrecognized output required');
end

% Check whether the keeptrials is correct for fooof
if startsWith(cfg.output, 'fooof')
  % ensure that Brainstorm is on the path: if the user uses their own
  % version of the code, assume that the paths are correctly set
  if keeprpt~=1
    ft_error('Keeping trials and/or tapers is not allowed when using fooof');
  end
  if ~isequal(cfg.method, 'mtmfft')
    ft_error('Fooof is only supported with cfg.method = ''mtmfft''');
  end
end

% prepare channel(cmb)
if ~isfield(cfg, 'channelcmb') && csdflg
  %set the default for the channelcombination
  cfg.channelcmb = {'all' 'all'};
elseif isfield(cfg, 'channelcmb') && ~csdflg
  % no cross-spectrum needs to be computed, hence remove the combinations from cfg
  cfg = rmfield(cfg, 'channelcmb');
end

if isfield(cfg, 'channelcmb')
  % the channels in the data are already the subset according to cfg.channel
  cfg.channelcmb = ft_channelcombination(cfg.channelcmb, data.label);
end

% determine the corresponding indices of all channels
chanind    = match_str(data.label, cfg.channel);
nchan      = numel(chanind);
if csdflg
  assert(nchan>1, 'CSD output requires multiple channels');
  % determine the corresponding indices of all channel combinations
  [dummy,chancmbind(:,1)] = match_str(cfg.channelcmb(:,1), data.label);
  [dummy,chancmbind(:,2)] = match_str(cfg.channelcmb(:,2), data.label);
  
  nchancmb   = size(chancmbind,1);
  chanind    = unique([chanind(:); chancmbind(:)]);
  nchan      = length(chanind);
  cutdatindcmb = zeros(size(chancmbind));
  for ichan = 1:nchan
    cutdatindcmb(chancmbind == chanind(ichan)) = ichan;
  end
end

% determine trial characteristics
ntrials = numel(data.trial);
trllength = zeros(1, ntrials);
for itrial = 1:ntrials
  trllength(itrial) = size(data.trial{itrial}, 2);
end
if strcmp(cfg.pad, 'maxperlen')
  padding = max(trllength);
  cfg.pad = padding/data.fsample;
elseif strcmp(cfg.pad, 'nextpow2')
  padding = 2^nextpow2(max(trllength));
  cfg.pad = padding/data.fsample;
else
  padding = cfg.pad*data.fsample;
  if padding<max(trllength)
    ft_error('the specified padding is too short');
  end
end

% correct foi and implement foilim 'backwards compatibility'
if ~isempty(cfg.foi) && ~isempty(cfg.foilim)
  ft_error('use either cfg.foi or cfg.foilim')
elseif ~isempty(cfg.foilim)
  % get the full foi in the current foilim and set it too be used as foilim
  fboilim = round(cfg.foilim .* cfg.pad) + 1;
  fboi    = fboilim(1):1:fboilim(2);
  cfg.foi = (fboi-1) ./ cfg.pad;
else
  % correct foi if foilim was empty and try to correct t_ftimwin (by detecting whether there is a constant factor between foi and t_ftimwin: cyclenum)
  oldfoi = cfg.foi;
  fboi   = round(cfg.foi .* cfg.pad) + 1;
  cfg.foi    = (fboi-1) ./ cfg.pad; % boi - 1 because 0 Hz is included in fourier output
  if strcmp(cfg.correctt_ftimwin, 'yes')
    cyclenum = oldfoi .* cfg.t_ftimwin;
    cfg.t_ftimwin = cyclenum ./ cfg.foi;
  end
end

% tapsmofrq compatibility between functions (make it into a vector if it's not)
if isfield(cfg, 'tapsmofrq')
  if strcmp(cfg.method, 'mtmconvol') && length(cfg.tapsmofrq) == 1 && length(cfg.foi) ~= 1
    cfg.tapsmofrq = ones(length(cfg.foi),1) * cfg.tapsmofrq;
  elseif strcmp(cfg.method, 'mtmfft') && length(cfg.tapsmofrq) ~= 1
    ft_warning('cfg.tapsmofrq should be a single number when cfg.method = mtmfft, now using only the first element')
    cfg.tapsmofrq = cfg.tapsmofrq(1);
  end
end

% options that don't change over trials
if isfield(cfg, 'tapsmofrq')
  options = {'pad', cfg.pad, 'padtype', cfg.padtype, 'freqoi', cfg.foi, 'tapsmofrq', cfg.tapsmofrq, 'polyorder', cfg.polyremoval, 'output', cfg.output};
else
  options = {'pad', cfg.pad, 'padtype', cfg.padtype, 'freqoi', cfg.foi, 'polyorder', cfg.polyremoval, 'output', cfg.output};
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Main loop over trials, inside fourierspectra are obtained and transformed into the appropriate outputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is done on trial basis to save memory

ft_progress('init', cfg.feedback, 'processing trials');
trlcnt = []; % only some methods need this variable, but it needs to be defined outside the trial loop
for itrial = 1:ntrials
  
  fbopt.i = itrial;
  fbopt.n = ntrials;
  
  dat = data.trial{itrial}; % chansel has already been performed
  time = data.time{itrial};
  
  clear spectrum % in case of very large trials, this lowers peak mem usage a bit
  
  % Perform specest call and set some specifics
  switch cfg.method
    
    case 'mtmconvol'
      [spectrum_mtmconvol,ntaper,foi,toi] = ft_specest_mtmconvol(dat, time, 'timeoi', cfg.toi, 'timwin', cfg.t_ftimwin, 'taper', ...
        cfg.taper, options{:}, 'dimord', 'chan_time_freqtap', 'feedback', fbopt);
      
      % the following variable is created to keep track of the number of
      % trials per time bin and is needed for proper normalization if
      % keeprpt==1 and the triallength is variable
      if itrial==1, trlcnt = zeros(1, numel(foi), numel(toi)); end
      
      hastime = true;
      % error for different number of tapers per trial
      if (keeprpt == 4) && any(ntaper(:) ~= ntaper(1))
        ft_error('currently you can only keep trials AND tapers, when using the number of tapers per frequency is equal across frequency')
      end
      % create tapfreqind for later indexing
      freqtapind = [];
      tempntaper = [0; cumsum(ntaper(:))];
      for iindfoi = 1:numel(foi)
        freqtapind{iindfoi} = tempntaper(iindfoi)+1:tempntaper(iindfoi+1);
      end
      
    case 'mtmfft'
      [spectrum,ntaper,foi] = ft_specest_mtmfft(dat, time, 'taper', cfg.taper, options{:}, 'feedback', fbopt);
      hastime = false;
      
    case 'irasa'
      [spectrum,ntaper,foi] = ft_specest_irasa(dat, time, options{:}, 'feedback', fbopt);
      hastime = false;
      
    case 'wavelet'
      [spectrum,foi,toi] = ft_specest_wavelet(dat, time, 'timeoi', cfg.toi, 'width', cfg.width, 'gwidth', cfg.gwidth, options{:}, 'feedback', fbopt);
      
      % the following variable is created to keep track of the number of
      % trials per time bin and is needed for proper normalization if
      % keeprpt==1 and the triallength is variable
      if itrial==1, trlcnt = zeros(1, numel(foi), numel(toi)); end
      
      hastime = true;
      % create FAKE ntaper (this requires very minimal code change below for compatibility with the other specest functions)
      ntaper = ones(1,numel(foi));
      % modify spectrum for same reason as fake ntaper
      spectrum = reshape(spectrum,[1 nchan numel(foi) numel(toi)]);
      
    case 'superlet'
      % calculate number of wavelets and respective cycle width dependent on superlet order
      % equivalent one-liners:
      %   multiplicative: cycles = arrayfun(@(order) arrayfun(@(wl_num) cfg.width*wl_num, 1:order), cfg.order,'uni',0)
      %   additive: cycles = arrayfun(@(order) arrayfun(@(wl_num) cfg.width+wl_num-1, 1:order), cfg.order,'uni',0)
      cycles = cell(length(cfg.foi),1);
      for i_f = 1:length(cfg.foi)
        frq_cyc = NaN(1,cfg.order(i_f));
        if strcmp(cfg.combine, 'multiplicative')
          for i_wl = 1:cfg.order(i_f)
            frq_cyc(i_wl) = cfg.width*i_wl;
          end
        elseif strcmp(cfg.combine, 'additive')
          for i_wl = 1:cfg.order(i_f)
            frq_cyc(i_wl) = cfg.width+i_wl-1;
          end
        end
        cycles{i_f} = frq_cyc;
      end
      
      % compute superlets
      spectrum = NaN(nchan,length(cfg.foi),length(cfg.toi));
      % index of 'freqoi' value in 'options'
      idx_freqoi = find(ismember(options(1:2:end), 'freqoi'))*2;
      foi = options{idx_freqoi};
      for i_f = 1:length(cfg.foi)
        % collext individual wavelets' responses per frequency
        spec_f = NaN(cfg.order(i_f), nchan, length(cfg.toi));
        opt = options;
        opt{idx_freqoi} = cfg.foi(i_f);
        % compute responses for individual wavelets
        for i_wl = 1:cfg.order(i_f)
          [spec_f(i_wl,:,:), dum, toi] = ft_specest_wavelet(dat, time, 'timeoi', cfg.toi, 'width', cycles{i_f}(i_wl), 'gwidth', cfg.gwidth, opt{:}, 'feedback', fbopt);
        end
        % geometric mean across individual wavelets
        spectrum(:,i_f,:) = prod(spec_f, 1).^(1/cfg.order(i_f));
      end
      clear spec_f
      
      % the following variable is created to keep track of the number of
      % trials per time bin and is needed for proper normalization if
      % keeprpt==1 and the triallength is variable
      if itrial==1, trlcnt = zeros(1, numel(foi), numel(toi)); end
      
      hastime = true;
      % create FAKE ntaper (this requires very minimal code change below for compatibility with the other specest functions)
      ntaper = ones(1,numel(foi));
      % modify spectrum for same reason as fake ntaper
      spectrum = reshape(spectrum,[1 nchan numel(foi) numel(toi)]);
      
    case 'tfr'
      [spectrum,foi,toi] = ft_specest_tfr(dat, time, 'timeoi', cfg.toi, 'width', cfg.width, 'gwidth', cfg.gwidth,options{:}, 'feedback', fbopt);
      
      % the following variable is created to keep track of the number of
      % trials per time bin and is needed for proper normalization if
      % keeprpt==1 and the triallength is variable
      if itrial==1, trlcnt = zeros(1, numel(foi), numel(toi)); end
      
      hastime = true;
      % create FAKE ntaper (this requires very minimal code change below for compatibility with the other specest functions)
      ntaper = ones(1,numel(foi));
      % modify spectrum for same reason as fake ntaper
      spectrum = reshape(spectrum,[1 nchan numel(foi) numel(toi)]);
      
    case 'hilbert'
      [spectrum,foi,toi] = ft_specest_hilbert(dat, time, 'timeoi', cfg.toi, 'width', cfg.width, bpfiltoptions{:}, options{:}, 'feedback', fbopt, 'edgeartnan', cfg.edgeartnan);
      hastime = true;
      % create FAKE ntaper (this requires very minimal code change below for compatibility with the other specest functions)
      ntaper = ones(1,numel(foi));
      % modify spectrum for same reason as fake ntaper
      spectrum = reshape(spectrum,[1 nchan numel(foi) numel(toi)]);
      
    case 'neuvar'
      [spectrum,foi] = ft_specest_neuvar(dat, time, options{:}, 'feedback', fbopt);
      hastime = false;
      % create FAKE ntaper (this requires very minimal code change below for compatibility with the other specest functions)
      ntaper = ones(1,numel(foi));
      
  end % switch
  
  % Set n's
  maxtap = max(ntaper);
  nfoi   = numel(foi);
  if hastime
    ntoi = numel(toi);
  else
    ntoi = 1; % this makes the same code compatible for hastime = false, as time is always the last dimension, and if singleton will disappear
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% Memory allocation
  if strcmp(cfg.method, 'mtmfft') && strcmp(cfg.taper, 'dpss')
    % memory allocation for mtmfft is slightly different because of the possiblity of
    % variable number of tapers over trials (when using dpss), the below exception is
    % made so memory can still be allocated fully (see bug #1025
    trllength = cellfun(@numel,data.time);
    % determine number of tapers per trial
    ntaptrl = sum(floor((2 .* (trllength./data.fsample) .* cfg.tapsmofrq) - 1)); % I floored it for now, because I don't know whether this formula is accurate in all cases, by flooring the memory allocated
    % will most likely be less than it should be, but this would still have the same effect of 'not-crashing-matlabs'.
    % I do have the feeling a round would be 100% accurate, but atm I cannot check this in Percival and Walden
    % - roevdmei
  else
    ntaptrl = ntrials .* maxtap; % the way it used to be in all cases (before bug #1025)
  end
  
  % by default, everything is has the time dimension, if not, some specifics are performed
  if itrial == 1
    % allocate memory to output variables
    if keeprpt == 1 % cfg.keeptrials, 'no' &&  cfg.keeptapers, 'no'
      if powflg, powspctrm     = zeros(nchan,nfoi,ntoi,cfg.precision);             end
      if csdflg, crsspctrm     = complex(zeros(nchancmb,nfoi,ntoi,cfg.precision)); end
      if fftflg, fourierspctrm = complex(zeros(nchan,nfoi,ntoi,cfg.precision));    end
      dimord    = 'chan_freq_time';
    elseif keeprpt == 2 % cfg.keeptrials, 'yes' &&  cfg.keeptapers, 'no'
      if powflg, powspctrm     = nan(ntrials,nchan,nfoi,ntoi,cfg.precision);                                                                 end
      if csdflg, crsspctrm     = complex(nan(ntrials,nchancmb,nfoi,ntoi,cfg.precision),nan(ntrials,nchancmb,nfoi,ntoi,cfg.precision)); end
      if fftflg, fourierspctrm = complex(nan(ntrials,nchan,nfoi,ntoi,cfg.precision),nan(ntrials,nchan,nfoi,ntoi,cfg.precision));       end
      dimord    = 'rpt_chan_freq_time';
    elseif keeprpt == 4 % cfg.keeptrials, 'yes' &&  cfg.keeptapers, 'yes'
      if powflg, powspctrm     = zeros(ntaptrl,nchan,nfoi,ntoi,cfg.precision);        end %
      if csdflg, crsspctrm     = complex(zeros(ntaptrl,nchancmb,nfoi,ntoi,cfg.precision)); end
      if fftflg, fourierspctrm = complex(zeros(ntaptrl,nchan,nfoi,ntoi,cfg.precision));    end
      dimord    = 'rpttap_chan_freq_time';
    end
    if ~hastime
      dimord = dimord(1:end-5); % cut _time
    end
    
    % prepare calcdof
    if strcmp(cfg.calcdof, 'yes')
      if hastime
        dof = zeros(nfoi,ntoi);
        %dof = zeros(ntrials,nfoi,ntoi);
      else
        dof = zeros(nfoi,1);
        %dof = zeros(ntrials,nfoi);
      end
    end
    
    % prepare cumtapcnt
    switch cfg.method %% IMPORTANT, SHOULD WE KEEP THIS SPLIT UP PER METHOD OR GO FOR A GENERAL SOLUTION NOW THAT WE HAVE SPECEST
      case 'mtmconvol'
        cumtapcnt = zeros(ntrials,nfoi);
      case 'mtmfft'
        cumtapcnt = zeros(ntrials,1);
    end
    
  end % itrial==1
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% Create output
  if keeprpt~=4
    
    % mtmconvol is a special case and needs special processing
    if strcmp(cfg.method, 'mtmconvol')
      foiind = ones(1,nfoi);
    else
      % by using this vector below for indexing, the below code does not need to be duplicated for mtmconvol
      foiind = 1:nfoi;
    end
    
    for ifoi = 1:nfoi
      if strcmp(cfg.method, 'mtmconvol')
        spectrum = reshape(permute(spectrum_mtmconvol(:,:,freqtapind{ifoi}),[3 1 2]),[ntaper(ifoi) nchan 1 ntoi]);
      end
      
      % set ingredients for below
      if ~hastime
        acttboi  = 1;
        nacttboi = 1;
      else
        acttboi   = ~all(isnan(spectrum(1,:,foiind(ifoi),:)), 2); % check over all channels, some channels might contain a NaN
        acttboi   = reshape(acttboi, [1 ntoi]);                   % size(spectrum) = [? nchan nfoi ntoi]
        nacttboi = sum(acttboi);
      end
      
      acttap = logical([ones(ntaper(ifoi),1);zeros(size(spectrum,1)-ntaper(ifoi),1)]);
      if powflg
        if strcmp(cfg.method, 'irasa') % ft_specest_irasa outputs power and not amplitude
          powdum = spectrum(acttap,:,foiind(ifoi),acttboi);
        else
          powdum = abs(spectrum(acttap,:,foiind(ifoi),acttboi)) .^2;
        end
        % sinetaper scaling is disabled, because it is not consistent with the other
        % tapers. if scaling is required, please specify cfg.taper =
        % 'sine_old'
        
        %         if isfield(cfg, 'taper') && strcmp(cfg.taper, 'sine')
        %             %sinetapscale = zeros(ntaper(ifoi),nfoi);  % assumes fixed number of tapers
        %             sinetapscale = zeros(ntaper(ifoi),1);  % assumes fixed number of tapers
        %             for isinetap = 1:ntaper(ifoi)  % assumes fixed number of tapers
        %               sinetapscale(isinetap,:) = (1 - (((isinetap - 1) ./ ntaper(ifoi)) .^ 2));
        %             end
        %             sinetapscale = reshape(repmat(sinetapscale,[1 1 nchan ntoi]),[ntaper(ifoi) nchan 1 ntoi]);
        %             powdum = powdum .* sinetapscale;
        %           end
      end
      if fftflg
        fourierdum = spectrum(acttap,:,foiind(ifoi),acttboi);
      end
      if csdflg
        csddum = spectrum(acttap,cutdatindcmb(:,1),foiind(ifoi),acttboi) .* conj(spectrum(acttap,cutdatindcmb(:,2),foiind(ifoi),acttboi));
      end
      
      % switch between keep's
      switch keeprpt
        
        case 1 % cfg.keeptrials, 'no' &&  cfg.keeptapers, 'no'
          if ~isempty(trlcnt)
            trlcnt(1, ifoi, :) = trlcnt(1, ifoi, :) + shiftdim(double(acttboi(:)'),-1);
          end
          
          if powflg
            powspctrm(:,ifoi,acttboi) = powspctrm(:,ifoi,acttboi) + (reshape(mean(powdum,1),[nchan 1 nacttboi]) ./ ntrials);
            %powspctrm(:,ifoi,~acttboi) = NaN;
          end
          if fftflg
            fourierspctrm(:,ifoi,acttboi) = fourierspctrm(:,ifoi,acttboi) + (reshape(mean(fourierdum,1),[nchan 1 nacttboi]) ./ ntrials);
            %fourierspctrm(:,ifoi,~acttboi) = NaN;
          end
          if csdflg
            crsspctrm(:,ifoi,acttboi) = crsspctrm(:,ifoi,acttboi) + (reshape(mean(csddum,1),[nchancmb 1 nacttboi]) ./ ntrials);
            %crsspctrm(:,ifoi,~acttboi) = NaN;
          end
          
        case 2 % cfg.keeptrials, 'yes' &&  cfg.keeptapers, 'no'
          if powflg
            powspctrm(itrial,:,ifoi,acttboi) = reshape(mean(powdum,1),[nchan 1 nacttboi]);
            powspctrm(itrial,:,ifoi,~acttboi) = NaN;
          end
          if fftflg
            fourierspctrm(itrial,:,ifoi,acttboi) = reshape(mean(fourierdum,1), [nchan 1 nacttboi]);
            fourierspctrm(itrial,:,ifoi,~acttboi) = NaN;
          end
          if csdflg
            crsspctrm(itrial,:,ifoi,acttboi) = reshape(mean(csddum,1), [nchancmb 1 nacttboi]);
            crsspctrm(itrial,:,ifoi,~acttboi) = NaN;
          end
          
      end % switch keeprpt
      
      % do calcdof  dof = zeros(numper,numfoi,numtoi);
      if strcmp(cfg.calcdof, 'yes')
        if hastime
          acttimboiind = ~all(isnan(spectrum(1,:,foiind(ifoi),:)), 2); % check over all channels, some channels might contain a NaN
          acttimboiind = reshape(acttimboiind, [1 ntoi]);
          dof(ifoi,acttimboiind) = ntaper(ifoi) + dof(ifoi,acttimboiind);
        else % hastime = false
          dof(ifoi) = ntaper(ifoi) + dof(ifoi);
        end
      end
    end %ifoi
    
  else
    % keep tapers
    if ~exist('tapcounter', 'var')
      tapcounter = 0;
    end
    
    if strcmp(cfg.method, 'mtmconvol')
      spectrum = permute(reshape(spectrum_mtmconvol,[nchan ntoi ntaper(1) nfoi]),[3 1 4 2]);
    end
    
    currrptind  = tapcounter + (1:maxtap);
    tapcounter  = currrptind(end);
    %rptind = reshape(1:ntrials .* maxtap,[maxtap ntrials]);
    %currrptind = rptind(:,itrial);
    if powflg
      if strcmp(cfg.method, 'irasa') % ft_specest_irasa outputs power and not amplitude
        powspctrm(currrptind,:,:) = spectrum;
      else
        powspctrm(currrptind,:,:) = abs(spectrum).^2;
      end
    end
    if fftflg
      fourierspctrm(currrptind,:,:,:) = spectrum;
    end
    if csdflg
      crsspctrm(currrptind,:,:,:) =          spectrum(cutdatindcmb(:,1),:,:) .* ...
        conj(spectrum(cutdatindcmb(:,2),:,:));
    end
    
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  % set cumptapcnt
  switch cfg.method %% IMPORTANT, SHOULD WE KEEP THIS SPLIT UP PER METHOD OR GO FOR A GENERAL SOLUTION NOW THAT WE HAVE SPECEST
    case {'mtmconvol' 'wavelet'}
      cumtapcnt(itrial,:) = ntaper;
    case 'mtmfft'
      cumtapcnt(itrial,1) = ntaper(1); % fixed number of tapers? for the moment, yes, as specest_mtmfft computes only one set of tapers
  end
  
end % for ntrials
ft_progress('close');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% END: Main loop over trials
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% re-normalise the TFRs if keeprpt==1
if (strcmp(cfg.method, 'mtmconvol') || strcmp(cfg.method, 'wavelet')) && keeprpt==1
  nanmask = trlcnt==0;
  if powflg
    powspctrm = powspctrm.*ntrials;
    powspctrm = powspctrm./trlcnt(ones(size(powspctrm,1),1),:,:);
    powspctrm(nanmask(ones(size(powspctrm,1),1),:,:)) = nan;
  end
  if fftflg
    fourierspctrm = fourierspctrm.*ntrials;
    fourierspctrm = fourierspctrm./trlcnt(ones(size(fourierspctrm,1),1),:,:);
    fourierspctrm(nanmask(ones(size(fourierspctrm,1),1),:,:)) = nan;
  end
  if csdflg
    crsspctrm = crsspctrm.*ntrials;
    crsspctrm = crsspctrm./trlcnt(ones(size(crsspctrm,1),1),:,:);
    crsspctrm(nanmask(ones(size(crsspctrm,1),1),:,:)) = nan;
  end
end

% set output variables
freq        = [];
freq.label  = data.label;
freq.dimord = dimord;
freq.freq   = foi;
hasdc       = find(foi==0);
hasnyq      = find(foi==data.fsample./2);
hasdc_nyq   = [hasdc hasnyq];
if exist('toi', 'var')
  freq.time = toi;
end
if powflg
  % correct the 0 Hz or Nyqist bin if present, scaling with a factor of 2 is only appropriate for ~0 Hz
  if ~isempty(hasdc_nyq)
    if keeprpt>1
      powspctrm(:,:,hasdc_nyq,:) = powspctrm(:,:,hasdc_nyq,:)./2;
    else
      powspctrm(:,hasdc_nyq,:) = powspctrm(:,hasdc_nyq,:)./2;
    end
  end
  
  if startsWith(cfg.output, 'fooof')
    % check for brainstorm functions on the path, and add if needed
    ft_hastoolbox('brainstorm', 1);
    
    TF(:,1,:) = powspctrm;
    Freqs     = freq.freq;
    Freqs(Freqs==0) = [];
    % This grabs the defaults from the brainstorm code
    opts_bst  = getfield(process_fooof('GetDescription'), 'options');
    
    % Fetch user settings, this is a chunk of code copied over from
    % process_fooof, to bypass the whole database etc handling.
    opt                     = ft_getopt(cfg, 'fooof', []);
    opt.freq_range          = ft_getopt(opt, 'freq_range', Freqs([1 end]));
    opt.peak_width_limits   = ft_getopt(opt, 'peak_width_limits', opts_bst.peakwidth.Value{1});
    opt.max_peaks           = ft_getopt(opt, 'max_peaks',         opts_bst.maxpeaks.Value{1});
    opt.min_peak_height     = ft_getopt(opt, 'min_peak_height',   opts_bst.minpeakheight.Value{1}/10); % convert from dB to B
    opt.aperiodic_mode      = ft_getopt(opt, 'aperiodic_mode',    opts_bst.apermode.Value);
    opt.peak_threshold      = ft_getopt(opt, 'peak_threshold',    2);   % 2 std dev: parameter for interface simplification
    opt.return_spectrum     = ft_getopt(opt, 'return_spectrum',   1);   % SPM/FT: set to 1
    opt.border_threshold    = ft_getopt(opt, 'border_threshold',  1);   % 1 std dev: proximity to edge of spectrum, static in Python 
    % Matlab-only options
    opt.power_line          = ft_getopt(opt, 'power_line',        '50'); % for some reason it should be a string, if you don't want a notch, use 'inf'. Brainstorm's default is '60'
    opt.peak_type           = ft_getopt(opt, 'peak_type',         opts_bst.peaktype.Value);
    opt.proximity_threshold = ft_getopt(opt, 'proximity_threshold', opts_bst.proxthresh.Value{1});
    opt.guess_weight        = ft_getopt(opt, 'guess_weight',      opts_bst.guessweight.Value);
    opt.thresh_after        = ft_getopt(opt, 'thresh_after',      true);   % Threshold after fitting always selected for Matlab (mirrors the Python FOOOF closest by removing peaks that do not satisfy a user's predetermined conditions)
    
    % Output options
    opt.sort_type  = opts_bst.sorttype.Value;
    opt.sort_param = opts_bst.sortparam.Value;
	  opt.sort_bands = opts_bst.sortbands.Value;

    % Check input frequency bounds
    if (any(opt.freq_range < 0) || opt.freq_range(1) >= opt.freq_range(2))
      bst_report('error','Invalid Frequency range');
      return
    end
    
    hasOptimTools = 0;
    if exist('fmincon', 'file')
      hasOptimTools = 1;
      disp('Using constrained optimization, Guess Weight ignored.')
    end
    
    [fs, fg] = process_fooof('FOOOF_matlab', TF, freq.freq, opt, hasOptimTools);
    
    % add the options back to the cfg
    cfg.fooof = opt;
    
    switch cfg.output
      case 'fooof'
        powspctrm_f = cat(1, fg.fooofed_spectrum);
      case 'fooof_peaks'
        powspctrm_f = cat(1, fg.peak_fit);
      case 'fooof_aperiodic'
        powspctrm_f = cat(1, fg.ap_fit);
    end
    fg = removefields(fg, {'fooofed_spectrum', 'peak_fit', 'ap_fit'});
    
    for k = 1:size(powspctrm_f,1)
      powspctrm(k,:) = interp1(fs, powspctrm_f(k,:), freq.freq, 'linear', nan);
      fg(k).label    = freq.label{k};
    end
    freq.powspctrm   = powspctrm;
    freq.fooofparams = fg(:);
    
  else
    freq.powspctrm = powspctrm;
  end
end
if fftflg
  % correct the 0 Hz or Nyqist bin if present, scaling with a factor of 2 is only appropriate for ~0 Hz
  if ~isempty(hasdc_nyq)
    if keeprpt>1
      fourierspctrm(:,:,hasdc_nyq,:) = fourierspctrm(:,:,hasdc_nyq,:)./sqrt(2);
    else
      fourierspctrm(:,hasdc_nyq,:) = fourierspctrm(:,hasdc_nyq,:)./sqrt(2);
    end
  end
  freq.fourierspctrm = fourierspctrm;
end
if csdflg
  % correct the 0 Hz or Nyqist bin if present, scaling with a factor of 2 is only appropriate for ~0 Hz
  if ~isempty(hasdc_nyq)
    if keeprpt>1
      crsspctrm(:,:,hasdc_nyq,:) = crsspctrm(:,:,hasdc_nyq,:)./2;
    else
      crsspctrm(:,hasdc_nyq,:) = crsspctrm(:,hasdc_nyq,:)./2;
    end
  end
  freq.labelcmb  = cfg.channelcmb;
  freq.crsspctrm = crsspctrm;
end
if strcmp(cfg.calcdof, 'yes')
  freq.dof = 2 .* dof;
end
if strcmp(cfg.method, 'mtmfft') && (keeprpt == 2 || keeprpt == 4)
  freq.cumsumcnt = trllength';
end
if exist('cumtapcnt', 'var') && (keeprpt == 2 || keeprpt == 4)
  freq.cumtapcnt = cumtapcnt;
end

% backwards compatability of foilim
if ~isempty(cfg.foilim)
  cfg = rmfield(cfg, 'foi');
else
  cfg = rmfield(cfg, 'foilim');
end

% some fields from the input should always be copied over in the output
freq = copyfields(data, freq, {'grad', 'elec', 'opto', 'topo', 'topolabel', 'unmixing'});

if isfield(data, 'trialinfo') && strcmp(cfg.keeptrials, 'yes')
  % copy the trialinfo into the output, but not the sampleinfo
  freq.trialinfo = data.trialinfo;
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous   data
ft_postamble provenance freq
ft_postamble history    freq
ft_postamble savevar    freq