cremez.m

matlabDigitalSigalProcess内有文件若干

  end
else
  % Domain is [0,+1] for user input, [0,.5] internally
  if any(edges < 0 | edges > 0.5),
    error(['Frequency band edges must be in the range [0,+1] for ' ...
           'designs with SYM=''' h_sym '''.']);
  end
end

% Generate frequency grid:
edge_pairs = reshape(edges, 2, num_bands)';
[tgrid, Lfft, vec_edges, indx_edges] = ...
               crmz_grid(edge_pairs, L, fdomain, grid_density);

%  Input: indx_edges = [a1 a2 a3 a4 ...]
% Output: IFGRD_CRMZ = [a1:a2 a3:a4 ...]
IFGRD_CRMZ = [];
for jj = 1:2:length(indx_edges),
  IFGRD_CRMZ = [IFGRD_CRMZ indx_edges(jj):indx_edges(jj+1)];
end

% Get just those points from the full grid (tgrid) which correspond
% to the frequency band intervals delimited by indx_edges (i.e., by
% edge_pairs, quantized to grid indices):
TGRID_CRMZ = tgrid(:);
GRID_CRMZ  = TGRID_CRMZ( IFGRD_CRMZ );
if (max(GRID_CRMZ) > edges(end)) | (min(GRID_CRMZ) < edges(1)),
  error('Internal error: Grid frequencies out of range.')
end

% Get desired frequency characteristics at the frequency points
% in the specified frequency band intervals:
%
% NOTE! The user needs to see normalized frequencies in the range
% [-1,+1], and not [-0.5,+0.5] as we use internally.
[DES_CRMZ, WT_CRMZ] = feval(filt_call, ...
                            M, 2*edges, 2*GRID_CRMZ, wgts, other_params{:}); 
% Cleanup the results, and check sizes:
DES_CRMZ = DES_CRMZ(:);
WT_CRMZ  = WT_CRMZ(:);
if ~isequal(size(DES_CRMZ), size(GRID_CRMZ)) | ...
   ~isequal(size(WT_CRMZ),  size(GRID_CRMZ)),
  str = ['Incorrect size of results from response function "' ...
         filt_call '".  Sizes must be the same size as the ' ...
         'frequency grid GF.'];
  error(str);
end

if strcmp(h_sym,'real') & all(edges >= 0)
   % We need to make DES and WT conjugate symmetric
   %error(['Frequency band edges must be specified over the entire ' ...
   %       'interval [-1,+1) for designs with SYM=''' h_sym '''.']);

   % crmz_grid moved the band edge grid points, so do the same
   % when constructing symmetric spectrum:
   len = length(TGRID_CRMZ);
   % If the DC term is included in the band edges, remove it:
   if any(indx_edges==len/2+1),
     indx_edges(1)=[]; % Throw away DC point
   end
   q = len+2-flipud(indx_edges(:));
   TGRID_CRMZ(flipud(q)) = -TGRID_CRMZ(indx_edges);
   % Adjust other grid vectors accordingly:
   indx_edges = [q; indx_edges(:)];
   IFGRD_CRMZ = [len+2-fliplr(IFGRD_CRMZ) IFGRD_CRMZ];
   GRID_CRMZ  = TGRID_CRMZ(IFGRD_CRMZ);
   % Now, impose conjugate symmetry:
   DES_CRMZ   = [conj(flipud(DES_CRMZ)); DES_CRMZ];
   WT_CRMZ    = [conj(flipud(WT_CRMZ));  WT_CRMZ];
end

% Complex Remez Stage:
[h,a,delta,not_optimal,iext,HH,EE,M_str,HH_str,h_str,Lf,Lb,Ls,Lc,A] = ...
             crmz( L, sym, Lfft, indx_edges, PLOTS, TRACE ); 
if PLOTS,
  plot_struct.L          = L;
  plot_struct.HH         = HH;
  plot_struct.EE         = EE;
  plot_struct.iext       = iext;
  plot_struct.indx_edges = indx_edges;
  plot_struct.sym        = sym;
  plot_struct.delta      = delta;
  plot_struct.plot       = 'plot-result';
  crmz( plot_struct );   % generate final plot
end

if not_optimal & allow_stage2,
  % Ascent-descent Stage:

  if TRACE,
     disp('           ***********************************************');
     disp('           *****    Invoking Second Ascent Stage      ****');
     disp('           ***********************************************');
  end

  [h,a,delta,HH,EE] = ...
       adesc( L, Lf, Lb, Ls, Lc, sym, indx_edges, iext, HH, EE, a,...
	      M_str, HH_str, h_str, A, delta, PLOTS, TRACE );
  if PLOTS,
    plot_struct.L          = L;
    plot_struct.HH         = HH;
    plot_struct.EE         = EE;
    plot_struct.iext       = iext;
    plot_struct.indx_edges = indx_edges;
    plot_struct.sym        = sym;
    plot_struct.delta      = delta;
    plot_struct.plot       = 'plot-result';
    crmz( plot_struct );   % generate final plot
  end
end

% Return a row-vector, and remove imag part if it's small:
h = h(:).';
if ~isreal(h) & (norm(imag(h)) < 1E-12*norm(real(h))),
  h = real(h);
end

if nargout>2,
  result.fgrid = 2*GRID_CRMZ;
  result.des   = DES_CRMZ;
  result.wt    = WT_CRMZ;
  result.H     = HH(IFGRD_CRMZ);
  result.error = EE;
  result.iextr = iext;
  result.fextr = 2*GRID_CRMZ(iext);
end

clear global DES_CRMZ WT_CRMZ GRID_CRMZ TGRID_CRMZ IFGRD_CRMZ

%========================================================

function  [grid,Ngrid,edge_vec,edge_idx] = crmz_grid(edge_pairs, L, fdomain, ...
                                                grid_density)
% CRMZ_GRID Generate grid for Remez exchange.
%
%   [Grid, Ngrid, V_edges, Indx_edges] = ...
%                    CRMZ_GRID(Edge_Pairs, L, fdomain, grid_density)
%           Grid: frequencies on the fine grid
%          Ngrid: number of grid pts if whole domain used [-1,+1]
%       edge_vec: specified edges reshaped into a vector of adjacent edge-pairs
%       edge_idx: index of specified band-edges within grid array
%
%     edge_pairs: specified band-edges, [Nbands X 2] array
%              L: filter length
%        fdomain: domain for frequency approximation
%                 default is [0,0.5] called 'half'
%                 'whole' means [-0.5,0.5), 
%   grid_density: density of grid

error(nargchk(3,4,nargin));

if (edge_pairs(1) == -0.5) & (edge_pairs(end) == 0.5),
   % -pi and +pi are the same point - move one a little:
   new_freq = 0.5 - 1/(50*L);
   if new_freq <= edge_pairs(end-1),
     % Last two freq points are too close to move - try first two:
     new_freq = -new_freq;
     if (new_freq >= edge_pairs(2)),
       error(['Both -1 and 1 have been specified as frequencies ' ...
              'in F, and the frequency spacing is too close to ' ...
              'move either of them toward its neighbor.']);
     else
       edge_pairs(1) = new_freq;
     end
   else
     edge_pairs(end) = new_freq;
   end
end

Ngrid = 2.^nextpow2(L*grid_density);
if (Ngrid/2) > 20*L,
   Ngrid = Ngrid/2;
end

del_f    = 1./Ngrid;
edge_vec = edge_pairs';  % M-by-2 to 2-by-M
edge_vec = edge_vec(:);  % single column of adjacent edge-pairs

switch fdomain
case 'whole'
  grid = (0:Ngrid-1)./Ngrid - 0.5;             % uniform grid points [-.5,.5)
  edge_idx = 1 + round((edge_vec+0.5)*Ngrid);  % closest indices in grid
case 'half'
  grid = (0:Ngrid/2)./Ngrid;                   % uniform grid points [0,.5]
  edge_idx = 1 + round(edge_vec*Ngrid);        % closest indices in grid
otherwise
  error('Internal error: domain must be "whole" or "half".');
end
edge_idx(end) = min(length(grid), edge_idx(end));  % Clip last index

% Fix repeated edges:
% This determines not only if
%         i(1)==i(2) and i(3)==i(4),
% but also if
%         i(2)==i(3), etc.
m = find(edge_idx(1:end-1) == edge_idx(2:end));
if ~isempty(m),
  % Replace REPEATED band edges with the uniform grid points
  % Could be a problem if [-1 -1] (if whole) or [0 0] (if half) specified
  edge_idx(m) = edge_idx(m) - 1;   % move 1 index lower
  edge_vec(m) = grid(edge_idx(m)); % change user's band edge accordingly
  m = m + 1;
  edge_idx(m) = edge_idx(m) + 1;
  edge_vec(m) = grid(edge_idx(m));
end

% Replace closest grid points with exact band edges:
grid(edge_idx) = edge_vec;

% end of crmz_grid