function [h,h2,rs] = cheblp(N,L,us,ls)
% h = cheblp(N,L,us,ls,g)
% A program for the design of symmetric lowpass FIR filters with flat 
% monotonically decreasing passbands and equiripple stopbands.
%
% output
%  h  : filter
% input
%  N  : length of filter
%  L  : degree of flatness at w=0
%  us : upper bound in stop band
%  ls : lower bound in stop band
%
% Author: Ivan W. Selesnick, Rice University
%
% % EXAMPLE
%     N  =  15;
%     L  =  6;
%     us =  0.01;
%     ls = -us;
%     [h,h2,rs] = cheblp(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

% require subprograms : localmax.m, refine2.m

if (rem(N,2) == 0) | (rem(L,2) == 1)
   disp('N must be odd and L must be even')
   return
else
   g = 2^ceil(log2(5*N));               % number of grid points
   SN = 1e-9;                           % SN : SMALL NUMBER
   q  = (N-L+1)/2;
   Y  = [ls*(1-(-1).^(1:q))/2 + us*((-1).^(1:q)+1)/2]';
   w  = [0:g]'*pi/g;
   wo = pi*(L/2)/(L/2+q);
   wo = wo - (us-ls);
   n  = 0:q-1;
   rs = n'*(pi-wo)/(q-1) + wo;
   Z   = zeros(2*(g-q)+1,1);
   A1  = sin(w/2).^L * (-1)^(L/2);
   A1r = sin(rs/2).^L * (-1)^(L/2);
   it = 0;
while 1 & (it < 10)
   a2 = cos(rs*n)\[(Y-1)./A1r];
   A2 = real(fft([a2(1);a2(2:q)/2;Z;a2(q:-1:2)/2])); A2 = A2(1:g+1);
   A  = 1 + A1.*A2;
   figure(1), plot(w/pi,A),
        hold on, plot(rs/pi,Y,'o'), hold off, pause(.1)
   figure(2), plot(w/pi,20*log10(abs(A))),
        hold on, plot(rs/pi,20*log10(abs(Y)),'o'), hold off, pause(.1)
   ri = sort([localmax(A); localmax(-A)]);
   ri(1:length(ri)-q) = [];
   rs = (ri-1)*pi/g;
   rs = refine2(a2,L/2,rs);
   A1r = sin(rs/2).^L * (-1)^(L/2);
   Ar = 1 + (cos(rs*n)*a2) .* A1r;
   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

h2 = [a2(q:-1:2)/2; a2(1); a2(2:q)/2];
h = h2;
for k = 1:L/2
   h = conv(h,[1 -2 1]')/4;
end
h((N+1)/2) = h((N+1)/2) + 1;

end