Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[8106] trunk/octave-forge/main/control/devel
|
From: <par...@us...> - 2011-02-11 22:16:22
|
Revision: 8106
http://octave.svn.sourceforge.net/octave/?rev=8106&view=rev
Author: paramaniac
Date: 2011-02-11 22:16:16 +0000 (Fri, 11 Feb 2011)
Log Message:
-----------
control: add first draft code for "multiplot" feature to /devel
Added Paths:
-----------
trunk/octave-forge/main/control/devel/__frequency_response_2__.m
trunk/octave-forge/main/control/devel/__frequency_vector_2__.m
trunk/octave-forge/main/control/devel/bode2.m
trunk/octave-forge/main/control/devel/bode2test.m
trunk/octave-forge/main/control/devel/zeromake.m
Added: trunk/octave-forge/main/control/devel/__frequency_response_2__.m
===================================================================
--- trunk/octave-forge/main/control/devel/__frequency_response_2__.m (rev 0)
+++ trunk/octave-forge/main/control/devel/__frequency_response_2__.m 2011-02-11 22:16:16 UTC (rev 8106)
@@ -0,0 +1,67 @@
+## Copyright (C) 2009, 2010 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope 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.
+##
+## LTI Syncope 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## Return frequency response H and frequency vector w.
+## If w is empty, it will be calculated by __frequency_vector__.
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: November 2009
+## Version: 0.2
+
+% function [H, w] = __frequency_response_2__ (sys, w = [], mimoflag = 0, resptype = 0, wbounds = "std", cellflag = false)
+function [H, w] = __frequency_response_2__ (mimoflag = 0, resptype = 0, wbounds = "std", cellflag = false, varargin)
+
+ sys_vec = cellfun (@(x) isa (x, "lti"), varargin) # true or false
+ sys_idx = find (sys_vec)
+
+ w_vec = cellfun (@is_real_vector, varargin)
+ w_idx = find (w_vec)
+
+% w_idx(end)
+
+ if (isempty (w_idx))
+ w = __frequency_vector_2__ (wbounds, varargin{sys_idx})
+ endif
+
+varargin{sys_idx}
+
+ ## check arguments
+ if(! isa (sys, "lti"))
+ error ("frequency_response: first argument sys must be an LTI system");
+ endif
+
+ if (! mimoflag && ! issiso (sys))
+ error ("frequency_response: require SISO system");
+ endif
+
+ if (isa (sys, "frd"))
+ if (isempty (w))
+ w = get (sys, "w");
+ else
+ warning ("frequency_response: second argument w is ignored");
+ endif
+ elseif (isempty (w)) # find interesting frequency range w if not specified
+ w = __frequency_vector__ (sys, wbounds);
+ elseif (! is_real_vector (w))
+ error ("frequency_response: second argument w must be a vector of frequencies");
+ endif
+
+ ## frequency response
+ H = __freqresp__ (sys, w, resptype, cellflag);
+
+endfunction
\ No newline at end of file
Added: trunk/octave-forge/main/control/devel/__frequency_vector_2__.m
===================================================================
--- trunk/octave-forge/main/control/devel/__frequency_vector_2__.m (rev 0)
+++ trunk/octave-forge/main/control/devel/__frequency_vector_2__.m 2011-02-11 22:16:16 UTC (rev 8106)
@@ -0,0 +1,167 @@
+## Copyright (C) 1996, 2000, 2004, 2005, 2006, 2007
+## Auburn University. All rights reserved.
+##
+##
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{w} =} __frequency_vector__ (@var{sys})
+## Get default range of frequencies based on cutoff frequencies of system
+## poles and zeros.
+## Frequency range is the interval
+## @iftex
+## @tex
+## $ [ 10^{w_{min}}, 10^{w_{max}} ] $
+## @end tex
+## @end iftex
+## @ifinfo
+## [10^@var{wmin}, 10^@var{wmax}]
+## @end ifinfo
+##
+## Used by @command{__frequency_response__}
+## @end deftypefn
+
+## Adapted-By: Lukas Reichlin <luk...@gm...>
+## Date: October 2009
+## Version: 0.2
+
+function w = __frequency_vector_2__ (wbounds = "std", varargin)
+
+ zer = cellfun ("@lti/zero", varargin, "uniformoutput", false)
+ pol = cellfun ("@lti/pole", varargin, "uniformoutput", false)
+
+% zer = cat (1, zer{:})
+% pol = cat (1, pol{:})
+
+ ct_idx = find (cellfun ("@lti/isct", varargin))
+ dt_idx = setdiff (1 : nargin-1, ct_idx)
+
+ %if (length (dt_idx) > 0)
+ tsam = cellfun (@(x) abs (get (x, "tsam")), varargin(dt_idx))
+ %endif
+
+% zer_ct = cat (1, zer{ct_idx})
+% pol_ct = cat (1, pol{ct_idx})
+
+% zer_dt = cat (1, zer{dt_idx})
+% pol_dt = cat (1, pol{dt_idx})
+
+ zer_ct = reshape (cat (1, zer{ct_idx}), 1, [])
+ pol_ct = reshape (cat (1, pol{ct_idx}), 1, [])
+
+ zer_dt = reshape (cat (1, zer{dt_idx}), 1, [])
+ pol_dt = reshape (cat (1, pol{dt_idx}), 1, [])
+
+
+ if (length (ct_idx) > 0) ## continuous
+ iip = find (abs(pol) > norm(pol)*eps);
+ iiz = find (abs(zer) > norm(zer)*eps);
+ czer = zer(iiz);
+ cpol = pol(iip);
+ endif
+
+% zer = cat (1, zer{:})
+% pol = cat (1, pol{:})
+
+
+% zer = zero (sys);
+% pol = pole (sys);
+% tsam = abs (get (sys, "tsam")); # tsam could be -1
+% discrete = ! isct (sys); # static gains (tsam = -2) are assumed continuous
+
+ ## make sure zer, pol are row vectors
+ pol = reshape (pol_ct, 1, [])
+ zer = reshape (zer, 1, [])
+
+ ## check for natural frequencies away from omega = 0
+ if (discrete)
+ ## The 2nd conditions prevents log(0) in the next log command
+ iiz = find (abs(zer-1) > norm(zer)*eps && abs(zer) > norm(zer)*eps);
+ iip = find (abs(pol-1) > norm(pol)*eps && abs(pol) > norm(pol)*eps);
+
+ ## avoid dividing empty matrices, it would work but looks nasty
+ if (! isempty (iiz))
+ czer = log (zer(iiz))/tsam;
+ else
+ czer = [];
+ endif
+
+ if (! isempty (iip))
+ cpol = log (pol(iip))/tsam;
+ else
+ cpol = [];
+ endif
+ else
+ ## continuous
+ iip = find (abs(pol) > norm(pol)*eps);
+ iiz = find (abs(zer) > norm(zer)*eps);
+
+ if (! isempty (zer))
+ czer = zer(iiz);
+ else
+ czer = [];
+ endif
+ if (! isempty (pol))
+ cpol = pol(iip);
+ else
+ cpol = [];
+ endif
+ endif
+
+ if (isempty (iip) && isempty (iiz))
+ ## no poles/zeros away from omega = 0; pick defaults
+ dec_min = 0; # -1
+ dec_max = 2; # 3
+ else
+ dec_min = floor (log10 (min (abs ([cpol, czer]))));
+ dec_max = ceil (log10 (max (abs ([cpol, czer]))));
+ endif
+
+ ## expand to show the entirety of the "interesting" portion of the plot
+ switch (wbounds)
+ case "std" # standard
+ if (dec_min == dec_max)
+ dec_min -= 2;
+ dec_max += 2;
+ else
+ dec_min--;
+ dec_max++;
+ endif
+ case "ext" # extended (for nyquist)
+ if (any (abs (pol) < sqrt (eps))) # look for integrators
+ ## dec_min -= 0.5;
+ dec_max += 2;
+ else
+ dec_min -= 2;
+ dec_max += 2;
+ endif
+ otherwise
+ error ("frequency_range: second argument invalid");
+ endswitch
+
+ ## run discrete frequency all the way to pi
+ if (discrete)
+ dec_max = log10 (pi/tsam);
+ endif
+
+ ## create frequency vector
+ zp = [abs(zer), abs(pol)];
+ idx = find (zp > 10^dec_min & zp < 10^dec_max);
+ zp = zp(idx);
+
+ w = logspace (dec_min, dec_max, 500);
+ w = unique ([w, zp]); # unique also sorts frequency vector
+
+endfunction
Added: trunk/octave-forge/main/control/devel/bode2.m
===================================================================
--- trunk/octave-forge/main/control/devel/bode2.m (rev 0)
+++ trunk/octave-forge/main/control/devel/bode2.m 2011-02-11 22:16:16 UTC (rev 8106)
@@ -0,0 +1,98 @@
+## Copyright (C) 2009, 2010 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope 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.
+##
+## LTI Syncope 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{mag}, @var{pha}, @var{w}] =} bode (@var{sys})
+## @deftypefnx {Function File} {[@var{mag}, @var{pha}, @var{w}] =} bode (@var{sys}, @var{w})
+## Bode diagram of frequency response. If no output arguments are given,
+## the response is printed on the screen.
+##
+## @strong{Inputs}
+## @table @var
+## @item sys
+## LTI system. Must be a single-input and single-output (SISO) system.
+## @item w
+## Optional vector of frequency values. If @var{w} is not specified,
+## it is calculated by the zeros and poles of the system.
+## @end table
+##
+## @strong{Outputs}
+## @table @var
+## @item mag
+## Vector of magnitude. Has length of frequency vector @var{w}.
+## @item pha
+## Vector of phase. Has length of frequency vector @var{w}.
+## @item w
+## Vector of frequency values used.
+## @end table
+##
+## @seealso{nichols, nyquist, sigma}
+## @end deftypefn
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: November 2009
+## Version: 0.2
+
+function [mag_r, pha_r, w_r] = bode2 (varargin)
+
+ ## TODO: multiplot feature: bode (sys1, "b", sys2, "r", ...)
+
+% if (nargin == 0 || nargin > 2)
+% print_usage ();
+% endif
+
+ [H, w] = __frequency_response_2__ (false, 0, "std", false, varargin{:});
+
+ H = reshape (H, [], 1);
+ mag = abs (H);
+ pha = unwrap (arg (H)) * 180 / pi;
+
+ if (! nargout)
+ mag_db = 20 * log10 (mag);
+
+ wv = [min(w), max(w)];
+ ax_vec_mag = __axis_limits__ ([reshape(w, [], 1), reshape(mag_db, [], 1)]);
+ ax_vec_mag(1:2) = wv;
+ ax_vec_pha = __axis_limits__ ([reshape(w, [], 1), reshape(pha, [], 1)]);
+ ax_vec_pha(1:2) = wv;
+
+ if (isct (sys))
+ xl_str = "Frequency [rad/s]";
+ else
+ xl_str = sprintf ("Frequency [rad/s] w_N = %g", pi / get (sys, "tsam"));
+ endif
+
+ subplot (2, 1, 1)
+ semilogx (w, mag_db)
+ axis (ax_vec_mag)
+ grid ("on")
+ title (["Bode Diagram of ", inputname(1)])
+ ylabel ("Magnitude [dB]")
+
+ subplot (2, 1, 2)
+ semilogx (w, pha)
+ axis (ax_vec_pha)
+ grid ("on")
+ xlabel (xl_str)
+ ylabel ("Phase [deg]")
+ else
+ mag_r = mag;
+ pha_r = pha;
+ w_r = w;
+ endif
+
+endfunction
\ No newline at end of file
Added: trunk/octave-forge/main/control/devel/bode2test.m
===================================================================
--- trunk/octave-forge/main/control/devel/bode2test.m (rev 0)
+++ trunk/octave-forge/main/control/devel/bode2test.m 2011-02-11 22:16:16 UTC (rev 8106)
@@ -0,0 +1,12 @@
+A = ss (-2, 3, 4, 5);
+B = ss (-1, 1, 1, 0);
+C = ss (-4, 7, 10, 12);
+w = [1 2 3 4 5 6 7];
+
+%bode2 (A, B, "b", w, C)
+%bode2 (A, B, "b", C)
+%bode2 (A, B, "b", C, WestlandLynx)
+%bode2 (WestlandLynx)
+%bode2 (A, B, "b", c2d (C, 0.375), WestlandLynx)
+bode2 (A, c2d (B, 0.223), "b", c2d (C, 0.375), WestlandLynx)
+
Added: trunk/octave-forge/main/control/devel/zeromake.m
===================================================================
--- trunk/octave-forge/main/control/devel/zeromake.m (rev 0)
+++ trunk/octave-forge/main/control/devel/zeromake.m 2011-02-11 22:16:16 UTC (rev 8106)
@@ -0,0 +1 @@
+mkoctfile TG04BX.f MB02RD.f MB02SD.f
\ No newline at end of file
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
