octave_packages/signal-1.1.3/cplxreal.m

   1 %% Copyright (C) 2005 Julius O. Smith III <jos@ccrma.stanford.edu>
   2 %%
   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
   6 %% version.
   7 %%
   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
  11 %% details.
  12 %%
  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/>.
  15
  16 %% -*- texinfo -*-
  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.
  20 %%
  21 %% INPUTS:@*
  22 %%   @itemize
  23 %%   @item
  24 %%   @var{z}      = row- or column-vector of complex numbers@*
  25 %%   @item
  26 %%   @var{thresh} = tolerance threshold for numerical comparisons (default = 100*eps)
  27 %%   @end itemize
  28 %%
  29 %% RETURNED:@*
  30 %%   @itemize
  31 %%   @item
  32 %% @var{zc} = elements of @var{z} having positive imaginary parts@*
  33 %%   @item
  34 %% @var{zr} = elements of @var{z} having zero imaginary part@*
  35 %%   @end itemize
  36 %%
  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}.
  40 %%
  41 %% @seealso{cplxpair}
  42 %% @end deftypefn
  43
  44 function [zc,zr] = cplxreal (z, thresh = 100*eps)
  45
  46   % interesting for testing: if nargin<2, thresh=1E-3; end
  47
  48   if isempty(z)
  49     zc=[];
  50     zr=[];
  51   else
  52     zcp = cplxpair(z); % sort complex pairs, real roots at end
  53     nz = length(z);
  54     nzrsec = 0;
  55     i=nz;
  56     while i && abs(imag(zcp(i)))<thresh % determine no. of real values
  57       zcp(i) = real(zcp(i));
  58       nzrsec = nzrsec+1;
  59       i=i-1;
  60     end
  61     nzsect2 = nz-nzrsec;
  62     if mod(nzsect2,2)~=0
  63       error('cplxreal: Odd number of complex values!');
  64     end
  65     nzsec = nzsect2/2;
  66     zc = zcp(2:2:nzsect2);
  67     zr = zcp(nzsect2+1:nz);
  68   end
  69 endfunction
  70
  71 %!test
  72 %! [zc,zr] = cplxreal(roots([1 0 0 1]));
  73 %! assert({zc,zr},{0.5+i*sin(pi/3),-1},10*eps);