function [h,h2,rs,del] = cheblp2(N,L,ws)
% h = cheblp2(N,L,us,ls)
% A second program for the design of symmetric lowpass FIR 
% filters with flat monotonically decreasing passbands and 
% equiripple stopbands. This program is similar to cheblp,
% but it uses a different set of input parameters.
%
% output
%  h  : filter
% input
%  N  : length of filter
%  L  : degree of flatness at w=0
%  ws : stopband edge (0-pi)
%
% Author: Ivan W. Selesnick, Rice University
%
% % EXAMPLE
%     N  =  25;
%     L  =  18;
%     ws =  0.6*pi;
%     [h,h2,rs,del] = cheblp2(N,L,ws);

% 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

%  calls subprograms: refine2.m, localmax.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(10*N));               	% number of grid points
   SN = 1e-10;                               	% SN : SMALL NUMBER
   q  = (N-L+1)/2;                              % number of filter parameters
   w = [0:g]'*pi/g;
   rs = [0:q]'*(pi-ws)/(q) + ws;
   Z = zeros(2*(g+1-q)-1,1);
   A1 = (-1)^(L/2) * (sin(w/2)).^L;
   si = (-1).^[0:q]';
   n = 0:q-1;
   A1r = (-1)^(L/2) * (sin(rs/2)).^L;
   it = 0;
while 1 & (it < 10)
   x = [cos(rs*n), si./A1r]\[1./A1r];
   a = -x(1:q);
   del = x(q+1);
   A2 = real(fft([a(1);a(2:q)/2;Z;a(q:-1:2)/2])); A2 = A2(1:g+1);
   A  = 1 + A1.*A2;
   Y  = si*del;
   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(a,L/2,rs);
   rs = [ws; rs];
   A1r = (-1)^(L/2) * (sin(rs/2)).^L;
   Ar = 1 + (cos(rs*n)*a).*A1r;
   Err = max([max(Ar)-del, -del-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 = [a(q:-1:2)/2; a(1); a(2:q)/2];
h = h2;
for k = 1:L/2
   h = conv(h,[1 -2 1]')/4;
	% fig(1), stem(h), pause
end
h((N+1)/2) = h((N+1)/2) + 1;

end