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