function [d1,d2,beta]=tf2ca(b,a)
%TF2CA Transfer function to coupled allpass conversion.

% References:[1] P.P. Vaidyanathan, ROBUST DIGITAL FILTER STRUCTURES, in
%               HANDBOOK FOR DIGITAL SIGNAL PROCESSING. S.K. Mitra and J.F.            
%               Kaiser Eds. Wiley-Interscience, N.Y., 1993, Chapter 7.
%            [2] P.P. Vaidyanathan, MULTIRATE SYSTEMS AND FILTER BANKS, 
%               Prentice Hall, N.Y., Englewood Cliffs, NJ, 1993, Chapter 3.

%   Copyright 1999-2017 The MathWorks, Inc.

narginchk(2,2);

[b,a] = parseinput(b,a);

% Compute the power complementary filter
try
    [bp,a] = iirpowcomp(b,a);
catch ME
    throw(ME);
end

% Sort numerators to call allpassdecomposition
[p,q] = sortnums(b,bp);

% Once you have sorted the numerators, compute the actual decomposition
[d1,d2,beta] = allpassdecomposition(p,q,a);

%------------------------------------------------------------------------
function [b,a] = parseinput(b,a)
%PARSEINPUT Make sure input args are real vectors. 
%   
%   Force them to be row vectors

% Check that b and a are vectors of the same length
if ~(any(size(b)==1) && any(size(a)==1))
   error(message('dsp:tf2ca:FilterErr1'));
end
if length(b)~=length(a)
   error(message('dsp:tf2ca:FilterErr2'));
end

% Make sure b and a are rows
b = b(:).';
a = a(:).';

if ~(isreal(b) && isreal(a))
    error(message('dsp:tf2ca:FilterErr3'));
end

% Make sure numerator is symmetric or antisymmetric
symstr = signalpolyutils('symmetrytest',b);

if strcmpi(symstr,'none')
    error(message('dsp:tf2ca:FilterErr4'));
end

%------------------------------------------------------------------------
function [p,q] = sortnums(b,bp)
%SORTNUMS Sort numerators prior to calling ALLPASSDECOMPOSITION.
%   ALLPASSDECOMPOSITION always requires the first argument to
%   be symmetric.  The second argument, can be
%   symmetric, or antisymmetric.
%
%   IIRPOWCOMP can be called in some cases with antisymmetric numerator
%   and return a symmetric power complementary numerator. In this
%   case, we must sort these arguments before proceeding

% If b is real, antisymmetric, make it the second arg
if max(abs(b + b(end:-1:1))) < eps^(2/3)
    p = bp;
    q = b;    
else
    p = b;
    q = bp;
end

%------------------------------------------------------------------------
function [d1,d2,beta] = allpassdecomposition(p,q,a)
%ALLPASSDECOMPOSITION  Compute the allpass decomposition.
%   Given an IIR filter P/A and its power complementary filter Q/A
%   find the allpass decomposition for the filter:
%
%           P(z)   1                                  (generalized
%   H(z) = ----- = - [conj(BETA)*H1(z) + BETA*H2(z)]  coupled allpass 
%           A(z)   2                                  decomposition).
%
%   NOTE: In this function, if P is real it must always be symmetric.
%         Make sure you sort the args correctly prior to calling this 
%         function.


% If q is real, antisymmetric, make it imaginary
if isreal(q) && (max(abs(q + q(end:-1:1))) < eps^(2/3))
    q = q*1i;
end

z = roots(p-1i*q);
   
% Initialize the allpass functions
d1 = 1;
d2 = 1;

% Separate the zeros inside the unit circle and the ones outside to form the allpass functions
for n=1:length(z)
   if abs(z(n)) < 1
      d2 = conv(d2,[1 -z(n)]);
   else
      d1 = conv(d1,[1 -1/conj(z(n))]);
   end
end

% Remove roundoff imaginary parts
d1 = signalpolyutils('imagprune',d1);
d2 = signalpolyutils('imagprune',d2);

beta = sum(d2)*(sum(p)+1i*sum(q))./sum(a)./sum(conj(d2));

% [EOF] - tf2ca.m

