|
0
|
1 ## Copyright (C) 2000 Paul Kienzle |
|
|
2 ## |
|
|
3 ## This program is free software; you can redistribute it and/or modify |
|
|
4 ## it under the terms of the GNU General Public License as published by |
|
|
5 ## the Free Software Foundation; either version 2 of the License, or |
|
|
6 ## (at your option) any later version. |
|
|
7 ## |
|
|
8 ## This program is distributed in the hope that it will be useful, |
|
|
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
11 ## GNU General Public License for more details. |
|
|
12 ## |
|
|
13 ## You should have received a copy of the GNU General Public License |
|
|
14 ## along with this program; if not, write to the Free Software |
|
|
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
|
|
16 |
|
|
17 ## usage: [z, n] = cplxpair(z [,tol]) |
|
|
18 ## |
|
|
19 ## Sort the numbers z into complex conjugate pairs ordered by increasing |
|
|
20 ## real part, or with identical real parts, increasing imaginary |
|
|
21 ## magnitude. Place the negative imaginary complex number first within |
|
|
22 ## each pair. Place all the real numbers after all the complex pairs |
|
|
23 ## (those with abs(imag(z)/z) < tol). Signal an error if some complex |
|
|
24 ## numbers could not be paired. Requires all complex numbers to be exact |
|
|
25 ## conjugates within tol, or signals an error. tol defaults to 100*eps. |
|
|
26 ## Note that there are no guarantees on the order of the returned pairs |
|
|
27 ## with identical real parts but differing imaginary parts. |
|
|
28 ## |
|
|
29 ## Returns the ordered list of values, plus the number of values which |
|
|
30 ## were complex. |
|
|
31 ## |
|
|
32 ## Example |
|
|
33 ## [ cplxpair (exp(2i*pi*[0:4]'/5)), exp(2i*pi*[3; 2; 4; 1; 0]/5) ] |
|
|
34 |
|
|
35 ## TODO: subsort returned pairs by imaginary magnitude |
|
|
36 ## TODO: Consider removing extra return value (the number of complex |
|
|
37 ## TODO: values in the list) for exact compatibility |
|
|
38 ## TODO: Why doesn't exp(2i*pi*[0:4]'/5) produce exact conjugates. Does |
|
|
39 ## TODO: it in Matlab? The reason is that complex pairs are supposed |
|
|
40 ## TODO: to be exact conjugates, and not rely on a tolerance test. |
|
|
41 function [y, n] = cplxpair(z, tol) |
|
|
42 |
|
|
43 if nargin < 1 || nargin > 2 |
|
|
44 usage ("[z, n] = cplxpair (z [, tol]);"); |
|
|
45 endif |
|
|
46 if nargin < 2, tol = 100*eps; endif |
|
|
47 |
|
|
48 y = zeros (size (z)); |
|
|
49 if (length (z) == 0), return; endif |
|
|
50 |
|
|
51 ## Sort the sequence in terms of increasing real values |
|
|
52 [q, idx] = sort (real (z)); |
|
|
53 z = z (idx); |
|
|
54 |
|
|
55 ## Put the purely real values at the end of the returned list |
|
|
56 idx = find ( abs(imag(z)) ./ (abs(z)+realmin) < tol ); |
|
|
57 n = length (z) - length (idx); |
|
|
58 if (length (idx) > 0) |
|
|
59 y (n+1 : length(z)) = z (idx); |
|
|
60 z (idx) = []; |
|
|
61 endif |
|
|
62 |
|
|
63 ## For each remaining z, place the value and its conjugate at the |
|
|
64 ## start of the returned list, and remove them from further |
|
|
65 ## consideration. |
|
|
66 for i=1:2:n |
|
|
67 if (i+1 > n) |
|
|
68 error ("cplxpair could not pair all complex numbers"); |
|
|
69 endif |
|
|
70 [v, idx] = min ( abs (z(i+1:n) - conj(z(i))) ); |
|
|
71 if (v > tol) |
|
|
72 error ("cplxpair could not pair all complex numbers"); |
|
|
73 endif |
|
|
74 if imag(z(i)) < 0 |
|
|
75 y ([i, i+1]) = z ([i, idx+i]); |
|
|
76 else |
|
|
77 y ([i, i+1]) = z ([idx+i, i]); |
|
|
78 endif |
|
|
79 z (idx+i) = z (i+1); |
|
|
80 endfor |
|
|
81 |
|
|
82 endfunction |
|
|
83 |
|
|
84 %!demo |
|
|
85 %! [ cplxpair (exp(2i*pi*[0:4]'/5)), exp (2i*pi*[3; 2; 4; 1; 0]/5) ] |
|
|
86 |
|
|
87 %!assert (isempty(cplxpair([]))); |
|
|
88 %!assert (cplxpair(1), 1) |
|
|
89 %!assert (cplxpair([1+1i, 1-1i]), [1-1i, 1+1i]) |
|
|
90 %!assert (cplxpair([1+1i, 1+1i, 1, 1-1i, 1-1i, 2]), \ |
|
|
91 %! [1-1i, 1+1i, 1-1i, 1+1i, 1, 2]) |
|
|
92 %!assert (cplxpair([1+1i; 1+1i; 1; 1-1i; 1-1i; 2]), \ |
|
|
93 %! [1-1i; 1+1i; 1-1i; 1+1i; 1; 2]) |
|
|
94 %!assert (cplxpair([0, 1, 2]), [0, 1, 2]); |
|
|
95 |
|
|
96 %!shared z |
|
|
97 %! z=exp(2i*pi*[4; 3; 5; 2; 6; 1; 0]/7); |
|
|
98 %!assert (cplxpair(z(randperm(7))), z); |
|
|
99 %!assert (cplxpair(z(randperm(7))), z); |
|
|
100 %!assert (cplxpair(z(randperm(7))), z); |
|
|
101 |
|
|
102 %!## tolerance test |
|
|
103 %!assert (cplxpair([1i, -1i, 1+(1i*eps)],2*eps), [-1i, 1i, 1+(1i*eps)]); |
