function [b,a,b1,b2] = maxflat2(L,M,N)
% MAXFLAT2: Generalized digital Butterworth filter (special values)
%
% [b,a,b1,b2] = maxflat2(L,M,N);
% input 
%   L : number of zeros at z=-1
%   M : number of zeros contributing to passband flatness
%   N : number of poles
%   need : L>N; M>=0; N>=0; N even
% output
%   b/a : IIR filter
%   b   : numerator coefficients
%   a   : denominator coefficients
%   b   : b = conv(b1,b2); b1 contains all zeros at z=-1,
%           b2 contains all other zeros.
%
% % Examples  
%	L = 10; M = 8; N = 4;
%	[b,a,b1,b2] = maxflat2(L,M,N);
%
%	L = 12; M = 12; N = 4;
%       [b,a,b1,b2] = maxflat2(L,M,N);

% Ivan W. Selesnick
% Rice University
% September 1995
% 
% see: Generalized Digital Butterworth Filter Design,
% by I. W. Selesnick and C. S. Burrus
%
%  required subprograms: choose.m

%  Check input
if rem(N,2)==1 
	disp('N must be even')
	b = []; a = []; b1 = []; b2 = [];
	break
end
if L<=N
	disp('L must be grater than N')
	b = []; a = []; b1 = []; b2 = [];
	break
end

k = M:-1:0;
S = choose(M+N-k,N).*choose(L-N+k-1,k);

k = N:-1:0;
Q = choose(M+N-k,M).*choose(M+L,k).*((-1).^k);

%%%%%%%%%%%%%%%%% Transform Roots %%%%%%%%%%%%%%%%%%%%%%

tmp = 1 - 2*roots(S);
br1 = tmp+sqrt(tmp.^2-1);
br2 = tmp-sqrt(tmp.^2-1);
br = sort([br1; br2]+eps*sqrt(-1));     % sort by absolute value
br = br(1:M);                           % take roots INSIDE unit circle
b2 = real(poly(br));

tmp = 1 - 2*roots(Q);
ar1 = tmp+sqrt(tmp.^2-1);
ar2 = tmp-sqrt(tmp.^2-1);
ar = sort([ar1; ar2]+eps*sqrt(-1));     % sort by absolute value
ar = ar(1:N);                           % take roots INSIDE unit circle
a = real(poly(ar));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b1 = 1;
for t = 1:L
   b1 = conv(b1,[1 1]); 
end
C  = sum(a)/(sum(b1)*sum(b2));
b2 = C*b2;
b  = conv(b1,b2);