Mercurial > hg > octave-nkf

comparison scripts/control/bodquist.m @ 3425:8625164a0a39

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2000-01-13 08:31:37 by jwe]
author jwe
date Thu, 13 Jan 2000 08:32:16 +0000
parents 1a8e2c0d627a
children
comparison
equal deleted inserted replaced
3424:61e40232a4e8 3425:8625164a0a39
1 ## Copyright (C) 1996, 1998 Auburn University. All rights reserved. 1 ## Copyright (C) 1996, 1998 Auburn University. All rights reserved.
2 ## 2 ##
3 ## This file is part of Octave. 3 ## This file is part of Octave.
4 ## 4 ##
5 ## Octave is free software; you can redistribute it and/or modify it 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 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 7 ## Free Software Foundation; either version 2, or (at your option) any
8 ## later version. 8 ## later version.
9 ## 9 ##
10 ## Octave is distributed in the hope that it will be useful, but WITHOUT 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 11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 ## for more details. 13 ## for more details.
14 ## 14 ##
15 ## You should have received a copy of the GNU General Public License 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 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. 17 ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File } { [@var{f}, @var{w}] =} bodquist (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx}) 20 ## @deftypefn {Function File} {[@var{f}, @var{w}] =} bodquist (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx})
21 ## used internally by bode, nyquist; compute system frequency response. 21 ## used internally by bode, nyquist; compute system frequency response.
22 ## 22 ##
23 ## @strong{Inputs} 23 ## @strong{Inputs}
24 ## @table @var 24 ## @table @var
25 ## @item sys 25 ## @item sys
26 ## input system structure 26 ## input system structure
27 ## @item w 27 ## @item w
34 ## name of routine that called bodquist ("bode" or "nyquist") 34 ## name of routine that called bodquist ("bode" or "nyquist")
35 ## @end table 35 ## @end table
36 ## @strong{Outputs} 36 ## @strong{Outputs}
37 ## @table @var 37 ## @table @var
38 ## @item w 38 ## @item w
39 ## list of frequencies 39 ## list of frequencies
40 ## @item f 40 ## @item f
41 ## frequency response of sys; @math{f(ii) = f(omega(ii))} 41 ## frequency response of sys; @math{f(ii) = f(omega(ii))}
42 ## @end table 42 ## @end table
43 ## @strong{Note} bodquist could easily be incorporated into a Nichols 43 ## @strong{Note} bodquist could easily be incorporated into a Nichols
44 ## plot function; this is in a "to do" list. 44 ## plot function; this is in a "to do" list.
45 ## 45 ##
46 ## Both bode and nyquist share the same introduction, so the common parts are 46 ## Both bode and nyquist share the same introduction, so the common parts are
47 ## in bodquist. It contains the part that finds the number of arguments, 47 ## in bodquist. It contains the part that finds the number of arguments,
48 ## determines whether or not the system is SISO, and computes the frequency 48 ## determines whether or not the system is SISO, and computes the frequency
49 ## response. Only the way the response is plotted is different between the 49 ## response. Only the way the response is plotted is different between the
50 ## two functions. 50 ## two functions.
51 ## @end deftypefn 51 ## @end deftypefn
52 52
53 function [f, w] = bodquist (sys, w, outputs, inputs, rname) 53 function [f, w] = bodquist (sys, w, outputs, inputs, rname)
54 54
55 ## check number of input arguments given 55 ## check number of input arguments given
56 if (nargin != 5) 56 if (nargin != 5)
57 usage("[f,w] = bodquist(sys,w,outputs,inputs,rname)"); 57 usage("[f,w] = bodquist(sys,w,outputs,inputs,rname)");
59 59
60 ## check each argument to see if it's in the correct form 60 ## check each argument to see if it's in the correct form
61 if (!is_struct(sys)) 61 if (!is_struct(sys))
62 error("sys must be a system data structure"); 62 error("sys must be a system data structure");
63 endif 63 endif
64 64
65 ## let freqresp determine w if it's not already given 65 ## let freqresp determine w if it's not already given
66 USEW = freqchkw(w); 66 USEW = freqchkw(w);
67 67
68 ## get initial dimensions (revised below if sysprune is called) 68 ## get initial dimensions (revised below if sysprune is called)
69 [nn,nz,mm,pp ] = sysdimensions(sys); 69 [nn,nz,mm,pp ] = sysdimensions(sys);
70 70
71 ## check for an output vector and to see whether it`s correct 71 ## check for an output vector and to see whether it`s correct
72 if (!isempty(outputs)) 72 if (!isempty(outputs))
73 if (isempty(inputs)) 73 if (isempty(inputs))
74 inputs = 1:mm; # use all inputs 74 inputs = 1:mm; # use all inputs
75 warning([rname,": outputs specified but not inputs"]); 75 warning([rname,": outputs specified but not inputs"]);
76 endif 76 endif
77 sys = sysprune(sys,outputs,inputs); 77 sys = sysprune(sys,outputs,inputs);
78 [nn,nz,mm,pp ] = sysdimensions(sys); 78 [nn,nz,mm,pp ] = sysdimensions(sys);
79 endif 79 endif
80 80
81 ## for speed in computation, convert local copy of 81 ## for speed in computation, convert local copy of
82 ## SISO state space systems to zero-pole form 82 ## SISO state space systems to zero-pole form
83 if( is_siso(sys) & strcmp( sysgettype(sys), "ss") ) 83 if( is_siso(sys) & strcmp( sysgettype(sys), "ss") )
84 [zer,pol,k,tsam,inname,outname] = sys2zp(sys); 84 [zer,pol,k,tsam,inname,outname] = sys2zp(sys);
85 sys = zp2sys(zer,pol,k,tsam,inname,outname); 85 sys = zp2sys(zer,pol,k,tsam,inname,outname);
86 endif 86 endif
87 87
88 ## get system frequency response 88 ## get system frequency response
89 [f,w] = freqresp(sys,USEW,w); 89 [f,w] = freqresp(sys,USEW,w);
90 90
91 phase = arg(f)*180.0/pi; 91 phase = arg(f)*180.0/pi;
92 92
93 if(!USEW) 93 if(!USEW)
94 ## smooth plots 94 ## smooth plots
95 pcnt = 5; # max number of refinement steps 95 pcnt = 5; # max number of refinement steps
96 dphase = 5; # desired max change in phase 96 dphase = 5; # desired max change in phase
97 dmag = 0.2; # desired max change in magnitude 97 dmag = 0.2; # desired max change in magnitude
98 while(pcnt) 98 while(pcnt)
99 pd = abs(diff(phase)); # phase variation 99 pd = abs(diff(phase)); # phase variation
100 pdbig = vec(find(pd > dphase)); 100 pdbig = vec(find(pd > dphase));
101 101
102 lp = length(f); lp1 = lp-1; # relative variation 102 lp = length(f); lp1 = lp-1; # relative variation
103 fd = abs(diff(f)); 103 fd = abs(diff(f));
104 fm = max(abs([f(1:lp1); f(2:lp)])); 104 fm = max(abs([f(1:lp1); f(2:lp)]));
105 fdbig = vec(find(fd > fm/10)); 105 fdbig = vec(find(fd > fm/10));
106 106
107 bigpts = union(fdbig, pdbig); 107 bigpts = union(fdbig, pdbig);
111 else 111 else
112 pcnt = pcnt - 1; 112 pcnt = pcnt - 1;
113 wnew = []; 113 wnew = [];
114 crossover_points = find ( phase(1:lp1).*phase(2:lp) < 0); 114 crossover_points = find ( phase(1:lp1).*phase(2:lp) < 0);
115 pd(crossover_points) = abs(359.99+dphase - pd(crossover_points)); 115 pd(crossover_points) = abs(359.99+dphase - pd(crossover_points));
116 np_pts = max(3,ceil(pd/dphase)+2); # phase points 116 np_pts = max(3,ceil(pd/dphase)+2); # phase points
117 nm_pts = max(3,ceil(log(fd./fm)/log(dmag))+2); # magnitude points 117 nm_pts = max(3,ceil(log(fd./fm)/log(dmag))+2); # magnitude points
118 npts = min(5,max(np_pts, nm_pts)); 118 npts = min(5,max(np_pts, nm_pts));
119 119
120 w1 = log10(w(1:lp1)); 120 w1 = log10(w(1:lp1));
121 w2 = log10(w(2:lp)); 121 w2 = log10(w(2:lp));
122 for ii=bigpts 122 for ii=bigpts
131 wnew = create_set(wnew); 131 wnew = create_set(wnew);
132 if(isempty(wnew)) # all small crossovers 132 if(isempty(wnew)) # all small crossovers
133 pcnt = 0; 133 pcnt = 0;
134 else 134 else
135 [fnew,wnew] = freqresp(sys,1,wnew); # get new freq resp points 135 [fnew,wnew] = freqresp(sys,1,wnew); # get new freq resp points
136 w = [w,wnew]; # combine with old freq resp 136 w = [w,wnew]; # combine with old freq resp
137 f = [f,fnew]; 137 f = [f,fnew];
138 [w,idx] = sort(w); # sort into order 138 [w,idx] = sort(w); # sort into order
139 f = f(idx); 139 f = f(idx);
140 phase = arg(f)*180.0/pi; 140 phase = arg(f)*180.0/pi;
141 endif 141 endif
142 endif 142 endif
143 endwhile 143 endwhile
150 w_diff = diff(w); 150 w_diff = diff(w);
151 w_dup = find(w_diff == 0); 151 w_dup = find(w_diff == 0);
152 w_idx = complement(w_dup,1:length(w)); 152 w_idx = complement(w_dup,1:length(w));
153 w = w(w_idx); 153 w = w(w_idx);
154 f = f(w_idx); 154 f = f(w_idx);
155 155
156 endfunction 156 endfunction