home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / POLYFUN.DI$ / RESIDUE.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  4.5 KB  |  136 lines
  1. function [coeffs,poles,k] = residue(u,v,k)
  2. %RESIDUE Partial-fraction expansion or residue computation.
  3. %    [R,P,K] = RESIDUE(B,A) finds the residues, poles and direct term of
  4. %    a partial fraction expansion of the ratio of two polynomials,
  5. %    B(s) and A(s).  If there are no multiple roots,
  6. %       B(s)       R(1)       R(2)             R(n)
  7. %       ----  =  -------- + -------- + ... + -------- + K(s)
  8. %       A(s)     s - P(1)   s - P(2)         s - P(n)
  9. %    Vectors B and A specify the coefficients of the polynomials in
  10. %    descending powers of s.  The residues are returned in the column
  11. %    vector R, the pole locations in column vector P, and the direct
  12. %    terms in row vector K.  The number of poles is
  13. %       n = length(A)-1 = length(R) = length(P)
  14. %    The direct term coefficient vector is empty if length(B) < length(A);
  15. %    otherwise
  16. %       length(K) = length(B)-length(A)+1
  17. %
  18. %    If P(j) = ... = P(j+m-1) is a pole of multplicity m, then the
  19. %    expansion includes terms of the form
  20. %                 R(j)        R(j+1)                R(j+m-1)
  21. %               -------- + ------------   + ... + ------------
  22. %               s - P(j)   (s - P(j))^2           (s - P(j))^m
  23. %    [B,A] = RESIDUE(R,P,K), with 3 input arguments and 2 output arguments,
  24. %    converts the partial fraction expansion back to the polynomials with
  25. %    coefficients in B and A.
  26. %
  27. %    Warning:
  28. %    Numerically, the partial fraction expansion of a ratio of polynomials
  29. %    represents an ill-posed problem.  If the denominator polynomial, A(s),
  30. %    is near a polynomial with multiple roots, then small changes in the
  31. %    data, including roundoff errors, can make arbitrarily large changes
  32. %    in the resulting poles and residues.  Problem formulations making use
  33. %    of state-space or zero-pole representations are preferable.
  34.  
  35. %    Reference:  A.V. Oppenheim and R.W. Schafer, Digital
  36. %    Signal Processing, Prentice-Hall, 1975, p. 56-58.
  37.  
  38. %    C.R. Denham and J.N. Little, MathWorks, 1986, 1989.
  39. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  40.  
  41. tol = 0.001;   % Repeated-root tolerance; adjust as needed.
  42.  
  43. % This section rebuilds the U/V polynomials from residues,
  44. %  poles, and remainder.  Multiple roots are those that
  45. %  differ by less than tol.
  46.  
  47. if nargin > 2
  48.    [p,i] = sort(-abs(v));   % Sort poles by magnitude.
  49.    p = v(i); r = u(i);   % Rearrange residues similarly.
  50.    n = length(p);
  51.    q = [p(:).' ; ones(1,n)];   % Poles and multiplicities.
  52.    for j = 2:n
  53.       if abs(p(j) - p(j-1)) < tol
  54.          q(2,j) = q(2,j-1) + 1;   % Multiplicity of pole.
  55.       end
  56.    end
  57.    v = poly(p); u = zeros(1,n);
  58.    for indx = 1:n
  59.       ptemp = q(1,:);
  60.       i = indx;
  61.       for j = 1:q(2,indx), ptemp(i) = nan; i = i-1; end
  62.       ptemp = ptemp(find(~isnan(ptemp))); temp = poly(ptemp);
  63.       j = length(temp);
  64.       if j < n, temp = [zeros(1,n-j) temp]; end
  65.       u = u + (r(indx) .* temp);
  66.    end
  67.    if ~isempty(k)
  68.       if any(k ~= 0)
  69.          u = [zeros(1,length(k)) u];
  70.          k = k(:).';
  71.          temp = conv(k,v);
  72.          u = u + temp;
  73.       end
  74.    end
  75.    coeffs = u; poles = v;    % Rename.
  76.    return
  77. end
  78.  
  79. % This section does the partial-fraction expansion
  80. %  for any order of pole.
  81.  
  82. u = u(:).'; v = v(:).'; k =[];
  83. f = find(u ~= 0);
  84. if length(f), u = u(f(1):length(u)); end
  85. f = find(v ~= 0);
  86. if length(f), v = v(f(1):length(v)); end
  87. u = u ./ v(1); v = v ./ v(1);   % Normalize.
  88. if length(u) >= length(v), [k,u] = deconv(u,v); end
  89. r = roots(v); [p,i] = sort(-abs(r)); p = r(i);
  90.  
  91. q = zeros(2,length(p));   % For poles and multiplicity.
  92. q(1,1) = p(1); q(2,1) = 1; j = 1;
  93. repeated = 0;
  94. for i = 2:length(p)
  95.    av = q(1,j) ./ q(2,j);
  96.    if abs(av - p(i)) <= tol    % Treat as repeated root.
  97.       q(1,j) = q(1,j) + p(i);   % Sum for average value.
  98.       q(2,j) = q(2,j) + 1;
  99.       repeated = 1;
  100.      else
  101.       j = j + 1; q(1,j) = p(i); q(2,j) = 1;
  102.    end
  103. end
  104. q(1,1:j) = q(1,1:j) ./ q(2,1:j);   % Multiple root average.
  105.  
  106. % Set desired = 1 if you want the output multiple poles
  107. %  to be averaged.
  108.  
  109. desired = 1;
  110. if repeated & desired
  111.    indx = 0;
  112.    for i = 1:j
  113.       for ii = 1:q(2,i), indx = indx+1; p(indx) = q(1,i); end
  114.    end
  115. end
  116.  
  117. poles = p(:);  % Rename.
  118.  
  119. coeffs = zeros(length(p),1); 
  120. if repeated   % Section for repeated root problem.
  121.    v = poly(p);
  122.    next = 0;
  123.    for i = 1:j
  124.       pole = q(1,i); n = q(2,i);
  125.       for indx = 1:n
  126.          next = next + 1;
  127.          coeffs(next) = resi2(u,v,pole,n,indx);
  128.       end
  129.    end
  130.   else   % No repeated roots.
  131.    for i = 1:j
  132.       temp = poly(p([1:i-1, i+1:j]));
  133.       coeffs(i) = polyval(u,p(i)) ./ polyval(temp,p(i));
  134.    end
  135. end
  136.