Mercurial > octave-nkf
comparison scripts/general/cplxpair.m @ 5820:27c966e4b2dc
[project @ 2006-05-17 21:00:54 by jwe]
author | jwe |
---|---|
date | Wed, 17 May 2006 21:00:54 +0000 |
parents | |
children | b0d4ff99a0c5 |
comparison
equal
deleted
inserted
replaced
5819:e54c11df0524 | 5820:27c966e4b2dc |
---|---|
1 ## Copyright (C) 2000 Paul Kienzle | |
2 ## | |
3 ## This file is part of Octave. | |
4 ## | |
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 2, or (at your option) | |
8 ## any later version. | |
9 ## | |
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. | |
14 ## | |
15 ## You should have received a copy of the GNU General Public License | |
16 ## along with Octave; see the file COPYING. If not, write to the Free | |
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA | |
18 ## 02110-1301, USA. | |
19 | |
20 ## -*- texinfo -*- | |
21 ## @deftypefn {Function File} {} cplxpair (@var{z}, @var{tol}, @var{dim}) | |
22 ## Sort the numbers @var{z} into complex conjugate pairs ordered by | |
23 ## increasing real part. With identical real parts, order by increasing | |
24 ## imaginary magnitude. Place the negative imaginary complex number | |
25 ## first within each pair. Place all the real numbers after all the | |
26 ## complex pairs (those with @code {abs ( imag (@var{z}) / @var{z}) < | |
27 ## @var{tol}}), where the default value of @var{tol} is @code{100 * | |
28 ## @var{eps}}. | |
29 ## | |
30 ## By default the complex pairs are sorted along the first non-singleton | |
31 ## dimension of @var{z}. If @var{dim} is specified, then the complex | |
32 ## pairs are sorted along this dimension. | |
33 ## | |
34 ## Signal an error if some complex numbers could not be paired. Requires | |
35 ## all complex numbers to be exact conjugates within tol, or signals an | |
36 ## error. Note that there are no guarantees on the order of the returned | |
37 ## pairs with identical real parts but differing imaginary parts. | |
38 ## | |
39 ## @example | |
40 ## cplxpair (exp(2i*pi*[0:4]'/5)) == exp(2i*pi*[3; 2; 4; 1; 0]/5) | |
41 ## @end example | |
42 ## @end deftypefn | |
43 | |
44 ## TODO: subsort returned pairs by imaginary magnitude | |
45 ## TODO: Why doesn't exp(2i*pi*[0:4]'/5) produce exact conjugates. Does | |
46 ## TODO: it in Matlab? The reason is that complex pairs are supposed | |
47 ## TODO: to be exact conjugates, and not rely on a tolerance test. | |
48 | |
49 ## 2006-05-12 David Bateman - Modified for NDArrays | |
50 | |
51 function y = cplxpair (z, tol, dim) | |
52 | |
53 if nargin < 1 || nargin > 3 | |
54 usage ("z = cplxpair (z, tol, dim);"); | |
55 endif | |
56 | |
57 if (length (z) == 0) | |
58 y = zeros (size (z)); | |
59 return; | |
60 endif | |
61 | |
62 if (nargin < 2 || isempty (tol)) | |
63 tol = 100*eps; | |
64 endif | |
65 | |
66 nd = ndims (z); | |
67 orig_dims = size (z); | |
68 if (nargin < 3) | |
69 ## Find the first singleton dimension. | |
70 dim = 0; | |
71 while (dim < nd && orig_dims(dim+1) == 1) | |
72 dim++; | |
73 endwhile | |
74 dim++; | |
75 if (dim > nd) | |
76 dim = 1; | |
77 endif | |
78 else | |
79 dim = floor(dim); | |
80 if (dim < 1 || dim > nd) | |
81 error ("cplxpair: invalid dimension along which to sort"); | |
82 endif | |
83 endif | |
84 | |
85 ## Move dimension to treat first, and convert to a 2-D matrix | |
86 perm = [dim:nd, 1:dim-1]; | |
87 z = permute (z, perm); | |
88 sz = size (z); | |
89 n = sz (1); | |
90 m = prod (sz) / n; | |
91 z = reshape (z, n, m); | |
92 | |
93 ## Sort the sequence in terms of increasing real values | |
94 [q, idx] = sort (real (z), 1); | |
95 z = z(idx + n * ones (n, 1) * [0:m-1]); | |
96 | |
97 ## Put the purely real values at the end of the returned list | |
98 [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin) < tol ); | |
99 q = sparse (idxi, idxj, 1, n, m); | |
100 nr = sum (q, 1); | |
101 [q, idx] = sort (q, 1); | |
102 z = z(idx); | |
103 y = z; | |
104 | |
105 ## For each remaining z, place the value and its conjugate at the | |
106 ## start of the returned list, and remove them from further | |
107 ## consideration. | |
108 for j = 1:m | |
109 p = n - nr(j); | |
110 for i=1:2:p | |
111 if (i+1 > p) | |
112 error ("cplxpair could not pair all complex numbers"); | |
113 endif | |
114 [v, idx] = min (abs (z(i+1:p) - conj (z(i)))); | |
115 if (v > tol) | |
116 error ("cplxpair could not pair all complex numbers"); | |
117 endif | |
118 if (imag (z(i)) < 0) | |
119 y([i, i+1]) = z([i, idx+i]); | |
120 else | |
121 y([i, i+1]) = z([idx+i, i]); | |
122 endif | |
123 z(idx+i) = z(i+1); | |
124 endfor | |
125 endfor | |
126 | |
127 ## Reshape the output matrix | |
128 y = ipermute (reshape (y, sz), perm); | |
129 | |
130 endfunction | |
131 | |
132 %!demo | |
133 %! [ cplxpair(exp(2i*pi*[0:4]'/5)), exp(2i*pi*[3; 2; 4; 1; 0]/5) ] | |
134 | |
135 %!assert (isempty(cplxpair([]))); | |
136 %!assert (cplxpair(1), 1) | |
137 %!assert (cplxpair([1+1i, 1-1i]), [1-1i, 1+1i]) | |
138 %!assert (cplxpair([1+1i, 1+1i, 1, 1-1i, 1-1i, 2]), \ | |
139 %! [1-1i, 1+1i, 1-1i, 1+1i, 1, 2]) | |
140 %!assert (cplxpair([1+1i; 1+1i; 1; 1-1i; 1-1i; 2]), \ | |
141 %! [1-1i; 1+1i; 1-1i; 1+1i; 1; 2]) | |
142 %!assert (cplxpair([0, 1, 2]), [0, 1, 2]); | |
143 | |
144 %!shared z | |
145 %! z=exp(2i*pi*[4; 3; 5; 2; 6; 1; 0]/7); | |
146 %!assert (cplxpair(z(randperm(7))), z); | |
147 %!assert (cplxpair(z(randperm(7))), z); | |
148 %!assert (cplxpair(z(randperm(7))), z); | |
149 %!assert (cplxpair([z(randperm(7)),z(randperm(7))]),[z,z]) | |
150 %!assert (cplxpair([z(randperm(7)),z(randperm(7))],[],1),[z,z]) | |
151 %!assert (cplxpair([z(randperm(7)).';z(randperm(7)).'],[],2),[z.';z.']) | |
152 | |
153 %!## tolerance test | |
154 %!assert (cplxpair([1i, -1i, 1+(1i*eps)],2*eps), [-1i, 1i, 1+(1i*eps)]); |