1 %% Copyright (C) 2005 Julius O. Smith III <jos@ccrma.stanford.edu>
3 %% This program is free software; you can redistribute it and/or modify it under
4 %% the terms of the GNU General Public License as published by the Free Software
5 %% Foundation; either version 3 of the License, or (at your option) any later
8 %% This program is distributed in the hope that it will be useful, but WITHOUT
9 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 %% You should have received a copy of the GNU General Public License along with
14 %% this program; if not, see <http://www.gnu.org/licenses/>.
17 %% @deftypefn {Function File} {[@var{zc}, @var{zr}] =} cplxreal (@var{z}, @var{thresh})
18 %% Split the vector z into its complex (@var{zc}) and real (@var{zr}) elements,
19 %% eliminating one of each complex-conjugate pair.
24 %% @var{z} = row- or column-vector of complex numbers@*
26 %% @var{thresh} = tolerance threshold for numerical comparisons (default = 100*eps)
32 %% @var{zc} = elements of @var{z} having positive imaginary parts@*
34 %% @var{zr} = elements of @var{z} having zero imaginary part@*
37 %% Each complex element of @var{z} is assumed to have a complex-conjugate
38 %% counterpart elsewhere in @var{z} as well. Elements are declared real
39 %% if their imaginary parts have magnitude less than @var{thresh}.
44 function [zc,zr] = cplxreal (z, thresh = 100*eps)
46 % interesting for testing: if nargin<2, thresh=1E-3; end
52 zcp = cplxpair(z); % sort complex pairs, real roots at end
56 while i && abs(imag(zcp(i)))<thresh % determine no. of real values
57 zcp(i) = real(zcp(i));
63 error('cplxreal: Odd number of complex values!');
66 zc = zcp(2:2:nzsect2);
67 zr = zcp(nzsect2+1:nz);
72 %! [zc,zr] = cplxreal(roots([1 0 0 1]));
73 %! assert({zc,zr},{0.5+i*sin(pi/3),-1},10*eps);