## Copyright (C) 2005 Julius O. Smith III ## Copyright (C) 2018-2019 Mike Miller ## ## This program is free software: you can redistribute it and/or modify ## it under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## This program is distributed in the hope that it will be useful, ## but WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with this program; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{zc}, @var{zr}] =} cplxreal (@var{z}) ## @deftypefnx {Function File} {[@var{zc}, @var{zr}] =} cplxreal (@var{z}, @var{tol}) ## @deftypefnx {Function File} {[@var{zc}, @var{zr}] =} cplxreal (@var{z}, @var{tol}, @var{dim}) ## Sort the numbers @var{z} into complex-conjugate-valued and real-valued ## elements. The positive imaginary complex numbers of each complex conjugate ## pair are returned in @var{zc} and the real numbers are returned in @var{zr}. ## ## @var{tol} is a weighting factor in the range [0, 1) which determines the ## tolerance of the matching. The default value is @code{100 * eps} and the ## resulting tolerance for a given complex pair is ## @code{@var{tol} * abs (@var{z}(i)))}. ## ## By default the complex pairs are sorted along the first non-singleton ## dimension of @var{z}. If @var{dim} is specified, then the complex pairs are ## sorted along this dimension. ## ## Signal an error if some complex numbers could not be paired. Signal an ## error if all complex numbers are not exact conjugates (to within @var{tol}). ## Note that there is no defined order for pairs with identical real parts but ## differing imaginary parts. ## @seealso{cplxpair} ## @end deftypefn function [zc, zr] = cplxreal (z, tol, dim) if (nargin < 1 || nargin > 3) print_usage (); endif if (isempty (z)) zc = zeros (size (z)); zr = zeros (size (z)); return; endif cls = ifelse (isa (z, "single"), "single", "double"); if (nargin < 2 || isempty (tol)) tol = 100 * eps (cls); endif args = cell (1, nargin); args{1} = z; args{2} = tol; if (nargin >= 3) args{3} = dim; endif zcp = cplxpair (args{:}); nz = length (z); idx = nz; # while ((idx > 0) && ((abs (imag (zcp(idx))) ./ abs (zcp(idx))) <= tol)) # not ok for real zero at 0 (handled as complex) while ((idx > 0) && (((abs (imag (zcp(idx))) ./ abs (zcp(idx))) <= tol) || (abs (zcp(idx)) <= tol))) # real zero at zero no more handled as complex zcp(idx) = real (zcp(idx)); idx--; endwhile if (mod (idx, 2) != 0) error ("cplxreal: odd number of complex values was returned from cplxpair"); endif zc = zcp(2:2:idx); zr = zcp(idx+1:nz); endfunction %!test %! [zc, zr] = cplxreal ([]); %! assert (isempty (zc)) %! assert (isempty (zr)) %!test %! [zc, zr] = cplxreal (1); %! assert (isempty (zc)) %! assert (zr, 1) %!test %! [zc, zr] = cplxreal ([1+1i, 1-1i]); %! assert (zc, 1+1i) %! assert (isempty (zr)) %!test %! [zc, zr] = cplxreal (roots ([1, 0, 0, 1])); %! assert (zc, complex (0.5, sin (pi/3)), 10*eps) %! assert (zr, -1, 2*eps) ## Test input validation %!error cplxreal () %!error cplxreal (1, 2, 3, 4) %!error cplxreal (1, ones (2, 3)) %!error cplxreal (1, -1) %!error cplxreal (1, [], 3)