function remezviz(makemovie)
% remezviz(makemovie)
%   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

if nargin < 1
  makemovie = 0;
end

% Initialize the movie
if makemovie
  moviename = 'remezviz.mov';
  disp(['writing movie ', moviename]);
  eval(['MakeQTMovie start ', moviename]);
  eval(['MakeQTMovie framerate 10']);
end


rtn = [char(13)];
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,[],makemovie);
if ~makemovie; pause; end

% 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,makemovie);
  if ~makemovie; pause; end
  
  e = polyval(p,r);
  d = remezplot([ii, 'Actual extrema'],p,w,h,r,e,d,makemovie);
  if ~makemovie; pause; end

  % 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,makemovie);
  if ~makemovie; pause; end

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

end

if makemovie
  MakeQTMovie finish
end

% ffmpeg -y -i remezviz.mov -vcodec mpeg4 -b 200kb -mbd 2 -flags +mv4 -trellis 2 -cmp 2 -subcmp 2 remezviz.mp4



function d = remezplot(t,p,w,h,r,e,D,makemovie)
% 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 = [];
end
if nargin < 8
  makemovie = 0;
end

if length(D) == 0;
  d.p = p;
  d.w = w;
  d.h = h;
  d.r = r;
  d.e = e;
else
  d = D;
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;

  % Save this image
  if makemovie
    MakeQTMovie addfigure
  else
    drawnow
  end

end

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