resid2wp.m

matlab有理分式多项式展开系数计算程序

function [coeffs,poles,rp,k] = resid2wp(u,v,k)
%   This routine 'fixes' the numerical problems associated with RESIDUE.M
%   Note that the problem is still numerically sensitive, depending on the 
%   user's choice of tol.
%   
%   This routine has been tested on Matlab 5.0 for Windows
%   on polynomials up to order 9.
%
%       Bradley T. Burchett, Asst. Prof.              	Ph     [phone]
%       Dept. of Civil and Mechanical Engr.		FAX [phone]
%       U.S. Military Academy
%      West Point, NY 10996 
%     [e-mail]
%   
%   Usage is essentially the same as RESIDUE.M with additional output 
%   argument 'rp', which is the pole exponent, available.
%
%   RESID2WP Partial-fraction expansion (residues).
%   [R,P,rp,K] = RESID2WP(B,A) finds the residues, poles and direct term of
%   a partial fraction expansion of the ratio of two polynomials B(s)/A(s).
%   If there are no multiple roots,
%      B(s)       R(1)       R(2)             R(n)
%      ----  =  -------- + -------- + ... + -------- + K(s)
%      A(s)     s - P(1)   s - P(2)         s - P(n)
%   Vectors B and A specify the coefficients of the numerator and
%   denominator polynomials in descending powers of s.  The residues
%   are returned in the column vector R, the pole locations in column
%   vector P, and the direct terms in row vector K.  The number of
%   poles is n = length(A)-1 = length(R) = length(P). The direct term
%   coefficient vector is empty if length(B) < length(A), otherwise
%   length(K) = length(B)-length(A)+1.
%
%   If P(j) = ... = P(j+m-1) is a pole of multplicity m, then the
%   expansion includes terms of the form
%                R(j)        R(j+1)                         R(j+m-1)
%              -------- + ------------------   + ... + --------------------
%              s - P(j)   (s - P(j))^rp(j+1)           (s - P(j))^rp(j+m-1)
%
%   [B,A] = RESID2WP(R,P,K), with 3 input arguments and 2 output arguments,
%   converts the partial fraction expansion back to the polynomials with
%   coefficients in B and A.
%
%   See also POLY, ROOTS, DECONV.


tol = 0.000000000001;   % Repeated-root tolerance; adjust as needed.

% This section rebuilds the U/V polynomials from residues,
%  poles, and remainder.  Multiple roots are those that
%  differ by less than tol.

if nargin >2
   
  p = v;
  r = u;
   
  n = length(p);
  eh = ones(n,n)*diag(p);
  bee = eh.';

  nail = abs(eh-bee) < tol;     % Matrix matching poles.
  thumb = sum(nail);            % Multiplicity of pole.
  
  
  v = poly(p); temp2 = zeros(1,n);
  q = [p(:).' ; zeros(1,n)];   % Poles and multiplicities.
  for indx = 1:n

       if q(2,indx)
       else
          
         % Index of repeated pole under consideration:
         fi = find(nail(:,indx) == 1);     
         
         % Exponents of repeated poles:
         for m = 1:length(fi)
            q(2,fi(m)) = m;
         end % for m
         
      end  % if q
      temp = r(indx);
      
      % Common Denominator with pole repetitions:
      for in = 1:thumb(indx)-q(2,indx)
         temp = conv(temp,[1 -p(indx)]);
      end
      
      % Common Denominator with other poles:
      for in = 1:n
         if nail(in,indx)
         else
            temp = conv(temp,[1 -p(in)]);
         end
      end
      
      % Match vector lengths (add leading zeros) for addition:
      if length(temp2) < length(temp)
         temp2 = [zeros(1,(length(temp)-length(temp2))),temp2];
      elseif length(temp) < length(temp2)
         temp = [zeros(1,(length(temp2)-length(temp))),temp];
      end
      
      temp2 = temp2 + temp;
   end  % for indx
   coeffs = temp2;
   poles = v;
   return
end

% This section does the partial-fraction expansion
%  for any order of pole.

n = length(v) - 1;
u = u(:).'; v = v(:).'; k =[];
f = find(u ~= 0);
if length(f), u = u(f(1):length(u)); end
f = find(v ~= 0);
if length(f), v = v(f(1):length(v)); end
u = u ./ v(1); v = v ./ v(1);   % Normalize.
if length(u) >= length(v), [k,u] = deconv(u,v); end
p = roots(v); 
if isempty(p), coeffs = []; poles = []; k = 1; , return, end

eh = ones(n,n)*diag(p);
bee = eh.';

nail = abs(eh-bee) < tol;  % Pole matching matrix
thumb = sum(nail);         % Pole multiplicities

if any(thumb > 1)
   rp = zeros(n,1);
   
   for nb = 1:n
      lam = p;
      
      if thumb(nb) > 1
         
         if rp(nb)     
         else
            
            % Index of repeated pole under consideration:
            fi = find(nail(:,nb) == 1);
            
            % Average repeated poles:
            pole = mean(p(fi));
            p(fi) = ones(size(fi))*pole;       
            
            lam(fi) = [];
            
            % Determine the exponent related to each pole repetition:
            for m = 1:length(fi)
               
               rp(fi(m)) = m;  % Exponents
               vee = poly(lam);
               u = 1;

        for j = 1:thumb(nb)-m, [u,vee] = polyder(u,vee); end
        c = 1; if m < thumb(nb), c = prod(1:thumb(n)-m); end
        nk(fi(m)) = (polyval(u,pole) ./ polyval(vee,pole)) ./ c;

      end %m

    end % if rp
  end % if thumb
end % for nb = 1:n
coeffs = nk;
poles = p;
   
else   % No repeated roots.
   for i = 1:n
      temp = poly(p([1:i-1, i+1:n]));
      coeffs(i) = polyval(u,p(i)) ./ polyval(temp,p(i));
   end
   poles = p;
   rp = ones(n,1);
end