1## Copyright (C) 1996, 2000, 2004, 2005, 2006, 2007
2##               Auburn University. All rights reserved.
3## Copyright (C) 2009-2016   Lukas F. Reichlin
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.4
39
40function w = __frequency_vector__ (sys_cell, wbounds = "std", wmin, wmax)
41
42  N = 1000;     # intervals within the w range
43  isc = iscell (sys_cell);
44
45  if (! isc)                                    # __sys2frd__ methods pass LTI models not in cells
46    sys_cell = {sys_cell};
47  endif
48
49  idx = cellfun (@(x) isa (x, "lti"), sys_cell);
50  sys_cell = sys_cell(idx);
51  len = numel (sys_cell);
52
53  [dec_min, dec_max, zp] = cellfun (@__frequency_range__, sys_cell, {wbounds}, "uniformoutput", false);
54
55  if (strcmpi (wbounds, "std"))                 # plots with explicit frequencies
56
57    if (nargin == 4)                        # w = {wmin, wmax}
58      dec_min = log10 (wmin);
59      dec_max = log10 (wmax);
60    else
61      dec_min = min (cell2mat (dec_min));
62      dec_max = max (cell2mat (dec_max));
63    endif
64
65    zp = horzcat (zp{:});
66
67    ## include zeros and poles for nice peaks in plots
68    idx = find (zp > 10^dec_min & zp < 10^dec_max);
69    zp = zp(idx);
70
71    w = logspace (dec_min, dec_max, N);
72    w = unique ([w, zp]);                       # unique also sorts frequency vector
73
74    w = repmat ({w}, 1, len);                   # return cell of frequency vectors
75
76  elseif (strcmpi (wbounds, "ext"))             # plots with implicit frequencies
77
78    if (nargin == 4)
79      dec_min = repmat ({log10(wmin)}, 1, len);
80      dec_max = repmat ({log10(wmax)}, 1, len);
81    endif
82
83    idx = cellfun (@(zp, dec_min, dec_max) find (zp > 10^dec_min & zp < 10^dec_max), ...
84                   zp, dec_min, dec_max, "uniformoutput", false);
85    zp = cellfun (@(zp, idx) zp(idx), zp, idx, "uniformoutput", false);
86
87    w = cellfun (@logspace, dec_min, dec_max, {N}, "uniformoutput", false);
88    w = cellfun (@(w, zp) unique ([w, zp]), w, zp, "uniformoutput", false);
89    ## unique also sorts frequency vector
90
91  else
92    error ("frequency_vector: invalid argument 'wbounds'");
93  endif
94
95  if (! isc)                                    # for __sys2frd__ methods
96    w = w{1};
97  endif
98
99endfunction
100
101
102function [dec_min, dec_max, zp] = __frequency_range__ (sys, wbounds = "std")
103
104  if (isa (sys, "frd"))
105    w = get (sys, "w");
106    dec_min = log10 (w(1));
107    dec_max = log10 (w(end));
108    zp = [];
109    return;
110  endif
111
112  zer = zero (sys);
113  pol = pole (sys);
114  tsam = abs (get (sys, "tsam"));               # tsam could be -1
115  discrete = ! isct (sys);                      # static gains (tsam = -2) are assumed continuous
116
117  ## make sure zer, pol are row vectors
118  pol = reshape (pol, 1, []);
119  zer = reshape (zer, 1, []);
120
121  ## check for natural frequencies away from omega = 0
122  if (discrete)
123    ## The 2nd conditions prevents log(0) in the next log command
124    iiz = find (abs(zer-1) > norm(zer)*eps && abs(zer) > norm(zer)*eps);
125    iip = find (abs(pol-1) > norm(pol)*eps && abs(pol) > norm(pol)*eps);
126
127    ## avoid dividing empty matrices, it would work but looks nasty
128    if (! isempty (iiz))
129      czer = log (zer(iiz))/tsam;
130    else
131      czer = [];
132    endif
133
134    if (! isempty (iip))
135      cpol = log (pol(iip))/tsam;
136    else
137      cpol = [];
138    endif
139  else
140    ## continuous
141    iip = find (abs(pol) > norm(pol)*eps);
142    iiz = find (abs(zer) > norm(zer)*eps);
143
144    if (! isempty (zer))
145      czer = zer(iiz);
146    else
147      czer = [];
148    endif
149    if (! isempty (pol))
150      cpol = pol(iip);
151    else
152      cpol = [];
153    endif
154  endif
155
156  if (isempty (iip) && isempty (iiz))
157    ## no poles/zeros away from omega = 0; pick defaults
158    dec_min = 0;                                # -1
159    dec_max = 2;                                # 3
160  else
161    dec_min = floor (log10 (min (abs ([cpol, czer]))));
162    dec_max = ceil (log10 (max (abs ([cpol, czer]))));
163  endif
164
165  ## expand to show the entirety of the "interesting" portion of the plot
166  switch (wbounds)
167    case "std"                                  # standard
168      if (dec_min == dec_max)
169        dec_min -= 2;
170        dec_max += 2;
171      else
172        dec_min--;
173        dec_max++;
174      endif
175    case "ext"                                  # extended (for nyquist)
176      if (any (abs (pol) < sqrt (eps)))         # look for integrators
177        dec_min -= 0.5;
178        dec_max += 2;
179      else
180        dec_min -= 2;
181        dec_max += 2;
182      endif
183    otherwise
184      error ("frequency_range: second argument invalid");
185  endswitch
186
187  ## run discrete frequency all the way to pi
188  if (discrete)
189    dec_max = log10 (pi/tsam);
190  endif
191
192  ## include zeros and poles for nice peaks in plots
193  zp = [abs(zer), abs(pol)];
194
195endfunction
196