home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / MPOLES.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.3 KB  |  45 lines
  1. function [mults, indx] = mpoles( p, mpoles_tol )
  2. %MPOLES    Identify repeated poles & their multiplicities.
  3. %     [MULTS, IDX] = mpoles( P, TOL )
  4. %        P:     the list of poles
  5. %        TOL:   tolerance for checking when two poles
  6. %                 are the "same"
  7. %        MULTS: list of pole multiplicities
  8. %        IDX:   indices used to sort P
  9. %     NOTE: this is a support function for RESIDUEZ.
  10.  
  11. %   EXAMPLE:
  12. %     input:  P = [1 2 3 2 1 1]
  13. %     output: MULTS = [1 2 3 1 2 1], IDX = [1 5 6 2 4 3]
  14. %            P[IDX] = [1 1 1 2 2 3]
  15. %   Thus, MULTS contains the exponent for each pole in a
  16. %     partial fraction expansion.
  17.  
  18. if( nargin == 1 )
  19.    mpoles_tol = 1.0e-03;   %--- DEFAULT, if not passed by user
  20. end
  21. Lp = length(p);
  22. mults = zeros(Lp,1);
  23. done = zeros(Lp,1);   %--- set to FALSE
  24. indx = [1:Lp]';
  25. ii = 1;
  26. for j = 1:Lp
  27.    if( ~done(j) )
  28.       mults(ii) = 1;
  29.       indx(ii) = j;  ii = ii+1;
  30.       done(j) = 1;
  31.       if( ii > Lp ) break; end
  32.       test = abs( p((j+1):Lp) - p(j) );
  33.       jkl = find( test < mpoles_tol*abs(p(j)) );
  34.       Lj = length(jkl);
  35.       if( Lj > 0 )
  36.          for k = 1:Lj
  37.             jj = jkl(k);
  38.             mults(ii) = 1+k;
  39.             indx(ii) = j+jj;  ii = ii+1;
  40.             done(j+jj) = 1;
  41.          end
  42.       end
  43.    end
  44. end
  45.