function [ r, p, k ] = residuez( b, a, t )
%RESIDUEZ Z-transform partial-fraction expansion.
%     [R,P,K] = RESIDUEZ(B,A)
%        finds the residues, poles and direct terms of a
%        partial-fraction expansion of B(z)/A(z)
%
%     B(z)      r(1)             r(n)
%     ---- = ---------- +...  ---------- + k(1) + k(2)z^-1 ...
%     A(z)   1-p(1)z^-1       1-p(n)z^-1
%
%          B: numerator polynomial coefficients
%          A: denominator coeffs (in ascending powers of z^-1)
%          R: the residues   (in a column vector)
%          P: the poles           (column vector)
%          K: the direct terms       (ROW vector)
%     [B,A] = RESIDUEZ(R,P,K)
%        convert partial-fraction expansion back to B/A form.
%     MULTIPLE POLES (order of residues):
%        residue for 1st power pole, then 2nd power, etc. 
%     see also RESIDUE, PRONY

%     J. McClellan 10-24-90
%    	Copyright (c) 1990 by the MathWorks, Inc.

%     REF: A.V. Oppenheim and R.W. Schafer, Discrete-Time
%          Signal Processing, Prentice-Hall, 1989, pp. 166-170.

mpoles_tol = 0.001;    %--- 0.1 percent
avg_roots = 0;         %--- FALSE, turns off averaging
if( nargin == 2 )
%================= convert B/A to partial fractions =========
  if( a(1) == 0 )
     error('first coefficient in a-vector must be non-zero')
  end
  b = b(:).'/a(1);   %--- b() & a() are now rows
  a = a(:).'/a(1);
  LB = length(b);  LA = length(a);
  if( LA == 1 )
     r = []; p = []; k = b;   %--- no poles!
     return
  end
  Nres = LA-1;
  p = zeros(Nres,1); r = p;
  LK = max([0;LB-Nres]);
  k = [];
  if( LK > 0 )
     [k,b] = deconv(fliplr(b),fliplr(a));
     k = fliplr(k);
     b(1:LK) = [];
     b = fliplr(b);
     LB = Nres;
  end

  p = roots(a);
  N = Nres;         %--- number of terms in r()
  xtra = 1;         %--- extra pts invoke least-squares approx
  imp = zeros(N+xtra,1);     %--- zero padding to make impulse
  imp(1) = 1;
  h = filter( b, a, imp );
  polmax = 1.0;max(abs(p));
  h = h.*filter(1,[1 -1/polmax],imp);
  [mults, idx] = mpoles( p, mpoles_tol );
  p = p(idx);
  S = zeros(N+xtra,N);
  for j = 1:Nres
     if( mults(j) > 1 )            %--- repeated pole
        S(:,j) = filter( 1, [1 -p(j)/polmax], S(:,j-1) );
     else
        S(:,j) = filter( 1, [1 -p(j)/polmax], imp );
     end
  end
  jdone = zeros(Nres,1);
  bflip = fliplr(b);
  for i=1:Nres
     ii = 0;
     if( i==Nres )
        ii = Nres;
     elseif( mults(i+1)==1 )
        ii = i;
     end
     if( ii > 0 )
        if( avg_roots )       %---  average repeated roots ??
          jkl = (i-mults(i)+1):i;
          p(jkl) = ones(mults(i),1)*mean( p(jkl) );
        end
        temp = fliplr(poly( p([1:(i-mults(i)), (i+1):Nres]) ));
        r(i) = polyval(bflip,1/p(i)) ./ polyval(temp,1/p(i));
        jdone(i) = 1;
     end
  end
%------ now do the repeated poles and the direct terms ------
  jjj = find(jdone==1);
  jkl = find(jdone~=1);
  if( any(jdone) )
     h = h - S(:,jjj)*r(jjj);
  end
  Nmultpoles = Nres - length(jjj);
  if( Nmultpoles > 0 )
     t = S(:,jkl)\h;
     r(find(jdone~=1)) = t(1:Nmultpoles);
  end
else
%============= partial fractions ---> B(z)/A(z) ============
     %----- NOW, B <--> R
     %-----      A <--> P
     %-----      T <--> K
  res = b(:);     %--- r is now a column
  pol = a(:);     %--- p is now a column
  LR = length(res);  LP = length(pol);  LK = length(t);
  if( LR ~= LP )
     error('length of r and p vectors must be the same')
  end
  if( LP == 0 )
     p = []; r = t(:).';     %--- no poles!
     return
  end
  N = LP+LK;            %--- number of terms in b() and a()

  [mults, idx] = mpoles( pol, mpoles_tol );
  pol = pol(idx);     %--- re-arrange poles & residues
  res = res(idx);
  p = poly(pol);  p = p(:);          %--- p is really A(z)

  if( LK > 0 )
     r = conv( p, t );   %--- A(z)K(z)
  else
     r = zeros(N,1);     %--- r is B(z), returned as COLUMN
  end
  for i=1:LP
     temp = poly( pol([1:(i-mults(i)), (i+1):LP]) );
     r = r + [res(i)*temp.'; zeros(N-length(temp),1)];
  end
end
