DFiltMPFIR (NCof, BSpec, SFreq, GridD)

function hM = DFiltMPFIR (NCof, BSpec, SFreq, GridD)
% h = DFiltMPFIR (NCof, BSpec, FType, SFreq, GridD)
% The program designs a minimum-phase filter. The filter specifications are
% in terms of a double order linear phase filter which is then factored
% to give the minimum-phase component. If the min-phase filter is of length
% N, the double-order linear phase filter is of length 2*N-1. The values,
% weights, and limits given apply to the double-order filter. The program
% imposes constraints to guarantee that any roots on the unit circle are
% double order to allow factorability of the filter into a min-phase part
% and a max-phase part. Each has the same amplitude response. The min-phase
% part is returned by this program. The max-phase part is the time-reversal
% of the min-phase filter. The program uses the locations of the unit
% circle zeros (known from the design procedure) to help in factoring the
% double-order filter. This strategy works well if the number of stopband
% zeros is large.
%
%   hM:    Filter coefficients (column vector)
%
%   NCof:  Number of filter coefficients (of the min-phase filter)
%   BSpec: Structure array containing the filter specifications of the
%          double-order filter, band by band. For each band, BSpec(i) has
%          the following fields.
%          BSpec(i).Freq: Frequency values (increasing order) for band i.
%            The frequencies within a band and for successive bands must
%            be increasing. The first and last frequencies of a band define
%            the extent of that band. Gaps between bands are transition
%            regions. Intermediate frequencies in a band allow for the
%            specification of values, weights, and limits at those
%            intermediate frequencies. Frequencies must lie in the
%            interval [0, 0.5*SFreq].
%          BSpec(i).Value: Values. This can be a single value that applies
%            to all frequencies in the band, or a value for each frequency
%            point in the band.
%          BSpec(i).Weight: Weights at the frequencies. This can be a
%            single weight that applies to all frequencies in the band, or
%            a weight for each frequency point in the band. The weights
%            must be positive. If a weight is specified for any band,
%            weights must be given for all bands.
%     The limit specifications below are optional. Note that this routine
%     will modify the limits to prevent the response from going negative.
%          BSpec(i).LLimit: Lower limits, either a single limit that
%            applies to all frequencies in the band or a limit for each
%            frequency point in the band. The default is 0.
%          BSpec(i).ULimit: Upper limits, either a single limit that
%            applies to all frequencies in the band or a limit for each
%            frequency point in the band. The default is Inf.
%     The following are optional.
%          BSpec(i).ValueInt: Interpolation domain. There is one
%             specification per band indicating the domain for
%             interpolation of the values. The choices are 'log', 'sqrt',
%             or 'linear'. The default is 'sqrt'.
%          BSpec(i).WeightInt: Interpolation domain. There is one
%             specification per band indicating the domain for
%             interpolation of the weights. The choices are 'log', 'sqrt',
%             or 'linear'. The default is 'linear'.
%          BSpec(i).LLimitInt: Interpolation domain. There is one
%             specification per band indicating the domain for
%             interpolation of the lower limits. The choices are 'log',
%            'sqrt', or 'linear'. The default is 'sqrt'.
%          BSpec(i).LLimitInt: Interpolation domain. There is one
%             specification per band indicating the domain for
%             interpolation of the upper limits. The choices are 'log',
%            'sqrt', or 'linear'. The default is 'sqrt'.
%   SFreq: Sampling frequency (used to interpret the frequencies in the
%           filter specifications). If missing, SFreq = 1.
%   GridD:  Grid density (number of grid points per extremum of the
%           response). The default is chosen to give at least 1000 grid
%           points or 32 points per extremum. If missing or equal to NaN,
%           the default value is used.
%
% Example: 73 coefficient minimum-phase filter, sampling frequency 48000
%   Hz.
%     [B(1:2).Freq]   = deal([0 3400], [4500 24000]);
%     [B(1:2).Value]  = deal(1, 0);
%     [B(1:2).Weight] = deal(1, 50);
%     hM = DFiltMPFIR(73, B, 48000);

%  $Id: DFiltMPFIR.m,v 1.11 2009/07/07 20:24:13 pkabal Exp $

% Defaults
if (~ exist ('SFreq', 'var'))
  SFreq = 1;
end
if (~ exist ('GridD', 'var'))
  GridD = NaN;
end

% Find the nominal number of grid points; the actual number will be
% higher since a higher density is used near band edges
NGRID_MIN = 1000;
GRIDD_DEF = 32;
if (isnan (GridD))
  NGrid = floor (max (NGRID_MIN, (NCof - 1) * GRIDD_DEF));
  GridD = NGrid / (NCof - 1);
end

FT = 1;

% Assemble the filter specifications into band-by-band structures
BSpecDef.ValueInt  = 'sqrt';
BSpecDef.WeightInt = 'sqrt';
BSpecDef.ULimitInt = 'sqrt';
BSpecDef.LLimitInt = 'sqrt';
B = MMBSpec (BSpec, SFreq, BSpecDef);

% Force the lower limits to be 0
NBand = length (B);
for (k = 1:NBand)
  B(k).LL = max (0, B(k).LL);
  B(k).UL = max (0, B(k).UL);
end

% Double-length filter 
NCofDL = 2 * NCof - 1;

% % Check and fill in missing band specifications
Gc = MMdSpec (FT, NCofDL, B, GridD);

% Filter design
[h, ier, Ex] = MMdesLPFIR (FT, NCofDL, Gc);
if (ier ~= 0)
  error ('DFiltFIR: Error during filter design');
end

% Print the double-order filter specifications
MMprintSpec (FT, NCofDL, Ex.dev, B);

% Factor the filter
fzeros = acos (Ex.x(Ex.EType == -2 & Ex.y == 0));
hM = MPfactorFilter (h, fzeros);
hM = hM(:);

% Print the min-phase filter specifications
MPprintSpec (NCof, Ex.dev, B);

return