function xrec = ifsst(sst,varargin)
%IFSST Inverse Fourier synchrosqueezed transform on GPU.
%
%   SST = IFSST(X)
%   SST = IFSST(X,WIN)
%   SST = IFSST(X,WIN,F,FREQRANGE)
%   SST = IFSST(X,WIN,...,IRIDGE)
%   SST = IFSST(X,WIN,...,IRIDGE,'NumFrequencyBins',NBINS)

%   Copyright 2022-2023 The MathWorks, Inc.

narginchk(1,5);

%Only want to perform computation on the GPU when either the SST matrix or
%the window is a gpuArray

%safe logic for single argument case due to logical shortcircuit when
%varagin is empty as this means SST has to be a gpuArray to call this
%overload
if  (~isgpuarray(sst) && (~isgpuarray(varargin{1}) || isscalar(varargin{1})))
    [sst,varargin{:}] = gather(sst,varargin{:});
    xrec = ifsst(sst,varargin{:});
    return
end

if ~isgpuarray(sst)
    sst=gpuArray(sst);
end

[iridge,f,freqrange,nfft,nbins,fReal,invScaledFreqInt] = signal.internal.fsst.ifsstParser(sst,varargin{:});

% Determine the last index to use in the second side of the spectrum.
% For an even length fft, the one-sided spectrum has (nfft/2)+1 elements,
% while for an odd length fft, the one-sided spectrum has (nfft+1)/2
% elements.

%Gathering theses arguments as there is no benefit for them being of the
%GPU
[f,freqrange] = gather(f,freqrange);

if signalwavelet.internal.isodd(nfft)
    iend = size(sst,1);
else
    iend = size(sst,1)-1;
end
% Extract ridges if provided
if ~isempty(iridge)
    %iridge must be a gpuArray
    iridge = gpuArray(iridge);

    sizIridge = size(iridge);

    xrec = zeros(sizIridge,'like',sst);

    % Do the inverse
    fIndices = transpose(gpuArray.colon(1,size(sst,1)));

    for j = 1:sizIridge(2)
        %For each ridge, apply the values of the SST nbins above and below
        %the ridge to capture similar frequencies near the ridge
        ridge = parallel.internal.gpu.extractRows(transpose(iridge),j);
        maskSST = arrayfun(@applyMask,sst,ridge,nbins,fIndices);

        if fReal
            xrecSingleRidge = invScaledFreqInt*real(sum(maskSST,1)+ ...
                conj(sum(parallel.internal.gpu.extractRows(maskSST, 2:iend),1)));

            xrec(:,j) = xrecSingleRidge;
        else
            xrec(:,j) = invScaledFreqInt*(sum(maskSST,1)).';
        end
    end

elseif ~isempty(freqrange)
    [~,i1] = min(abs(f - freqrange(1)));
    [~,i2] = min(abs(f - freqrange(2)));
    % Do the inverse

    if fReal
        xrec = invScaledFreqInt*real(sum(parallel.internal.gpu.extractRows(sst, i1:i2),1) + ...
            conj(sum(parallel.internal.gpu.extractRows(sst, max(i1,2):min(i2,iend)),1)));
    else
        xrec = invScaledFreqInt*(sum(parallel.internal.gpu.extractRows(sst, i1:i2),1));
    end

    xrec = colonize(xrec);
else
    % Do the inverse
    if fReal

        xrec = invScaledFreqInt*real(sum(sst,1) + conj(sum(parallel.internal.gpu.extractRows(sst, 2:iend),1)));
    else
        xrec = invScaledFreqInt*(sum(sst,1));
    end

    xrec = colonize(xrec);

end
end

function value = applyMask(sst,freqRidgeIdx,numBins,fIdx)

%Checking if the frequency falls between freqRidge +- numbins
tf = fIdx >= freqRidgeIdx-numBins && fIdx <= freqRidgeIdx+numBins;

value = zeros('like',sst);
if tf
    value = sst;
end

end
