%%%%%%%%%% crmz.m (Complex Remez) %%%%%%%%%% 
% crmz.m
% 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. 
%
% @(#) crmz.m: version 5.0 created on 02/95                        
% (c) 1995 by Lina Karam and James McClellan
%
% The author can be reached at:
%
% Lina Karam
% School of Electrical and Computer Engineering
% Goergia Institute of Technology 
% Atlanta, GA 30332
% (404) 894-2910
% karam@eedsp.gatech.edu
%main program
disp(' ');
disp('             ***********************************************');
disp('             *****           Complex  Remez            *****');
disp('             ***********************************************');
if ~exist('sym'), sym = 0; end;
if ~exist('N1'), N1 = 0; end; 
N2 = N1 + L - 1;   % Last index of impulse response 
is_odd = mod(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),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)'];
end;
Lc = Lf;
Ls = L - Lf;
if (~sym)
  Lb = L;
  M_str = ['[cos(W)  sin(' Ws ')]'];
  h_str = [hr_str '+J*' hi_str];
  HH_str = ['HH = fftshift( fft(hc,Ng) ).*' ph_str ';'];
elseif (sym == 1) 
  Lb = Lc;
  M_str = ['[cos(W)]'];
  h_str = [hr_str];
  HH_str = ['HH=fft(hc,Ng); HH((Ng/2+2):Ng)=[]; HH=HH.* ' ph_str ';'];
elseif (sym == 2) 
  Lb = Ls;
  Lc = 0;
  M_str = ['[sin(' Ws ')]'];
  h_str = ['J*' hi_str];
  HH_str = ['HH=fft(hc,Ng); HH((Ng/2+2):Ng)=[]; HH=HH.* ' ph_str ';'];
end;
flag = ~sym; 
A = Des .* exp(J*2*pi*grd*(N1+N2)/2);
clf;
plot(grd,abs(Des),'w--')
hold on
plot(grd, angle(Des),'w-')
hold off 
title('Desired Filter Frequency Response')
xlabel('normalized frequency') 
ylabel('magnitude (dashed) and phase (solid)')
disp('Press <return> to start ...')
pause
Ng = round(1.0/delf);
ifgrid = round((grd+0.5*flag)/delf) + 1; 
if (max(ifgrid) > Ng) 
  ifgrid(length(ifgrid)) = Ng;
end;
disp('Calculating initial solution .....')
[fext, iext] = crmz_guess(edges, grd, Lb);
it = 0; delta = 0.0; delta_old = -1;
%while ( abs(delta) > delta_old )
no_stp = 1;
while (no_stp)
  it = it + 1;
  delta_old = abs(delta);
  fext = grd(iext);
  W = 2*pi*fext*([0:(Lf-1)]+(~is_odd)*0.5);
  Mb = eval(M_str);
  M = [ Mb ((-1).^[0:Lb]')./WT(iext) ];
  a = M \ A(iext);
  delta = a(Lb+1);
  h = eval(h_str);
  hc = crmz_rotate(zeropad(h,Ng-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 .* (A - HH(ifgrid));
  EE(iext) = delta*( (-1) .^ [2:length(iext)+1] ); 
  [jext,EEj] = crmz_find(EE,iext);  
  if (prg_dspl < 2)
    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)));
     if (~prg_dspl)
      crmz_plerror;
      keyboard
    end;
  end;
  if (all(iext==jext'))
    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))
  not_optimal = 0; 
  fprintf('Optimal Solution Obtained !\n')
  if (prg_dspl == 2)
    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
  crmz_chckopt;
  not_optimal = 1; 
  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 Complex Remez Stage  %%%%%%%%%%%%%%%%%%%%%%%%
