Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[11011] trunk/octave-forge/main/control
|
From: <par...@us...> - 2012-09-14 11:40:49
|
Revision: 11011
http://octave.svn.sourceforge.net/octave/?rev=11011&view=rev
Author: paramaniac
Date: 2012-09-14 11:40:39 +0000 (Fri, 14 Sep 2012)
Log Message:
-----------
control: move multiplot functions into place
Modified Paths:
--------------
trunk/octave-forge/main/control/devel/bode2.m
trunk/octave-forge/main/control/inst/__frequency_response__.m
trunk/octave-forge/main/control/inst/__frequency_vector__.m
trunk/octave-forge/main/control/inst/bode.m
trunk/octave-forge/main/control/inst/nyquist.m
Modified: trunk/octave-forge/main/control/devel/bode2.m
===================================================================
--- trunk/octave-forge/main/control/devel/bode2.m 2012-09-14 11:33:41 UTC (rev 11010)
+++ trunk/octave-forge/main/control/devel/bode2.m 2012-09-14 11:40:39 UTC (rev 11011)
@@ -109,4 +109,4 @@
w_r = w{1};
endif
-endfunction
\ No newline at end of file
+endfunction
Modified: trunk/octave-forge/main/control/inst/__frequency_response__.m
===================================================================
--- trunk/octave-forge/main/control/inst/__frequency_response__.m 2012-09-14 11:33:41 UTC (rev 11010)
+++ trunk/octave-forge/main/control/inst/__frequency_response__.m 2012-09-14 11:40:39 UTC (rev 11011)
@@ -21,35 +21,47 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: November 2009
-## Version: 0.4
+## Version: 0.5
-function [H, w] = __frequency_response__ (sys, w = [], mimoflag = 0, resptype = 0, wbounds = "std", cellflag = false)
+function [H, w] = __frequency_response__ (args, mimoflag = 0, resptype = 0, wbounds = "std", cellflag = false)
+ %if (! iscell (args))
+ % args = {args};
+ %endif
+
+ sys_idx = cellfun (@isa, args, {"lti"}); # look for LTI models
+ w_idx = cellfun (@is_real_vector, args); # look for frequency vectors
+ r_idx = cellfun (@iscell, args); # look for frequency ranges {wmin, wmax}
+
+ sys_cell = args(sys_idx); # extract LTI models
+ frd_idx = cellfun (@isa, sys_cell, {"frd"}); # look for FRD models
+
## check arguments
- if(! isa (sys, "lti"))
- error ("frequency_response: first argument sys must be an LTI system");
+ if (! mimoflag && ! all (cellfun (@issiso, sys_cell)))
+ error ("frequency_response: require SISO systems");
endif
- if (! mimoflag && ! issiso (sys))
- error ("frequency_response: require SISO system");
- endif
-
- if (isa (sys, "frd"))
- if (! isempty (w))
- warning ("frequency_response: second argument w is ignored");
+ ## determine frequencies
+ if (any (r_idx)) # if there are frequency ranges
+ r = args(r_idx){end}; # take the last one
+ if (numel (r) == 2 && issample (r{1}) && issample (r{2}))
+ w = __frequency_vector__ (sys_cell, wbounds, r{1}, r{2});
+ else
+ error ("frequency_response: invalid cell");
endif
- w = get (sys, "w");
- H = __freqresp__ (sys, [], resptype, cellflag);
- elseif (isempty (w)) # find interesting frequency range w if not specified
- w = __frequency_vector__ (sys, wbounds);
- H = __freqresp__ (sys, w, resptype, cellflag);
- elseif (iscell (w) && numel (w) == 2 && issample (w{1}) && issample (w{2}))
- w = __frequency_vector__ (sys, wbounds, w{1}, w{2});
- H = __freqresp__ (sys, w, resptype, cellflag);
- elseif (! is_real_vector (w))
- error ("frequency_response: second argument w must be a vector of frequencies");
- else
- H = __freqresp__ (sys, w, resptype, cellflag);
+ elseif (any (w_idx)) # are there any frequency vectors?
+ w = args(w_idx){end};
+ else # there are neither frequency ranges nor vectors
+ w = __frequency_vector__ (sys_cell, wbounds);
endif
+ w_frd = w(frd_idx); # temporarily save frequency vectors of FRD models
+ w(frd_idx) = {[]}; # freqresp returns all frequencies of FRD models for w=[]
+
+ ## compute frequency response H for all LTI models
+ H = cellfun (@__freqresp__, sys_cell, w, {resptype}, {cellflag}, "uniformoutput", false);
+
+ ## restore frequency vectors of FRD models in w
+ w(frd_idx) = w_frd;
+
endfunction
Modified: trunk/octave-forge/main/control/inst/__frequency_vector__.m
===================================================================
--- trunk/octave-forge/main/control/inst/__frequency_vector__.m 2012-09-14 11:33:41 UTC (rev 11010)
+++ trunk/octave-forge/main/control/inst/__frequency_vector__.m 2012-09-14 11:40:39 UTC (rev 11011)
@@ -1,7 +1,7 @@
## Copyright (C) 1996, 2000, 2004, 2005, 2006, 2007
## Auburn University. All rights reserved.
+## Copyright (C) 2009 - 2012 Lukas F. Reichlin
##
-##
## 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
@@ -35,14 +35,85 @@
## Adapted-By: Lukas Reichlin <luk...@gm...>
## Date: October 2009
-## Version: 0.3
+## Version: 0.4
-function w = __frequency_vector__ (sys, wbounds = "std", wmin, wmax)
+function w = __frequency_vector__ (sys_cell, wbounds = "std", wmin, wmax)
+ isc = iscell (sys_cell);
+
+ if (! isc) # __sys2frd__ methods pass LTI models not in cells
+ sys_cell = {sys_cell}
+ endif
+
+ idx = cellfun (@(x) isa (x, "lti"), sys_cell);
+ sys_cell = sys_cell(idx);
+ len = numel (sys_cell);
+
+ [dec_min, dec_max, zp] = cellfun (@__frequency_range__, sys_cell, {wbounds}, "uniformoutput", false);
+
+ if (strcmpi (wbounds, "std")) # plots with explicit frequencies
+
+ if (nargin == 2)
+ dec_min = min (cell2mat (dec_min));
+ dec_max = max (cell2mat (dec_max));
+ elseif (nargin == 4) # w = {wmin, wmax}
+ dec_min = log10 (wmin);
+ dec_max = log10 (wmax);
+ else
+ print_usage ();
+ endif
+
+ zp = horzcat (zp{:});
+
+ ## include zeros and poles for nice peaks in plots
+ 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
+
+ w = repmat ({w}, 1, len); # return cell of frequency vectors
+
+ elseif (strcmpi (wbounds, "ext")) # plots with implicit frequencies
+
+ if (nargin == 4)
+ dec_min = repmat ({log10 (wmin)}, 1, len);
+ dec_max = repmat ({log10 (wmax)}, 1, len);
+ endif
+
+ idx = cellfun (@(zp, dec_min, dec_max) find (zp > 10^dec_min & zp < 10^dec_max), \
+ zp, dec_min, dec_max, "uniformoutput", false);
+ zp = cellfun (@(zp, idx) zp(idx), zp, idx, "uniformoutput", false);
+
+ w = cellfun (@logspace, dec_min, dec_max, {500}, "uniformoutput", false);
+ w = cellfun (@(w, zp) unique ([w, zp]), w, zp, "uniformoutput", false);
+ ## unique also sorts frequency vector
+
+ else
+ error ("frequency_vector: invalid argument 'wbounds'");
+ endif
+
+ if (! isc) # for __sys2frd__ methods
+ w = w{1};
+ endif
+
+endfunction
+
+
+function [dec_min, dec_max, zp] = __frequency_range__ (sys, wbounds = "std")
+
+ if (isa (sys, "frd"))
+ w = get (sys, "w");
+ dec_min = log10 (w(1));
+ dec_max = log10 (w(end));
+ zp = [];
+ return;
+ endif
+
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
+ 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, 1, []);
@@ -85,8 +156,8 @@
if (isempty (iip) && isempty (iiz))
## no poles/zeros away from omega = 0; pick defaults
- dec_min = 0; # -1
- dec_max = 2; # 3
+ 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]))));
@@ -94,7 +165,7 @@
## expand to show the entirety of the "interesting" portion of the plot
switch (wbounds)
- case "std" # standard
+ case "std" # standard
if (dec_min == dec_max)
dec_min -= 2;
dec_max += 2;
@@ -102,8 +173,8 @@
dec_min--;
dec_max++;
endif
- case "ext" # extended (for nyquist)
- if (any (abs (pol) < sqrt (eps))) # look for integrators
+ case "ext" # extended (for nyquist)
+ if (any (abs (pol) < sqrt (eps))) # look for integrators
## dec_min -= 0.5;
dec_max += 2;
else
@@ -119,17 +190,7 @@
dec_max = log10 (pi/tsam);
endif
- if (nargin == 4) # w = {wmin, wmax}
- dec_min = log10 (wmin);
- dec_max = log10 (wmax);
- endif
-
- ## create frequency vector
+ ## include zeros and poles for nice peaks in plots
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
Modified: trunk/octave-forge/main/control/inst/bode.m
===================================================================
--- trunk/octave-forge/main/control/inst/bode.m 2012-09-14 11:33:41 UTC (rev 11010)
+++ trunk/octave-forge/main/control/inst/bode.m 2012-09-14 11:40:39 UTC (rev 11011)
@@ -48,50 +48,65 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: November 2009
-## Version: 0.4
+## Version: 0.5
-function [mag_r, pha_r, w_r] = bode (sys, w = [])
+function [mag_r, pha_r, w_r] = bode (varargin)
- ## TODO: multiplot feature: bode (sys1, "b", sys2, "r", ...)
-
- if (nargin == 0 || nargin > 2)
+ if (nargin == 0)
print_usage ();
endif
- [H, w] = __frequency_response__ (sys, w, false, 0, "std");
+ [H, w] = __frequency_response__ (varargin, false, 0, "std", false);
+
+ H = cellfun (@reshape, H, {[]}, {1}, "uniformoutput", false);
+ mag = cellfun (@abs, H, "uniformoutput", false);
+ pha = cellfun (@(H) unwrap (arg (H)) * 180 / pi, H, "uniformoutput", false);
- H = reshape (H, [], 1);
- mag = abs (H);
- pha = unwrap (arg (H)) * 180 / pi;
-
if (! nargout)
- mag_db = 20 * log10 (mag);
+ mag_db = cellfun (@(mag) 20 * log10 (mag), mag, "uniformoutput", false);
- if (isct (sys))
- xl_str = "Frequency [rad/s]";
- else
- xl_str = sprintf ("Frequency [rad/s] w_N = %g", pi / get (sys, "tsam"));
- endif
+ tmp = cellfun (@isa, varargin, {"lti"});
+ sys_idx = find (tmp);
+ tmp = cellfun (@ischar, varargin);
+ style_idx = find (tmp);
+ len = numel (H);
+ mag_args = {};
+ pha_args = {};
+ legend_args = cell (len, 1);
+
+ for k = 1:len
+ if (k == len)
+ lim = nargin;
+ else
+ lim = sys_idx(k+1);
+ endif
+ style = varargin(style_idx(style_idx > sys_idx(k) & style_idx <= lim));
+ mag_args = cat (2, mag_args, w(k), mag_db(k), style);
+ pha_args = cat (2, pha_args, w(k), pha(k), style);
+ legend_args{k} = inputname(sys_idx(k));
+ endfor
+
subplot (2, 1, 1)
- semilogx (w, mag_db)
+ semilogx (mag_args{:})
axis ("tight")
ylim (__axis_margin__ (ylim))
grid ("on")
- title (["Bode Diagram of ", inputname(1)])
+ title ("Bode Diagram")
ylabel ("Magnitude [dB]")
subplot (2, 1, 2)
- semilogx (w, pha)
+ semilogx (pha_args{:})
axis ("tight")
ylim (__axis_margin__ (ylim))
grid ("on")
- xlabel (xl_str)
+ xlabel ("Frequency [rad/s]")
ylabel ("Phase [deg]")
+ legend (legend_args)
else
- mag_r = mag;
- pha_r = pha;
- w_r = w;
+ mag_r = mag{1};
+ pha_r = pha{1};
+ w_r = w{1};
endif
-endfunction
\ No newline at end of file
+endfunction
Modified: trunk/octave-forge/main/control/inst/nyquist.m
===================================================================
--- trunk/octave-forge/main/control/inst/nyquist.m 2012-09-14 11:33:41 UTC (rev 11010)
+++ trunk/octave-forge/main/control/inst/nyquist.m 2012-09-14 11:40:39 UTC (rev 11011)
@@ -48,35 +48,68 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: November 2009
-## Version: 0.3
+## Version: 0.4
-function [re_r, im_r, w_r] = nyquist (sys, w = [])
+function [re_r, im_r, w_r] = nyquist2 (varargin)
- ## TODO: multiplot feature: nyquist (sys1, "b", sys2, "r", ...)
-
- if (nargin == 0 || nargin > 2)
+ if (nargin == 0)
print_usage ();
endif
- [H, w] = __frequency_response__ (sys, w, false, 0, "ext");
+ [H, w] = __frequency_response__ (varargin, false, 0, "ext");
- H = reshape (H, [], 1);
- re = real (H);
- im = imag (H);
+ H = cellfun (@reshape, H, {[]}, {1}, "uniformoutput", false);
+ re = cellfun (@real, H, "uniformoutput", false);
+ im = cellfun (@imag, H, "uniformoutput", false);
if (! nargout)
- plot (re, im, "b", re, -im, "r")
+ tmp = cellfun (@isa, varargin, {"lti"});
+ sys_idx = find (tmp);
+ tmp = cellfun (@ischar, varargin);
+ style_idx = find (tmp);
+
+ len = numel (H);
+ pos_args = {};
+ neg_args = {};
+ legend_args = cell (len, 1);
+ colororder = get (gca, "colororder");
+ rc = rows (colororder);
+
+ for k = 1:len
+ col = colororder(1+rem (k-1, rc), :);
+ if (k == len)
+ lim = nargin;
+ else
+ lim = sys_idx(k+1);
+ endif
+ style = varargin(style_idx(style_idx > sys_idx(k) & style_idx <= lim));
+ if (isempty (style))
+ pos_args = cat (2, pos_args, re{k}, im{k}, {"-", "color", col});
+ neg_args = cat (2, neg_args, re{k}, -im{k}, {"-.", "color", col});
+ else
+ pos_args = cat (2, pos_args, re{k}, im{k}, style);
+ neg_args = cat (2, neg_args, re{k}, -im{k}, style);
+ endif
+ legend_args{k} = inputname(sys_idx(k));
+ endfor
+
+ ## FIXME: pos_args = cat (2, pos_args, re{k}, im{k}, {"-", "color", col}, style);
+ ## doesn't work! it would be nice to have default arguments that can be
+ ## (partially) overwritten by user-specified plot styles.
+
+ h = plot (pos_args{:}, neg_args{:});
axis ("tight")
xlim (__axis_margin__ (xlim))
ylim (__axis_margin__ (ylim))
grid ("on")
- title (["Nyquist Diagram of ", inputname(1)])
+ title ("Nyquist Diagram")
xlabel ("Real Axis")
ylabel ("Imaginary Axis")
+ legend (h(1:len), legend_args)
else
- re_r = re;
- im_r = im;
- w_r = w;
+ re_r = re{1};
+ im_r = im{1};
+ w_r = w{1};
endif
-endfunction
\ No newline at end of file
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
