function varargout = resamp(varargin)

% RESAMP   Time signal resampling
%   [R1,R2,...,TR] = RESAMP(X1,T1,X2,T2,...[,T][,MODE]); RESAMP resamples the
%   signal-time pairs XI,TI on a common time base TR. XI can be a matrix of
%   signals, one per column. TI can be a column vector or a matrix of times, one
%   per column. The time base TR is either explicitly given by T, a time vector
%   or a sampling period, and/or a three letter mode argument (default 'fel') :
%   - <f>astest or <s>lowest period of all signals (if T absent)
%   - <r>estricted, NaN padded <e>xtended or <p>rolonged with 1st/last sample
%   - <n>earest sample, <l>inear or <c>ubic or <s>pline interpolation
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% default arguments
ngroup = floor(nargin/2);
mode = varargin{end}; % the last input argument
if ischar(mode)
 mode(4) = mode(2); if mode(4) == 'p', mode(2) = 'e'; end
 if rem(nargin,2)
  imposed = []; % for clarty
 else
  ngroup = ngroup - 1;
  imposed = varargin{end-1}; % the one before the last
  mode(2-(length(imposed)==1)) = 'i';
 end
else
 if rem(nargin,2)
  imposed = mode; mode = 'feln';
  mode(2-(length(imposed)==1)) = 'i';
 else
  imposed = []; % for clarty
  mode = 'feln'; % default default
 end
end

% evaluate sampling characteristics of each signal
tf = []; tl = []; Tf = Inf; Ts = -Inf;
for group = 1:ngroup
 time = varargin{2*group};
 tf = [tf min(time)];
 tl = [tl max(time)];
 time = diff(time);
 Tf = min([Tf;time]);
 Ts = max([Ts;time]);
end

% choose sampling characteristics
switch mode(1)
 case 'f',  T = Tf;
 case 's',  T = Ts;
 otherwise, T = imposed;
end
switch mode(2)
 case 'r',  time = (max(tf):T:min(tl))';
 case 'e',  time = (min(tf):T:max(tl))';
 otherwise, time = imposed(:);
end

% interpolation
for group = 1:ngroup
 [nsamps,nsignal] = size(varargin{2*group-1});
 [nsampt,ntime]   = size(varargin{2*group});
 if nsamps ~= nsampt | (ntime ~= 1 & ntime ~= nsignal)
  error(sprintf('Size mismatch in a signal-time pair #%d.',group))
 end
 varargout{group} = zeros(length(time),nsignal);
 nsignal = 1:nsignal;
 if ntime == 1 & mode(3) ~= 'c', nsignal = nsignal'; end % do all in one step

 for signal = nsignal
  X = varargin{2*group-1}(:,signal);
  T = varargin{2*group}  (:,min(signal));
  if isequal(time,T) % the wanted time may be identical to the signal time
   varargout{group}(:,signal) = X;
  else
   varargout{group}(:,signal) = interp1(T,X,time,mode(3));
   if mode(4) == 'p'
    k = find(time < T(1));
    varargout{group}(k,signal) = repmat(X(1,:),  length(k),1);
    k = find(time > T(end));
    varargout{group}(k,signal) = repmat(X(end,:),length(k),1);
   end
  end
 end
end

% place the time in the output argument
varargout{ngroup+1} = time;
