crmz.m

有关matlab的电子书籍有一定的帮助希望有用

function [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 )
%CRMZ Complex Remez multiple-exchange filter design algorithm.
%   Designs FIR filters with arbitary magnitude and 
%   phase specifications; reduces to the Remez (Parks-McClellan)
%   algorithm in the real case.

%   Authors: L. Karam, J. McClellan
%   Revised: 22-Oct-96, D. Orofino
%
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.8 $  $Date: 1997/12/02 20:47:50 $

global DES_CRMZ WT_CRMZ GRID_CRMZ TGRID_CRMZ IFGRD_CRMZ

J = 1i;
N1 = 0;

if nargin==1 & isstruct(L),
 if isstr(L.plot) & strcmp(L.plot,'plot-result'),
    N2 = N1 + L.L - 1;
    H = L.HH .* exp(-J*2*pi*TGRID_CRMZ*(N1+N2)/2);
    crmz_plresult( H, L.EE, L.iext, L.iext, L.indx_edges, N1, N2, L.delta, L.sym );
    drawnow
 end
 return
end

if TRACE,
  disp(' ');
  disp('             ***********************************************');
  disp('             *****           Complex  Remez            *****');
  disp('             ***********************************************');
end

if nargin<2, sym = 0; end;

N2 = N1 + L - 1;   % Last index of impulse response 
is_odd = rem(L,2); 
if (is_odd)
  Lf = (L+1)/2;
  Ws = ['W(:,2:Lf)'];
  hr_str = ['[a([Lc:-1:2])/2; a(1); a([2:Lc])/2]'];
  hi_str = ['[-a([Lb:-1:(Lc+1)])/2;    0; a([(Lc+1):Lb])/2]'];
  ph_str = ['ones(length(TGRID_CRMZ),1)'];
else
  Lf = L / 2;
  Ws = ['W'];
  hr_str = ['[a([Lc:-1:1])/2; a([1:Lc])/2]'];
  hi_str = ['[-a([Lb:-1:(Lc+1)])/2; a([(Lc+1):Lb])/2]'];
  ph_str = ['exp(-J*pi*TGRID_CRMZ)'];
end;
Lc = Lf;
Ls = L - Lf;

switch sym
case 0,
  Lb = L;
  M_str = ['[cos(W)  sin(' Ws ')]'];
  h_str = [hr_str '+J*' hi_str];
  HH_str = ['HH = fftshift( fft(hc,Lfft) ).*' ph_str ';'];
case 1,
  % Even-symmetry:
  Lb = Lc;
  M_str = ['cos(W)'];
  h_str = [hr_str];
  HH_str = ['HH=fft(hc,Lfft); HH((Lfft/2+2):Lfft)=[]; HH=HH.* ' ph_str ';'];
case 2,
  % Odd-symmetry:
  Lb = Ls;
  Lc = 0;
  M_str = ['sin(' Ws ')'];
  h_str = ['J*' hi_str];
  HH_str = ['HH=fft(hc,Lfft); HH((Lfft/2+2):Lfft)=[]; HH=HH.* ' ph_str ';'];
end

%------------------------
A = DES_CRMZ .* exp(J*2*pi*GRID_CRMZ*(N1+N2)/2);
if PLOTS
  % clf;
  plot(GRID_CRMZ*2,abs(DES_CRMZ),'--')
  hold on
  plot(GRID_CRMZ*2, angle(DES_CRMZ),'-')
  hold off 
  title('Desired Filter Frequency Response')
  xlabel('normalized frequency') 
  ylabel('magnitude (dashed) and phase (solid)')
  drawnow
end

%-----------------------
if TRACE, disp('Calculating initial solution .....'), end

vec_edges    = TGRID_CRMZ(indx_edges);
edges        = reshape(vec_edges,2,length(vec_edges)/2)';
[fext, iext] = crmz_guess(edges, GRID_CRMZ, Lb);

it = 0; delta = 0.0; delta_old = -1; no_stp = 1;

exactTol = 1e-15;  % tolerance for exactness (arbitrary)
% if the maximum error relative to the absolute maximum of the
%  desired function is less than exactTol, the filter designed
%  is considered 'exact' and the iteration is terminated successfully.

while (no_stp)
  it = it + 1;
  delta_old = abs(delta);
  fext = GRID_CRMZ(iext);
  W    = 2*pi*fext*([0:(Lf-1)]+(~is_odd)*0.5);
  Mb   = eval(M_str);
  M    = [ Mb ((-1).^[0:Lb]')./WT_CRMZ(iext) ];
  a    = M \ A(iext);
  delta = a(Lb+1);
  h  = eval(h_str);
  hc = crmz_rotate(zeropad(h,Lfft-L), -Ls);
  eval(HH_str);
  W  = 2*pi*vec_edges*([0:(Lf-1)]+(~is_odd)*0.5);
  Mb = eval(M_str);

  HH(indx_edges) = Mb * a(1:Lb);
  EE = WT_CRMZ .* (A - HH(IFGRD_CRMZ));
  EE(iext)   = delta*( (-1) .^ [2:length(iext)+1] ); 
  [jext,EEj] = crmz_find(EE,iext);  
  if TRACE,
    fprintf('(%3.0f): ', it) 
    fprintf('delta = %8.5f+j%8.5f, ', real(delta), imag(delta)) 
    fprintf('|delta| = %5.2e, ', abs(delta));
    fprintf('e_max = %5.2e\n', max(abs(EE)));
  end
  if PLOTS,
    crmz_plerror(EE,HH,iext,jext,indx_edges,GRID_CRMZ,delta,sym);
    drawnow
  end
  if all(iext == jext') | (max(abs(EE))/max(abs(A))<exactTol)
    no_stp = 0;
  end
  iext = jext';
end
%%% end of Complex Remez algorithm %%%%

%%%%% Assessing Optimality of the Solution %%%%
e_max = max(abs(EE));
tlr = abs(delta)/100;  %% Needed due to limited machine accuracy %% 
if (e_max <= (abs(delta)+tlr)) | (e_max/max(abs(A))<exactTol)
  not_optimal = 0; 
  if TRACE,
    fprintf('Optimal Solution Obtained !\n')
    fprintf('(%3.0f): ', it) 
    fprintf('delta = %8.5f+j%8.5f, ', real(delta), imag(delta)) 
    fprintf('|delta| = %5.2e, ', abs(delta));
    fprintf('e_max = %5.2e\n',max(abs(EE)));
  end
else
  not_optimal = 1; 
  if TRACE
    crmz_chckopt( EE, delta, tlr, Lfft )
    fprintf('(%3.0f): ', it) 
    fprintf('delta = %8.5f+j%8.5f, ', real(delta), imag(delta)) 
    fprintf('|delta| = %5.2e, ', abs(delta));
    fprintf('e_max = %5.2e\n',max(abs(EE)));
  end
end

%%%%%%%%%%%%%  End Complex Remez Stage  %%%%%%%%%%%%%%%%%%%%%%%%

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

function crmz_chckopt( EE, delta, tlr, Lfft )
% Checks degree of suboptimality of solution not optimal 
% on desired frequency bands. 
%
subint = find( abs(EE) <= (abs(delta)+tlr) );
rel_size = round((length(subint)/Lfft)*100);
disp(['Optimal solution found on a frequency subset comprising ' ...
      sprintf('%3.0f percent \n',rel_size) 'of the original frequency set'])

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

function  [fnew,enew] = crmz_find( error, fold )
%  Construct new subset of points using exchange rules

%    Authors: LJ Karam and JH McClellan
%  Reference: "Complex Chebyshev approximation for FIR filter design",
%             IEEE Trans. on Circuits and Systems II, March 1995. 

fold      = fold(:);
Nx        = length( fold );
Ngrid     = length( error );
if error(fold(1))~=0
    sgn_error = error(fold(1))/abs(error(fold(1)));
else
    sgn_error = 1;
end
error     = real( conj(sgn_error)*error );
delta     = min( abs( error(fold) ) ); %--- present value of delta
up        = sign(error(fold(1))) > 0;
if up,
  fence = [ [1;fold(1:2:Nx)]  [fold(1:2:Nx);Ngrid] ];
else
  fence = [ [1;fold(2:2:Nx)]  [fold(2:2:Nx);Ngrid] ];
end
Lf  = length(fence(:,1));
emn = zeros(Lf,1); imn = emn; emx = emn; imx = emn;
for i=1:Lf,
   [emn(i),imn(i)] = min( error(fence(i,1):fence(i,2) ));
   imn(i) = imn(i) + fence(i,1) -1;
end
imn(find((emn > -delta))) = [];
fence = [ [1;imn]  [imn;Ngrid] ];
Lf = length(fence(:,1));
for i=1:Lf
   [emx(i),imx(i)] = max( error(fence(i,1):fence(i,2) ));
   imx(i) = imx(i) + fence(i,1) -1;
end
imx(find((emx < delta))) = [];
fnew = sort([imx;imn]);
Nf   = length(fnew);
if ( Nf > Nx)
 if ( abs(error(fnew(1))) >= abs(error(fnew(Nf))) )
   fnew = fnew(1:Nx);
 else
   fnew = fnew(Nf-Nx+1:Nf);
 end
end
fnew = fnew';
enew = error(fnew);

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

function  [fext,iext] = crmz_guess( edges, grid, nfcns )
%CRMZ_GUESS   generate initial guess of "EXTREMAL FREQS"
%               for remez exchange algorithm
%   usage:
%        [Fext,Iext] = crmz_guess( Edges, Grid, Nfcns )
%
%     Edges:  band-edges moved onto the grid (Revised edges)
%     Grid:   frequencies on the fine grid
%     Nfcns:  number of basis functions in the approx
%     Fext:   initial guess of "extremal frequencies"
%     Iext:   indices for the "extremal frequencies"
%
TOL    = 5*eps;
next   = nfcns+1;
Nbands = length(edges(:,2));
tt     = edges;
merged = tt;
if Nbands > 1,
  jkl = find( abs(tt(1:(Nbands-1),2)-tt(2:Nbands,1)) > TOL );
  merged = [ [tt(1,1); tt(jkl+1,1)] , [tt(jkl,2); tt(Nbands,2)] ];
end
Nbands = length(merged(:,1));
bw = merged(:,2)-merged(:,1);
if any(bw < 0),
   edges
   error('Internal error: obtained a negative bandwidth');
end
percent_bw = bw / sum(bw);
fext = zeros(next,1);
n = 0;  i = 1;
while (n < next) & (i <= Nbands),
   nfreqs_i = min( next-n, ceil( percent_bw(i)*next ) );
   if( nfreqs_i == 0 )
      n=n+1;   fext(n) = merged(i,1);
   else
      fext(n+[1:nfreqs_i]) = ...
           linspace(merged(i,1),merged(i,2),nfreqs_i);
      n = n+nfreqs_i;
   end