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