function [h,h1,h2] = fir_orthog(L,N,del)
%
%  [h,h1,h2] = fir_orthog(L,N,del)
%  Design of a minimum phase lowpass filter for a two 
%  channel orthogonal filter bank (wavelet system).
%  by : Ivan Selesnick, Rice University.
%  input
%    L   : length of filter (must be even)
%    N   : degree of flatness
%    del : ripple size (in magnitude squared)
%  need L/2-N even and L/2 >= N
%  output
%    h : minimum phase wavelet filter
%    h = conv(h1,h2); h1 contains the roots at z=-1,
%                     h2 contains all the other roots.
%  subprograms: 
%    oref.m, local_max.m, choose.m
%  % EXAMPLE
%	L = 14;
%	N = 5;
%	del = 0.01;
%	[h,h1,h2] = fir_orthog(L,N,del);

if rem(L,2) == 1
	disp('L must be even')
	break
end
if rem(L/2-N,2) == 1
	disp('L/2-N must be even, thanks!')
	break
end
if L/2 < N 
	disp('L/2 must be greater than or equal to N')
	break
end

SN = 1e-10;		% SN : SMALL NUMBER
g  = round(20*L);
m  = L/2-N;
Y  = del*(1/2+(-1).^[1:L/2-N]'/2);
x  = [0:g]'/g;
cs = 1/2 + [1:L/2-N]'/(2*(L/2-N+1));
pn = choose(N-1+[N-1:-1:0],N-1:-1:0)';

while 1	

	% use Daubechies parameterization:

	rhs = (Y./(1-cs).^N-polyval(pn,cs))./(cs.^N.*(1/2-cs));

	s = vander((1/2-cs).^2)\rhs;
	% can use a better polynomial interpolation algorithm here

	f = (1-x).^N.*(polyval(pn,x)+x.^N.*(1/2-x).*polyval(s,(1/2-x).^2));

	% plot stuff
	figure(1)	
	plot(x,f)
	hold on
	plot(cs,Y,'o')
	hold off

	ri = sort([local_max(f); local_max(-f)]);
	ri = ri(x(ri)>1/2);
	ri = ri(1:L/2-N);
	cs = x(ri);
	cs = oref(s,pn,N,cs);
	fcs = (1-cs).^N.*(polyval(pn,cs)+cs.^N.*(1/2-cs).*polyval(s,(1/2-cs).^2));
	Err = max([abs(fcs - Y);0]);
	fprintf(1,'    Err = %18.15f\n',Err);
	if Err < SN
		fprintf(1,'\n    I have converged \n'), break
	end

end

si = roots(s);
ti = [1/2+sqrt(si); 1/2-sqrt(si)];
if length(s) > 0
   t  = real(poly(ti)*s(1));
end
t  = t(:);
q = [conv([-1 1/2],t); pn];

% g = (1-x).^N .* polyval(q,x); % g is f.
% fig(2)
% plot(x,g)

%%%%%%%%%%%%%%%%%%%% convert to z-domain %%%%%%%%%%%%%%%%%%%%%%%%%%

hr1 = [exp(acos(1-2*cs(1:2:L/2-N))*sqrt(-1)); ...
	exp(-acos(1-2*cs(1:2:L/2-N))*sqrt(-1))];
qrts = roots(q);
hr2 = [1-2*qrts+sqrt((1-2*qrts).^2-1); ...
	1-2*qrts-sqrt((1-2*qrts).^2-1)];
hr2 = hr2(abs(hr2)<0.99);
h2 = real(poly([hr1;hr2]));

h1 = 1;
for t = 1:N
   h1 = conv(h1,[1 1]);
end

%%%%%%%%%%%%%%%%%%%% normalize %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h2 = h2/(sum(h1)*sum(h2));
h = conv(h1,h2);

% [mag,w] = h2mag(h);
% plot(w/pi,mag)