cremez.m

matlabDigitalSigalProcess内有文件若干

function [h,delta,result] = cremez(M, edges, filt_str, varargin)
%CREMEZ Complex and nonlinear phase equiripple FIR filter design.
%  CREMEZ allows arbitrary frequency-domain constraints to be
%  specified for the design of a possibly complex FIR filter.  The 
%  Chebyshev (or minimax) filter error is optimized, producing
%  equiripple FIR filter designs.
%
%  B=CREMEZ(N,F,'fresp',W) returns a length N+1 FIR filter which
%  has the best approximation to the desired frequency response
%  as returned by function 'fresp'.  The function is called
%  from within CREMEZ using the syntax:
%
%                    [DH,DW] = fresp(N,F,GF,W);
%  where:
%  N is the filter order.
%  F is the vector of frequency band edges which must appear
%    monotonically between -1 and +1, where 1 is the Nyquist
%    frequency.  The frequency bands span F(k) to F(k+1) for k odd;
%    the intervals  F(k+1) to F(k+2) for k odd are "transition bands"
%    or "don't care" regions during optimization.
%  GF is a vector of grid points which have been linearly interpolated
%    over each specified frequency band by CREMEZ, and determines the
%    frequency grid at which the response function will be evaluated.
%  W is a vector of real, positive weights, one per band, for use
%    during optimization.  W is optional; if not specified, it is set
%    to unity weighting before being passed to 'fresp'.
%  DH and DW are the desired complex frequency response and
%    optimization weight vectors, respectively, evaluated at each
%    frequency in grid GF.
%
%  Predefined frequency response functions for 'fresp' include:
%       'lowpass'  'bandpass' 'multiband'      'hilbfilt'
%       'highpass' 'bandstop' 'differentiator' 'invsinc'
%  See the help for PRIVATE/LOWPASS, etc., for more information.
%
%  B=CREMEZ(N,F,{'fresp',P1,P2,...},W) specifies optional arguments
%  P1, P2, etc., to be passed to the response function 'fresp'.
%
%  B=CREMEZ(N,F,A,W) is a synonym for B=CREMEZ(N,F,{'multiband',A},W),
%  where A is a vector of response amplitudes at each band edge in F.
%
%  B=CREMEZ(..., SYM) imposes a symmetry constraint on the impulse
%  response of the design, where SYM may be one of the following:
%          'none' - Default if any negative band edge frequencies are
%                   passed, or if 'fresp' does not supply a default.
%          'even' - Impulse response will be real and even.  This is
%                   the default for highpass, lowpass, bandpass,
%                   bandstop, and multiband designs.
%           'odd' - Impulse response will be real and odd.  This is
%                   the default for Hilbert and differentiator designs.
%          'real' - Impose conjugate symmetry on frequency response.
%  Each frequency response function 'fresp' provides a default value
%  for SYM; see help on private/lowpass, etc., for more information.
%  If any SYM option other than 'none' is specified, the band edges
%  should only be specified over positive frequencies; the negative
%  frequency region will be filled in from symmetry.
%
%  B=CREMEZ(..., 'skip_stage2') disables the second-stage optimization
%  algorithm, which executes only when CREMEZ determines that an
%  optimal solution has not been reached by the standard Remez
%  error-exchange.  Disabling this algorithm may increase the speed of
%  computation, but with a reduction in accuracy.  By default, the
%  second-stage optimization is enabled.
%
%  B=CREMEZ(..., DEBUG) enables the display of intermediate results
%  during the filter design, where DEBUG may be one of 'trace',
%  'plots', 'both', or 'off'.  By default, DEBUG is set to'off'.
%
%  B=CREMEZ(...,{LGRID}), where {LGRID} is a one-by-one cell array containing
%  an integer, controls the density of the frequency grid. The frequency grid 
%  size is roughly 2^nextpow2(LGRID*N).  LGRID defaults to 25.
%
%  Note that any combination of the SYM, DEBUG, 'skip_stage2', and {LGRID}
%  options may be specified.
%
%  [B,ERR]=CREMEZ(...) returns the maximum ripple height ERR.
%
%  [B,ERR,RES]=CREMEZ(...) returns a structure RES of optional results
%  computed by CREMEZ, and contains the following fields:
%
%     RES.fgrid: vector containing the frequency grid used in
%                the filter design optimization
%       RES.des: desired response on fgrid
%        RES.wt: weights on fgrid
%         RES.H: actual frequency response on the grid
%     RES.error: error at each point on the frequency grid
%     RES.iextr: vector of indices into fgrid of extremal frequencies
%     RES.fextr: vector of extremal frequencies
%
%  See also REMEZ, FIR1, FIRLS, FILTER,  PRIVATE/LOWPASS,
%  PRIVATE/HIGHPASS,   PRIVATE/BANDPASS, PRIVATE/BANDSTOP,
%  PRIVATE/MULTIBAND,  PRIVATE/INVSINC,  PRIVATE/HILBFILT,
%  PRIVATE/DIFFERENTIATOR.

%   Authors: L. Karam, J. McClellan
%   Revised: October 1996, D. Orofino
%
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.10 $  $Date: 1997/12/02 18:36:01 $

% NOTE: This algorithm is equivalent to Remez for real B
%       when the filter specs are exactly linear phase.

% Declare globals:
global DES_CRMZ WT_CRMZ GRID_CRMZ TGRID_CRMZ IFGRD_CRMZ

if nargin<3, error('Not enough input arguments.'); end

L         = M+1;             % convert ORDER to LENGTH
edges     = edges/2;         % compensate for [-1,1) frequency normalization
num_bands = length(edges)/2; % # of frequency bands specified

% Some quick checks on band edge vector:
if (num_bands ~= floor(num_bands)),
  error('F must contain an even number of band edge entries.');
end
if any(diff(edges) <= 0),
   error('Band edges must be monotonically increasing.')
end

% Assign default parameter values:
wgts         = ones(1, num_bands);
h_sym        = 'unspecified';
plot_flag    = 'off';
allow_stage2 = 1;
grid_density = 25;

plot_opts = {'plots','trace','both','off'};
sym_opts  = {'even','odd','real','none'};
for ii=1:length(varargin),
  if iscell(varargin{ii})
     grid_density = varargin{ii}{1};
  elseif ~isstr(varargin{ii}),
     if ii==1,
       wgts = varargin{ii};
     else
       error('Expecting a string argument.');
     end
  else
    ppv = lower(varargin{ii});
    switch ppv
    case plot_opts
      plot_flag = ppv;      % plot option
    case sym_opts
      h_sym = ppv;          % symmetry option
    case 'skip_stage2'
      allow_stage2 = 0;     % skip 2nd stage optimization
    otherwise
      error(['Invalid argument "' ppv '" specified.']);
    end
  end
end

PLOTS = any(strcmp(plot_flag,{'plots','both'}));
TRACE = any(strcmp(plot_flag,{'trace','both'}));

% Force function name into a cell-array:
if ~iscell(filt_str)
   if isstr(filt_str)
      % Function name passed:
      filt_str = {filt_str};
   else
      % User passed a non-string, non-cell -- assume it to be a
      % magnitude vector.  Pass this to multiband:
      filt_str = {'multiband', filt_str};
   end
end
filt_call    = lower(filt_str{1});  % string name of function
other_params = filt_str(2:end);     % cell_array of additional args

% Determine symmetry options:
h_sym = lower(h_sym);
if strcmp(h_sym,'unspecified'),

   % If symmetry was not specified, AND there are negative freqencies
   % passed in the edge vector, use 'none':
   if any(edges < 0),
     h_sym = 'none';
   else
     % 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.
     h_sym = feval(filt_call, 'defaults', ...
                  {M, 2*edges, [], wgts, other_params{:}} );
     if ~any(strcmp(h_sym,sym_opts)),
       error(['Invalid default symmetry option "' h_sym '" returned ' ...
              'from response function "' filt_call '".  Must be one ' ...
              'of ''none'',''real'',''even'', or ''odd''']);
     end
  end
end

if TRACE,
  disp(['   Symmetry option: ' h_sym]);
end

if any(strcmp(h_sym,{'real','none'})),
  fdomain = 'whole';     % [-0.5,0.5)
  sym = 0;               % h_sym = 'real' or 'none'
else
  fdomain = 'half';      % [0,0.5]
  if strcmp(h_sym,'even')
    sym = 1;           % h_sym = 'even'
  else
    sym = 2;           % h_sym = 'odd'
  end
end

% Check domain before generating frequency grid:
if strcmp(fdomain, 'whole'),
  % Domain is [-1,+1) for user input, [-0.5,0.5) internally
  if any(edges < -0.5 | edges > 0.5),
    error(['Frequency band edges must be in the range [-1,+1] for ' ...
           'designs with SYM=''' h_sym '''.']);