X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fm%2Fgeneral%2Fcplxpair.m;fp=octave_packages%2Fm%2Fgeneral%2Fcplxpair.m;h=dff9aecaee8d64e0943157c81f7c923c4c3e9547;hb=1c0469ada9531828709108a4882a751d2816994a;hp=0000000000000000000000000000000000000000;hpb=63de9f36673d49121015e3695f2c336ea92bc278;p=CreaPhase.git

diff --git a/octave_packages/m/general/cplxpair.m b/octave_packages/m/general/cplxpair.m
new file mode 100644
index 0000000..dff9aec
--- /dev/null
+++ b/octave_packages/m/general/cplxpair.m
@@ -0,0 +1,164 @@
+## Copyright (C) 2000-2012 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {} cplxpair (@var{z})
+## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol})
+## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol}, @var{dim})
+## Sort the numbers @var{z} into complex conjugate pairs ordered by
+## increasing real part.  Place the negative imaginary complex number
+## first within each pair.  Place all the real numbers (those with
+## @code{abs (imag (@var{z}) / @var{z}) < @var{tol})}) after the
+## complex pairs.
+##
+## If @var{tol} is unspecified the default value is 100*@code{eps}.
+##
+## By default the complex pairs are sorted along the first non-singleton
+## dimension of @var{z}.  If @var{dim} is specified, then the complex
+## pairs are sorted along this dimension.
+##
+## Signal an error if some complex numbers could not be paired.  Signal an
+## error if all complex numbers are not exact conjugates (to within
+## @var{tol}).  Note that there is no defined order for pairs with identical
+## real parts but differing imaginary parts.
+## @c Set example in small font to prevent overfull line
+##
+## @smallexample
+## cplxpair (exp(2i*pi*[0:4]'/5)) == exp(2i*pi*[3; 2; 4; 1; 0]/5)
+## @end smallexample
+## @end deftypefn
+
+## FIXME: subsort returned pairs by imaginary magnitude
+## FIXME: Why doesn't exp(2i*pi*[0:4]'/5) produce exact conjugates. Does
+## FIXME:    it in Matlab?  The reason is that complex pairs are supposed
+## FIXME:    to be exact conjugates, and not rely on a tolerance test.
+
+## 2006-05-12 David Bateman - Modified for NDArrays
+
+function y = cplxpair (z, tol, dim)
+
+  if nargin < 1 || nargin > 3
+    print_usage ();
+  endif
+
+  if (length (z) == 0)
+    y = zeros (size (z));
+    return;
+  endif
+
+  if (nargin < 2 || isempty (tol))
+    if (isa (z, "single"))
+      tol = 100 * eps("single");
+    else
+      tol = 100*eps;
+    endif
+  endif
+
+  nd = ndims (z);
+  orig_dims = size (z);
+  if (nargin < 3)
+    ## Find the first singleton dimension.
+    dim = 0;
+    while (dim < nd && orig_dims(dim+1) == 1)
+      dim++;
+    endwhile
+    dim++;
+    if (dim > nd)
+      dim = 1;
+    endif
+  else
+    dim = floor(dim);
+    if (dim < 1 || dim > nd)
+      error ("cplxpair: invalid dimension along which to sort");
+    endif
+  endif
+
+  ## Move dimension to treat first, and convert to a 2-D matrix.
+  perm = [dim:nd, 1:dim-1];
+  z = permute (z, perm);
+  sz = size (z);
+  n = sz (1);
+  m = prod (sz) / n;
+  z = reshape (z, n, m);
+
+  ## Sort the sequence in terms of increasing real values.
+  [q, idx] = sort (real (z), 1);
+  z = z(idx + n * ones (n, 1) * [0:m-1]);
+
+  ## Put the purely real values at the end of the returned list.
+  cls = "double";
+  if (isa (z, "single"))
+    cls = "single";
+  endif
+  [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin(cls)) < tol);
+  q = sparse (idxi, idxj, 1, n, m);
+  nr = sum (q, 1);
+  [q, idx] = sort (q, 1);
+  z = z(idx);
+  y = z;
+
+  ## For each remaining z, place the value and its conjugate at the
+  ## start of the returned list, and remove them from further
+  ## consideration.
+  for j = 1:m
+    p = n - nr(j);
+    for i = 1:2:p
+      if (i+1 > p)
+        error ("cplxpair: could not pair all complex numbers");
+      endif
+      [v, idx] = min (abs (z(i+1:p) - conj (z(i))));
+      if (v > tol)
+        error ("cplxpair: could not pair all complex numbers");
+      endif
+      if (imag (z(i)) < 0)
+        y([i, i+1]) = z([i, idx+i]);
+      else
+        y([i, i+1]) = z([idx+i, i]);
+      endif
+      z(idx+i) = z(i+1);
+    endfor
+  endfor
+
+  ## Reshape the output matrix.
+  y = ipermute (reshape (y, sz), perm);
+
+endfunction
+
+%!demo
+%! [ cplxpair(exp(2i*pi*[0:4]'/5)), exp(2i*pi*[3; 2; 4; 1; 0]/5) ]
+
+%!assert (isempty(cplxpair([])));
+%!assert (cplxpair(1), 1)
+%!assert (cplxpair([1+1i, 1-1i]), [1-1i, 1+1i])
+%!assert (cplxpair([1+1i, 1+1i, 1, 1-1i, 1-1i, 2]), \
+%!        [1-1i, 1+1i, 1-1i, 1+1i, 1, 2])
+%!assert (cplxpair([1+1i; 1+1i; 1; 1-1i; 1-1i; 2]), \
+%!        [1-1i; 1+1i; 1-1i; 1+1i; 1; 2])
+%!assert (cplxpair([0, 1, 2]), [0, 1, 2]);
+
+%!shared z
+%! z=exp(2i*pi*[4; 3; 5; 2; 6; 1; 0]/7);
+%!assert (cplxpair(z(randperm(7))), z);
+%!assert (cplxpair(z(randperm(7))), z);
+%!assert (cplxpair(z(randperm(7))), z);
+%!assert (cplxpair([z(randperm(7)),z(randperm(7))]),[z,z])
+%!assert (cplxpair([z(randperm(7)),z(randperm(7))],[],1),[z,z])
+%!assert (cplxpair([z(randperm(7)).';z(randperm(7)).'],[],2),[z.';z.'])
+
+%!## tolerance test
+%!assert (cplxpair([1i, -1i, 1+(1i*eps)],2*eps), [-1i, 1i, 1+(1i*eps)]);