function [h,h2,rs] = chebdiff(N,L,us,ls)
% h = chebdiff(N,L,us,ls,g)
% A program for the design of symmetric FIR lowpass
% differentiators with flat passbands and equiripple stopbands.
%
% output
%  h  : filter
% input 
%  us : upper bound in stop band
%  ls : lower bound in stop band
%
% % EXAMPLE
%    N  =  16;
%    L  =  6;
%    us =  0.01;
%    ls = -us;
%    [h,h2,rs] = chebdiff(N,L,us,ls);

% Reference:
% "Exchange Algorithms for the Design of Linear Phase FIR Filters 
% and Differentiators Having Flat Monotonic Passbands and Equiripple 
% Stopbands" by I. W. Selesnick and C. S. Burrus, 
% IEEE Trans. on Cicuits and Systems II, 43(9):671-675, Sept 1996


% required subprogram: localmax.m

if (rem(N,2) == 1) | (rem(L,2) == 1)
   disp('N and L must both be even')
   return
else

g = 2^ceil(log2(4*N));               % number of grid points

q  = (N-L)/2;			% number of filter parameters

% compute "d" coefficients of structure
d = 2;
for k = 1:L/2
	d(k+1) = prod(1:2:2*k-1)/(k*(2*k+1)*prod(1:(k-1))*(2^(2*k-1)));
end

SN  = 1e-9;                     % SN : SMALL NUMBER
Y   = [ls*(1-(-1).^(1:q))/2 + us*((-1).^(1:q)+1)/2]';
w   = [0:g]'*pi/g;		% w : frequency grid
wo  = pi*(L/2)/((L/2)+(q-1));
wo  = wo - (us-ls);		% initial estimate of stopband edge
m   = 0:q-1;
rs  = m'*(pi-wo)/(q-1) + wo;	% rs: reference set.
Z   = zeros(2*(g+2-q)-3,1);

C  = 1-cos(w);
Cr = 1-cos(rs);
S  = sin(w/2);
Sr = sin(rs/2);
it = 0;
while 1 & (it < 10)

   rhs = Y./Sr;
   for k = 1:L/2
      rhs = (rhs-d(k))./Cr;
   end
   a3 = cos(rs*m)\rhs;		% solve for coefficients
   A3 = real(fft([a3(1);a3(2:q)/2;Z;a3(q:-1:2)/2])); A3 = A3(1:g+1);
				% A3 : Amplitude
   A = A3;
   for k = L/2:-1:1
      A = A.*C+d(k);
   end
   A = A.*S;			% A : Amplitude of total filter

   figure(1), plot(w/pi,A),
        hold on, plot(rs/pi,Y,'o'), hold off, 

   % update the reference set
   ri = sort([localmax(A); localmax(-A)]);
   ri(1:length(ri)-q) = [];
   rs = (ri-1)*pi/g;
   % rs = refine3(a,L/2,rs);
   Cr = 1-cos(rs);
   Sr  = sin(rs/2);
  
   % evaluate Amplitude over the reference set rs 
   Ar = cos(rs*m)*a3;
   for k = L/2:-1:1
      Ar = Ar.*Cr+d(k); 
   end
   Ar = Ar.*Sr;

   % check stopping criterion
   Err = max([max(Ar)-us, ls-min(Ar)]);
   fprintf(1,'    Err = %18.15f\n',Err);
   if Err < SN, fprintf(1,'\n    I have converged \n'), break, end
   it = it + 1;
end

% construct symmetric filter from cosine coefficients
h3 = [a3(q:-1:2)/2; a3(1); a3(2:q)/2];
h = h3;
% use the parameterization 
for k = L/2:-1:1
   h = conv(h,[-1 2 -1]')/2;
   h((length(h)+1)/2) = h((length(h)+1)/2) + d(k);
end
h = conv(h,[1 -1]/2);

end