home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / DATAFUN.DI$ / CPLXPAIR.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.2 KB  |  100 lines
  1. function z = cplxpair(x,tol)
  2. %CPLXPAIR Sort numbers into complex conjugate pairs.
  3. %    Y = CPLXPAIR(X) rearranges the elements of vector X so that
  4. %    complex numbers are collected into matched pairs of complex
  5. %    conjugates.  The pairs are ordered by increasing real part.
  6. %    Any purely real elements are placed after all the complex pairs.
  7. %    Y = CPLXPAIR(X,TOL) uses a relative tolerance of TOL for
  8. %    comparison purposes.  The default is TOL = 100*EPS.
  9.  
  10. %    L. Shure 1-27-88
  11. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  12.  
  13. if nargin < 2
  14.     tol = 100*eps;    % need to give a tolerance more generous than eps
  15. end
  16. if isempty(x)
  17.     z = x;
  18.     return;
  19. end
  20. r = [];
  21. j = sqrt(-1);
  22. [m,n] = size(x);
  23. z  = zeros(m,n); % return number in same shape vector as input
  24. x = x(:);     % make sure x is a column vector
  25. nz = max(size(z));
  26.  
  27. % first segregate reals
  28. ir = find(abs(imag(x)) <= tol*abs(x));
  29. nr = max(size(ir));
  30. if ~isempty(ir)
  31.     r  = sort(real(x(ir)));    % return sorted reals first
  32.     x(ir) = [];    % remove these values from input array
  33. end
  34. clear ir
  35.  
  36. nc = max(size(x));
  37. if nc == 0    % no complex numbers
  38.     z = r;
  39.     return
  40. end
  41. if rem(nc,2)
  42.     error('Complex numbers can''t be paired');
  43. end
  44.  
  45. % copy complex x to a work vector c
  46. c = x;
  47. np = nc/2;    % number of pairs supposedly
  48.  
  49. % get values to upper half plane
  50. c = real(c) + j*abs(imag(c));
  51.  
  52. % ordered by real part (since sort is very sensitive to magnitudes)
  53. [cc,cind] = sort(real(c));
  54. c = cc +j*imag(c(cind));
  55. % check to see if real parts are at least in pair
  56. if any(abs(cc(1:2:nc)-cc(2:2:nc)) > tol*abs(c(1:2:nc)))
  57.     error('Complex numbers can''t be paired')
  58. end
  59. x = x(cind);    % reorder x the same way c has been
  60. clear cc
  61.  
  62. % check real part pairs to see if imag parts are paired by conjugates
  63. % be careful with multiple roots!
  64. ip = 1;    % initialize pair counter
  65. while (ip <= np)
  66. % find indices for same real parts - but be careful because a real part can be 0
  67.     ii = find(abs(real(c(1:2:nc))-real(c(2*ip))) <= tol*abs(c(2*ip)));
  68.     if isempty(ii)
  69.         error('Complex numbers can''t be paired')
  70.     end
  71.     if max(size(ii)) > 1
  72. % multiple pairs with same real part - sort on imag(c(ii)) for all with real part
  73. % ij below are indices with same real part
  74.         ij = find(abs(real(c(1:nc)) - real(c(2*ip))) <= tol*abs(c(2*ip)));
  75.         nn = max(size(ij));
  76.         xtemp = x(ij);
  77.         [xi, xind] = sort(imag(xtemp));    % sort on imag parts - should be paired
  78.         xtemp = xtemp(xind);
  79. % check pairing
  80.         if any((abs(xtemp-conj(xtemp(nn:-1:1))) > tol*abs(xtemp)))
  81.             error('Complex numbers can''t be paired')
  82.         else
  83.             x(ij(1):2:ij(nn-1)) = xtemp(1:nn/2);
  84.             x(ij(2):2:ij(nn)) = conj(xtemp(1:nn/2));
  85.         end
  86.     else    % only one pair with that real part
  87.         if x(2*ip)-conj(x(2*ip-1)) > tol*abs(x(2*ip))
  88.             error('Complex numbers can''t be paired')
  89.         else
  90.             xtmp = real(x(2*ip))-sqrt(-1)*abs(imag(x(2*ip)));
  91.             x(2*ip-1:2*ip) = [xtmp conj(xtmp)];
  92.         end
  93.     end
  94.     ip = ip + max(size(ii));    % increment pair counter
  95. end
  96. % copy complex pairs into return vector
  97. z(1:nc) = x(1:nc);
  98. % append reals to this
  99. z(nc+1:nc+nr) = r;
  100.