|
7017
|
1 ## Copyright (C) 1996, 1998, 2000, 2002, 2003, 2004, 2005, 2006, 2007 |
|
|
2 ## Auburn University. All rights reserved. |
|
3437
|
3 ## |
|
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
7016
|
7 ## under the terms of the GNU General Public License as published by |
|
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
9 ## your option) any later version. |
|
3437
|
10 ## |
|
7016
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
14 ## General Public License for more details. |
|
3437
|
15 ## |
|
|
16 ## You should have received a copy of the GNU General Public License |
|
7016
|
17 ## along with Octave; see the file COPYING. If not, see |
|
|
18 ## <http://www.gnu.org/licenses/>. |
|
3437
|
19 |
|
6945
|
20 ## Undocumented internal function. |
|
|
21 |
|
3437
|
22 ## -*- texinfo -*- |
|
3500
|
23 ## @deftypefn {Function File} {[@var{f}, @var{w}, @var{rsys}] =} __bodquist__ (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx}) |
|
5016
|
24 ## Used internally by @command{bode}, @command{nyquist}; compute system frequency response. |
|
3437
|
25 ## |
|
|
26 ## @strong{Inputs} |
|
|
27 ## @table @var |
|
|
28 ## @item sys |
|
|
29 ## input system structure |
|
|
30 ## @item w |
|
|
31 ## range of frequencies; empty if user wants default |
|
|
32 ## @item out_idx |
|
3462
|
33 ## @itemx in_idx |
|
|
34 ## names or indices of output/input signal names; empty if user wants all |
|
3437
|
35 ## @item rname |
|
3455
|
36 ## name of routine that called __bodquist__ ("bode", "nyquist", or "nichols") |
|
3437
|
37 ## @end table |
|
|
38 ## @strong{Outputs} |
|
|
39 ## @table @var |
|
|
40 ## @item w |
|
|
41 ## list of frequencies |
|
|
42 ## @item f |
|
|
43 ## frequency response of sys; @math{f(ii) = f(omega(ii))} |
|
3455
|
44 ## @item rsys |
|
|
45 ## system with selected inputs and outputs |
|
3437
|
46 ## @end table |
|
|
47 ## |
|
3455
|
48 ## @code{bode}, @code{nichols}, and @code{nyquist} share the same |
|
|
49 ## introduction, so the common parts are |
|
3438
|
50 ## in __bodquist__. It contains the part that finds the number of arguments, |
|
5016
|
51 ## determines whether or not the system is @acronym{SISO}, and computes the frequency |
|
3437
|
52 ## response. Only the way the response is plotted is different between the |
|
3455
|
53 ## these functions. |
|
3437
|
54 ## @end deftypefn |
|
|
55 |
|
3455
|
56 function [f, w, rsys] = __bodquist__ (sys, w, outputs, inputs, rname) |
|
3437
|
57 |
|
|
58 ## check number of input arguments given |
|
|
59 if (nargin != 5) |
|
6046
|
60 print_usage (); |
|
3437
|
61 endif |
|
|
62 |
|
|
63 ## check each argument to see if it's in the correct form |
|
7126
|
64 if (! isstruct (sys)) |
|
|
65 error ("sys must be a system data structure"); |
|
3437
|
66 endif |
|
|
67 |
|
3438
|
68 ## let __freqresp__ determine w if it's not already given |
|
7126
|
69 USEW = freqchkw (w); |
|
3437
|
70 |
|
|
71 ## get initial dimensions (revised below if sysprune is called) |
|
7126
|
72 [nn, nz, mm, pp] = sysdimensions (sys); |
|
3437
|
73 |
|
|
74 ## check for an output vector and to see whether it`s correct |
|
7126
|
75 if (! isempty (outputs)) |
|
|
76 if (isempty (inputs)) |
|
3437
|
77 inputs = 1:mm; # use all inputs |
|
7126
|
78 warning ("%s: outputs specified but not inputs", rname); |
|
|
79 elseif (is_signal_list (inputs) || ischar (inputs)) |
|
|
80 inputs = sysidx (sys, "in", inputs); |
|
3437
|
81 endif |
|
7126
|
82 if (is_signal_list (outputs) || ischar (outputs)) |
|
|
83 outputs = sysidx (sys, "out", outputs); |
|
3455
|
84 end |
|
7126
|
85 sys = sysprune (sys, outputs, inputs); |
|
|
86 [nn, nz, mm, pp] = sysdimensions (sys); |
|
3437
|
87 endif |
|
|
88 |
|
|
89 ## for speed in computation, convert local copy of |
|
|
90 ## SISO state space systems to zero-pole form |
|
7126
|
91 if (is_siso (sys) && strcmp (sysgettype (sys), "ss")) |
|
|
92 [zer, pol, k, tsam, inname, outname] = sys2zp (sys); |
|
|
93 sys = zp (zer, pol, k, tsam, inname, outname); |
|
3437
|
94 endif |
|
|
95 |
|
|
96 ## get system frequency response |
|
3438
|
97 [f, w] = __freqresp__ (sys, USEW, w); |
|
3437
|
98 |
|
|
99 phase = arg(f)*180.0/pi; |
|
|
100 |
|
7126
|
101 if (! USEW) |
|
3437
|
102 ## smooth plots |
|
|
103 pcnt = 5; # max number of refinement steps |
|
|
104 dphase = 5; # desired max change in phase |
|
|
105 dmag = 0.2; # desired max change in magnitude |
|
7126
|
106 while (pcnt) |
|
|
107 pd = abs (diff (phase)); # phase variation |
|
|
108 pdbig = find (pd > dphase); |
|
|
109 |
|
|
110 ## relative variation |
|
|
111 lp = length (f); |
|
|
112 lp1 = lp-1; |
|
3437
|
113 |
|
7126
|
114 fd = abs (diff (f)); |
|
|
115 fm = max (abs ([f(1:lp1); f(2:lp)])); |
|
|
116 fdbig = find (fd > fm/10); |
|
3437
|
117 |
|
7126
|
118 bigpts = union (fdbig, pdbig); |
|
3437
|
119 |
|
7126
|
120 if (isempty (bigpts)) |
|
3437
|
121 pcnt = 0; |
|
|
122 else |
|
|
123 pcnt = pcnt - 1; |
|
|
124 wnew = []; |
|
7126
|
125 crossover_points = find (phase(1:lp1).*phase(2:lp) < 0); |
|
|
126 pd(crossover_points) = abs (359.99+dphase - pd(crossover_points)); |
|
|
127 np_pts = max (3, ceil(pd/dphase)+2); # phase points |
|
|
128 nm_pts = max (3, ceil(log(fd./fm)/log(dmag))+2); # magnitude points |
|
|
129 npts = min (5, max(np_pts, nm_pts)); |
|
3437
|
130 |
|
7126
|
131 w1 = log10 (w(1:lp1)); |
|
|
132 w2 = log10 (w(2:lp)); |
|
|
133 for ii = bigpts |
|
|
134 if (npts(ii)) |
|
|
135 wtmp = logspace (w1(ii), w2(ii), npts(ii)); |
|
3437
|
136 wseg(ii,1:(npts(ii)-2)) = wtmp(2:(npts(ii)-1)); |
|
|
137 endif |
|
|
138 endfor |
|
|
139 wnew = vec(wseg)'; # make a row vector |
|
7126
|
140 wnew = wnew(find (wnew != 0)); |
|
|
141 wnew = sort (wnew); |
|
|
142 wnew = create_set (wnew); |
|
|
143 if (isempty (wnew)) # all small crossovers |
|
3437
|
144 pcnt = 0; |
|
|
145 else |
|
3438
|
146 ## get new freq resp points, combine with old, and sort. |
|
|
147 [fnew, wnew] = __freqresp__ (sys, 1, wnew); |
|
|
148 w = [w, wnew]; |
|
|
149 f = [f, fnew]; |
|
|
150 [w, idx] = sort (w); |
|
|
151 f = f (idx); |
|
3437
|
152 phase = arg(f)*180.0/pi; |
|
|
153 endif |
|
|
154 endif |
|
|
155 endwhile |
|
|
156 endif |
|
|
157 |
|
|
158 ## ensure unique frequency values |
|
7126
|
159 [w, idx] = sort (w); |
|
3437
|
160 f = f(idx); |
|
|
161 |
|
7126
|
162 w_diff = diff (w); |
|
|
163 w_dup = find (w_diff == 0); |
|
|
164 w_idx = complement (w_dup, 1:length(w)); |
|
3437
|
165 w = w(w_idx); |
|
|
166 f = f(w_idx); |
|
|
167 |
|
3455
|
168 ## set rsys to pruned system |
|
|
169 rsys = sys; |
|
|
170 |
|
3437
|
171 endfunction |
