octave_packages/control-2.3.52/__frequency_vector__.m

   1 ## Copyright (C) 1996, 2000, 2004, 2005, 2006, 2007
   2 ##               Auburn University. All rights reserved.
   3 ##
   4 ##
   5 ## This program is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## This program is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with this program; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {@var{w} =} __frequency_vector__ (@var{sys})
  21 ## Get default range of frequencies based on cutoff frequencies of system
  22 ## poles and zeros.
  23 ## Frequency range is the interval
  24 ## @iftex
  25 ## @tex
  26 ## $ [ 10^{w_{min}}, 10^{w_{max}} ] $
  27 ## @end tex
  28 ## @end iftex
  29 ## @ifnottex
  30 ## [10^@var{wmin}, 10^@var{wmax}]
  31 ## @end ifnottex
  32 ##
  33 ## Used by @command{__frequency_response__}
  34 ## @end deftypefn
  35
  36 ## Adapted-By: Lukas Reichlin <lukas.reichlin@gmail.com>
  37 ## Date: October 2009
  38 ## Version: 0.3
  39
  40 function w = __frequency_vector__ (sys, wbounds = "std", wmin, wmax)
  41
  42   zer = zero (sys);
  43   pol = pole (sys);
  44   tsam = abs (get (sys, "tsam"));        # tsam could be -1
  45   discrete = ! isct (sys);               # static gains (tsam = -2) are assumed continuous
  46
  47   ## make sure zer, pol are row vectors
  48   pol = reshape (pol, 1, []);
  49   zer = reshape (zer, 1, []);
  50
  51   ## check for natural frequencies away from omega = 0
  52   if (discrete)
  53     ## The 2nd conditions prevents log(0) in the next log command
  54     iiz = find (abs(zer-1) > norm(zer)*eps && abs(zer) > norm(zer)*eps);
  55     iip = find (abs(pol-1) > norm(pol)*eps && abs(pol) > norm(pol)*eps);
  56
  57     ## avoid dividing empty matrices, it would work but looks nasty
  58     if (! isempty (iiz))
  59       czer = log (zer(iiz))/tsam;
  60     else
  61       czer = [];
  62     endif
  63
  64     if (! isempty (iip))
  65       cpol = log (pol(iip))/tsam;
  66     else
  67       cpol = [];
  68     endif
  69   else
  70     ## continuous
  71     iip = find (abs(pol) > norm(pol)*eps);
  72     iiz = find (abs(zer) > norm(zer)*eps);
  73
  74     if (! isempty (zer))
  75       czer = zer(iiz);
  76     else
  77       czer = [];
  78     endif
  79     if (! isempty (pol))
  80       cpol = pol(iip);
  81     else
  82       cpol = [];
  83     endif
  84   endif
  85
  86   if (isempty (iip) && isempty (iiz))
  87     ## no poles/zeros away from omega = 0; pick defaults
  88     dec_min = 0;                         # -1
  89     dec_max = 2;                         # 3
  90   else
  91     dec_min = floor (log10 (min (abs ([cpol, czer]))));
  92     dec_max = ceil (log10 (max (abs ([cpol, czer]))));
  93   endif
  94
  95   ## expand to show the entirety of the "interesting" portion of the plot
  96   switch (wbounds)
  97     case "std"                           # standard
  98       if (dec_min == dec_max)
  99         dec_min -= 2;
 100         dec_max += 2;
 101       else
 102         dec_min--;
 103         dec_max++;
 104       endif
 105     case "ext"                           # extended (for nyquist)
 106       if (any (abs (pol) < sqrt (eps)))  # look for integrators
 107         ## dec_min -= 0.5;
 108         dec_max += 2;
 109       else
 110         dec_min -= 2;
 111         dec_max += 2;
 112       endif
 113     otherwise
 114       error ("frequency_range: second argument invalid");
 115   endswitch
 116
 117   ## run discrete frequency all the way to pi
 118   if (discrete)
 119     dec_max = log10 (pi/tsam);
 120   endif
 121
 122   if (nargin == 4)                       # w = {wmin, wmax}
 123     dec_min = log10 (wmin);
 124     dec_max = log10 (wmax);
 125   endif
 126
 127   ## create frequency vector
 128   zp = [abs(zer), abs(pol)];
 129   idx = find (zp > 10^dec_min & zp < 10^dec_max);
 130   zp = zp(idx);
 131
 132   w = logspace (dec_min, dec_max, 500);
 133   w = unique ([w, zp]);                  # unique also sorts frequency vector
 134
 135 endfunction