function [a1,a2,p,q] = mfaps3(K,L,d,wc)
% 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] = mfaps3(K,L,d,wc)
% 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
% wc : cut-off frequency in [0,pi]
% a1, a2 : the denominators of the allpass filters A1(z), A2(z)
% p/q : overall transfer function
%
% % example
% K = 3; L = 5; d = 8; wc = 0.4554*pi;
% [a1,a2,p,q] = mfaps3(K,L,d,wc);

% Calls the function: nrst.m

% Ivan W. Selesnick
% Polytechnic University
% June, 1997

% 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

N = K+L+1;

[tmp,d1] = flatdelay(K,L+1,(d-N)/2);
[tmp,d0] = flatdelay(K+1,L,(d-N)/2);

D0c = polyval(d0(K+L+2:-1:1),exp(-i*wc));
D1c = polyval(d1(K+L+2:-1:1),exp(-i*wc));

r0 = abs(D0c); b0 = angle(D0c);
r1 = abs(D1c); b1 = angle(D1c);

bc = pi/3 - wc/2 * (N-d);
alpha1 = (r0*sin(b0-b1)-r0*cos(b0-b1)*tan(bc-b1)) / ...
        (r0*sin(b0-b1)+(r1-r0*cos(b0-b1))*tan(bc-b1));

bc = -pi/3 - wc/2 * (N-d);
alpha2 = (r0*sin(b0-b1)-r0*cos(b0-b1)*tan(bc-b1)) / ...
        (r0*sin(b0-b1)+(r1-r0*cos(b0-b1))*tan(bc-b1));

v = [alpha1 alpha2];
alpha = nrst(v, 1/2)

a = alpha*d1 + (1-alpha)*d0;
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