function   [outnum, outden] = iirftransf(orignum, origden, ftfnum, ftfden)
%IIRFTRANSF IIR frequency transformation of the digital filter.
%   [OutNum,OutDen] = IIRFTRANSF(OrigNum,OrigDen,FTFNum,FTFDen) returns
%   numerator and denominator coefficients, OUTNUM and OUTDEN of the target
%   filter, which is the result of transforming the prototype filter
%   specified by the numerator, ORIGNUM, and denominator, ORIGDEN, with the
%   mapping filter given by the numerator, FTFNUM, and the denominator,
%   FTFDEN.
%
%   Inputs:
%     OrigNum - Numerator of the prototype lowpass filter provided as a
%               vector, or a as a matrix of coefficients for filter
%               implemented using second or higher order sections.
%     OrigDen - Denominator of the prototype lowpass filter provided as a
%               vector, or a as a matrix of coefficients for filter
%               implemented using second or higher order sections.
%     FTFNum  - Numerator of the mapping filter
%     FTFDen  - Denominator of the mapping filter
%   Outputs:
%     OutNum  - Numerator of the target filter
%     OutDen  - Denominator of the target filter
%
%   Example:
%        [b, a]           = ellip(3, 0.1, 30, 0.409);      
%        [alpnum, alpden] = allpasslp2lp(0.5, 0.25);
%        [num, den]       = iirftransf(b, a, alpnum, alpden);
%        ss = tf2sos(b,a);
%        [alpnum2, alpden2] = allpasslp2lp(0.5, 0.125);
%        [num2, den2] = iirftransf(ss(:,1:3),ss(:,4:6), alpnum2, alpden2);
%        fvtool(b, a, num, den,[num2 den2]);
%
%   See also ZPKFTRANSF.

%   Author(s): Dr. Artur Krukowski, University of Westminster, London, UK.
%   Copyright 1999-2021 The MathWorks, Inc.

%#codegen
% --------------------------------------------------------------------
% Verify the number of input arguments
narginchk(2,4);
% Verify the coefficients of mapping filter
switch nargin
    case 4
        ftransfargchk(ftfnum,'Numerator of the mapping filter',  'vector');
        ftransfargchk(ftfden,'Denominator of the mapping filter','vector');
        % Verify that Numerator and Denominator of the mapping filter have 
        % the same length
        coder.internal.errorIf(length(ftfnum)~=length(ftfden),...
            'dsp:iirftransf:InvalidMappingFilter');

   % FIR case -> denominator defaults to unity
    case 3
        ftfden = fliplr(conj(ftfnum));
        ftransfargchk(ftfnum,'Numerator of the mapping filter',  'vector');

   % Transformation filter defaults to constant equal one
    case 2
        ftfnum = 1;
        ftfden = 1;
end

checkVector = isvector(orignum) && isvector(origden);

% Check if numerator and denominator are vectors or if the filter is
% implemented using sections
if checkVector
    [outnum, outden] = transformFilter(orignum, origden, ftfnum, ftfden);
else
    % Check validity of input arguments
    coder.internal.errorIf(~isnumeric(orignum) || ~ismatrix(orignum),...
        'dsp:iirftransf:InvalidDimension','Numerator');
    coder.internal.errorIf(~isnumeric(origden) || ~ismatrix(origden),...
        'dsp:iirftransf:InvalidDimension','Denominator');
    % Check if number of sections are equal
    coder.internal.errorIf(size(orignum,1)~=size(origden,1),...
        'dsp:iirftransf:SectionsUnequal');    
    
    % Preallocate numerator and denominator matrices for max possible order
    nSections= size(orignum,1);
    maxOrder = (max(size(orignum,2),max(size(origden,2)))-1)*(length(ftfnum)-1);
    maxLength = maxOrder+1;
    % Obtain datatype for Numerator
    if isreal(orignum) && isreal(origden) && isreal(ftfnum) && isreal(ftfden)
        pNum = 1;
    else
        pNum = 1i;
    end         
    if isa(orignum,'single') || isa(origden,'single') || isa(ftfnum,'single') || isa(ftfden,'single')
        pNumType = cast(pNum,'single');
    else
        pNumType = pNum;
    end
    % Obtain dataType for Denominator
    if isreal(origden) && isreal(ftfnum) && isreal(ftfden)
        pDen = 1;
    else
        pDen = 1i;
    end         
    if isa(origden,'single') || isa(ftfnum,'single') || isa(ftfden,'single')
        pDenType = cast(pDen,'single');
    else
        pDenType = pDen;
    end

    outnum = zeros(nSections,maxLength,'like',pNumType);
    outden = zeros(nSections,maxLength,'like',pDenType);
    lenSectionsMax = 0;
    % Iterate and perform transformation on each section
    for i = 1:nSections
        [tempNum,tempDen]  = transformFilter(orignum(i,:), origden(i,:), ftfnum, ftfden);
        lenSection = length(tempNum);
        outnum(i,1:lenSection) = tempNum;
        outden(i,1:lenSection) = tempDen;
        lenSectionsMax = max(lenSectionsMax,lenSection);
    end    
    outnum = outnum(:,1:lenSectionsMax);
    outden = outden(:,1:lenSectionsMax);
end
end


function [outnum, outden] = transformFilter(orignum, origden, ftfnum, ftfden)
    % Function to perform filter transformation
    % Check input arguments
    ftransfargchk(orignum, 'Numerator of the original filter',     'vector');
    ftransfargchk(origden, 'Denominator of the original filter',   'vector');
    
    if length(ftfnum) == 1 && length(ftfden) == 1 && ftfnum == 1 && ftfden == 1
        [outnum,outden] = eqtflength(orignum,origden);
    else
        [outnum, outden] = polyallpasssub(orignum, origden, ftfnum, ftfden);
    
        % Force stability
        s = checkstable(outden);
        if s == 0
            if isreal(outden)
                outden = real(polystable(outden));
            else
                outden = polystable(outden);
            end
        end
    end
end