function [datout, tim, Fnew] = ft_preproc_resample(dat, Fold, Fnew, method)

% FT_PREPROC_RESAMPLE resamples all channels in the data matrix
%
% Use as
%   dat = ft_preproc_resample(dat, Fold, Fnew, method)
% where
%   dat    = matrix with the input data (Nchans X Nsamples)
%   Fold   = scalar, original sampling frequency in Hz
%   Fnew   = scalar, desired sampling frequency in Hz
%   method = string, can be 'resample', 'decimate', 'downsample', 'fft'
%
% The resample method applies an anti-aliasing (lowpass) FIR filter to
% the data during the resampling process, and compensates for the filter's
% delay. For the other two methods you should apply an anti-aliassing
% filter prior to calling this function.
%
% If the data contains NaNs, these are ignored for the computation, but
% retained in the output.
%
% See also PREPROC, FT_PREPROC_LOWPASSFILTER

% Copyright (C) 2006-2012, Robert Oostenveld
%
% 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$

[nchans, nsamples] = size(dat);

if nargout>1
  tim = 1:nsamples;
  tim = ft_preproc_resample(tim, Fold, Fnew, method);
end

if Fold==Fnew
  ft_warning('requested sampling rate after resampling is the same as the sampling rate of the data, no resampling is performed');
  datout = dat;
  return
end

% resample and decimate require double formatted input
if ~strcmp(method, 'downsample')
  typ = class(dat);
  dat = cast(dat, 'double');
end

% preprocessing fails on channels that contain NaN
nanchan = any(isnan(dat),2);
if any(nanchan)
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
  if ft_platform_supports('matlabversion', '2020a', inf)
    % temporarily replace with zero, this is not needed for older MATLAB versions
    dat(nanchan,:) = 0;
  end
end

switch method
  case 'resample'
    [fold, fnew] = rat(Fold./Fnew); %account for non-integer fs
    Fnew         = Fold.*(fnew./fold); %get new fs exact
    
    % the actual implementation resamples along columns
    datout = resample(dat', fnew, fold)';
    
  case 'decimate'
    fac         = round(Fold/Fnew);
    % this only works one channel at the time
    nresampled  = ceil(nsamples/fac);
    datout      = zeros(nchans, nresampled);
    for i=1:nchans
      datout(i,:) = decimate(dat(i,:), fac);
    end
    
  case 'downsample'
    fac = round(Fold/Fnew);
    % the actual implementation resamples along columns
    datout = downsample(dat', fac)';
    
  case 'fft'
    % Code written for SPM by Jean Daunizeau
    fac         = Fnew/Fold;
    nresampled  = floor(nsamples*fac);
    fac         = nresampled/nsamples;
    datfft      = fftshift(fft(dat,[],2),2);
    middle = floor(size(datfft,2)./2)+1;
    if fac>1 % upsample
      npad = floor((nresampled-nsamples)./2);
      
      if nsamples/2 == floor(nsamples/2)
        datfft(:,1) = []; % throw away non symmetric DFT coef
      end
      
      datfft  = [zeros(size(datfft,1),npad), datfft,zeros(size(datfft,1),npad)];
    else % downsample
      ncut    = floor(nresampled./2);
      datfft  = datfft(:,middle-ncut:middle+ncut);
    end
    datout      = fac*ifft(ifftshift(datfft,2),[],2);
  otherwise
    ft_error('unsupported resampling method');
end

if any(nanchan) && ft_platform_supports('matlabversion', '2020a', inf)
  % replace the zeros back to nan, this is not needed for older MATLAB versions
  datout(nanchan,:) = nan;
end

if ~strcmp(method, 'downsample')
  % convert back into the original input format
  datout = cast(datout, typ);
end
