1 ## Copyright (C) 2000-2012 Paul Kienzle
3 ## This file is part of Octave.
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 3 of the License, or (at
8 ## your option) any later version.
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
20 ## @deftypefn {Function File} {} cplxpair (@var{z})
21 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol})
22 ## @deftypefnx {Function File} {} cplxpair (@var{z}, @var{tol}, @var{dim})
23 ## Sort the numbers @var{z} into complex conjugate pairs ordered by
24 ## increasing real part. Place the negative imaginary complex number
25 ## first within each pair. Place all the real numbers (those with
26 ## @code{abs (imag (@var{z}) / @var{z}) < @var{tol})}) after the
29 ## If @var{tol} is unspecified the default value is 100*@code{eps}.
31 ## By default the complex pairs are sorted along the first non-singleton
32 ## dimension of @var{z}. If @var{dim} is specified, then the complex
33 ## pairs are sorted along this dimension.
35 ## Signal an error if some complex numbers could not be paired. Signal an
36 ## error if all complex numbers are not exact conjugates (to within
37 ## @var{tol}). Note that there is no defined order for pairs with identical
38 ## real parts but differing imaginary parts.
39 ## @c Set example in small font to prevent overfull line
42 ## cplxpair (exp(2i*pi*[0:4]'/5)) == exp(2i*pi*[3; 2; 4; 1; 0]/5)
46 ## FIXME: subsort returned pairs by imaginary magnitude
47 ## FIXME: Why doesn't exp(2i*pi*[0:4]'/5) produce exact conjugates. Does
48 ## FIXME: it in Matlab? The reason is that complex pairs are supposed
49 ## FIXME: to be exact conjugates, and not rely on a tolerance test.
51 ## 2006-05-12 David Bateman - Modified for NDArrays
53 function y = cplxpair (z, tol, dim)
55 if nargin < 1 || nargin > 3
64 if (nargin < 2 || isempty (tol))
65 if (isa (z, "single"))
66 tol = 100 * eps("single");
75 ## Find the first singleton dimension.
77 while (dim < nd && orig_dims(dim+1) == 1)
86 if (dim < 1 || dim > nd)
87 error ("cplxpair: invalid dimension along which to sort");
91 ## Move dimension to treat first, and convert to a 2-D matrix.
92 perm = [dim:nd, 1:dim-1];
93 z = permute (z, perm);
97 z = reshape (z, n, m);
99 ## Sort the sequence in terms of increasing real values.
100 [q, idx] = sort (real (z), 1);
101 z = z(idx + n * ones (n, 1) * [0:m-1]);
103 ## Put the purely real values at the end of the returned list.
105 if (isa (z, "single"))
108 [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin(cls)) < tol);
109 q = sparse (idxi, idxj, 1, n, m);
111 [q, idx] = sort (q, 1);
115 ## For each remaining z, place the value and its conjugate at the
116 ## start of the returned list, and remove them from further
122 error ("cplxpair: could not pair all complex numbers");
124 [v, idx] = min (abs (z(i+1:p) - conj (z(i))));
126 error ("cplxpair: could not pair all complex numbers");
129 y([i, i+1]) = z([i, idx+i]);
131 y([i, i+1]) = z([idx+i, i]);
137 ## Reshape the output matrix.
138 y = ipermute (reshape (y, sz), perm);
143 %! [ cplxpair(exp(2i*pi*[0:4]'/5)), exp(2i*pi*[3; 2; 4; 1; 0]/5) ]
145 %!assert (isempty(cplxpair([])));
146 %!assert (cplxpair(1), 1)
147 %!assert (cplxpair([1+1i, 1-1i]), [1-1i, 1+1i])
148 %!assert (cplxpair([1+1i, 1+1i, 1, 1-1i, 1-1i, 2]), \
149 %! [1-1i, 1+1i, 1-1i, 1+1i, 1, 2])
150 %!assert (cplxpair([1+1i; 1+1i; 1; 1-1i; 1-1i; 2]), \
151 %! [1-1i; 1+1i; 1-1i; 1+1i; 1; 2])
152 %!assert (cplxpair([0, 1, 2]), [0, 1, 2]);
155 %! z=exp(2i*pi*[4; 3; 5; 2; 6; 1; 0]/7);
156 %!assert (cplxpair(z(randperm(7))), z);
157 %!assert (cplxpair(z(randperm(7))), z);
158 %!assert (cplxpair(z(randperm(7))), z);
159 %!assert (cplxpair([z(randperm(7)),z(randperm(7))]),[z,z])
160 %!assert (cplxpair([z(randperm(7)),z(randperm(7))],[],1),[z,z])
161 %!assert (cplxpair([z(randperm(7)).';z(randperm(7)).'],[],2),[z.';z.'])
164 %!assert (cplxpair([1i, -1i, 1+(1i*eps)],2*eps), [-1i, 1i, 1+(1i*eps)]);