%%%%%%%%%% crmz_find.m %%%%%%%%%%
function  [fnew,enew] = crmz_find( error, fold )
%  Construct new subset of points using the exchange
%  rules derived in Chapter 3 of Thesis. (see also 
%  "Complex Chebyshev approximation for FIR filter design" 
%  by Lina J. Karam and James H. McClellan in the IEEE 
%  Trans. on Circuits and Systems II - March 1995. 
%
fold = fold(:);
Nx = length( fold );
Ngrid = length( error );
sgn_error = error(fold(1))/abs(error(fold(1)));
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';
%%%%%%%%%% End crmz_find.m %%%%%%%%%%
