|
7017
|
1 ## Copyright (C) 1996, 1998, 2000, 2003, 2004, 2005, 2006, 2007 |
|
|
2 ## Auburn University. All rights reserved. |
|
3431
|
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. |
|
3431
|
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. |
|
3431
|
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/>. |
|
3431
|
19 |
|
|
20 ## -*- texinfo -*- |
|
3500
|
21 ## @deftypefn {Function File} {[@var{realp}, @var{imagp}, @var{w}] =} nyquist (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx}, @var{atol}) |
|
|
22 ## @deftypefnx {Function File} {} nyquist (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx}, @var{atol}) |
|
3431
|
23 ## Produce Nyquist plots of a system; if no output arguments are given, Nyquist |
|
|
24 ## plot is printed to the screen. |
|
|
25 ## |
|
|
26 ## Compute the frequency response of a system. |
|
5016
|
27 ## |
|
3431
|
28 ## @strong{Inputs} (pass as empty to get default values) |
|
|
29 ## @table @var |
|
|
30 ## @item sys |
|
|
31 ## system data structure (must be either purely continuous or discrete; |
|
5016
|
32 ## see @code{is_digital}) |
|
3431
|
33 ## @item w |
|
|
34 ## frequency values for evaluation. |
|
5016
|
35 ## If sys is continuous, then bode evaluates @math{G(@var{jw})}; |
|
|
36 ## if sys is discrete, then bode evaluates @math{G(exp(@var{jwT}))}, |
|
|
37 ## where @var{T} is the system sampling time. |
|
3431
|
38 ## @item default |
|
|
39 ## the default frequency range is selected as follows: (These |
|
5016
|
40 ## steps are @strong{not} performed if @var{w} is specified) |
|
3431
|
41 ## @enumerate |
|
5016
|
42 ## @item via routine @command{__bodquist__}, isolate all poles and zeros away from |
|
3431
|
43 ## @var{w}=0 (@var{jw}=0 or @math{exp(@var{jwT})=1}) and select the frequency |
|
|
44 ## range based on the breakpoint locations of the frequencies. |
|
|
45 ## @item if @var{sys} is discrete time, the frequency range is limited |
|
|
46 ## to @var{jwT} in |
|
|
47 ## @ifinfo |
|
|
48 ## [0,2p*pi] |
|
|
49 ## @end ifinfo |
|
|
50 ## @iftex |
|
5016
|
51 ## @tex |
|
|
52 ## $ [ 0,2 p \pi ] $ |
|
|
53 ## @end tex |
|
3431
|
54 ## @end iftex |
|
5016
|
55 ## @item A ``smoothing'' routine is used to ensure that the plot phase does |
|
3431
|
56 ## not change excessively from point to point and that singular |
|
|
57 ## points (e.g., crossovers from +/- 180) are accurately shown. |
|
|
58 ## @end enumerate |
|
|
59 ## @item atol |
|
|
60 ## for interactive nyquist plots: atol is a change-in-slope tolerance |
|
|
61 ## for the of asymptotes (default = 0; 1e-2 is a good choice). This allows |
|
|
62 ## the user to ``zoom in'' on portions of the Nyquist plot too small to be |
|
|
63 ## seen with large asymptotes. |
|
|
64 ## @end table |
|
|
65 ## @strong{Outputs} |
|
|
66 ## @table @var |
|
|
67 ## @item realp |
|
|
68 ## @itemx imagp |
|
|
69 ## the real and imaginary parts of the frequency response |
|
|
70 ## @math{G(jw)} or @math{G(exp(jwT))} at the selected frequency values. |
|
5016
|
71 ## @item w |
|
3431
|
72 ## the vector of frequency values used |
|
|
73 ## @end table |
|
|
74 ## |
|
|
75 ## If no output arguments are given, nyquist plots the results to the screen. |
|
|
76 ## If @var{atol} != 0 and asymptotes are detected then the user is asked |
|
|
77 ## interactively if they wish to zoom in (remove asymptotes) |
|
|
78 ## Descriptive labels are automatically placed. |
|
|
79 ## |
|
5016
|
80 ## Note: if the requested plot is for an @acronym{MIMO} system, a warning message is |
|
3431
|
81 ## presented; the returned information is of the magnitude |
|
5016
|
82 ## @iftex |
|
|
83 ## @tex |
|
|
84 ## $ \Vert G(jw) \Vert $ or $ \Vert G( {\rm exp}(jwT) \Vert $ |
|
|
85 ## @end tex |
|
|
86 ## @end iftex |
|
|
87 ## @ifinfo |
|
|
88 ## ||G(jw)|| or ||G(exp(jwT))|| |
|
|
89 ## @end ifinfo |
|
|
90 ## only; phase information is not computed. |
|
3431
|
91 ## @end deftypefn |
|
|
92 |
|
|
93 ## Author: R. Bruce Tenison <btenison@eng.auburn.edu> |
|
|
94 ## Created: July 13, 1994 |
|
|
95 ## A. S. Hodel July 1995 (adaptive frequency spacing, |
|
|
96 ## remove acura parameter, etc.) |
|
|
97 ## Revised by John Ingram July 1996 for system format |
|
|
98 |
|
|
99 function [realp, imagp, w] = nyquist (sys, w, outputs, inputs, atol) |
|
|
100 |
|
|
101 ## Both bode and nyquist share the same introduction, so the common |
|
3438
|
102 ## parts are in a file called __bodquist__.m. It contains the part that |
|
3431
|
103 ## finds the number of arguments, determines whether or not the system |
|
|
104 ## is SISO, andd computes the frequency response. Only the way the |
|
|
105 ## response is plotted is different between the two functions. |
|
|
106 |
|
|
107 ## check number of input arguments given |
|
6448
|
108 if (nargin < 1 || nargin > 5) |
|
6046
|
109 print_usage (); |
|
3431
|
110 endif |
|
6448
|
111 if (nargin < 2) |
|
3431
|
112 w = []; |
|
|
113 endif |
|
6448
|
114 if (nargin < 3) |
|
3431
|
115 outputs = []; |
|
|
116 endif |
|
6448
|
117 if (nargin < 4) |
|
3431
|
118 inputs = []; |
|
|
119 endif |
|
6448
|
120 if (nargin < 5) |
|
3431
|
121 atol = 0; |
|
6448
|
122 elseif (! (is_sample (atol) || atol == 0)) |
|
|
123 error ("nyquist: atol must be a nonnegative scalar") |
|
3431
|
124 endif |
|
|
125 |
|
3438
|
126 ## signal to __bodquist__ who's calling |
|
3431
|
127 |
|
3455
|
128 [f, w, sys] = __bodquist__ (sys, w, outputs, inputs, "nyquist"); |
|
3431
|
129 |
|
|
130 ## Get the real and imaginary part of f. |
|
6448
|
131 realp = real (f); |
|
|
132 imagp = imag (f); |
|
3431
|
133 |
|
|
134 ## No output arguments, then display plot, otherwise return data. |
|
|
135 if (nargout == 0) |
|
|
136 dnplot = 0; |
|
6448
|
137 while (! dnplot) |
|
|
138 plot (realp, imagp, "- ;+w;", realp, -imagp, "-@ ;-w;"); |
|
|
139 |
|
3431
|
140 grid ("on"); |
|
|
141 |
|
6448
|
142 if (is_digital (sys)) |
|
3431
|
143 tstr = " G(e^{jw}) "; |
|
|
144 else |
|
|
145 tstr = " G(jw) "; |
|
|
146 endif |
|
6448
|
147 xlabel (sprintf ("Re(%s)", tstr)); |
|
|
148 ylabel (sprintf ("Im(%s)", tstr)); |
|
3431
|
149 |
|
6448
|
150 [stn, inn, outn] = sysgetsignals (sys); |
|
|
151 if (is_siso (sys)) |
|
|
152 title (sprintf ("Nyquist plot from %s to %s, w (rad/s) in [%e, %e]", |
|
|
153 inn{1}, outn{1}, w(1), w(end))); |
|
3431
|
154 endif |
|
|
155 |
|
6448
|
156 axis (axis2dlim ([[realp(:), imagp(:)]; [realp(:), -imagp(:)]])); |
|
3431
|
157 |
|
|
158 ## check for interactive plots |
|
|
159 dnplot = 1; # assume done; will change later if atol is satisfied |
|
6448
|
160 if (atol > 0 && length (f) > 2) |
|
3431
|
161 |
|
|
162 ## check for asymptotes |
|
6448
|
163 fmax = max (abs (f)); |
|
|
164 fi = find (abs (f) == fmax, 1, "last"); |
|
3431
|
165 |
|
|
166 ## compute angles from point to point |
|
6448
|
167 df = diff (f); |
|
|
168 th = atan2 (real (df), imag (df)) * 180 / pi; |
|
3431
|
169 |
|
|
170 ## get angle at fmax |
|
7126
|
171 if (fi == length (f)) |
|
6448
|
172 fi = fi-1; |
|
|
173 endif |
|
3431
|
174 thm = th(fi); |
|
|
175 |
|
|
176 ## now locate consecutive angles within atol of thm |
|
6448
|
177 ith_same = find (abs (th - thm) < atol); |
|
|
178 ichk = union (fi, find (diff (ith_same) == 1)); |
|
3431
|
179 |
|
|
180 ## locate max, min consecutive indices in ichk |
|
6448
|
181 loval = max (complement (ichk, 1:fi)); |
|
|
182 if (isempty (loval)) |
|
|
183 loval = fi; |
|
|
184 else |
|
|
185 loval = loval + 1; |
|
|
186 endif |
|
3431
|
187 |
|
6448
|
188 hival = min (complement (ichk, fi:length(th))); |
|
|
189 if (isempty (hival)) |
|
|
190 hival = fi+1; |
|
|
191 endif |
|
3431
|
192 |
|
6448
|
193 keep_idx = complement (loval:hival, 1:length(w)); |
|
3431
|
194 |
|
6448
|
195 if (length (keep_idx)) |
|
|
196 resp = input ("Remove asymptotes and zoom in (y or n): ", 1); |
|
|
197 if (resp(1) == "y") |
|
3431
|
198 dnplot = 0; # plot again |
|
|
199 w = w(keep_idx); |
|
|
200 f = f(keep_idx); |
|
6448
|
201 realp = real (f); |
|
|
202 imagp = imag (f); |
|
3431
|
203 endif |
|
|
204 endif |
|
|
205 |
|
|
206 endif |
|
6448
|
207 endwhile |
|
|
208 w = realp = imagp = []; |
|
|
209 endif |
|
3431
|
210 |
|
|
211 endfunction |
