function X = sosfilt(SOS, X, dim)
%SOSFILT Second-order (biquadratic) IIR filtering of gpuArray data.
%   Y = SOSFILT(SOS,X) converts the direct-form II coefficients in the
%   matrix SOS into parallel second-order subsections and applies it to the
%   input signal X. The coefficients of the SOS matrix must be expressed
%   using an L-by-6 second-order section matrix, where L is the number of
%   second-order sections. If X is a matrix, SOSFILT operates along the
%   columns of X.
%
%   Y = SOSFILT(SOS,X,DIM) operates along the dimension DIM.
%
%   Usage notes and limitations:
%   - Input SOS filter must be stable if it has one or more IIR
%   subsections. Use isstable(SOS) to check for filter stability.
%   - gpuArray/sosfilt uses a parallel algorithm that is different from the
%   MATLAB version. The algorithms give different results for
%   complex-valued SOS inputs with Nan or Inf values. If the input SOS
%   filter has one or more IIR subsections, then:
%      - In the MATLAB version, then NaNs and Infs propagate only in the
%      real part.
%      - In the gpuArray version, the NaNs and Infs propagate in both the
%      real part and the imaginary part.
%
%  See also SOSFILT, ISSTABLE, GPUARRAY.

%  References:
%  Bank, Balazs. "Converting Infinite Impulse Response Filters to Parallel
%  Form". IEEE Signal Processing Magazine. Vol. 35, Number 3, May 2018, pp.
%  124-130.


% Copyright 2020-2022 The MathWorks, Inc.

narginchk(2,3);
if nargin>2
    [SOS,dim] = gather(SOS,dim);
else
    SOS = gather(SOS);
    dim = [];
end

% If X not on GPU, redirect to in-memory version
if ~isgpuarray(X)
    X = parallel.internal.array.callToolboxFunction("signal/signal/sosfilt", ...
        SOS, X, dim);
    return
end

% We duplicate the basic input checking to ensure correct error messages
if ~isfloat(X)
    error(message('signal:sosfilt:VectMustBeFloat', 'X'));
end
if ~isfloat(SOS)
    error(message('signal:sosfilt:MatrixMustBeFloat', 'SOS'));
end
[m,n] = size(SOS);
if (m<1) || (n~=6)
    error(message('signal:sosfilt:InvalidDimensions'));
end

% Parsosfilt has an additional limitation that all subsections must be IIR.
% If any subsection is FIR, use filter instead.
isIIR = isIIRSubsection(SOS);
if ~all(isIIR)
    % Mixed FIR and IIR subsections or only FIR subsections. Process FIR
    % and IIR subsections separately.
    FIRSOS = SOS(~isIIR, :);
    % Process each FIR second-order subsection at a time.
    for k = 1:size(FIRSOS,1)
        if isempty(dim)
            X = filter(FIRSOS(k,1:3),1,X);
        else
            X = filter(FIRSOS(k,1:3),1,X,[],dim);
        end
    end
    SOS = SOS(isIIR,:);
end

if isempty(SOS)
    % All subsections were FIR and have already been processed.
    return
end

% Convert Direct Form II SOS coefficients to Parallel SOS coefficients
[B, A, FIR] = signal.internal.gpu.sos2parsos(SOS);

% Make sure X is 2D and we are operating down the first dimension
szX = size(X);
[X,perm,nshifts] = shiftdata(X,dim);
shiftSzX = size(X); % Shifted size of X
X = reshape(X,size(X,1),[]); % Force into 2-D

% Perform Parallel SOS filter.
X = parallel.internal.gpu.filter.parsosfilt(B, A, FIR, X);

% Put back to original size/shape
X = reshape(X,shiftSzX); % Back to N-D array
X = unshiftdata(X,perm,nshifts);
X = reshape(X,szX);
end

function isIIR = isIIRSubsection(SOS)
% Checks to see if each subsection of the input SOS matrix is IIR (i.e.
% every row has a non-zero in position 5 or 6 or both).
isIIR = SOS(:,5) | SOS(:,6);
end