|
3437
|
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 |
|
7016
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
3437
|
9 ## |
|
7016
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
3437
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
3437
|
18 |
|
6945
|
19 ## Undocumented internal function. |
|
|
20 |
|
3437
|
21 ## -*- texinfo -*- |
|
3500
|
22 ## @deftypefn {Function File} {[@var{f}, @var{w}, @var{rsys}] =} __bodquist__ (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx}) |
|
5016
|
23 ## Used internally by @command{bode}, @command{nyquist}; compute system frequency response. |
|
3437
|
24 ## |
|
|
25 ## @strong{Inputs} |
|
|
26 ## @table @var |
|
|
27 ## @item sys |
|
|
28 ## input system structure |
|
|
29 ## @item w |
|
|
30 ## range of frequencies; empty if user wants default |
|
|
31 ## @item out_idx |
|
3462
|
32 ## @itemx in_idx |
|
|
33 ## names or indices of output/input signal names; empty if user wants all |
|
3437
|
34 ## @item rname |
|
3455
|
35 ## name of routine that called __bodquist__ ("bode", "nyquist", or "nichols") |
|
3437
|
36 ## @end table |
|
|
37 ## @strong{Outputs} |
|
|
38 ## @table @var |
|
|
39 ## @item w |
|
|
40 ## list of frequencies |
|
|
41 ## @item f |
|
|
42 ## frequency response of sys; @math{f(ii) = f(omega(ii))} |
|
3455
|
43 ## @item rsys |
|
|
44 ## system with selected inputs and outputs |
|
3437
|
45 ## @end table |
|
|
46 ## |
|
3455
|
47 ## @code{bode}, @code{nichols}, and @code{nyquist} share the same |
|
|
48 ## introduction, so the common parts are |
|
3438
|
49 ## in __bodquist__. It contains the part that finds the number of arguments, |
|
5016
|
50 ## determines whether or not the system is @acronym{SISO}, and computes the frequency |
|
3437
|
51 ## response. Only the way the response is plotted is different between the |
|
3455
|
52 ## these functions. |
|
3437
|
53 ## @end deftypefn |
|
|
54 |
|
3455
|
55 function [f, w, rsys] = __bodquist__ (sys, w, outputs, inputs, rname) |
|
3437
|
56 |
|
|
57 ## check number of input arguments given |
|
|
58 if (nargin != 5) |
|
6046
|
59 print_usage (); |
|
3437
|
60 endif |
|
|
61 |
|
|
62 ## check each argument to see if it's in the correct form |
|
4030
|
63 if (!isstruct(sys)) |
|
3437
|
64 error("sys must be a system data structure"); |
|
|
65 endif |
|
|
66 |
|
3438
|
67 ## let __freqresp__ determine w if it's not already given |
|
3437
|
68 USEW = freqchkw(w); |
|
|
69 |
|
|
70 ## get initial dimensions (revised below if sysprune is called) |
|
|
71 [nn,nz,mm,pp ] = sysdimensions(sys); |
|
|
72 |
|
|
73 ## check for an output vector and to see whether it`s correct |
|
|
74 if (!isempty(outputs)) |
|
|
75 if (isempty(inputs)) |
|
|
76 inputs = 1:mm; # use all inputs |
|
|
77 warning([rname,": outputs specified but not inputs"]); |
|
5443
|
78 elseif(is_signal_list(inputs) | ischar(inputs)) |
|
3455
|
79 inputs = sysidx(sys,"in",inputs); |
|
3437
|
80 endif |
|
5443
|
81 if(is_signal_list(outputs) | ischar(outputs)) |
|
3455
|
82 outputs = sysidx(sys,"out",outputs); |
|
|
83 end |
|
3437
|
84 sys = sysprune(sys,outputs,inputs); |
|
|
85 [nn,nz,mm,pp ] = sysdimensions(sys); |
|
|
86 endif |
|
|
87 |
|
|
88 ## for speed in computation, convert local copy of |
|
|
89 ## SISO state space systems to zero-pole form |
|
|
90 if( is_siso(sys) & strcmp( sysgettype(sys), "ss") ) |
|
|
91 [zer,pol,k,tsam,inname,outname] = sys2zp(sys); |
|
4771
|
92 sys = zp(zer,pol,k,tsam,inname,outname); |
|
3437
|
93 endif |
|
|
94 |
|
|
95 ## get system frequency response |
|
3438
|
96 [f, w] = __freqresp__ (sys, USEW, w); |
|
3437
|
97 |
|
|
98 phase = arg(f)*180.0/pi; |
|
|
99 |
|
|
100 if(!USEW) |
|
|
101 ## smooth plots |
|
|
102 pcnt = 5; # max number of refinement steps |
|
|
103 dphase = 5; # desired max change in phase |
|
|
104 dmag = 0.2; # desired max change in magnitude |
|
|
105 while(pcnt) |
|
|
106 pd = abs(diff(phase)); # phase variation |
|
4373
|
107 pdbig = find(pd > dphase); |
|
3437
|
108 |
|
|
109 lp = length(f); lp1 = lp-1; # relative variation |
|
|
110 fd = abs(diff(f)); |
|
|
111 fm = max(abs([f(1:lp1); f(2:lp)])); |
|
4373
|
112 fdbig = find(fd > fm/10); |
|
3437
|
113 |
|
|
114 bigpts = union(fdbig, pdbig); |
|
|
115 |
|
|
116 if(isempty(bigpts) ) |
|
|
117 pcnt = 0; |
|
|
118 else |
|
|
119 pcnt = pcnt - 1; |
|
|
120 wnew = []; |
|
|
121 crossover_points = find ( phase(1:lp1).*phase(2:lp) < 0); |
|
|
122 pd(crossover_points) = abs(359.99+dphase - pd(crossover_points)); |
|
|
123 np_pts = max(3,ceil(pd/dphase)+2); # phase points |
|
|
124 nm_pts = max(3,ceil(log(fd./fm)/log(dmag))+2); # magnitude points |
|
|
125 npts = min(5,max(np_pts, nm_pts)); |
|
|
126 |
|
|
127 w1 = log10(w(1:lp1)); |
|
|
128 w2 = log10(w(2:lp)); |
|
|
129 for ii=bigpts |
|
|
130 if(npts(ii)) |
|
|
131 wtmp = logspace(w1(ii),w2(ii),npts(ii)); |
|
|
132 wseg(ii,1:(npts(ii)-2)) = wtmp(2:(npts(ii)-1)); |
|
|
133 endif |
|
|
134 endfor |
|
|
135 wnew = vec(wseg)'; # make a row vector |
|
|
136 wnew = wnew(find(wnew != 0)); |
|
|
137 wnew = sort(wnew); |
|
|
138 wnew = create_set(wnew); |
|
|
139 if(isempty(wnew)) # all small crossovers |
|
|
140 pcnt = 0; |
|
|
141 else |
|
3438
|
142 ## get new freq resp points, combine with old, and sort. |
|
|
143 [fnew, wnew] = __freqresp__ (sys, 1, wnew); |
|
|
144 w = [w, wnew]; |
|
|
145 f = [f, fnew]; |
|
|
146 [w, idx] = sort (w); |
|
|
147 f = f (idx); |
|
3437
|
148 phase = arg(f)*180.0/pi; |
|
|
149 endif |
|
|
150 endif |
|
|
151 endwhile |
|
|
152 endif |
|
|
153 |
|
|
154 ## ensure unique frequency values |
|
|
155 [w,idx] = sort(w); |
|
|
156 f = f(idx); |
|
|
157 |
|
|
158 w_diff = diff(w); |
|
|
159 w_dup = find(w_diff == 0); |
|
|
160 w_idx = complement(w_dup,1:length(w)); |
|
|
161 w = w(w_idx); |
|
|
162 f = f(w_idx); |
|
|
163 |
|
3455
|
164 ## set rsys to pruned system |
|
|
165 rsys = sys; |
|
|
166 |
|
3437
|
167 endfunction |
