function [h, Ex] = DFiltFIR (NCof, BSpec, FType, SFreq, GridD)
% [h, Ex] = DFiltFIR (NCof, BSpec, FType, SFreq, GridD)
% This program designs four types of linear phase finite impulse response
% filters.
% 1: Equiripple Linear Phase Filter (multiple passbands and stopbands).
% The filter is specified in terms of the desired response in bands.
% The band specifications include desired values, relative weighting
% and limits on the allowed response values in the band. The resulting
% filter is a weighted minimax approximation (with constraints) to the
% given specifications. The filter coefficients have an even symmetry
% about the middle of the filter.
% 2: Equiripple Linear Phase Filter (multiple passbands and stopbands)
% with sin(x)/x compensation. The filter is specified in terms of the
% desired response in bands. The band specifications include desired
% values, relative weighting and limits on the allowed response
% values in the band. The resulting filter is a weighted minimax
% approximation (with constraints) to the given specifications. The
% filter coefficients have an even symmetry about the middle of the
% filter.
% 3: Equiripple Differentiator. The filter is specified in terms of the
% desired slopes in bands. The band specifications include desired
% slope, relative weighting and limits on the allowed slopes in the
% band. The resulting filter is a weighted minimax approximation (with
% constraints) to the given specifications. The filter coefficients
% have an odd symmetry about the middle of the filter.
% 4: Equiripple Hilbert Transform Filter. The filter is specified in
% terms of the desired response in bands. The band specifications
% include desired values, relative weighting and limits on the allowed
% values in the band. The resulting filter is a weighted minimax
% approximation (with constraints) to the given specifications. The
% filter coefficients have an odd symmetry about the middle of the
% filter.
%
% 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.
%
% The filter will have weighted deviations from the desired values which
% have equal peak values. The weighting is specified as a function of
% frequency. Higher (relative) weights mean smaller deviations.
%
% Multiple desired values, weights, and limits can be specified for a
% band. Monotonic cubic interpolation is used between the values. The
% interpolation can be done in the linear, square root, or log domain.
% If a single value, weight, or limit is given for a band, the value,
% weight or limit is constant for the band. If two values, weights, or
% limits are given for a band, the interpolation is linear (in the
% linear, square root, or log domain).
%
% h: Filter coefficients (column vector)
% Ex: Structure containing the positions and values of the extrema.
% The positions of the extrema are in x units, where x=cos(omega).
%
% NCof: Number of filter coefficients
% BSpec: Structure array containing the filter specifications, 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.
% 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).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 'linear'.
% 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'.
% FType: Filter type (optional)
% 'bpf' - bandpass filter (default)
% 'rec' - receive filter (sin(x) / x compensation)
% 'dif' - differentiator
% 'hil' - Hilbert transform filter
% SFreq: Sampling frequency (used to interpret the frequencies in the
% filter specifications). If not specified, 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 16 points per extremum. If this value is missing or
% equal to NaN, the default value is used.
%
% Example: 32 tap constrained bandpass filter:
% The stopbands are from 0 to 0.1 and 0.425 to 0.5, and the passband is
% from 0.2 to 0.35. The relative weight for the first stopband is 10;
% the weight for the passband is 1, and the weight for the second
% stopband goes linearly from 10 to 20. The response for the first
% stopband is constrained to be positive.
% [B(1:3).Freq] = deal ([0 0.1], [0.2 0.3], [0.425 0.5]);
% [B(1:3).Value] = deal (0, 1, 0);
% [B(1:3).Weight] = deal (10, 1, [10 20]);
% B(1).LLimit = 0;
% h = DFiltFIR (32, B, 'bpf');
%
% Notes:
% Limits are signed values which apply to the zero-phase response. The
% zero-phase response is the real response after the linear phase term
% has been removed. The zero-phase response can go negative. For example,
% consider a lowpass response. The passband approximates a constant
% (usually set to be positive). The stopband oscillates around zero. The
% limits apply to this response. For instance setting the lower limit
% to be zero in the stopband will prevent the response from going
% negative in the stopband.
% $Id: DFiltFIR.m,v 1.12 2009/07/07 20:24:49 pkabal Exp $
% Defaults
if (~ exist ('FType', 'var'))
FType = 'bpf';
end
if (~ exist ('SFreq', 'var'))
SFreq = 1;
end
if (~ exist ('GridD', 'var'))
GridD = NaN;
end
% Decode the filter type
FT = strmatch (lower (FType), {'bpf', 'rec', 'dif', 'hil'});
if (isempty (FT))
error ('DFiltFIR: Invalid filter type');
end
% Check and fill in missing band specifications
B = MMBSpec (BSpec, SFreq);
% Convert to a specification on a dense grid of frequencies
Gc = MMdSpec (FT, NCof, B, GridD);
% Filter design
[h, ier, Ex] = MMdesLPFIR (FT, NCof, Gc);
if (ier ~= 0)
error ('DFiltFIR: Error during filter design');
end
% Print the filter specifications
MMprintSpec (FT, NCof, Ex.dev, B);
return
