DFiltNyquist (M, NCof, BSpec, SFreq)

function h = DFiltNyquist (M, NCof, BSpec, SFreq)
% This program designs an equiripple Nyquist filter. The design is
% specified in terms of a zero crossing interval and a stopband weighting.
% The resulting filter will be a minimax approximation to the stopband
% specification with the constraint of the regular zero crossings. The
% filter is linear phase, with even symmetry. The total number of
% coefficients must be odd. The coefficients are normalized such that the
% middle coefficient is 1. This means that the passband gain approximates
% M, where M is the zero crossing interval.
%
%   h:     Filter coefficients (column vector)
%
%   M:     Zero crossing interval in samples
%   NCof:  Number of filter coefficients (must be odd)
%   BSpec: Stopband specifications. Normally this is a single band, with
%          the lower band edge greater than 0.5 SFreq/M and upper band edge
%          equal to 0.5 SFreq. 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.5*SFreq/M, 0.5*SFreq].
%          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.
%          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 -Inf.
%          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).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 'linear'.
%          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 'linear'.
%   SFreq: Sampling frequency (used to interpret the frequencies in the
%           filter specifications). If not specified, SFreq = 1.
% Example: Nyquist filter with 15 coefficients, M = 2
%   The sampling frequency is 16 kHz and stopband starts at 4500 Hz. The
%   stopband weights vary from 1 to 10 to give better attenuation at the
%   top of the stopband.
%     B.Freq   = [4500 8000];
%     B.Weight = [1, 10];
%     h = DFiltNyquist(2, 15, B, 16000);

% $Id: DFiltNyquist.m,v 1.5 2011/07/20 15:32:16 pkabal Exp $

% Check for DFiltFIR
if (~ exist('DFiltFIR', 'file'))
  error('Cannot access routines from the DFiltFIR package');
end

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

% Filter length
if (mod(NCof, 2) == 0)
  error ('DFiltNyquist: Number of coefficients must be odd');
end

% Let H(z) = H0(z) H1(z), where H0(z) and H1(z) are linear phase FIR
% filters. H1(z) controls the stopband. H0(z) forces the zero crossings
% into the response.
% Set up the specifications for H1(z)
B = NYQ_BSpec(BSpec, M, SFreq);

% Nyquist design
[h, Ex] = NYQdesNyquist(M, NCof, B, GridD);

% Print the specifications
NYQprintSpec(M, NCof, Ex.dev, B);

return

% ----- -----
function B = NYQ_BSpec (BSpec, M, SFreq)

% Check frequency points
FreqSB = [BSpec(:).Freq];
if (min(FreqSB) < 0.5*SFreq/M)
  error ('DFiltNyquist: Invalid lower stopband frequency edge');
end
if (max(FreqSB) < 0.5*SFreq)
  warning ('DFiltNyquist - Upper stopband edge below 0.5*SFreq');
end

% Add a passband point (add at end and then move to beginning)
NB = length(BSpec) + 1;
BSpec(NB).Freq = 0;
BSpec = [BSpec(NB) BSpec(1:NB-1)];

BSpec(1).Value = 1;
if (isfield (BSpec, 'Weight'))
  BSpec(1).Weight = 1;
end
BSpec(1).Value = 1;
BSpec(1).LLimit = 1;
BSpec(1).ULimit = 1;
BSpec(2:end).Value = 0;

% Check and fill in missing band specifications
B = MMBSpec(BSpec, SFreq);

return