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