function [h,err,res] = firgr(order, ff1, aa, varargin)
%FIRGR Generalized Remez FIR filter design.

%   References:
%     Shpak, D. and A. Antoniou, A Generalized Remez Method for the Design
%     of FIR Digital filters, IEEE Trans. Circ. & Sys, vol 37, No 2, 1990.

%   Copyright 1999-2023 The MathWorks, Inc.

coder.internal.errorIf(nargin < 3, 'dsp:firgr:NARGCHK');

minOrder = 0;
estOrder = false;

s = parseOrder(order,estOrder,minOrder);

s = parse_inputs(s, ff1, aa, varargin{:});

% Estimate the filter order with firpmord?
if s.estOrder
    s.nfilt = estimateOrder(s);
end

if s.params.minOrder == 4
    h = s.nfilt - 1;
    err = 0;       
    res = 0;
    return;
end

% For minorder designs, enforce the boundary conditions for
% the filter types
if s.params.minOrder == 1 && ~iscell(s.aa) && ~isequal(s.frf, @genlp) 
    s = enforceType(s);
end

if (s.params.minOrder == 1 && isequal(s.frf, @genhp))
  % Minimum order designs using genhp need to be forced to even order if no
  % order constraint was specified. 
 s.ftype = 1;  
 s.params.minOrder = 2;
end
 
 
% For minOrder designs the minimum increment of the order is 1 unless
% the filter is specified to be even-order, odd-order, or minphase
s.minOrderInc = 1;
if s.params.minOrder
    s = determineIncrement(s);
    doMinOrder(true);
end

% Try up to 50 different orders if we are finding the minimum order
[h,err,s] = design(s);

% Validate and generate filter
[h,err,s] = genfilter(h,err,s);

if nargout > 2
    res = genresstruct(h,s);
end

if s.params.check
    k = find(s.checks == 0, 1);
    if ~isempty(k)
        warning(message('dsp:firgr:transitionAnomaly'));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [h,err,s] = design(s)

% Try up to 50 different orders if we are finding the minimum order
done = false;
s.valid = false;
saveit = [];
s.lastNfilt = s.nfilt;
for tries=1:50
    if done, break; end

    % Design the filter
    [h,err,s] = firgrso(s);

    if s.ret_code == -3 || s.ret_code == -4 || s.ret_code == -5
        break;
    end

    if s.params.minOrder

        % Check if we met specs
        [done,saveit,restoreIt,s] = checkDesign(h,err,saveit,s);

        if restoreIt
            h=saveit{1};		err=saveit{2};		s.iext=saveit{3};
            s.checks = saveit{4};	s.iters=saveit{5};
            s.zero=saveit{6};     s.ftype = saveit{7};	s.nfilt = saveit{8};
            s.grid=saveit{9};		s.des=saveit{10};		s.wt=saveit{11};
            s.ret_code = 0;
        end
    else
        done = true;
        if s.ret_code == 0
            s.valid = true;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hRet,errRet,s] = genfilter(h,err,s)
% Genfilter - generate the final filter when design is valid

hRet = h;
minPhase = s.params.minPhase;
if minPhase
    errRet = sqrt(abs(err));
else
    errRet = abs(err);
end

s.Warnings = {};

if s.valid == 0
    if err(1) > 1.0
        if s.SuppressWarnings
            s.Warnings{end+1}= message('dsp:firgr:finalFilterOrder1', s.lastNfilt);
        else
            warning(message('dsp:firgr:finalFilterOrder1', s.lastNfilt));
        end
    else
        if s.SuppressWarnings
            s.Warnings{end+1} = message('dsp:firgr:finalFilterOrder2', s.lastNfilt);
        else
            warning(message('dsp:firgr:finalFilterOrder2', s.lastNfilt));
        end
    end
end

if ~s.params.minOrder && ~s.valid
    s=handleErrorCode(s.ret_code, s.iters(1), s);
end

if s.valid && any(abs(err) < 100*max(abs(s.des))*eps)
    warning(message('dsp:firgr:machineAccuracy'));
end

% For minimum-phase filters, we now find the required half-order polynomial
if minPhase
    hRet = firminphase(hRet, s.zero, 'angles');
    % Handle maximum phase
    if s.params.minPhase == 2
        hRet = fliplr(hRet);
    end
end

% Correct the sign for a differentiator
hRet = hRet * s.sign_val;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Handle error codes
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s=handleErrorCode(ret_code, iters,s)
% if ret_code == -1
    %error('%s\n%s\n%s\n%s%s', ...
    %    'Approximation error was not computable:', ...
    %    '1) Check the specifications.', ...
    %    '2) Transition region could be much too narrow or too wide for filter order.', ...
    %    '3) For multiband filters, try making the transition regions', ...
    %    ' more similar in width.');
if ret_code == -2
    if s.SuppressWarnings
        s.Warnings{end+1}=message('dsp:firgr:notConverging', iters);
    else
        warning(message('dsp:firgr:notConverging', iters));
    end
elseif ret_code == -3
    error(message('dsp:firgr:InvalidDimensions'));
elseif ret_code == -4
    error(message('dsp:firgr:FilterErr1'));
elseif ret_code == -5
    error(message('dsp:firgr:FilterErr2'));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res = genresstruct(h,s)
% arrange 'results' structure

res.order = s.nfilt - 1;
res.fgrid = s.grid(:);
res.H = freqz(h,1,res.fgrid*pi);
if s.neg  % asymmetric impulse response
    linphase = exp(sqrt(-1)*(res.fgrid*pi*(res.order/2) - pi/2));
else
    linphase = exp(sqrt(-1)*res.fgrid*pi*(res.order/2));
end
if s.hilbert == 1  % hilbert
    res.error = real(s.des(:) + res.H.*linphase);
elseif ~s.params.minPhase
    res.error = real(s.des(:) - res.H.*linphase);
else
    res.error = s.des(:) - abs(res.H);
end
if s.params.minPhase
    res.des = sqrt(s.des(:));
    res.wt = sqrt(s.wt(:));
else
    res.des = s.des(:);
    res.wt = s.wt(:);
end
res.iextr = s.iext(1:end);
res.iextr(res.iextr > length(s.grid)) = length(s.grid); % CLIP MAX INDEX
res.fextr = s.grid(res.iextr);  % extremal frequencies
res.fextr = res.fextr(:);
res.iterations = s.iters(1);
res.evals = s.iters(2);
res.edgeCheck = s.checks;
res.returnCode = s.ret_code;
res.Warnings = s.Warnings;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [done,saveit,restoreIt,s] = checkDesign(h,err,saveit,s)
% Check if we met specs

if s.ret_code == 0
    ok = metSpecs(h,s,err);
else
    ok = false;
end

if ok
    % Save what might be the desired solution
    saveit = {h, err, s.iext, s.checks, s.iters, s.zero, s.ftype, s.nfilt, s.grid, s.des,s.wt};
end
[done,restoreIt,s] = doMinOrder(false,err,saveit,ok,s);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Check if the filter meets the specifications when we
% are finding the lowest-order filter
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ok = metSpecs (h,s,err)
ok = true;
if s.ret_code ~= 0
    ok = false;
else
    for m=1:s.nbands
        if ~s.constr(m) && abs(err(s.banderr(m)))/s.wtx1(m) > s.devs(m)
            ok = false;
        end
    end
end

% For filter of short length, the check above may give false positive.
% Test more thoroughly: 
if s.nfilt <= 5 && ok == true    
    grididx = zeros(size(s.ff));
    for k = 1:length(grididx)
        [~,grididx(k)] = min(abs(s.grid-s.ff(k)));        
    end
    H    = freqz(h,1,s.grid,2);
    Hmag = abs(H);
    
    desIdx = 1;
    for l = 1:2:length(grididx)
        measuredError = norm(Hmag(grididx(l):grididx(l+1)) - ...
            s.des(grididx(l):grididx(l+1)),inf);
        if measuredError > s.wtx(desIdx)
            ok = false;
            break;
        end
        desIdx = desIdx + 1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Check if a minorder design is successful
% Compute the next filter order to try
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [done,restoreIt,s] = doMinOrder(init,err,saveit,ok,s)
persistent orderFactor orderDirection orderHigh orderLow lowErr highErr;
persistent tried orderSlow lowNfilt highNfilt;
if init
    % If the minimum order design is not yet bracketed by two designs,
    % we want to increase (or decrease) the order by a multiplicative
    % factor
    orderFactor = [0.3 0.1];
    orderDirection = 1;  % 1=increase, -1=decrease
    % Both orderHigh and orderLow must become true to bracket the solution
    orderHigh = false;  lowErr = 0;
    orderLow = false;   highErr = 0;
    saveit = []; %#ok
    tried  = [];
    orderSlow = false;
    lowNfilt = -1; % Make sure that lower-order filter is attempted
    highNfilt = 0;
    return;
else
    restoreIt = false;           done = false;

    if ok
        % Found a filter with sufficient order
        orderHigh = true;
        orderDirection = -1;
        s.valid = 1;
        done = true;
        lowErr = err;    highNfilt = s.nfilt;
    elseif abs(err(1)) > 1
        % Found a filter with insufficient order
        orderLow = true;
        orderDirection = 1;
        highErr = err;    lowNfilt = s.nfilt;
    elseif orderLow && orderHigh && highNfilt-s.nfilt == s.minOrderInc
        orderLow = true;  lowNfilt = s.nfilt;
    end

    if ~ok && lowNfilt >= 0 && ~isempty(saveit) ...
            && highNfilt-lowNfilt == s.minOrderInc
        % Restore the previous design
        restoreIt = true;
        s.valid = 1;
        done = true;
    end
    s.lastNfilt = s.nfilt;
    % Must have solution bracketed by the minimum order increment
    if ~orderLow || ~orderHigh || highNfilt-lowNfilt > s.minOrderInc
        done = false;
    elseif done
        return;
    end

    if any(tried == s.nfilt2), orderSlow = true; end
    if orderSlow
        % Avoid a limit cycle in the order
        s.nfilt = s.nfilt + orderDirection * s.minOrderInc;
    elseif orderLow && orderHigh && isfinite(err)
        % The solution is bracketed. Interpolate the approximation
        % errors to get the next order estimate
        lowErr = saveit{2};
        dOrder = (abs(err(1)) - 1) / (abs(highErr(1)) - abs(lowErr(1))) ...
            * (saveit{8} - lowNfilt);
        signum = sign(dOrder);
        dOrder = round(dOrder);
        if abs(dOrder) < s.minOrderInc
            dOrder = s.minOrderInc * signum;
        elseif s.minOrderInc == 2 && mod(dOrder,2)
            dOrder = dOrder + 1;
        end
        if abs(dOrder) > s.nfilt/2, dOrder = round(sign(dOrder)*s.nfilt/2); end
        s.nfilt = s.nfilt + dOrder;
    else
        % Solution is not bracketed.  Extend the interval.
        if abs(err(1)) < 0.5 || abs(err(1)) > 2.0
            dOrder = round(s.nfilt * orderFactor(1));
        else
            dOrder = round(s.nfilt * orderFactor(2));
        end
        if dOrder > s.minOrderInc
            if s.minOrderInc == 2 && mod(dOrder,2)
                dOrder = dOrder + 1;
            end
        else
            dOrder = s.minOrderInc;
        end
        s.nfilt = s.nfilt + orderDirection*dOrder;
    end

    if s.nfilt < 4
        if orderHigh
            s.nfilt = highNfilt;
            return
        else
            error(message('dsp:firgr:orderTooSmall1'));
        end
    end
    tried = [tried s.nfilt2];
    s.ftype = new_ftype(s.ftype, s.nfilt, s.ff, s.aa);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% determine Increment
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = determineIncrement(s)
% Determine order increment

if s.params.minPhase
    s.minOrderInc = 2; % Because order is doubled
    % If even or odd order is requested, make it so.
elseif s.params.minOrder == 2  % Even-order (odd-length) design
    if rem(s.nfilt,2) == 0, s.nfilt = s.nfilt + 1; end
    s.minOrderInc = 2;
elseif s.params.minOrder == 3 % Odd-order design
    if rem(s.nfilt,2) == 1, s.nfilt = s.nfilt + 1; end
    s.minOrderInc = 2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% enforce filter type
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = enforceType(s)
% Enforce the correct filter type

if (s.ftype == 1 || s.ftype == 2) && s.ff(end) == 1 && s.aa(end) ~=0
    s.ftype = 1;    s.params.minOrder = 2;
elseif (s.ftype == 3 || s.ftype == 4) && s.ff(1) == 0 && s.aa(1) ~= 0
    error(message('dsp:firgr:oddSymmetry'));
elseif s.ftype == 3 && s.ff(end) == 1 && s.aa(end) ~= 0
    s.ftype = 4;    s.params.minOrder = 3;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function nfilt = estimateOrder(s)

% We can only estimate the order using firpmord if it's
% a typical multiband filter.
% Otherwise, we'll use start with the order that was
% specified in the call to firpm
if ~isequal(s.frf, @firpmfrf2)
    error(message('dsp:firgr:orderEstimation1'));
end
if any(diff(s.aa(1:2:end)-s.aa(2:2:end))~= 0.0)
    error(message('dsp:firgr:orderEstimation2'));
end
if s.ff(1) ~= 0 || s.ff(end) ~= 1
    error(message('dsp:firgr:orderEstimation3'));
end

% If we are designing with an estimated minimum order, need to call the frf
% to pre-process the specifications
grid = s.ff(1):s.ff(2);
[des,wt,wtx1,s.ff,s.free_edges,devs,constr,banderr,lb,ub,errindex] = ...
    s.frf(0, s.ff, s.free_edges, grid, s.wtx, s.wtx_prop, ...
    s.params.minPhase, s.params.minOrder, s.frf_params{:}); %#ok

nfilt = firpmord(s.ff(2:end-1), s.aa(1:2:end), devs, 2);
if (s.params.minPhase && rem(nfilt,2) == 0)
    nfilt = nfilt + 1;
end
if nfilt < 3
    % disp('Estimated filter order was < 3.  Trying order=3.');
    nfilt = 3;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Process the specified order, which may be an integer, string
% or cell array
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = parseOrder(order,estOrder,minOrder)
% If we have a character string for the order, we will
% need to call firpmord to estimate the order
ordstring = [];
if dsp.util.isChar(order)
    estOrder = true;
    ordstring = char(order);
    order = 0;
    %params.newFeatures = 1;
end
if iscell(order)
    %params.newFeatures = 1;
    ordstring = order{1};
    if length(order) == 2
        % An estimate of the order has been specified
        order = order{2};
    else
        estOrder = true;
        order = 0;
    end
end

% If we have a string for the order, parse it
if ~isempty(ordstring)
    if strncmp(ordstring, 'minorder', 5)
        minOrder = 1;
    elseif strncmp(ordstring, 'mineven', 5)
        minOrder = 2;
    elseif strncmp(ordstring, 'minodd', 5)
        minOrder = 3;
    elseif strncmp(ordstring, 'estorder', 5)
        minOrder = 4;
        estOrder = true;
    else
        error(message('dsp:firgr:wrongOrder'));
    end
end

s.order = order; s.estOrder = estOrder; s.params.minOrder = minOrder;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = parse_inputs(s, ff1, aa, varargin)

% Define default values for input arguments
ftype = 'f';
lgrid = 16;   % Grid density (should be at least 16)

% parse inputs and alter defaults
%
%  First find cell array and remove it if present
%  This cell array is used for changing lgrid
%
%  Also, check for any features that are new in the
%  C-mex version.  If no new features are used,
%  we will use the new Remez in "legacy" mode for
%  maximum compatibility with the FORTRAN version.

s.params.newFeatures = 0;  % Legacy mode

suppressWarnings = false;
for i=1:length(varargin)
    if ischar(varargin{i}) || isstring(varargin{i})
        if strcmpi(varargin{i},'SuppressWarnings')
            suppressWarnings = varargin{i+1};
            varargin(i:i+1) = [];
            break;
        end
    end
end

s.SuppressWarnings = suppressWarnings;

% If we have a cell array argument in varargin
% containing one scalar, it is the grid density specification.
% Use this value and remove it from varargin
for i=1:length(varargin)
    if (iscell(varargin{i}) && length(varargin{i}) == 1 && isnumeric(varargin{i}{:}))
        lgrid = varargin{i}{:};
        if lgrid < 16
            warning(message('dsp:firgr:lowDensity'));
        end
        if lgrid < 1
            error(message('dsp:firgr:negativeGridDensity'));
        end
        varargin(i) = [];
        break
    end
end

% If the first varargin is a cell array, then it must be
% a cell array of characters specifying band edge properties
edge_prop = [];
if ~isempty(varargin) && (iscell(varargin{1}) || isstring(varargin{1}))
    edge_prop = cellstr(varargin{1});
    varargin(1) = [];
    s.params.newFeatures = 1;
end

% The default is a vanilla linear-phase filter
% and no checking for transition-region anomalies
s.params.minPhase = 0;
s.params.check = 0;
wtx = [];
wtx_prop = [];
for i=1:length(varargin)
    if dsp.util.isChar(varargin{i})
        % We have a string argument, so look for text-based options
        charVal = char(varargin{i});
        str = lower(charVal);
        if length(charVal) == 1
            ftype = charVal;
        elseif strncmp(str, 'minphase', 5)
            s.params.minPhase = 1;
        elseif strncmp(str, 'maxphase', 5)
            s.params.minPhase = 2;
        elseif strncmpi(str, 'hilbert',length(str))
            ftype = 'h';
        elseif strncmpi(str, 'differentiator',length(str))
            ftype = 'd';
        elseif strcmp (str, 'check')
            s.params.check = 1;
        end
    elseif iscell(varargin{i}) || isstring(varargin{i})
        % If there is a cell array remaining,
        % it must be the weight properties
        wtx_prop = cellstr(varargin{i});
        s.params.newFeatures = 1;
        if any(strcmpi(wtx_prop{1}, {'w', 'c'}))
            warning(message('dsp:firgr:DeprecatedFeature'));
        end
    else
        % Otherwise it is the error weights for the bands
        wtx = varargin{i};
    end
end

if s.params.check ~= 0 || s.params.minPhase ~= 0
    s.params.newFeatures = 1;
end

% To get accurate location for stopband zeros, minPhase filter
% need more grid points
if s.params.minPhase
    if lgrid < 64; lgrid = 64; end
end

if isempty(ftype), ftype = 'f'; end

if length(ftype) == 1 && ftype(1)=='m'
    error(message('dsp:firgr:noMFile'));
end

if ~s.estOrder && (fix(s.order) ~= s.order || s.order < 3)
    error(message('dsp:firgr:orderTooSmall2'))
end
%
% Error checking
%
% Need to frequencies and band-edge properties
[ff,aa,free_edges] = firpmedge(ff1, aa, edge_prop);
if isempty(wtx)
    if s.params.minOrder ~= 0
        error(message('dsp:firgr:missingDeviations'));
    end
    wtx = ones(fix((1+length(ff))/2),1);
end
wtx1 = wtx;

if rem(length(ff),2)
    error(message('dsp:firgr:oddNumberOfFreqPoints'));
end
% Error check the band-edge frequencies
if any(ff < 0) || any(ff > 1)
    error(message('dsp:firgr:freqOutOfRange'));
end
df = diff(ff);
if (any(df < 0))
    error(message('dsp:firgr:decreasingFreqs'));
end
if length(wtx) ~= fix((1+length(ff))/2)
    if s.params.minOrder
        error(message('dsp:firgr:wrongNumberOfDevs'));
    else
        error(message('dsp:firgr:wrongNumberOfWeights'));
    end
end

if (any(sign(wtx) == 1) && any(sign(wtx) == -1)) || any(sign(wtx) == 0)
    if(s.params.minOrder)
        error(message('dsp:firgr:nonPositiveRipples'));
    else
        error(message('dsp:firgr:nonPositiveWeights'));
    end
end

nbands = length(ff)/2;

%
% Determine "Frequency Response Function" (frf)
% by inspecting the amplitude specification
%

% The following comment, MATLAB compiler pragma, is necessary to allow the
% Compiler to find the private functions. Do not remove.
%#function firpmfrf
%#function firpmfrf2

if dsp.util.isChar(aa)
    % If still using a string, convert to function handle
    frf = str2func(char(aa));
    frf_params = {};
elseif isa(aa, 'function_handle')
    % We have a user-specified frf
    frf = aa;
    frf_params = {};
elseif iscell(aa)
    frf = aa{1};
    if dsp.util.isChar(frf)
        frf = str2func(char(frf));
    end
    frf_params = aa(2:end);
else
    % We have a standard bandwise filter design
    frf = @firpmfrf2;
    % A differentiator gets a special weighting function
    frf_params = { aa, strcmpi(ftype(1),'d') };
end

% We need the following code to generate fcn handles to private functions.
% Otherwise the syntax h=firgr(n,f,{@firpmfrf2,m},w); will not work
fcns = functions(frf);
if isempty(fcns.file)
    frf = str2func(func2str(frf));
end


% Work with filter length
nfilt = s.order + 1;        % filter length

% Create a structure of outputs to return
s.nfilt = nfilt; s.lgrid = lgrid; s.ff=ff; s.aa = aa;
s.ftype=ftype; s.free_edges=free_edges; s.nbands = nbands; s.wtx=wtx;
s.wtx1=wtx1; s.wtx_prop = wtx_prop; s.frf=frf; s.frf_params=frf_params; 


% Determine symmetry of filter
s = doFtype(s);

% Make sure that we can find the mex file wrapper
if (exist('filterdesign.internal.gremezmex', 'file') ~= 0)
    error(message('dsp:firgr:noExecutable'));
end

if s.params.minPhase && rem(s.nfilt, 2) == 0
    error(message('dsp:firgr:evenOrder'));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Determine and verify the filter type
% (1 = even-order symmetric)
% (2 = odd-order symmetric)
% (3 = even-order antisymmetric)
% (4 = odd-order antisymmetric)
% Also, the special cases of Hilbert transformer or differentiator
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = doFtype(s)

nodd = rem(s.nfilt,2);      % nodd == 1 ==> filter length is odd
% nodd == 0 ==> filter length is even
sign_val = 1; % Signs of coefficients get changed for a Differentiator

% Support the existing filter types 'hilbert' and 'differentiator'
% but also let the user specify types 1, 2, 3, or 4
hilbert = 0;
switch (lower(s.ftype(1)))
    case 'h'
        s.ftype = 3;  % Hilbert transformer
        hilbert = 1;
        if (rem(s.order, 2) == 1)
            s.ftype = 4;
        end
    case 'd'
        s.ftype = 4;  % Differentiator
        sign_val = -1;
        if (rem(s.order, 2) == 0)
            s.ftype = 3;
        end
    case '1'
        s.ftype = 1;
        if (nodd == 0)
            error(message('dsp:firgr:incompatibleSpecs1'));
        end
    case '2'
        s.ftype = 2;
        if (nodd == 1)
            error(message('dsp:firgr:incompatibleSpecs2'));
        end
    case '3'
        s.ftype = 3;
        if (nodd == 0)
            error(message('dsp:firgr:incompatibleSpecs3'));
        end
    case '4'
        s.ftype = 4;
        if (nodd == 1)
            error(message('dsp:firgr:incompatibleSpecs4'));
        end
    otherwise
        % If symmetry was not specified, call the fresp function
        % with 'defaults' string and a cell-array of the actual
        % function call arguments to query the default value.
        if isequal(s.frf, @firpmfrf2)
            h_sym = 'even';
        else
            try
                h_sym = s.frf('defaults', [{s.order, s.ff, [], s.wtx}, s.frf_params(:).']);
            catch
                h_sym = 'even';
            end
        end

        if ~any(strcmp(h_sym,{'even' 'odd'}))
            error(message('dsp:firgr:invalidSymmetry', h_sym, s.frf));
        end

        switch h_sym
            case 'even' % Symmetric filter
                s.ftype = new_ftype(1, s.nfilt, s.ff, s.aa);
            case 'odd'		% Odd (antisymmetric) filter
                s.ftype = new_ftype(3, s.nfilt, s.ff, s.aa);
        end
end

if (s.ftype == 3 || s.ftype == 4)
    neg = 1;  % neg == 1 ==> antisymmetric imp resp,
else
    neg = 0;  % neg == 0 ==> symmetric imp resp
end

% Add properties to structure
s.sign_val = sign_val; s.neg=neg; s.hilbert=hilbert;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Since the order can change when finding the lowest filter
% order, the filter type can change too.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function newtype = new_ftype(ftype, nfilt, ff, aa)
switch ftype
    case {1, 2}
        if rem(nfilt, 2)
            newtype = 1;
        else
            newtype = 2;
            if iscell(aa), return; end
            if (ff(end) == 1 && aa(end) ~= 0)
                error(message('dsp:firgr:incompatibleSpecs5'));
            end
        end
    case {3, 4}
        if rem(nfilt, 2)
            newtype = 3;
            if iscell(aa), return; end
            if (ff(1) == 0 && aa(1) ~= 0) || (ff(end) == 1 && aa(end) ~= 0)
                error(message('dsp:firgr:incompatibleSpecs6'));
            end

        else
            newtype = 4;
            if iscell(aa), return; end
            if (ff(1) == 0 && aa(1) ~= 0)
                error(message('dsp:firgr:incompatibleSpecs7'));
            end
        end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Specify order design
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [h,err,s] = firgrso(s)
%FIRGRSO   Specify-order FIRGR.

% Core design function. Generate frequency grid and call mex-file for
% design

% Create grid of frequencies on which to perform firpm exchange iteration
%
% We use the new grid function even in legacy mode since
% the results are more accurate
s.grid = firpmgrid(s.nfilt,s.lgrid,s.ff,s.ftype,s.free_edges);
while length(s.grid)<=s.nfilt
    s.lgrid = s.lgrid*4;  % need more grid points
    s.grid = firpmgrid(s.nfilt,s.lgrid,s.ff,s.ftype,s.free_edges);
end

%
% Get desired frequency characteristics at the frequency points
% in the specified frequency band intervals.
%
% NOTE! The frf needs to see normalized frequencies in the range [0,1].
%
% Unconstrained designs have empty constraints lb and ub
lb = [];
ub = [];
errindex = ones(length(s.grid),1);
s.banderr = ones(s.nbands, 1);

try
    % First we will see if the "frf" returns the maximum nargout
    [s.des,s.wt,s.wtx1,s.ff,s.free_edges,s.devs,s.constr,s.banderr,lb,ub,errindex] = ...
        s.frf(s.order, s.ff, s.free_edges, s.grid, s.wtx, s.wtx_prop, ...
        s.params.minPhase, s.params.minOrder, s.frf_params{:}); 
catch
    if isequal(s.frf, @firpmfrf2)
        error(message('dsp:firgr:cannotEvalSpecs'));
    end
    % Otherwise we must be using a user-supplied legacy frf
    [s.des,s.wt] = s.frf(s.order, s.ff, s.grid, s.wtx, s.frf_params{:}); 
end

s.ftype = new_ftype(s.ftype, s.nfilt, s.ff, s.aa); 

if (s.nfilt < 100)
    s.params.fast = 1;  % Use fast Vandermonde solver
else
    s.params.fast = 0;
end

s.params.dof = max(errindex);

if s.params.minPhase
    % Minimum-phase designs start with a filter of twice the order
    s.nfilt2 = 2*(s.nfilt-1) + 1;
    % We also need to square the weights, constraints, and the desired values
    k = find(s.des == 0);
    s.wt(k) = s.wt(k).^2;
    k = find(s.des > 0);
    s.wt(k) = s.des(k).*s.wt(k)/2;
    %s.wt(k) = sqrt(s.wt(k));
    lb = lb.^2;
    ub = ub.^2;
    s.des = s.des.^2;
else
    s.nfilt2 = s.nfilt;
end

s.params.ftype = s.ftype;

% Call MEX-file
% Note that the mex file does not support single precision. If for any
% reason a single precision variable was able to sneak in, cast it here to a double
[h,err,s.iext,s.ret_code,s.checks,s.iters,s.zero] = filterdesign.internal.gremezmex(...
    double(s.nfilt2), ...
    double(s.ff), ...
    double(s.grid), ...
    double(s.des), ...
    double(s.wt), ...
    double(s.wtx1), ...
    double(s.free_edges), ...
    double(lb), ...
    double(ub), ...
    double(errindex), ...
    s.params);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set up the grid of frequency points for the approximation
% on a compact subset of the interval [0, pi]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function grid = firpmgrid(nfilt,lgrid,ff,ftype,free_edges)
nfcns = fix(nfilt/2);
if ftype == 1
    nfcns = nfcns + 1;
end
nbands = length(ff) / 2;
ngrid = lgrid*nfcns;

% Sum the parts of the frequency axis that are
% used (i.e., exclude "don't care" regions
fsum = 0.0;
for i=1:nbands
    fsum = fsum + (ff(2*i) - ff(2*i-1));
end

if (fsum == 0)
    error(message('dsp:firgr:finiteWidthBand'));
end
grid = [];
delf = fsum / ngrid;

if (nbands == 1)
    grid = ff(1) + (0:ngrid-2)*delf;
    grid(ngrid) = ff(2);
else
    % Double the number of grid points near internal band edges
    % This results in more accurate designs with a small increase
    % in computational effort
    for i=1:nbands
        if free_edges(2*i-1) == 1
            if i == 1
                left = 0;
            else
                left = ff(2*i-2) + delf;
            end
        else
            left = ff(2*i - 1);
        end
        if free_edges(2*i) == 1
            if i == nbands
                right = 1;
            else
                right = ff(2*i + 1);
            end
        else
            right = ff(2*i);
        end
        % Do not duplicate grid points
        if (i > 1)
            if (grid(end) == left)
                if (length(grid) > 1)
                    grid = grid(1:end-1);
                else
                    grid = [];
                end
            end
        end
        ngrid = 2*round((right - left) / delf);
        delf1 = (right - left) / (ngrid-1);
        if (ngrid <= 1 && (right-left) < eps)
            % One grid point
            newgrid = [grid 0.5*(left+right)];
            grid = newgrid;
        elseif (ngrid < 3*lgrid)
            % Not many grid points in the band, so
            % increase density for accuracy and robustness
            ngrid = 3*lgrid;
            delf1 = (right - left) / (ngrid-1);
            newgrid = [grid left:delf1:right];
            grid = newgrid;
            grid(end) = right;
        elseif (i == 1)
            % Normal spacing
            stop = right - lgrid*(2*delf1);
            newgrid = [grid left:2*delf1:stop];
            grid = newgrid;
            % Narrow spacing
            grid = [grid(1:end-1) grid(end):delf1:right+eps];
            grid(end) = right;
        elseif (i == nbands)
            % Narrow spacing
            stop = left + lgrid*(2*delf1);
            newgrid = [grid left:delf1:stop;];
            grid = newgrid;
            % Normal spacing
            grid = [grid(1:end-1) grid(end):2*delf1:right+eps];
            grid(end) = right;
        else
            % Narrow spacing
            stop = left + lgrid*(1.3*delf1);
            newgrid = [grid left:delf1:stop];
            grid = newgrid;
            % Normal spacing
            stop = right - lgrid*(1.3*delf1);
            grid = [grid(1:end-1) grid(end):2*delf1:stop];
            % Narrow spacing
            grid = [grid(1:end-1) grid(end):delf1:right+eps];
            grid(end) = right;
        end
    end
end

% If the filter has odd symmetry, the value at f=0
% is constrained to be zero and the first grid point
% must not be at zero.
if (ftype == 3 || ftype == 4) && (grid(1) < delf)
    grid = grid(2:end);
end
% If the filter has even length, the value at f=1
% is constrained to be zero and the last grid
% point must not be at 1.0.
if (ftype == 2 || ftype == 3) && grid(end) > (1 - delf)
    grid = grid(1:end-1);
end

% [EOF]
