function [a1,a2,p,q] = mfaps2(K,L,d,alpha)
% Design of a maximally flat lowpass filter H(z) as the
% sum of two allpass filters: H(z) = z^(-d) A2(z) + A1(z).
%
% [a1,a2,p,q] = mfaps2(K,L,d,alpha)
% K : number of conditions at w=0
% L : number of conditions at w=pi
% d : degree of delay
% Note: two conditions must be satisfied
%   (1) abs(K-L) <= d <= K+L+2
%   (2) d must be same parity as K+L
% alpha : alpha in [0,1]
% a1, a2 : the denominators of the allpass filters A1(z), A2(z)
% p/q : overall transfer function
%
% % example
% K = 3; L = 5; d = 8; alpha = 0.5;
% [a1,a2,p,q] = mfaps2(K,L,d,alpha);

% Ivan W. Selesnick
% Rice University
% December, 1996

% check input for validity:
b1 = (abs(K-L) <= d) & (d <= K+L+2);
b2 = rem(K+L+2-d,2)==0;
if ~(b1 & b2)
   disp('For this K and L, d must be one of the following:');
   disp(abs(K-L):2:(K+L+2));
   break
end

[tmp,at1] = flatdelay(K,L+1,(d-K-L-1)/2);
[tmp,at2] = flatdelay(K+1,L,(d-K-L-1)/2);
a = alpha*at1 + (1-alpha)*at2;
rts = roots(a);
v = abs(rts)<1;
r1 = rts(v);                 	% roots inside unit circle
r2 = rts(~v);                	% roots outside unit circle
a1 = real(poly(r1));
a2 = real(poly(1./r2));

% compute overall transfer function p/q
p = [zeros(1,d), a] + [a(K+L+2:-1:1) zeros(1,d)];
q = conv(a1,a2);
p = p*sum(q)/sum(p);		% normalize