1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
function zsort = cplxpair (z, tol, dim)
// Description
// Sort the numbers z into complex conjugate pairs ordered by increasing real part.
// The negative imaginary complex numbers are placed first within each pair. All real numbers (those with abs (imag (z) / z) < tol) are placed after the complex pairs.
// and the resulting tolerance for a given complex pair is 100 * eps (abs (z(i))).
// Signal an error if some complex numbers could not be paired. Signal an error if all complex numbers are not exact conjugates (to within tol).
// Note that there is no defined order for pairs with identical real parts but differing imaginary parts
// Calling Sequence
// cplxpair (z)
// cplxpair (z, tol)
// cplxpair (z, tol, dim)
// Parameters:
// z is a matrix or vector.
// tol is a weighting factor which determines the tolerance of matching. The default value is 100.
// By default the complex pairs are sorted along the first non-singleton dimension of z. If dim is specified, then the complex pairs are sorted along this dimension.
// Dependencies: ipermute
// Example:
// The following code
// [ cplxpair(exp(2*%i*%pi*[0:4]'/5)), exp(2*%i*%pi*[3; 2; 4; 1; 0]/5) ]
// Produces the following output
// ans =
// -0.80902 - 0.58779i -0.80902 - 0.58779i
// -0.80902 + 0.58779i -0.80902 + 0.58779i
// 0.30902 - 0.95106i 0.30902 - 0.95106i
// 0.30902 + 0.95106i 0.30902 + 0.95106i
// 1.00000 + 0.00000i 1.00000 + 0.00000i
if (nargin < 1)
error("Invalid inputs");
end
// default double
realmin = 2.2251e-308
if (isempty (z))
zsort = zeros (size (z,1) , size (z,2));
return;
end
if (nargin < 2 || isempty (tol))
tol = 100* %eps;
elseif (~ isscalar (tol) || tol < 0 || tol >= 1)
error ("cplxpair: TOL must be a scalar number in the range 0 <= TOL < 1");
end
nd = ndims (z);
if (nargin < 3)
// Find the first singleton dimension.
sz = size (z);
dim = find (sz > 1, 1);
if isempty(dim)
dim = 1;
end
else
dim = floor (dim);
if (dim < 1 || dim > nd)
error ("cplxpair: invalid dimension DIM");
end
end
// Move dimension to analyze to first position, 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 = matrix (z, n, m);
// Sort the sequence in terms of increasing real values.
[temp, idx] = gsort (real (z), 1 , "i");
z = z(idx + n * ones (n, 1) * [0:m-1]);
// Put the purely real values at the end of the returned list.
[idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin) <= tol);
// Force values detected to be real within tolerance to actually be real.
z(idxi + n*(idxj-1)) = real (z(idxi + n*(idxj-1)));
//if idxi and idxj are not scalers
if ~isscalar(idxi) then
v = ones(size(idxi,1),size(idxi,2));
else
v = 1 ;
end
q = sparse ([idxi' idxj'], v, [n m]);
nr = sum (q, 1);
[temp, idx] = gsort (q, 'r','i');
midx = idx + size (idx,1) * ones (size (idx,1), 1) * [0:size(idx,2)-1];
z = z(midx);
zsort = 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");
end
[v, idx] = min (abs (z(i+1:p,j) - conj (z(i,j))));
if (v >= tol * abs (z(i,j)))
error ("cplxpair: could not pair all complex numbers");
end
// For pairs, select the one with positive imaginary part and use it and
// it's conjugate, but list the negative imaginary pair first.
if (imag (z(i,j)) > 0)
zsort([i, i+1],j) = [conj(z(i,j)), z(i,j)];
else
zsort([i, i+1],j) = [conj(z(idx+i,j)), z(idx+i,j)];
end
z(idx+i,j) = z(i+1,j);
end
end
// Reshape the output matrix.
zsort = ipermute (matrix (zsort, sz), perm);
endfunction
|
