function varargout = resample(varargin)
%RESAMPLE Resample uniformly spaced data to a new fixed rate on a GPU.
%   Y = RESAMPLE(X,P,Q)
%   Y = RESAMPLE(X,P,Q,N)
%   Y = RESAMPLE(X,P,Q,N,BTA)
%   Y = RESAMPLE(X,P,Q,N,B)
%   [Y,B] = RESAMPLE(X,...)
%   Y = RESAMPLE(X, TX)
%   Y = RESAMPLE(X, TX, FS)
%   Y = RESAMPLE(X, TX, FS, P, Q)
%   Y = RESAMPLE(X, TX, ..., method)
%   [Y,TY] = RESAMPLE(X, TX, ..., method)
%   [Y,TY,B] = RESAMPLE(X,TX,...)
%   [...] = RESAMPLE(...,'Dimension',dim)
%
%   See also UPFIRDN, FIRLS, GPUARRAY.

%   Copyright 2020-2023 The MathWorks, Inc.

% Validate the number of I/O arguements.
narginchk(2,10);
nargoutchk(0,3);

% Determine if we should be using the GPU or not.
[varargin,useGPU,isUniform,isMultichannel] = gpuParser(varargin{:});

%----------------------------------------------------------------------
% Dispatch to in-memory RESAMPLE if useGPU flag is false.
if ~useGPU
    [varargout{1:max(1,nargout)}] = resample(varargin{:});
    return
end

%----------------------------------------------------------------------
% Validate input args and define defaults.
% Initial parse of inputs.
[dim,numArgsExceptDim,isDimValSet,method,numNumericArgs] = ...
    signal.internal.resample.sharedParser(varargin{:});

% Extract numeric args here. This is not shared with in-memory function due
% to in-memory codegen support.
numericArgs = cell(1,numNumericArgs);
idx = 0;

for kk = 1:length(varargin)
    if isnumeric(varargin{kk}) || isdatetime(varargin{kk}) ...
            || istimetable(varargin{kk})
        idx = idx+1;
        numericArgs{idx} = varargin{kk};
    end
end

% No GPU support for PCHIP interpolation.
if strcmp(method,'pchip')
    error(message('signal:resample:NoGpuPchipInterpolation'))
end

% Permute dimensions of signal array, so first dimension is to be
% resampled.
[x, dimIn] = signal.internal.resample.orientToNdim (varargin{1}, dim);

% If we are resampling a singleton dimension for multichannel arrays,
% flow-through in-memory instead for performance, rather than looping over
% channels on GPU.
if (size(x,1) == 1) && isMultichannel
    % Gather all inputs.
    [varargin{1:nargin}] = gather(varargin{:});
    
    % Call in-memory resample.
    [varargout{1:max(1,nargout)}] = resample(varargin{:});
    
    % Transfer any numeric output back to GPU.
    [varargout{1:max(1,nargout)}] = sendToGpu(varargout{:});
    
    return
end

%----------------------------------------------------------------------
% Perform the computation on the GPU.

if isUniform
    % [...] = RESAMPLE(X,P,Q,...)
    if numArgsExceptDim < 3
        error(message('signal:resample:MissingQ'));
    end
    
    % interpolation is not performed when using uniformly sampled data
    if strlength(method) > 0
        error(message('signal:resample:UnexpectedInterpolation',method));
    end
    
    [varargout{1:max(1,nargout)}] = uniformResampleOnGpu(x,size(x),...
        isDimValSet,dimIn,numericArgs{2:end});
    
else
    % [...] = RESAMPLE(X,Tx,...)
    
    if isempty(method)
        % use linear method by default
        method = 'linear';
    end
    
    [varargout{1:max(1,nargout)}] = ...
        nonUniformResampleOnGpu(isDimValSet,dim, numArgsExceptDim,...
        dimIn, method, x,numericArgs{:});
end

end


function [y,hout] = uniformResampleOnGpu(x,sz,isDimValSet,dimIn,p,q,varargin)

% parse uniform inputs
[p,q,isN] = signal.internal.resample.uniformParser(isDimValSet,p,q,varargin{:});

xrDt = real(x(1)); % data class of input

if p == 1 && q == 1
    % Reshape dimensions to input (after permutation of dimensions).
    y = reshape(x,sz);
    % Permute back to original input order.
    y = ipermute(y,dimIn);
    hout = gpuArray(cast(1,'like',xrDt));
    return
end

% Ensure signal is a vector or matrix, if not already.
if isvector(x)
    Lx = length(x);
else
    Lx = sz(1);
    x = reshape(x,Lx,[]);
end

if isN
    % parse N and Beta and design FIR filter
    [h,L] = signal.internal.resample.uniformParserNBetaAndDesignFilter(isDimValSet,xrDt,p,q,varargin{:});
else
    h = varargin{1};
    signal.internal.resample.validateFIRFilterCoeff(h);
    L = length(h);
    [x,h] = signal.internal.utilities.castInputToSingle(x,h);
end

% find FIR filter group delay and zero-pad FIR filter
[h,hout,delay] = signal.internal.resample.findDelayAndZeroPadFilter(h,L,Lx,p,q);

% Get rid of trailing and leading data so input and output signals line up
% temporally:
Ly = ceil(Lx*p/q);  % output length

idx = delay + (1:Ly).';

% Call vectorised upfirdn.
y = upfirdn(x,h,p,q);

% index the result to account for delay.
if isvector(y)
    y = y(idx);
else
    y = y(idx,:);
end

if ~isvector(y)
    % Update output size and reshape.
    sz(1) = Ly;
    y = reshape(y,sz);
    
    % Permute back to original input order.
    y = ipermute(y,dimIn);
end

end

function [y,tyOut,h] = nonUniformResampleOnGpu(...
    isDimValSet,Dim,numArgsExceptDim,dimIn,method,xIn,varargin)

% Parse non-uniform inputs.
[x,tx,tGrid,fs,p,q] = signal.internal.resample.nonUniformParser(...
    isDimValSet,Dim,numArgsExceptDim,xIn,varargin{:});

% Ensure time points are columnar.
if iscolumn(tx)
    tGrid = reshape(tGrid,length(tGrid),1);
end

% Ensure scalars are gathered.
[fs,p,q] = gather(fs,p,q);

% Reshape to single-channel (vector) or multichannel (matrix) if required.
sz = size(x);
x = reshape(x,sz(1),[]);

% Perform the interpolation.
xGrid = interpOnGpu(tx,x,tGrid,method);

% Update resample dimension to new grid size.
sz(1) = size(xGrid,1);

% recover the desired sampling rate by resampling the grid at the
% specified ratio

% If the dimension to be resampled, has been interpolated to 1, flow
% through in-memory code rather than looping over channels on GPU.
if (size(xGrid,1) == 1) && (size(x,1) ~= 1)
    % Gather any GPU arrays to memory.
    [xGrid,varargin{1:length(varargin)}] = gather(xGrid,varargin{:});
    
    % Call in-memory uniform resample.
    [y,h] = signal.internal.resample.uniformResample(xGrid,isDimValSet,...
        Dim,dimIn,varargin{1},p,q,varargin{6:end});
    
    % Transfer outputs back to GPU.
    [y,h] = sendToGpu(y,h);
else
    % Peform uniform resampling on GPU.
    [y, h] = uniformResampleOnGpu(xGrid,sz,isDimValSet,...
        dimIn,p,q,varargin{6:end});
end

% create output time vector
ny = size(y,Dim);

if isdatetime(tx)
    ty = tx(1) + seconds((0:ny-1)/fs);
else
    % tx is a gpuArray
    ty = tx(1) + (0:ny-1)/fs;
end

% match dimensionality of output time vector to input
if iscolumn(tx)
    tyOut = reshape(ty,[],1);
else
    tyOut = ty;
end

end

function xGrid = interpOnGpu(tx,x,tGrid,method)

% If time instances are datetimes, convert to numeric double vectors for
% GPU interpolation.
if isdatetime(tx)
    tx = datenum(tx);
    tGrid = datenum(tGrid);
end

if any(isnan(x),'all')
    % we have to loop over each channel to remove NaNs.
    xGrid = signal.internal.resample.nDInterp1(tx,x,tGrid,method);
else
    % Check for duplicated time points.
    [tx,~,idx] = unique(tx);
    
    % Count the duplicated occurances. num is a vector with length equal to
    % (now duplicate-free) tx.
    num = accumarray(idx,1);
    
    % If there are duplicate time points, we need to remove these, and
    % average signal data at these points.
    if any(num > 1)
        if isrow(x)
            % accumarray will return a column vector, so flag that ouptut
            % should be transposed.
            isRow = true;
        else
            isRow = false;
            numChannels = size(x,2);
            
            % create an array of linear indices to support accumarray with
            % multichannel (columns) input. This utilises an outer sum to
            % create an array from two vectors.
            idx = idx + length(tx)*(0:numChannels-1);
            
            % reshape from array to vector.
            idx = reshape(idx,[],1);
            
            % replicate the duplicate number count for each channel.
            num = repmat(num,numChannels,1);
        end
        
        % Average duplicated data points.
        x = accumarray(idx,x)./num;
        
        % If input signal was row vector, ensure it is transposed from
        % column back to row vector.
        if isRow
            x = x.';
        else
            % Reshape from vector to matrix for multichannel interpolation.
            x = reshape(x,[],numChannels);
        end
        
    end
    
    % vectorised call to interp1.
    xGrid = callInterp1(tx,x,tGrid,method);
    
    % Ensure the interpolated signal data has same orientation if vector.
    if iscolumn(x)
        xGrid = reshape(xGrid,[],1);
    elseif isrow(x)
        xGrid = reshape(xGrid,1,[]);
    end
    
end

end

function xGrid = callInterp1(tx,x,tGrid,method)

% call interp1 with supported methods.
switch method
    case 'spline'
        xGrid = interp1(tx, x, tGrid, 'spline', 'extrap');
    otherwise
        xGrid = interp1(tx, x, tGrid, 'linear', 'extrap');
end

end

function [varargin,useGPU,isUniform,isMultichannel] = gpuParser(varargin)
% GPUPARSER determines whether we should be performing the
% RESAMPLE operations on the GPU or in-memory, and if the operation is
% uniform or non-uniform.

% Check if signal data is a gpuArray.
useGPU = isgpuarray(varargin{1});

% Check if signal data is multichannel (matrix or nD-array).
isMultichannel = ~isvector(varargin{1});

% Check second arguement to determine uniform/non-uniform branch.
isUniform = isscalar(varargin{2});

% Determine index for filter coefficients (if they exist).
if isUniform
    % Define varargin index for filter (if exists).
    filterIdx = 4;
else
    % Define varargin index for filter (if exists).
    filterIdx = 6;
    
    if isgpuarray(varargin{2})
        % Non-uniform case and time-instants are gpuArray.
        useGPU = true;
    end
    
end

% Check if we are using a user-supplied filter.
useFilter = (filterIdx <= nargin) && ...
    isnumeric(varargin{filterIdx}) && ...
    isvector(varargin{filterIdx}) && ...
    ~isscalar(varargin{filterIdx});

% If we are using a user-supplied, check if it is on the GPU.
if useFilter && isgpuarray(varargin{filterIdx})
    useGPU = true;
end

% Ensure all data is transferred to the GPU, if required.
if useGPU
    % Transfer signal X to GPU (if not already there).
    varargin{1} = gpuArray(varargin{1});
    
    % Transfer numerical time instants (if exists) to GPU.
    if ~isUniform && isnumeric(varargin{2})
        varargin{2} = gpuArray(varargin{2});
    else
        % Ensure it is gathered to CPU.
        varargin{2} = gather(varargin{2});
    end
    
    % Transfer filter coefficients (if exists) to GPU.
    if useFilter
        varargin{filterIdx} = gpuArray(varargin{filterIdx});
    end
    
    % If required, gather other variables.
    if nargin > 2
        % Define varargin indices greater than 2.
        gIdx = 3:nargin;
        
        % if a user supplied filter, remove its index from the gather.
        if useFilter
            gIdx(filterIdx - 2) = [];
        end
        
        % Peform the gather operation.
        [varargin{gIdx}] = gather(varargin{gIdx});
    end
    
else
    % Not using GPU, gather all variables and return.
    [varargin{1:nargin}] = gather(varargin{:});
end

end

function [varargout] = sendToGpu(varargin)
% Send to GPU if numeric.
for n = 1:max(1,nargout)
    if isnumeric(varargin{n})
        varargout{n} = gpuArray(varargin{n}); %#ok<AGROW>
    end
end

end
