|
3279
|
1 # Copyright (C) 1996,1998 Auburn University. All Rights Reserved |
|
3213
|
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, 675 Mass Ave, Cambridge, MA 02139, USA. |
|
|
18 |
|
3282
|
19 function [mag,phase,w] = bode(sys,w,outputs,inputs,plot_style) |
|
|
20 # [mag,phase,w] = bode(sys[,w,outputs,inputs,plot_style]) |
|
3213
|
21 # Produce Bode plots of a system |
|
|
22 # |
|
|
23 # Compute the frequency response of a system. |
|
|
24 # inputs: |
|
|
25 # sys: system data structure (must be either purely continuous or discrete; |
|
|
26 # see is_digital) |
|
|
27 # w: frequency values for evaluation. |
|
|
28 # if sys is continuous, then bode evaluates G(jw) |
|
|
29 # if sys is discrete, then bode evaluates G(exp(jwT)), where T=sys.tsam |
|
|
30 # (the system sampling time) |
|
|
31 # default: the default frequency range is selected as follows: (These |
|
|
32 # steps are NOT performed if w is specified) |
|
|
33 # (1) via routine bodquist, isolate all poles and zeros away from |
|
|
34 # w=0 (jw=0 or exp(jwT)=1) and select the frequency |
|
|
35 # range based on the breakpoint locations of the frequencies. |
|
|
36 # (2) if sys is discrete time, the frequency range is limited |
|
|
37 # to jwT in [0,2p*pi] |
|
|
38 # (3) A "smoothing" routine is used to ensure that the plot phase does |
|
|
39 # not change excessively from point to point and that singular |
|
|
40 # points (e.g., crossovers from +/- 180) are accurately shown. |
|
|
41 # outputs, inputs: the indices of the output(s) and input(s) to be used in |
|
|
42 # the frequency response; see sysprune. |
|
3282
|
43 # plot_style: An optional argument specifying the type of plot to |
|
|
44 # produce (if plotting is being done). Valid values are |
|
|
45 # "dB" or "mag". If omitted, "dB" is assumed. |
|
|
46 # |
|
3213
|
47 # outputs: |
|
|
48 # mag, phase: the magnitude and phase of the frequency response |
|
|
49 # G(jw) or G(exp(jwT)) at the selected frequency values. |
|
|
50 # w: the vector of frequency values used |
|
|
51 # If no output arguments are given, bode plots the results to the screen. |
|
|
52 # Descriptive labels are automatically placed. See xlabel, ylable, title, |
|
|
53 # and replot. |
|
|
54 # |
|
|
55 # Note: if the requested plot is for an MIMO system, mag is set to |
|
|
56 # ||G(jw)|| or ||G(exp(jwT))|| and phase information is not computed. |
|
|
57 |
|
|
58 # Written by John Ingram July 10th, 1996 |
|
|
59 # Based on previous code |
|
|
60 # By R. Bruce Tenison, July 13, 1994 |
|
|
61 # Modified by David Clem November 13, 1994 |
|
|
62 # again by A. S. Hodel July 1995 (smart plot range, etc.) |
|
|
63 # Modified by Kai P. Mueller September 28, 1997 (multiplot mode) |
|
|
64 |
|
|
65 # check number of input arguments given |
|
3282
|
66 if (nargin < 1 | nargin > 5) |
|
|
67 usage("[mag,phase,w] = bode(sys[,w,outputs,inputs,plot_style])"); |
|
3213
|
68 endif |
|
|
69 if(nargin < 2) |
|
|
70 w = []; |
|
|
71 endif |
|
|
72 if(nargin < 3) |
|
|
73 outputs = []; |
|
|
74 endif |
|
|
75 if(nargin < 4) |
|
|
76 inputs = []; |
|
|
77 endif |
|
3282
|
78 if(nargin < 5) |
|
|
79 plot_style = "dB"; |
|
|
80 endif |
|
|
81 |
|
|
82 if (strcmp (plot_style, "dB")) |
|
|
83 do_db_plot = 1; |
|
|
84 elseif (strcmp (plot_style, "mag")) |
|
|
85 do_db_plot = 0; |
|
|
86 else |
|
|
87 error ("bode: invalid value of plot_style specified"); |
|
|
88 endif |
|
3213
|
89 |
|
|
90 [f, w] = bodquist(sys,w,outputs,inputs,"bode"); |
|
|
91 |
|
|
92 [stname,inname,outname] = sysgetsignals(sys); |
|
|
93 systsam = sysgettsam(sys); |
|
|
94 |
|
|
95 # Get the magnitude and phase of f. |
|
|
96 mag = abs(f); |
|
|
97 phase = arg(f)*180.0/pi; |
|
|
98 |
|
|
99 if (nargout < 1), |
|
|
100 # Plot the information |
|
|
101 if(gnuplot_has_multiplot) |
|
|
102 oneplot(); |
|
|
103 endif |
|
|
104 gset autoscale; |
|
|
105 if(gnuplot_has_multiplot) |
|
|
106 gset nokey; |
|
|
107 endif |
|
|
108 clearplot(); |
|
|
109 gset data style lines; |
|
|
110 if(is_digital(sys)) |
|
|
111 xlstr = ["Digital frequency w=rad/sec. pi/T=",num2str(pi/systsam)]; |
|
|
112 tistr = "(exp(jwT)) "; |
|
|
113 else |
|
|
114 xlstr = "Frequency in rad/sec"; |
|
|
115 tistr = "(jw)"; |
|
|
116 endif |
|
|
117 xlabel(xlstr); |
|
|
118 if(is_siso(sys)) |
|
|
119 if (gnuplot_has_multiplot) |
|
|
120 subplot(2,1,1); |
|
|
121 endif |
|
3228
|
122 title(["|[Y/U]",tistr,"|, u=", nth(inname,1),", y=",nth(outname,1)]); |
|
3213
|
123 else |
|
|
124 title([ "||Y(", tistr, ")/U(", tistr, ")||"]); |
|
|
125 disp("MIMO plot from") |
|
3228
|
126 disp(outlist(inname," ")); |
|
3213
|
127 disp("to") |
|
3228
|
128 disp(outlist(outname," ")); |
|
3213
|
129 endif |
|
|
130 wv = [min(w), max(w)]; |
|
3282
|
131 if(do_db_plot && max(mag) > 0) |
|
3236
|
132 ylabel("Gain in dB"); |
|
|
133 md = 20*log10(mag); |
|
|
134 else |
|
|
135 ylabel("Gain |Y/U|") |
|
|
136 md = mag; |
|
|
137 endif |
|
3213
|
138 |
|
|
139 axvec = axis2dlim([vec(w),vec(md)]); |
|
|
140 axvec(1:2) = wv; |
|
|
141 axis(axvec); |
|
|
142 grid("on"); |
|
3282
|
143 if (do_db_plot) |
|
|
144 semilogx(w,md); |
|
|
145 else |
|
|
146 loglog(w,md); |
|
|
147 endif |
|
3213
|
148 if (is_siso(sys)) |
|
|
149 if (gnuplot_has_multiplot) |
|
|
150 subplot(2,1,2); |
|
|
151 else |
|
|
152 prompt('Press any key for phase plot'); |
|
|
153 endif |
|
|
154 axvec = axis2dlim([vec(w),vec(phase)]); |
|
|
155 axvec(1:2) = wv; |
|
|
156 axis(axvec); |
|
|
157 xlabel(xlstr); |
|
|
158 ylabel("Phase in deg"); |
|
|
159 title([ "phase([Y/U]", tistr, ... |
|
3228
|
160 "), u=", nth(inname,1),", y=",nth(outname,1)]); |
|
3213
|
161 grid("on"); |
|
|
162 semilogx(w,phase); |
|
|
163 # This should be the default for subsequent plot commands. |
|
|
164 if(gnuplot_has_multiplot) |
|
|
165 oneplot(); |
|
|
166 endif |
|
|
167 endif |
|
|
168 mag = phase = w = []; |
|
|
169 endif |
|
|
170 endfunction |
