function [h,h2,rs] = chebbp(N,L,wp,us,ls)
% h = chebbp(N,L,wp,us,ls,g)
% A program for the design of symmetric bandpass FIR filters with 
% flat monotonically decreasing passbands (on either side of wp)
% and equiripple stopbands.
%
% output
%  h  : filter
% input
%  N  : length of total filter
%  L  : degree of flatness
%  wp : passband frequency of flatness
%  us : upper bound in stop band
%  ls : lower boudn in stop band
%
% Author: Ivan W. Selesnick, Rice University
%
% % EXAMPLE
%    N  =  25;
%    L  =  8;
%    wp  = .4*pi;
%    us =  0.01;
%    ls = -us;
%    [h,h2,rs] = chebbp(N,L,wp,us,ls);
% or
%    N = 31; L = 16; 
%    chebbp(N,L,wp,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



if (rem(N,2) == 0) | (rem(L,4) ~= 0)
   disp('N must be odd and L must be divisible by 4')
   return
else
   g = 2^ceil(log2(4*N));		% number of grid points
   SN = 1e-9;                           % SN : SMALL NUMBER
   q  = (N-L+1)/2;
   w  = [0:g]'*pi/g;
   ws1 = wp*0.8;
   a = 5; b = 1;
   ws2 = (a*wp+b*pi)/(a+b);

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

   Y1 = [ls*(1-(-1).^(q1:-1:1))/2 + us*((-1).^(q1:-1:1)+1)/2]';
   Y2 = [ls*(1-(-1).^(1:q2))/2 + us*((-1).^(1:q2)+1)/2]';
   Y = [Y1; Y2];

   n  = 0:q-1;
   Z   = zeros(2*(g-q)+1,1);
   % A1  = (-1)^(L/2) * (sin(w/2-wp/2).*sin(w/2+wp/2)).^(L/2);
   % A1r = (-1)^(L/2) * (sin(rs/2-wp/2).*sin(rs/2+wp/2)).^(L/2);
   A1  = (-1)^(L/2) * ((cos(wp)-cos(w))/2).^(L/2);
   A1r = (-1)^(L/2) * ((cos(wp)-cos(rs))/2).^(L/2);
   it = 0;
while 1 & (it < 20)
	if length(rs) ~= length(n)
		keyboard
	end
   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)
	pause(.5)
   ri = sort([localmax(A); localmax(-A)]);
   lri = length(ri);
   % ri(1:length(ri)-q) = [];
        if lri ~= 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) = [];
   q1 = sum(rs < wp);
   q2 = sum(rs > wp);

   Y1 = [ls*(1-(-1).^(q1:-1:1))/2 + us*((-1).^(q1:-1:1)+1)/2]';
   Y2 = [ls*(1-(-1).^(1:q2))/2 + us*((-1).^(1:q2)+1)/2]';
   Y = [Y1; Y2];

   % rs = refine2(a2,L/2,rs);
   % A1r = (-1)^(L/2) * (sin(rs/2-wp/2).*sin(rs/2+wp/2)).^(L/2);
   A1r = (-1)^(L/2) * ((cos(wp)-cos(rs))/2).^(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*cos(wp) 1]')/4;
end
h((N+1)/2) = h((N+1)/2) + 1;

end