function y = fftfilt(b,x,nfft)
%FFTFILT Overlap-add method of FIR filtering using FFT on the GPU.
%
%   Y = FFTFILT(B,X)
%   Y = FFTFILT(B,X,N)
%
%   See also FILTFILT, GPUARRAY.

%   Copyright 2012-2023 The MathWorks, Inc.


narginchk(2,3); %check input argument count

if ~(isgpuarray(x) || isgpuarray(b))
    %if only nfft is a gpuArray, dispatch to cpu version
    nfft = gather(nfft);
    y = fftfilt(b,x,nfft);
    return;
else
    %put data in the right place
    if nargin > 2
        nfft = gather(nfft);
    end
    x = gpuArray(x);
    b = gpuArray(b);
end

m = size(x, 1);
if m == 1
    x = reshape(x,[],1);    % turn row into a column
end

nx = size(x,1);

if min(size(b))>1
    if (size(b,2)~=size(x,2))&&(size(x,2)>1)
        error(message('signal:fftfilt:InvalidDimensions'));
    end
else
    b = reshape(b, [],1);   % make input a column
end
nb = size(b,1);

% This GPU algorithm handles empty arrays differently than the CPU.
% To avoid differences, explicitly handle the empties here and return.
if ( isempty(x) || isempty(b) )
    y = gpuArray.zeros(size(x));
    return;
end



if nargin < 3
    % figure out which nfft and L to use
    if nb >= nx || nb > 2^20    % take a single FFT in this case
        nfft = 2^nextpow2(nb+nx-1);
        L = nx;
    else
        fftflops = [ 18 59 138 303 660 1441 3150 6875 14952 32373 69762 ...
            149647 319644 680105 1441974 3047619 6422736 13500637 28311786 59244791];
        n = 2.^(1:20);
        validset = find(n>(nb-1));   % must have nfft > (nb-1)
        n = n(validset);
        fftflops = fftflops(validset);
        % minimize (number of blocks) * (number of flops per fft)
        L = n - (nb - 1);
        [dum,ind] = min( ceil(nx./L) .* fftflops ); %#ok
        nfft = n(ind);
        L = L(ind);
    end

else  % nfft is given
    if nfft < nb
        nfft = nb;
    end
    nfft = 2.^(ceil(log(nfft)/log(2))); % force this to a power of 2 for speed
    L = nfft - nb + 1;
end



% Check the input data type. Single precision is not supported.
if ~(checkDoubles(b) && checkDoubles(x) && checkDoubles(nfft))
    error(message('signal:gpu:OnlyDoubles'));
end

if isreal(b) && isreal(x)
    ifftflag  = 'symmetric';
else
    ifftflag = 'nonsymmetric';
end


B = fft(b,nfft,1);


%Buffer x into sections of length L by inserting zeros and pad with an
%extra column of zeros. FFT this buffered x, multiply in frequency
%domain and then IFFT. The bottom L:nfft rows are the overlap. Using
%indexing, add the overlap from the bottom of one column to the top of the
%next. Then remove the overlap rows. Reshape to get back to the proper
%output size.


xbuf = privBuffer(x, L);
X = fft(xbuf, nfft,1);

[Xr, ~] = size(X);
colsOut = max(size(b,2), size(x,2));

X = reshape(X, Xr, [], size(x,2));
B = reshape(B, size(B,1), 1, size(b,2));
Y = X.*B;
Y = reshape(Y, size(Y,1), []);

y = ifft(Y, nfft, ifftflag);

overlap = y((L+1):nfft, 1:(size(y,2)-1));
tmp = y(1:size(overlap,1), 2:size(y,2)) + overlap;
y(1:size(overlap,1), 2:size(y,2)) = tmp;

%now remove overlap
y = y(1:L, 1:size(y,2));

y = reshape(y, [], colsOut);

%remove tail
y = y(1:size(x,1), :);

if (m == 1)&&(size(y,2) == 1)
    y = reshape(y,1, []);    % turn column back into a row
end

end


function y = privBuffer(x,n)
% Add a matrix of zeros to the bottom of x,
% then n more zeros per column. Reshape to the appropriate size.
% At the end the matrix y should look like:
%   y = [ x(a:b,1)  x(b+1:c, 1) ... 0  x(a:b,2) x(b+1:c,2) ...0
%            0           0       ...0      0        0      ...0
%                   ....
%             0          0          0      0        0         0
%

rem = mod(size(x,1),n);

if (rem ~= 0)
    pad = n-rem;
    y = [x ;...
        gpuArray.zeros(pad,size(x,2))];
else
    y = x;
end

%one more column worth so the overlap doesn't cross from one frame to the
%next.
y = [y; gpuArray.zeros(n, size(x,2))];
y = reshape(y, n, []);

end

function tf = checkDoubles(z)
tf = isUnderlyingType(z,'double');
end

