octave_packages/control-2.3.52/sigma.m

   1 ## Copyright (C) 2009, 2010, 2011, 2012   Lukas F. Reichlin
   2 ##
   3 ## This file is part of LTI Syncope.
   4 ##
   5 ## LTI Syncope is free software: you can redistribute it and/or modify
   6 ## it under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation, either version 3 of the License, or
   8 ## (at your option) any later version.
   9 ##
  10 ## LTI Syncope is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with LTI Syncope.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn{Function File} {[@var{sv}, @var{w}] =} sigma (@var{sys})
  20 ## @deftypefnx{Function File} {[@var{sv}, @var{w}] =} sigma (@var{sys}, @var{w})
  21 ## @deftypefnx{Function File} {[@var{sv}, @var{w}] =} sigma (@var{sys}, @var{[]}, @var{ptype})
  22 ## @deftypefnx{Function File} {[@var{sv}, @var{w}] =} sigma (@var{sys}, @var{w}, @var{ptype})
  23 ## Singular values of frequency response.  If no output arguments are given,
  24 ## the singular value plot is printed on the screen;
  25 ##
  26 ## @strong{Inputs}
  27 ## @table @var
  28 ## @item sys
  29 ## LTI system.  Multiple inputs and/or outputs (MIMO systems) make practical sense.
  30 ## @item w
  31 ## Optional vector of frequency values.  If @var{w} is not specified,
  32 ## it is calculated by the zeros and poles of the system.
  33 ## Alternatively, the cell @code{@{wmin, wmax@}} specifies a frequency range,
  34 ## where @var{wmin} and @var{wmax} denote minimum and maximum frequencies
  35 ## in rad/s.
  36 ## @item ptype = 0
  37 ## Singular values of the frequency response @var{H} of system @var{sys}.  Default Value.
  38 ## @item ptype = 1
  39 ## Singular values of the frequency response @code{inv(H)}; i.e. inversed system.
  40 ## @item ptype = 2
  41 ## Singular values of the frequency response @code{I + H}; i.e. inversed sensitivity
  42 ## (or return difference) if @code{H = P * C}.
  43 ## @item ptype = 3
  44 ## Singular values of the frequency response @code{I + inv(H)}; i.e. inversed complementary
  45 ## sensitivity if @code{H = P * C}.
  46 ## @end table
  47 ##
  48 ## @strong{Outputs}
  49 ## @table @var
  50 ## @item sv
  51 ## Array of singular values.  For a system with m inputs and p outputs, the array sv
  52 ## has @code{min (m, p)} rows and as many columns as frequency points @code{length (w)}.
  53 ## The singular values at the frequency @code{w(k)} are given by @code{sv(:,k)}.
  54 ## @item w
  55 ## Vector of frequency values used.
  56 ## @end table
  57 ##
  58 ## @seealso{bodemag, svd}
  59 ## @end deftypefn
  60
  61 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  62 ## Created: May 2009
  63 ## Version: 0.6
  64
  65 function [sv_r, w_r] = sigma (sys, w = [], resptype = 0)
  66
  67   ## TODO: multiplot feature:   sigma (sys1, "b", sys2, "r", ...)
  68
  69   if (nargin == 0 || nargin > 3)
  70     print_usage ();
  71   endif
  72
  73   [H, w] = __frequency_response__ (sys, w, true, resptype, "std", true);
  74
  75   sv = cellfun (@svd, H, "uniformoutput", false);
  76   sv = horzcat (sv{:});
  77
  78   if (! nargout)  # plot the information
  79
  80     ## convert to dB for plotting
  81     sv_db = 20 * log10 (sv);
  82
  83     ## determine xlabel
  84     if (isct (sys))
  85       xl_str = "Frequency [rad/s]";
  86     else
  87       xl_str = sprintf ("Frequency [rad/s]     w_N = %g", pi / get (sys, "tsam"));
  88     endif
  89
  90     ## plot results
  91     semilogx (w, sv_db, "b")
  92     axis ("tight")
  93     ylim (__axis_margin__ (ylim))
  94     grid ("on")
  95     title (["Singular Values of ", inputname(1)])
  96     xlabel (xl_str)
  97     ylabel ("Singular Values [dB]")
  98   else            # return values
  99     sv_r = sv;
 100     w_r = reshape (w, [], 1);
 101   endif
 102
 103 endfunction
 104
 105
 106 %!shared sv_exp, w_exp, sv_obs, w_obs
 107 %! A = [1, 2; 3, 4];
 108 %! B = [5, 6; 7, 8];
 109 %! C = [4, 3; 2, 1];
 110 %! D = [8, 7; 6, 5];
 111 %! w = [2, 3, 4];
 112 %! sv_exp = [7.9176, 8.6275, 9.4393;
 113 %!           0.6985, 0.6086, 0.5195];
 114 %! w_exp = [2; 3; 4];
 115 %! [sv_obs, w_obs] = sigma (ss (A, B, C, D), w);
 116 %!assert (sv_obs, sv_exp, 1e-4);
 117 %!assert (w_obs, w_exp, 1e-4);