ft_specest_mtmfft.m

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

% FT_SPECEST_MTMFFT computes a fast Fourier transform using multitapering with
% multiple tapers from the DPSS sequence or using a variety of single tapers.
%
% Use as
%   [spectrum, ntaper, freqoi, se] = ft_specest_mtmfft(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 ntaper*nchan*nfreq of fourier coefficients
%   ntaper     = vector containing number of tapers per element of freqoi
%   freqoi     = vector of frequencies in spectrum
%   se         = (only with dpss-tapers) effective degrees of freedom
%
% Optional arguments should be specified in key-value pairs and can include
%   freqoi     = vector, containing frequencies of interest
%   taper      = 'dpss', 'hanning' or many others, see WINDOW (default = 'dpss')
%   taperopt   = additional taper options to be used in the WINDOW function, see WINDOW
%   weightopt  = 'mean', 'eig', 'adapt' (default = 'mean'), taper weights for dpss tapers, see ADAPTSPEC_DPSS
%   tapsmofrq  = the amount of spectral smoothing through multi-tapering. Note: 4 Hz smoothing means plus-minus 4 Hz, i.e. a 8 Hz smoothing box
%   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')
%   dimord     = 'tap_chan_freq' (default) or 'chan_time_freqtap' for memory efficiency (only when using variable number of slepian tapers)
%   polyorder  = number, the order of the polynomial to fitted to and removed from the data prior to the fourier transform (default = 0 -> remove DC-component)
%   verbose    = output progress to console (0 or 1, default 1)
%
% See also FT_FREQANALYSIS, FT_SPECEST_MTMCONVOL, FT_SPECEST_TFR, FT_SPECEST_HILBERT, FT_SPECEST_WAVELET

% Copyright (C) 2010, Donders Institute for Brain, Cognition and Behaviour
%
% 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 previous_w

% get the optional input arguments
taper     = ft_getopt(varargin, 'taper'); if isempty(taper), ft_error('You must specify a taper'); end
pad       = ft_getopt(varargin, 'pad');
padtype   = ft_getopt(varargin, 'padtype', 'zero');
freqoi    = ft_getopt(varargin, 'freqoi', 'all');
tapsmofrq = ft_getopt(varargin, 'tapsmofrq');
dimord    = ft_getopt(varargin, 'dimord', 'tap_chan_freq');
fbopt     = ft_getopt(varargin, 'feedback');
verbose   = ft_getopt(varargin, 'verbose', true);
polyorder = ft_getopt(varargin, 'polyorder', 0);
tapopt    = ft_getopt(varargin, 'taperopt');
weightopt = ft_getopt(varargin, 'weightopt', 'mean');

if isempty(fbopt)
  fbopt.i = 1;
  fbopt.n = 1;
end

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

% throw errors for required input
if isempty(tapsmofrq) && (strcmp(taper, 'dpss') || strcmp(taper, 'sine'))
  ft_error('you need to specify tapsmofrq when using dpss or sine tapers')
end

if strcmp(taper, 'dpss')
  switch weightopt
    case 'mean'
      adaptflag = 0;
    case 'eig'
      adaptflag = 1;
    case 'adapt'
      adaptflag = 2;
    otherwise
      ft_error('unknown weightopt specified for dpss-tapering');
  end
elseif ~strcmp(weightopt, 'mean')
  ft_warning('no dpss tapers are specified, so the value %s of weightopt does not have an effect', weightopt);
end

% this does not work on integer data
dat = cast(dat, 'double');

% Set n's
[nchan,ndatsample] = size(dat);

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

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

% Determine fsample and set total time-length of data
fsample = 1./mean(diff(time));
dattime = ndatsample / fsample; % total time in seconds of input data

% Zero padding
if round(pad * fsample) < ndatsample
  ft_error('the padding that you specified is shorter than the data');
end
if isempty(pad) % if no padding is specified padding is equal to current data length
  pad = dattime;
end
endnsample = round(pad * fsample);    % total number of samples of padded data
postpad    = endnsample - ndatsample; % number of samples for zero padding
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);
if (strcmp(taper, 'dpss') || strcmp(taper, 'sine')) && numel(tapsmofrq)~=1 && (numel(tapsmofrq)~=nfreqoi)
  ft_error('tapsmofrq needs to contain a smoothing parameter for every frequency when requesting variable number of slepian tapers')
end

% 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*1e9)
      ft_warning('output frequencies are different from input frequencies');
    end
  end
end

% determine whether tapers need to be recomputed
current_argin = {time, postpad, taper, tapsmofrq, freqoi, tapopt, dimord}; % reasoning: if time and postpad 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;
  w   = previous_w;
else
  % recompute tapers
  w = [];
  switch taper
    
    case 'dpss'
      if isscalar(tapsmofrq)
        % create a sequence of DPSS tapers, ensure that the input arguments are double precision
        [tap, w] = double_dpss(ndatsample,ndatsample*(tapsmofrq./fsample));

        % remove the last taper because the last slepian taper has poor spectral concentration properties
        tap = tap(:, 1:(end-1))';
        w   = w(1:(end-1));

        % give error/warning about number of tapers
        if isempty(tap)
          ft_error('datalength to short for specified smoothing\ndatalength: %.3f s, smoothing: %.3f Hz, minimum smoothing: %.3f Hz',ndatsample/fsample,tapsmofrq,fsample/ndatsample);
        elseif size(tap,1) == 1
          ft_warning('using only the first Slepian taper for the requested smoothing');
        end
      elseif numel(tapsmofrq)>1
        tap = cell(1,nfreqoi);
        for ifreqoi = 1:nfreqoi
          % create a sequence of DPSS tapers, ensure that the input arguments are double precision
          currtap = double_dpss(ndatsample, ndatsample .* (tapsmofrq(ifreqoi) ./ fsample))';
          % remove the last taper because the last slepian taper is always messy
          currtap = currtap(1:(end-1), :);
          
          % give error/warning about number of tapers
          if isempty(currtap)
            ft_error('%.3f Hz: datalength to short for specified smoothing\ndatalength: %.3f s, smoothing: %.3f Hz, minimum smoothing: %.3f Hz',freqoi(ifreqoi), ndatsample/fsample,tapsmofrq(ifreqoi),fsample/ndatsample(ifreqoi));
          elseif size(currtap,1) == 1
            disp([num2str(freqoi(ifreqoi)) ' Hz: WARNING: using only the first Slepian taper for the requested smoothing'])
          end
          tap{ifreqoi} = currtap;
        end
      end
      
    case 'sine'
      if isscalar(tapsmofrq)
        % create a sequence of sine tapers, 
        tap = sine_taper(ndatsample, ndatsample*(tapsmofrq./fsample))';
        % remove the last taper 
        tap = tap(1:(end-1), :);
        
        % give error/warning about number of tapers
        if isempty(tap)
          ft_error('datalength to short for specified smoothing\ndatalength: %.3f s, smoothing: %.3f Hz, minimum smoothing: %.3f Hz',ndatsample/fsample,tapsmofrq,fsample/ndatsample);
        elseif size(tap,1) == 1
          ft_warning('using only one taper for specified smoothing');
        end
      elseif numel(tapsmofrq)>1
        tap = cell(1,nfreqoi);
        for ifreqoi = 1:nfreqoi
          % create a sequence of sine tapers
          currtap = sine_taper(ndatsample, ndatsample .* (tapsmofrq(ifreqoi) ./ fsample))';
          % remove the last taper because the last slepian taper is always messy
          currtap = currtap(1:(end-1), :);
          
          % give error/warning about number of tapers
          if isempty(currtap)
            ft_error('%.3f Hz: datalength to short for specified smoothing\ndatalength: %.3f s, smoothing: %.3f Hz, minimum smoothing: %.3f Hz',freqoi(ifreqoi), ndatsample/fsample,tapsmofrq(ifreqoi),fsample/ndatsample(ifreqoi));
          elseif size(currtap,1) == 1
            disp([num2str(freqoi(ifreqoi)) ' Hz: WARNING: using only one taper for specified smoothing'])
          end
          tap{ifreqoi} = currtap;
        end
      end          
      
    case 'sine_old'
      % to provide compatibility with the tapers being scaled (which was default
      % behavior prior to 29apr2011) yet this gave different magnitude of power
      % when comparing with slepian multi tapers
      tap = sine_taper_scaled(ndatsample, ndatsample*(tapsmofrq./fsample))';
      tap = tap(1:(end-1), :); % remove the last taper
      
    case 'alpha'
      tap = cell(1,nfreqoi);
      for ifreqoi = 1:nfreqoi
        tap{ifreqoi} = alpha_taper(ndatsample, freqoi(ifreqoi)./ fsample)';
        tap{ifreqoi} = tap{ifreqoi}./norm(tap{ifreqoi}, 'fro');
      end

    case 'hanning'
      tap = hanning(ndatsample)';
      tap = tap./norm(tap, 'fro');
      
    otherwise
      % create the taper and ensure that it is normalized
      if isempty(tapopt) % some windowing functions don't support nargin>1, and window.m doesn't check it
        tap = window(taper, ndatsample)';
      else
        tap = window(taper, ndatsample, tapopt)';
      end
      tap = tap ./ norm(tap,'fro');
      
  end % switch taper
end % isequal currargin

% set ntaper
if ~((strcmp(taper,'dpss') || strcmp(taper,'sine')) && numel(tapsmofrq)>1) % variable number of slepian tapers not requested
  ntaper = repmat(size(tap,1),1,nfreqoi);
else % variable number of slepian tapers requested
  ntaper = cellfun(@size,tap,repmat({1},[1 nfreqoi]));
end

% determine phase-shift so that for all frequencies angle(t=0) = 0
timedelay = time(1);
if timedelay ~= 0
  angletransform = complex(zeros(1,nfreqoi));
  for ifreqoi = 1:nfreqoi
    missedsamples = round(timedelay * fsample);
    % determine angle of freqoi if oscillation started at 0
    % the angle of wavelet(cos,sin) = 0 at the first point of a cycle, with sine being in upgoing flank, which is the same convention as in mtmconvol
    anglein = (missedsamples) .* ((2.*pi./fsample) .* freqoi(ifreqoi));
    coswav  = cos(anglein);
    sinwav  = sin(anglein);
    angletransform(ifreqoi) = atan2(sinwav, coswav);
  end
  angletransform = repmat(angletransform,[nchan,1]);
end

% compute fft
st = dbstack;
if ~(ismember(taper, {'dpss' 'sine'}) && numel(tapsmofrq)>1)  && ~isequal(taper, 'alpha') % variable number of slepian/sine tapers not requested, or alpha taper
  str = sprintf('nfft: %d samples, datalength: %d samples, %d tapers',endnsample,ndatsample,ntaper(1));
  if length(st)>1 && strcmp(st(2).name, 'ft_freqanalysis')
    % ft_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
  
  siz = [size(dat,1), size(dat,2)+postpad, ntaper(1)];
  dum = complex(zeros(siz), zeros(siz));
  for itap = 1:ntaper(1)
    % fft of zero-padded tapered data segment
    dum(:,:,itap) = fft(ft_preproc_padding(bsxfun(@times,dat,tap(itap,:)), padtype, 0, postpad),[], 2);
  end 
  
  if strcmp(taper, 'dpss')
    % compute taper weights
    [spec, se, wt] = adaptspec_dpss(dum, w, adaptflag);
    wt  = wt(:,freqboi,:);
    se  = se(:,freqboi);
    if all(se==se(1))
      se = se(1,:);
    end
  else
    wt = 1;
    se = 2;
  end
  dum = dum(:,freqboi,:);
  
  % weight the spectra
  dum = dum .* (wt./sqrt(mean(wt.^2,3)));

  % phase-shift according to above angles
  if timedelay ~= 0
    dum = dum .* exp(-1i*angletransform(:,:,ones(1,ntaper(1))));
  end
  dum = dum .* sqrt(2 ./ endnsample);
  spectrum = permute(dum, [3 1 2]);
  wt       = permute(wt,  [3 1 2]);
  
  %spectrum{itap} = dum;
  %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]);
  
else % variable number of slepian tapers requested
  switch dimord
    
    case 'tap_chan_freq' % default
      % start fft'ing
      spectrum = complex(nan([max(ntaper) nchan nfreqoi]));
      for ifreqoi = 1:nfreqoi
        str = sprintf('nfft: %d samples, datalength: %d samples, frequency %d (%.2f Hz), %d tapers',endnsample,ndatsample,ifreqoi,freqoi(ifreqoi),ntaper(ifreqoi));
        if length(st)>1 && strcmp(st(2).name, 'ft_freqanalysis') && verbose
          % ft_specest_mtmfft has been called by ft_freqanalysis, meaning that ft_progress has been initialised
          ft_progress(fbopt.i./fbopt.n, ['processing trial %d, ',str,'\n'], fbopt.i);
        elseif verbose
          fprintf([str, '\n']);
        end

        siz = [size(dat,1), size(dat,2)+postpad, ntaper(1)];
        dum = complex(zeros(siz), zeros(siz));
        for itap = 1:ntaper(ifreqoi)
          dum(:,:,itap) = fft(ft_preproc_padding(bsxfun(@times,dat,tap{ifreqoi}(itap,:)), padtype, 0, postpad), [], 2);
        end

        if strcmp(taper, 'dpss')
          if adaptflag~=0
            ft_warning('adaptively weighted multitapering with variable numbers of tapers across frequencies is performed at your own risk');
          end
          [spec, se, wt] = adaptspec_dpss(dum, w, adaptflag);
          wt  = wt(:,freqboi(ifreqoi),:);
          se  = se(:,freqboi(ifreqoi));
        else
          wt = 1;
          se = 2;
        end
        dum = dum(:,freqboi(ifreqoi),:);

        % weight the spectra
        dum = dum .* (wt./sqrt(mean(wt.^2,3)));

        % phase-shift according to above angles
        if timedelay ~= 0
          dum = dum .* exp(-1i*angletransform(:,ifreqoi,ones(1,ntaper(ifreqoi))));
        end
        dum = dum .* sqrt(2 ./ endnsample);

        spectrum(:,:,ifreqoi) = permute(dum, [3 1 2]); 

      end % for nfreqoi
      
    case 'chan_freqtap' % memory efficient representation
      if strcmp(taper, 'dpss') && adaptflag~=0
        ft_error('adaptively weighted multitapering is not possible with the specified requested dimord');
      end
      
      % create tapfreqind
      freqtapind = cell(1,nfreqoi);
      tempntaper = [0; cumsum(ntaper(:))];
      for ifreqoi = 1:nfreqoi
        freqtapind{ifreqoi} = tempntaper(ifreqoi)+1:tempntaper(ifreqoi+1);
      end
      
      % start fft'ing
      spectrum = complex(zeros([nchan sum(ntaper)]));
      for ifreqoi = 1:nfreqoi
        str = sprintf('nfft: %d samples, datalength: %d samples, frequency %d (%.2f Hz), %d tapers',endnsample,ndatsample,ifreqoi,freqoi(ifreqoi),ntaper(ifreqoi));
        if length(st)>1 && strcmp(st(2).name, 'ft_freqanalysis') && verbose
          % specest_mtmconvol has been called by ft_freqanalysis, meaning that ft_progress has been initialised
          ft_progress(fbopt.i./fbopt.n, ['processing trial %d, ',str,'\n'], fbopt.i);
        elseif verbose
          fprintf([str, '\n']);
        end
        for itap = 1:ntaper(ifreqoi)
          
          dum = fft(ft_preproc_padding(bsxfun(@times,dat,tap{ifreqoi}(itap,:)), padtype, 0, postpad), [], 2);
          
          dum = dum(:,freqboi(ifreqoi));
          % phase-shift according to above angles
          if timedelay ~= 0
            dum = dum .* exp(-1i*angletransform(:,ifreqoi));
          end
          dum = dum .* sqrt(2 ./ endnsample);
          spectrum(:,freqtapind{ifreqoi}(itap)) = dum;
        end
      end % for nfreqoi
  end % switch dimord
end

% 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;
previous_w     = w;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION ensure that the first two input arguments are of double
% precision this prevents an instability (bug) in the computation of the
% tapers for MATLAB 6.5 and 7.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tap, w] = double_dpss(a, b, varargin)
[tap, w] = dpss(double(a), double(b), varargin{:});