function remezviz()
% remezviz()
%   Visualize development of an equiripple polynomial by something 
%   that behaves a little like the remez exchange algorithm
%   (but without controlling epsilon in different ranges)
% 2006-11-21 Dan Ellis dpwe@ee.columbia.edu

rtn = [char(13),char(10)];
disp(['This animation attempts to give some insight into the Remez ',rtn, ...
      'exchange algorithm, which is able to find the minimax ',rtn, ...
      '(equiripple) polynomial fitting certain initial constraints. ',rtn, ...
      'The essence is that we find a polynomial that passes ',rtn, ...
      'through a set of frequencies where the ripple is equal, ',rtn, ...
      'but then at iteration we shift those frequencies to ',rtn, ...
      'match the extrema of the resulting polynomial, which ',rtn, ...
      'in general will correspond to errors outside our ',rtn, ...
      'equiripple bound.  The target points rapidly converge ',rtn, ...
      'to the actual extrema of the polynomial, giving the ',rtn, ...
      'equiripple solution.',rtn,rtn,'Press any key to continue']);


% w is the set of candidate extremal freqs - intialize uniformly
w = 0:.1:1;

% intialize record of last iteration's frequencies
wm1 = 0*w;

% h is the desired magnitude response at the frequencies.
% Set it to alternate between 0 and 1 (equiripple)
h = ones(1,length(w));
h(2:2:end) = 0;

N = length(h);

zz = zeros(N-2,1);
d = remezplot('Initial target points', zeros(1,N),w,h,zz,zz);
pause;

% x and y locations of actual extrema from last iteration
rm1 = zz;
em1 = zz;

for i = 1:7

  ii = (['Iteration ', num2str(i),': ']);
  
  [p,r] = lagrangepoly(w,h);

  d = remezplot([ii, 'Lagrange fit'],p,w,h,rm1,em1,d);
  pause;

  e = polyval(p,r);
  d = remezplot([ii, 'Actual extrema'],p,w,h,r,e,d);
  
  pause;

  % replace control points with nearest extrema
  deltas = repmat(w', 1, length(r)) - repmat(r', length(w), 1);
  [vv,ix] = min(abs(deltas));
  wm1 = w;
  w(ix) = r';

  d = remezplot([ii, 'Shift target points to extremal freqs'],p,w,h,r,e,d);
  pause;

  rm1 = r;
  em1 = e;
  
  disp([ii,'Mean square frequency shift: ', num2str(norm(w - wm1)/N)]);

end


function d = remezplot(t,p,w,h,r,e,d)
% Subfunction to plot intermediate stages of remez
% t is the title string
% p are the polynomial coefficients
% w, h are the x and y co-ordinates of the (equiripple) target points
% r, e are the x and y co-ordinates of the true extrema of the poly
% d is a state structure, pass in the return from the last call to 
%   get continuous animation

if nargin < 7
  d.p = p;
  d.w = w;
  d.h = h;
  d.r = r;
  d.e = e;
end

steps = 10;

xx = -.05:.001:1.05;
aa =  [-.1 1.1 -3 4];

for i = 1:steps;
  
  u = i/steps; v = 1-u;
  
  pp = u*p + v*d.p;
  ww = u*w + v*d.w;
  hh = u*h + v*d.h;
  rr = u*r + v*d.r;
  ee = u*e + v*d.e;

  plot([rr,rr]',repmat([aa(3) aa(4)],length(r),1)','--g'); 
  hold on; 
  plot(xx, polyval(pp,xx),'-b',rr,ee,'ob',ww,hh,'or');
  hold off
  title(t);
  axis(aa); grid; drawnow;

end

d.p = p;
d.w = w;
d.h = h;
d.r = r;
d.e = e;