3381
|
1 ## Copyright (C) 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. |
3213
|
18 |
3381
|
19 ## [mag,phase,w] = nichols(sys[,w,outputs,inputs]) |
|
20 ## Produce Nichols plot of a system |
|
21 ## |
|
22 ## Compute the frequency response of a system. |
|
23 ## inputs: |
|
24 ## sys: system data structure (must be either purely continuous or discrete; |
3213
|
25 # see is_digital) |
3381
|
26 ## w: frequency values for evaluation. |
|
27 ## if sys is continuous, then nichols evaluates G(jw) |
|
28 ## if sys is discrete, then nichols evaluates G(exp(jwT)), where T=sys.tsam |
|
29 ## (the system sampling time) |
|
30 ## default: the default frequency range is selected as follows: (These |
|
31 ## steps are NOT performed if w is specified) |
|
32 ## (1) via routine bodquist, isolate all poles and zeros away from |
|
33 ## w=0 (jw=0 or exp(jwT)=1) and select the frequency |
|
34 ## range based on the breakpoint locations of the frequencies. |
|
35 ## (2) if sys is discrete time, the frequency range is limited |
|
36 ## to jwT in [0,2p*pi] |
|
37 ## (3) A "smoothing" routine is used to ensure that the plot phase does |
|
38 ## not change excessively from point to point and that singular |
|
39 ## points (e.g., crossovers from +/- 180) are accurately shown. |
|
40 ## outputs, inputs: the indices of the output(s) and input(s) to be used in |
|
41 ## the frequency response; see sysprune. |
|
42 ## outputs: |
|
43 ## mag, phase: the magnitude and phase of the frequency response |
|
44 ## G(jw) or G(exp(jwT)) at the selected frequency values. |
|
45 ## w: the vector of frequency values used |
|
46 ## If no output arguments are given, nichols plots the results to the screen. |
|
47 ## Descriptive labels are automatically placed. See xlabel, ylable, title, |
|
48 ## and replot. |
|
49 ## |
|
50 ## Note: if the requested plot is for an MIMO system, mag is set to |
|
51 ## ||G(jw)|| or ||G(exp(jwT))|| and phase information is not computed. |
3213
|
52 |
3381
|
53 function [mag,phase,w] = nichols(sys,w,outputs,inputs) |
|
54 |
|
55 ## check number of input arguments given |
3213
|
56 if (nargin < 1 | nargin > 4) |
|
57 usage("[mag,phase,w] = nichols(sys[,w,outputs,inputs])"); |
|
58 endif |
|
59 if(nargin < 2) |
|
60 w = []; |
|
61 endif |
|
62 if(nargin < 3) |
|
63 outputs = []; |
|
64 endif |
|
65 if(nargin < 4) |
|
66 inputs = []; |
|
67 endif |
|
68 |
|
69 [f, w] = bodquist(sys,w,outputs,inputs,"nichols"); |
|
70 |
|
71 [stname,inname,outname] = sysgetsignals(sys); |
|
72 systsam = sysgettsam(sys); |
|
73 |
3381
|
74 ## Get the magnitude and phase of f. |
3213
|
75 mag = abs(f); |
|
76 phase = arg(f)*180.0/pi; |
|
77 |
|
78 if (nargout < 1), |
3381
|
79 ## Plot the information |
3213
|
80 if(gnuplot_has_multiplot) |
|
81 oneplot(); |
|
82 endif |
|
83 gset autoscale; |
|
84 if(gnuplot_has_multiplot) |
|
85 gset nokey; |
|
86 endif |
|
87 clearplot(); |
|
88 grid("on"); |
|
89 gset data style lines; |
|
90 if(is_digital(sys)) |
|
91 tistr = "(exp(jwT)) "; |
|
92 else |
|
93 tistr = "(jw)"; |
|
94 endif |
|
95 xlabel("Phase (deg)"); |
|
96 if(is_siso(sys)) |
3229
|
97 title(["Nichols plot of |[Y/U]",tistr,"|, u=", ... |
|
98 sysgetsignals(sys,"in",1,1), ", y=",sysgetsignals(sys,"out",1,1)]); |
3213
|
99 else |
|
100 title([ "||Y(", tistr, ")/U(", tistr, ")||"]); |
3229
|
101 printf("MIMO plot from\n%s\nto\n%s\n",outlist(inname," "), ... |
|
102 outlist(outname," ")); |
3213
|
103 endif |
3236
|
104 if(max(mag) > 0) |
|
105 ylabel("Gain in dB"); |
|
106 md = 20*log10(mag); |
|
107 else |
|
108 ylabel("Gain |Y/U|") |
|
109 md = mag; |
|
110 endif |
|
111 |
3213
|
112 axvec = axis2dlim([vec(phase),vec(md)]); |
|
113 axis(axvec); |
|
114 plot(phase,md); |
|
115 mag = phase = w = []; |
|
116 endif |
|
117 endfunction |