MMdesLPFIR (FType, NCof, Gc)

function [h, ier, Ex] = MMdesLPFIR (FType, NCof, Gc)
%  This routine designs multiple band bandpass filters, differentiators,
%  or Hilbert transform filters using the Remez exchange algorithm. This
%  program implements the filter design procedure outlined by McClellan,
%  Parks and Rabiner. This procedure has been extended to include
%  constraints as described by Grenez.
%  References:
%    J. H. McClellan, T.W. Parks and L.R. Rabiner,
%    "A Computer Program for Designing Optimum FIR Linear Phase Digital
%    Filters", IEEE Trans. Audio and Electroacoustics, vol. AU-21,
%    pp. 506-526, December 1973.
%    F. Grenez,
%    "Design of Linear or Minimum-Phase FIR filters by Constrained
%    Chebyshev Approximation", Signal Processing, vol. 5, pp. 325-332,
%    July 1983.
%
% h: Resultant filter coefficients (column vector)
% ier: Error parameter coded as follows.
%      0  - No error
%      1  - Remez algorithm failed to converge, the deviation is not
%           monotonically increasing
%      2  - Remez algorithm failed to converge, too many iterations
%      For these cases, this routine returns coefficient values. For all
%      other errors, an error message is printed and execution is halted.
% Ex:  Structure with the extremal frequencies and deviation.
%
% FType: Filter type.
%      1 - Multiple passband/stopband filter
%      2 - Multiple passband/stopband filter (sin(x)/x compensation)
%      3 - Differentiator
%      4 - Hilbert transform filter
% NCof:  Number of filter coefficients
% Gc:  Structure containing the specifications on a dense grid of points

% Resolve the filter symmetry
[FCase, NPar] = MM_FCase (FType, NCof);

% Find the extremal values for the minimax approximation
NGrid = length (Gc.x);
NExt = min (NGrid, NPar + 1);
[Ex, ier] = MMdesRemez (NExt, Gc);

% Find the coefficients of the approximation
xs = Gc.x(1);
xf = Gc.x(end);
alpha = MMcosCof (xs, xf, Ex);
alpha(end+1:NPar) = 0;

% Calculate the impulse response
h = MMfiltCof (FCase, alpha);
h = h(:);

return

%-----------------
function [FCase, NPar] = MM_FCase (FType, NCof)

%             Sym   Asym
%     FType:  1  2  3  4    bpf rec dif hil
FCaseData = [ 1  1  3  3;   % NCof odd
              2  2  4  4];  % NCof even

Parity = (mod (NCof, 2) == 0) + 1;
FCase = FCaseData (Parity, FType);

switch FCase
  case 1
    NPar = (NCof + 1) / 2;
  case {2, 4}
    NPar = NCof / 2;
  case 3
    NPar = (NCof - 1) / 2;
end

return