function [spectrum, ntaper, freqoi] = ft_specest_irasa(dat, time, varargin)

% FT_SPECEST_IRASA separates the fractal components from the orginal power spectrum
% using Irregular-Resampling Auto-Spectral Analysis (IRASA)
%
% Use as
%   [spectrum, ntaper, freqoi] = ft_specest_irasa(dat, time, ...)
% where the input arguments are
%   dat        = matrix of chan*sample
%   time       = vector, containing time in seconds for each sample
% and the output arguments are
%   spectrum   = matrix of taper*chan*freqoi of fourier coefficients
%   ntaper     = vector containing number of tapers per element of freqoi
%   freqoi     = vector of frequencies in spectrum
%
% Optional arguments should be specified in key-value pairs and can include
%   freqoi     = vector, containing frequencies of interest
%   output     = string, indicating type of output ('fractal' or 'original', default 'fractal')
%   pad        = number, total length of data after zero padding (in seconds)
%   padtype    = string, indicating type of padding to be used, can be 'zero', 'mean', 'localmean', 'edge', or 'mirror' (default = 'zero')
%   polyorder  = number, the order of the polynomial to fitted to and removed from the data prior to the Fourier transform (default = 0, which removes the DC-component)
%   verbose    = boolean, output progress to console (0 or 1, default 1)
%
% This implements: Wen.H. & Liu.Z.(2016), Separating fractal and oscillatory components in the power
% spectrum of neurophysiological signal. Brain Topogr. 29(1):13-26. The source code accompanying the
% original paper is avaible from https://purr.purdue.edu/publications/1987/1
%
% For more information about the difference between the current and previous version and how to use this
% function, please see https://www.fieldtriptoolbox.org/example/irasa/
%
% See also FT_FREQANALYSIS, FT_SPECEST_MTMFFT, FT_SPECEST_MTMCONVOL, FT_SPECEST_TFR, FT_SPECEST_HILBERT, FT_SPECEST_WAVELET

% Copyright (C) 2019-2020, Rui Liu, Arjen Stolk
%
% 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/>.
%
% $Id$

% these are for speeding up computation of tapers on subsequent calls
persistent previous_argin previous_tap

% get the optional input arguments
freqoi    = ft_getopt(varargin, 'freqoi', 'all');
output    = ft_getopt(varargin, 'output', 'fractal');
pad       = ft_getopt(varargin, 'pad');
padtype   = ft_getopt(varargin, 'padtype', 'zero');
polyorder = ft_getopt(varargin, 'polyorder', 1);
verbose   = ft_getopt(varargin, 'verbose', true);
fbopt     = ft_getopt(varargin, 'feedback');

verbose = istrue(verbose); % if the calling function has 'yes'/'no'/etc

% check output option
if ~strcmp(output, {'fractal','original'})
  ft_error('The current version ft_specest_irasa outputs ''fractal'' or ''original'' power only. For more information about the update, see https://www.fieldtriptoolbox.org/example/irasa/');
end

% this does not work on integer data
dat = cast(dat, 'double');
if ~isa(dat, 'double') && ~isa(dat, 'single')
  dat = cast(dat, 'double');
end

% dat param
[nchan,ndatsample] = size(dat);
fsample = 1./mean(diff(time));

% subset param
subset_nsample = 2^floor(log2(ndatsample*0.9)); % the number of sub-segments in samples is the power of 2 that does not exceed 90% of input data
subset_num     = 10;                            % the number of sub-segments
subset_dist    = floor((ndatsample  - subset_nsample)/(subset_num - 1)); % distance between sub-segements in sample, to evenly distribute the sub-subsegments within the total length of input data

% resampling ratio
hset  = 1.1:0.05:1.9;
nhset = length(hset);

% remove polynomial fit from the data -> default is demeaning
if polyorder >= 0
  dat = ft_preproc_polyremoval(dat, polyorder, 1, ndatsample);
end

% zero padding
if isempty(pad)
  % if no padding is specified this is set equal to the current data length
  pad = subset_nsample/fsample;
end
if round(pad * fsample) < subset_nsample
  ft_error('the padding that you specified is shorter than the data');
end
endnsample = round(pad * fsample);  % total number of samples of padded data
endtime    = pad;                   % total time in seconds of padded data

% set freqboi and freqoi
freqoiinput = freqoi;
if isnumeric(freqoi) % if input is a vector
  freqboi   = round(freqoi ./ (fsample ./ endnsample)) + 1; % is equivalent to: round(freqoi .* endtime) + 1;
  freqboi   = unique(freqboi);
  freqoi    = (freqboi-1) ./ endtime; % boi - 1 because 0 Hz is included in fourier output
elseif strcmp(freqoi, 'all') % if input was 'all'
  freqboilim = round([0 fsample/2] ./ (fsample ./ endnsample)) + 1;
  freqboi    = freqboilim(1):1:freqboilim(2);
  freqoi     = (freqboi-1) ./ endtime;
end
nfreqboi = length(freqboi);
nfreqoi  = length(freqoi);

% throw a warning if input freqoi is different from output freqoi
if isnumeric(freqoiinput)
  % check whether padding is appropriate for the requested frequency resolution
  rayl = 1/endtime;
  if any(rem(freqoiinput,rayl)) % not always the case when they mismatch
    ft_warning('padding not sufficient for requested frequency resolution, for more information please see the FAQs on www.ru.nl/neuroimaging/fieldtrip');
  end
  if numel(freqoiinput) ~= numel(freqoi) % freqoi will not contain double frequency bins when requested
    ft_warning('output frequencies are different from input frequencies, multiples of the same bin were requested but not given');
  else
    if any(abs(freqoiinput-freqoi) >= eps*1e6)
      ft_warning('output frequencies are different from input frequencies');
    end
  end
end

% determine whether tapers need to be recomputed
current_argin = {output, time, endnsample, freqoi}; % reasoning: if time and endnsample are equal, it's the same length trial, if the rest is equal then the requested output is equal
if isequal(current_argin, previous_argin)
  % don't recompute tapers
  tap = previous_tap;
else
  if strcmp(output, 'fractal')
    % (re)compute tapers, 1:mid are upsample tapers, mid+1:end are downsample tapers
    tap = cell(nhset*2,1);
    for ih = 1:nhset
      [n, d] = rat(hset(ih)); % n > d
      tmp = hanning(size(resample(zeros(subset_nsample,nchan), n, d),1))';
      tap{ih,1} = tmp./norm(tmp, 'fro');% for upsampled subsets
      tmp = hanning(size(resample(zeros(subset_nsample,nchan), d, n),1))';
      tap{ih+nhset,1} = tmp./norm(tmp, 'fro');% for downsampled subsets
    end
  elseif strcmp(output, 'original')
    tap = hanning(subset_nsample)';
    tap = tap./norm(tap, 'fro');
  end
end

% set ntaper
if strcmp(output, 'fractal')
  ntaper = repmat(size(tap,2),nfreqoi,1); % pretend there's only one taper
elseif strcmp(output, 'original')
  ntaper = repmat(size(tap,1),nfreqoi,1);
end

% feedback of computing progress
if isempty(fbopt)
  fbopt.i = 1;
  fbopt.n = 1;
end
str = sprintf('nfft: %d samples, datalength: %d samples, %d tapers',endnsample,ndatsample,ntaper(1));
st  = dbstack;
if length(st)>1 && strcmp(st(2).name, 'ft_freqanalysis')
  % specest_mtmfft has been called by ft_freqanalysis, meaning that ft_progress has been initialised
  ft_progress(fbopt.i./fbopt.n, ['processing trial %d/%d ',str,'\n'], fbopt.i, fbopt.n);
elseif verbose
  fprintf([str, '\n']);
end

% compute irasa or fft
spectrum = cell(ntaper(1),1);
for itap = 1:ntaper(1)
  %%%%%%%% IRASA %%%%%%%%%%
  if strcmp(output,'fractal')
    pow  = zeros(nchan,nfreqboi,nhset);
    for ih = 1:nhset % loop resampling ratios hset
      upow = zeros(nchan,nfreqboi);
      dpow = zeros(nchan,nfreqboi);
      [n, d] = rat(hset(ih)); % n > d
      for k = 0 : subset_num-1 % loop #subset
        % pair-wised resampling
        subset_start = subset_dist*k + 1;
        subset_end = subset_start+subset_nsample-1;
        subdat = dat(:,subset_start:subset_end);
        udat = resample(subdat', n, d)'; % upsample
        ddat = resample(subdat', d, n)'; % downsample
        
        % compute auto-power of resampled data
        % upsampled
        postpad = ceil(round(pad * fsample) - size(udat,2));
        tmp = fft(ft_preproc_padding(bsxfun(@times,udat,tap{ih,1}), padtype, 0, postpad), endnsample, 2);
        tmp = tmp(:,freqboi);
        tmp = tmp .* sqrt(2 ./ endnsample);
        tmp = abs(tmp).^2;
        upow = upow + tmp;
        % downsampled
        postpad = ceil(round(pad * fsample) - size(ddat,2));
        tmp = fft(ft_preproc_padding(bsxfun(@times,ddat,tap{ih+nhset,1}), padtype, 0, postpad), endnsample, 2);
        tmp = tmp(:,freqboi);
        tmp = tmp .* sqrt(2 ./ endnsample);
        tmp = abs(tmp).^2;
        dpow = dpow + tmp;
      end % loop #subset
      
      % average across subsets for noise filtering
      upow = upow / subset_num;
      dpow = dpow / subset_num;
      
      % geometric mean for relocating oscillatory component
      pow(:,:,ih) = sqrt(upow.*dpow);
    end % loop resampling ratios hset
    
    % median across resampling factors
    spectrum{itap} = median(pow,3);
    
    %%%%%%%% FFT %%%%%%%%%%
  elseif strcmp(output,'original')
    pow  = zeros(nchan,nfreqboi);
    for k = 0 : subset_num-1 % loop #subset
      % pair-wised resampling
      subset_start = subset_dist*k + 1;
      subset_end = subset_start+subset_nsample-1;
      subdat = dat(:,subset_start:subset_end);
      % compute auto-power
      postpad = ceil(round(pad * fsample) - size(subdat,2));
      tmp = fft(ft_preproc_padding(bsxfun(@times,subdat,tap(itap,:)), padtype, 0, postpad), endnsample, 2);
      tmp = tmp(:,freqboi);
      tmp = tmp .* sqrt(2 ./ endnsample);
      tmp = abs(tmp).^2;
      pow = pow + tmp;
    end % loop #subset
    
    % average across subsets for noise filtering
    spectrum{itap} = pow / subset_num;
    
  end % if output
end % for taper

spectrum = reshape(vertcat(spectrum{:}),[nchan ntaper(1) nfreqboi]); % collecting in a cell-array and later reshaping provides significant speedups
spectrum = permute(spectrum, [2 1 3]);

% remember the current input arguments, so that they can be
% reused on a subsequent call in case the same input argument is given
previous_argin = current_argin;
previous_tap   = tap;
