home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / RESIDUEZ.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.1 KB  |  134 lines
  1. function [ r, p, k ] = residuez( b, a, t )
  2. %RESIDUEZ Z-transform partial-fraction expansion.
  3. %     [R,P,K] = RESIDUEZ(B,A)
  4. %        finds the residues, poles and direct terms of a
  5. %        partial-fraction expansion of B(z)/A(z)
  6. %
  7. %     B(z)      r(1)             r(n)
  8. %     ---- = ---------- +...  ---------- + k(1) + k(2)z^-1 ...
  9. %     A(z)   1-p(1)z^-1       1-p(n)z^-1
  10. %
  11. %          B: numerator polynomial coefficients
  12. %          A: denominator coeffs (in ascending powers of z^-1)
  13. %          R: the residues   (in a column vector)
  14. %          P: the poles           (column vector)
  15. %          K: the direct terms       (ROW vector)
  16. %     [B,A] = RESIDUEZ(R,P,K)
  17. %        convert partial-fraction expansion back to B/A form.
  18. %     MULTIPLE POLES (order of residues):
  19. %        residue for 1st power pole, then 2nd power, etc. 
  20. %     see also RESIDUE, PRONY
  21.  
  22. %     J. McClellan 10-24-90
  23. %        Copyright (c) 1990 by the MathWorks, Inc.
  24.  
  25. %     REF: A.V. Oppenheim and R.W. Schafer, Discrete-Time
  26. %          Signal Processing, Prentice-Hall, 1989, pp. 166-170.
  27.  
  28. mpoles_tol = 0.001;    %--- 0.1 percent
  29. avg_roots = 0;         %--- FALSE, turns off averaging
  30. if( nargin == 2 )
  31. %================= convert B/A to partial fractions =========
  32.   if( a(1) == 0 )
  33.      error('first coefficient in a-vector must be non-zero')
  34.   end
  35.   b = b(:).'/a(1);   %--- b() & a() are now rows
  36.   a = a(:).'/a(1);
  37.   LB = length(b);  LA = length(a);
  38.   if( LA == 1 )
  39.      r = []; p = []; k = b;   %--- no poles!
  40.      return
  41.   end
  42.   Nres = LA-1;
  43.   p = zeros(Nres,1); r = p;
  44.   LK = max([0;LB-Nres]);
  45.   k = [];
  46.   if( LK > 0 )
  47.      [k,b] = deconv(fliplr(b),fliplr(a));
  48.      k = fliplr(k);
  49.      b(1:LK) = [];
  50.      b = fliplr(b);
  51.      LB = Nres;
  52.   end
  53.  
  54.   p = roots(a);
  55.   N = Nres;         %--- number of terms in r()
  56.   xtra = 1;         %--- extra pts invoke least-squares approx
  57.   imp = zeros(N+xtra,1);     %--- zero padding to make impulse
  58.   imp(1) = 1;
  59.   h = filter( b, a, imp );
  60.   polmax = 1.0;max(abs(p));
  61.   h = h.*filter(1,[1 -1/polmax],imp);
  62.   [mults, idx] = mpoles( p, mpoles_tol );
  63.   p = p(idx);
  64.   S = zeros(N+xtra,N);
  65.   for j = 1:Nres
  66.      if( mults(j) > 1 )            %--- repeated pole
  67.         S(:,j) = filter( 1, [1 -p(j)/polmax], S(:,j-1) );
  68.      else
  69.         S(:,j) = filter( 1, [1 -p(j)/polmax], imp );
  70.      end
  71.   end
  72.   jdone = zeros(Nres,1);
  73.   bflip = fliplr(b);
  74.   for i=1:Nres
  75.      ii = 0;
  76.      if( i==Nres )
  77.         ii = Nres;
  78.      elseif( mults(i+1)==1 )
  79.         ii = i;
  80.      end
  81.      if( ii > 0 )
  82.         if( avg_roots )       %---  average repeated roots ??
  83.           jkl = (i-mults(i)+1):i;
  84.           p(jkl) = ones(mults(i),1)*mean( p(jkl) );
  85.         end
  86.         temp = fliplr(poly( p([1:(i-mults(i)), (i+1):Nres]) ));
  87.         r(i) = polyval(bflip,1/p(i)) ./ polyval(temp,1/p(i));
  88.         jdone(i) = 1;
  89.      end
  90.   end
  91. %------ now do the repeated poles and the direct terms ------
  92.   jjj = find(jdone==1);
  93.   jkl = find(jdone~=1);
  94.   if( any(jdone) )
  95.      h = h - S(:,jjj)*r(jjj);
  96.   end
  97.   Nmultpoles = Nres - length(jjj);
  98.   if( Nmultpoles > 0 )
  99.      t = S(:,jkl)\h;
  100.      r(find(jdone~=1)) = t(1:Nmultpoles);
  101.   end
  102. else
  103. %============= partial fractions ---> B(z)/A(z) ============
  104.      %----- NOW, B <--> R
  105.      %-----      A <--> P
  106.      %-----      T <--> K
  107.   res = b(:);     %--- r is now a column
  108.   pol = a(:);     %--- p is now a column
  109.   LR = length(res);  LP = length(pol);  LK = length(t);
  110.   if( LR ~= LP )
  111.      error('length of r and p vectors must be the same')
  112.   end
  113.   if( LP == 0 )
  114.      p = []; r = t(:).';     %--- no poles!
  115.      return
  116.   end
  117.   N = LP+LK;            %--- number of terms in b() and a()
  118.  
  119.   [mults, idx] = mpoles( pol, mpoles_tol );
  120.   pol = pol(idx);     %--- re-arrange poles & residues
  121.   res = res(idx);
  122.   p = poly(pol);  p = p(:);          %--- p is really A(z)
  123.  
  124.   if( LK > 0 )
  125.      r = conv( p, t );   %--- A(z)K(z)
  126.   else
  127.      r = zeros(N,1);     %--- r is B(z), returned as COLUMN
  128.   end
  129.   for i=1:LP
  130.      temp = poly( pol([1:(i-mults(i)), (i+1):LP]) );
  131.      r = r + [res(i)*temp.'; zeros(N-length(temp),1)];
  132.   end
  133. end
  134.