|
3431
|
1 ## Copyright (C) 1996, 1998 Auburn University. All rights reserved. |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by the |
|
|
7 ## Free Software Foundation; either version 2, or (at your option) any |
|
|
8 ## later version. |
|
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
13 ## for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
5307
|
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
18 ## 02110-1301 USA. |
|
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 |
|
|
108 if (nargin < 1 | nargin > 5) |
|
|
109 usage("[realp,imagp,w] = nyquist(sys[,w,outputs,inputs,atol])"); |
|
|
110 endif |
|
|
111 if(nargin < 2) |
|
|
112 w = []; |
|
|
113 endif |
|
|
114 if(nargin < 3) |
|
|
115 outputs = []; |
|
|
116 endif |
|
|
117 if(nargin < 4) |
|
|
118 inputs = []; |
|
|
119 endif |
|
|
120 if(nargin < 5) |
|
|
121 atol = 0; |
|
|
122 elseif(!(is_sample(atol) | atol == 0)) |
|
|
123 error("atol must be a nonnegative scalar.") |
|
|
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. |
|
|
131 realp = real(f); |
|
|
132 imagp = imag(f); |
|
|
133 |
|
|
134 ## No output arguments, then display plot, otherwise return data. |
|
|
135 if (nargout == 0) |
|
|
136 dnplot = 0; |
|
|
137 while(!dnplot) |
|
4422
|
138 oneplot(); |
|
5215
|
139 __gnuplot_set__ key; |
|
3431
|
140 clearplot(); |
|
|
141 grid ("on"); |
|
5215
|
142 __gnuplot_set__ data style lines; |
|
3431
|
143 |
|
|
144 if(is_digital(sys)) |
|
|
145 tstr = " G(e^{jw}) "; |
|
|
146 else |
|
|
147 tstr = " G(jw) "; |
|
|
148 endif |
|
|
149 xlabel(["Re(",tstr,")"]); |
|
|
150 ylabel(["Im(",tstr,")"]); |
|
|
151 |
|
|
152 [stn, inn, outn] = sysgetsignals(sys); |
|
|
153 if(is_siso(sys)) |
|
|
154 title(sprintf("Nyquist plot from %s to %s, w (rad/s) in [%e, %e]", ... |
|
4771
|
155 inn{1}, outn{1}, w(1), w(length(w))) ) |
|
3431
|
156 endif |
|
|
157 |
|
5215
|
158 __gnuplot_set__ nologscale xy; |
|
3431
|
159 |
|
|
160 axis(axis2dlim([[vec(realp),vec(imagp)];[vec(realp),-vec(imagp)]])); |
|
|
161 plot(realp,imagp,"- ;+w;",realp,-imagp,"-@ ;-w;"); |
|
|
162 |
|
|
163 ## check for interactive plots |
|
|
164 dnplot = 1; # assume done; will change later if atol is satisfied |
|
|
165 if(atol > 0 & length(f) > 2) |
|
|
166 |
|
|
167 ## check for asymptotes |
|
|
168 fmax = max(abs(f)); |
|
|
169 fi = max(find(abs(f) == fmax)); |
|
|
170 |
|
|
171 ## compute angles from point to point |
|
|
172 df = diff(f); |
|
|
173 th = atan2(real(df),imag(df))*180/pi; |
|
|
174 |
|
|
175 ## get angle at fmax |
|
|
176 if(fi == length(f)) fi = fi-1; endif |
|
|
177 thm = th(fi); |
|
|
178 |
|
|
179 ## now locate consecutive angles within atol of thm |
|
|
180 ith_same = find(abs(th - thm) < atol); |
|
|
181 ichk = union(fi,find(diff(ith_same) == 1)); |
|
|
182 |
|
|
183 ## locate max, min consecutive indices in ichk |
|
|
184 loval = max(complement(ichk,1:fi)); |
|
|
185 if(isempty(loval)) loval = fi; |
|
|
186 else loval = loval + 1; endif |
|
|
187 |
|
|
188 hival = min(complement(ichk,fi:length(th))); |
|
|
189 if(isempty(hival)) hival = fi+1; endif |
|
|
190 |
|
|
191 keep_idx = complement(loval:hival,1:length(w)); |
|
|
192 |
|
|
193 if(length(keep_idx)) |
|
|
194 resp = input("Remove asymptotes and zoom in (y or n): ",1); |
|
|
195 if(resp(1) == "y") |
|
|
196 dnplot = 0; # plot again |
|
|
197 w = w(keep_idx); |
|
|
198 f = f(keep_idx); |
|
|
199 realp = real(f); |
|
|
200 imagp = imag(f); |
|
|
201 endif |
|
|
202 endif |
|
|
203 |
|
|
204 endif |
|
|
205 endwhile |
|
|
206 w = []; |
|
|
207 realp=[]; |
|
|
208 imagp=[]; |
|
|
209 endif |
|
|
210 |
|
|
211 endfunction |
