%RESIDUED computes the partial fraction form from transfer-function form B/A.
%   [R,P,F] = residued(B,A) converts the transfer function B(z)/A(z) to
%   the partial fraction form:
%
%       B(z)              N       R(n)           NF
%       ---- = z^(-NF) * SUM ---------------  + SUM F(i)*z^(-i+1)
%       A(z)             n=1 (1-P(n)/z)^M(n)    i=1
%
%   where NF = nb-na+1 (length of FIR part), N=na-1 (number of poles),
%   nb=length(B), na = length(A).  When nb<na, F is empty.
%
%   RETURNED:
%    R = column vector of filter-pole residues
%    P = column vector of filter poles
%    F = row vector containing the FIR part, or an empty vector
%
%   The difference between RESIDUED and RESIDUEZ is the z^(-NF) term
%   above when there an FIR part F exists.  That is, RESIDUED delays
%   the sum of one-pole terms so it doesn't overlap with the FIR part. See
%   the companion paper for details:
%
%   Balazs Bank and Julius O. Smith, III, "A delayed parallel filter structure
%   with an FIR part having improved numerical properties", 136th AES
%   Convention, Preprint No. 9084, Berlin, April 2014.
%   Home page: http://www.mit.bme.hu/~bank/delparf
%   Initial version: http://ccrma.stanford.edu/~jos/filters/residued.html
%
%   SEE ALSO: filter, residue, residuez, tf2delparf, rpk2parf, delparfilt
%
%   C. Julius Smith, 2014.

function [r, p, f] = residued(b, a)
NUM = b(:)';
DEN = a(:)';
nb = length(NUM);
na = length(DEN);
f = [];
if na<=nb
  f = filter(NUM,DEN,[1 zeros(1,nb-na)]);
  NUM = NUM - conv(DEN,f);
  NUM = NUM(nb-na+2:end);
end;
[r,p,f2] = residuez(NUM,DEN);