function p = remezviz2(makemovie) % p = remezviz2(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 = -1:.2: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 % scale the w's so largest is at one. Arrangement is symmetric, % so smallest will fall to -1 scale = max(w); w = w*(1/scale); % move green lines too r = w(ix)'; p = shiftpoly(p,0,scale); d = remezplot([ii, 'Re-stretch to full range'],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 = -1.1:.001:1.1; aa = [-1.2 1.2 -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;