function [h,h2,rs,del] = chebbp2(N,L,wp,ws1,ws2)
% h = chebbp2(N,L,wp,ws1,ws2)
% A second program for the design of symmetric bandpass FIR filters with 
% flat monotonically decreasing passbands (on either side of wp)
% and equiripple stopbands. This program is similar to chebbp.m,
% but it uses a different set of input parameters.
%
% output
%  h  : filter
% input
%  N   : length of total filter
%  L   : degree of flatness
%  wp  : passband frequency of flatness 
%  ws1 : first stopband edge
%  ws2 : second stopband edge
%  Need: 0 < ws1 < wp < ws2 < pi
%
% Author: Ivan W. Selesnick, Rice University
%
% % EXAMPLE
%    N  =  55;
%    L  =  8;
%    wp  = 0.4*pi;
%    ws1 = 0.3*pi;
%    ws2 = 0.6*pi;
%    [h,h2,rs,del] = chebbp2(N,L,wp,ws1,ws2);

% 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


if (rem(N,2) == 0) | (rem(L,4) ~= 0)
   disp('N must be odd and L must be divisible by 4')
   return
elseif (0 >= ws1) | (ws1 >= wp) | (wp >= ws2) | (ws2 >= pi)
   disp('need : 0 < ws1 < wp < ws2 < pi')
elseif (L<1)
   disp('L must be positive')
   return
else

   q  = (N-L+1)/2; 	% number of filter parameters
			% num. of ref. set frequencies = q + 1

   g = 2^ceil(log2(8*N));       % number of grid points
   SN = 1e-8;                   % SN : SMALL NUMBER
   w = [0:g]'*pi/g;

   d = ws1/(pi-ws2);		% q1 : number of ref. set freq. in 1st stopband
   q1 = round((q+1)/(1+1/d));	% q2 : number of ref. set freq. in 2nd stopband
   if q1 == 0
      q1 = 1;
   elseif q1 == q+1
      q1 = q;
   end
   q2 = q + 1 - q1;	

   if q1 == 1;
	rs1 = ws1;
   else
	rs1 = [0:q1-1]'*(ws1/(q1-1));
   end
   if q2 == 1
	rs2 = ws2;
   else
	rs2 = [0:q2-1]'*(pi-ws2)/(q2-1)+ws2;
   end
   rs = [rs1; rs2];


   Z = zeros(2*(g+1-q)-1,1);
   A1 = (-1)^(L/2) * (sin(w/2-wp/2).*sin(w/2+wp/2)).^(L/2);
   si = (-1).^[0:q]';
   n = 0:q-1;
   A1r = (-1)^(L/2) * (sin(rs/2-wp/2).*sin(rs/2+wp/2)).^(L/2);
   it = 0;
while 1 & (it < 15)
   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, 
   figure(2), plot(w/pi,20*log10(abs(A))),
        hold on, plot(rs/pi,20*log10(abs(Y)),'o'), hold off, 
	pause(.5)
   ri = sort([localmax(A); localmax(-A)]);
   lri = length(ri);
	if length(ri) ~= q+1 
		if abs(A(ri(lri))-A(ri(lri-1))) < abs(A(ri(1))-A(ri(2)))
			ri(lri) = [];
		else
			ri(1) = [];
		end
	end
   rs = (ri-1)*pi/g;
   [temp, k] = min(abs(rs - wp)); rs(k) = [];
   Aws1 = 1 + (cos(ws1*n)*a).*((-1)^(L/2) * (sin(ws1/2-wp/2).*sin(ws1/2+wp/2)).^(L/2));
   Aws2 = 1 + (cos(ws2*n)*a).*((-1)^(L/2) * (sin(ws2/2-wp/2).*sin(ws2/2+wp/2)).^(L/2));
   if (Aws1 > Aws2) | any((rs > wp) & (rs < ws2))
      rs = sort([ws1; rs]);
   else
      rs = sort([ws2; rs]);
   end
   A1r = (-1)^(L/2) * (sin(rs/2-wp/2).*sin(rs/2+wp/2)).^(L/2);
   Ar = 1 + (cos(rs*n)*a).*A1r;
   Err = max([max(Ar)-abs(del), abs(del)-abs(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*cos(wp) 1]')/4;
	% fig(1), stem(h), pause
end
h((N+1)/2) = h((N+1)/2) + 1;

end