Mercurial > octave-nkf

comparison scripts/control/base/nyquist.m @ 6448:2110cc251779

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-03-24 02:47:36 by jwe]
author jwe
date Sat, 24 Mar 2007 02:47:36 +0000
parents 34f96dd5441b
children 93c65f2a5668
comparison
equal deleted inserted replaced
6447:3f79532415b5 6448:2110cc251779
103 ## finds the number of arguments, determines whether or not the system 103 ## finds the number of arguments, determines whether or not the system
104 ## is SISO, andd computes the frequency response. Only the way the 104 ## is SISO, andd computes the frequency response. Only the way the
105 ## response is plotted is different between the two functions. 105 ## response is plotted is different between the two functions.
106 106
107 ## check number of input arguments given 107 ## check number of input arguments given
108 if (nargin < 1 | nargin > 5) 108 if (nargin < 1 || nargin > 5)
109 print_usage (); 109 print_usage ();
110 endif 110 endif
111 if(nargin < 2) 111 if (nargin < 2)
112 w = []; 112 w = [];
113 endif 113 endif
114 if(nargin < 3) 114 if (nargin < 3)
115 outputs = []; 115 outputs = [];
116 endif 116 endif
117 if(nargin < 4) 117 if (nargin < 4)
118 inputs = []; 118 inputs = [];
119 endif 119 endif
120 if(nargin < 5) 120 if (nargin < 5)
121 atol = 0; 121 atol = 0;
122 elseif(!(is_sample(atol) | atol == 0)) 122 elseif (! (is_sample (atol) || atol == 0))
123 error("atol must be a nonnegative scalar.") 123 error ("nyquist: atol must be a nonnegative scalar")
124 endif 124 endif
125 125
126 ## signal to __bodquist__ who's calling 126 ## signal to __bodquist__ who's calling
127 127
128 [f, w, sys] = __bodquist__ (sys, w, outputs, inputs, "nyquist"); 128 [f, w, sys] = __bodquist__ (sys, w, outputs, inputs, "nyquist");
129 129
130 ## Get the real and imaginary part of f. 130 ## Get the real and imaginary part of f.
131 realp = real(f); 131 realp = real (f);
132 imagp = imag(f); 132 imagp = imag (f);
133 133
134 ## No output arguments, then display plot, otherwise return data. 134 ## No output arguments, then display plot, otherwise return data.
135 if (nargout == 0) 135 if (nargout == 0)
136 dnplot = 0; 136 dnplot = 0;
137 while(!dnplot) 137 while (! dnplot)
138 oneplot(); 138 plot (realp, imagp, "- ;+w;", realp, -imagp, "-@ ;-w;");
139 __gnuplot_set__ key; 139
140 clearplot();
141 grid ("on"); 140 grid ("on");
142 __gnuplot_set__ data style lines; 141
143 142 if (is_digital (sys))
144 if(is_digital(sys))
145 tstr = " G(e^{jw}) "; 143 tstr = " G(e^{jw}) ";
146 else 144 else
147 tstr = " G(jw) "; 145 tstr = " G(jw) ";
148 endif 146 endif
149 xlabel(["Re(",tstr,")"]); 147 xlabel (sprintf ("Re(%s)", tstr));
150 ylabel(["Im(",tstr,")"]); 148 ylabel (sprintf ("Im(%s)", tstr));
151 149
152 [stn, inn, outn] = sysgetsignals(sys); 150 [stn, inn, outn] = sysgetsignals (sys);
153 if(is_siso(sys)) 151 if (is_siso (sys))
154 title(sprintf("Nyquist plot from %s to %s, w (rad/s) in [%e, %e]", ... 152 title (sprintf ("Nyquist plot from %s to %s, w (rad/s) in [%e, %e]",
155 inn{1}, outn{1}, w(1), w(length(w))) ) 153 inn{1}, outn{1}, w(1), w(end)));
156 endif 154 endif
157 155
158 __gnuplot_set__ nologscale xy; 156 axis (axis2dlim ([[realp(:), imagp(:)]; [realp(:), -imagp(:)]]));
159
160 axis(axis2dlim([[vec(realp),vec(imagp)];[vec(realp),-vec(imagp)]]));
161 plot(realp,imagp,"- ;+w;",realp,-imagp,"-@ ;-w;");
162 157
163 ## check for interactive plots 158 ## check for interactive plots
164 dnplot = 1; # assume done; will change later if atol is satisfied 159 dnplot = 1; # assume done; will change later if atol is satisfied
165 if(atol > 0 & length(f) > 2) 160 if (atol > 0 && length (f) > 2)
166 161
167 ## check for asymptotes 162 ## check for asymptotes
168 fmax = max(abs(f)); 163 fmax = max (abs (f));
169 fi = find(abs(f) == fmax, 1, "last"); 164 fi = find (abs (f) == fmax, 1, "last");
170 165
171 ## compute angles from point to point 166 ## compute angles from point to point
172 df = diff(f); 167 df = diff (f);
173 th = atan2(real(df),imag(df))*180/pi; 168 th = atan2 (real (df), imag (df)) * 180 / pi;
174 169
175 ## get angle at fmax 170 ## get angle at fmax
176 if(fi == length(f)) fi = fi-1; endif 171 if (fi == length(f))
172 fi = fi-1;
173 endif
177 thm = th(fi); 174 thm = th(fi);
178 175
179 ## now locate consecutive angles within atol of thm 176 ## now locate consecutive angles within atol of thm
180 ith_same = find(abs(th - thm) < atol); 177 ith_same = find (abs (th - thm) < atol);
181 ichk = union(fi,find(diff(ith_same) == 1)); 178 ichk = union (fi, find (diff (ith_same) == 1));
182 179
183 ## locate max, min consecutive indices in ichk 180 ## locate max, min consecutive indices in ichk
184 loval = max(complement(ichk,1:fi)); 181 loval = max (complement (ichk, 1:fi));
185 if(isempty(loval)) loval = fi; 182 if (isempty (loval))
186 else loval = loval + 1; endif 183 loval = fi;
187 184 else
188 hival = min(complement(ichk,fi:length(th))); 185 loval = loval + 1;
189 if(isempty(hival)) hival = fi+1; endif 186 endif
190 187
191 keep_idx = complement(loval:hival,1:length(w)); 188 hival = min (complement (ichk, fi:length(th)));
192 189 if (isempty (hival))
193 if(length(keep_idx)) 190 hival = fi+1;
194 resp = input("Remove asymptotes and zoom in (y or n): ",1); 191 endif
195 if(resp(1) == "y") 192
193 keep_idx = complement (loval:hival, 1:length(w));
194
195 if (length (keep_idx))
196 resp = input ("Remove asymptotes and zoom in (y or n): ", 1);
197 if (resp(1) == "y")
196 dnplot = 0; # plot again 198 dnplot = 0; # plot again
197 w = w(keep_idx); 199 w = w(keep_idx);
198 f = f(keep_idx); 200 f = f(keep_idx);
199 realp = real(f); 201 realp = real (f);
200 imagp = imag(f); 202 imagp = imag (f);
201 endif 203 endif
202 endif 204 endif
203 205
204 endif 206 endif
205 endwhile 207 endwhile
206 w = []; 208 w = realp = imagp = [];
207 realp=[]; 209 endif
208 imagp=[];
209 endif
210 210
211 endfunction 211 endfunction