function [b,a,b1,b2] = maxflat2(L,M,N)
%  [b,a,b1,b2] = maxflat2(L,M,N);
%  Generalized Digital Butterworth Filter.
%  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   : length L+M+1 vector of polynomial coefficients
%    a   : length N+1 vector of polynomial 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);


