1## Copyright (C) 1996, 2000, 2004, 2005, 2006, 2007 2## Auburn University. All rights reserved. 3## Copyright (C) 2009-2016 Lukas F. Reichlin 4## 5## This program is free software; you can redistribute it and/or modify it 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. 9## 10## This program 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. 14## 15## You should have received a copy of the GNU General Public License 16## along with this program; see the file COPYING. If not, see 17## <http://www.gnu.org/licenses/>. 18 19## -*- texinfo -*- 20## @deftypefn {Function File} {@var{w} =} __frequency_vector__ (@var{sys}) 21## Get default range of frequencies based on cutoff frequencies of system 22## poles and zeros. 23## Frequency range is the interval 24## @iftex 25## @tex 26## $ [ 10^{w_{min}}, 10^{w_{max}} ] $ 27## @end tex 28## @end iftex 29## @ifnottex 30## [10^@var{wmin}, 10^@var{wmax}] 31## @end ifnottex 32## 33## Used by @command{__frequency_response__} 34## @end deftypefn 35 36## Adapted-By: Lukas Reichlin <lukas.reichlin@gmail.com> 37## Date: October 2009 38## Version: 0.4 39 40function w = __frequency_vector__ (sys_cell, wbounds = "std", wmin, wmax) 41 42 N = 1000; # intervals within the w range 43 isc = iscell (sys_cell); 44 45 if (! isc) # __sys2frd__ methods pass LTI models not in cells 46 sys_cell = {sys_cell}; 47 endif 48 49 idx = cellfun (@(x) isa (x, "lti"), sys_cell); 50 sys_cell = sys_cell(idx); 51 len = numel (sys_cell); 52 53 [dec_min, dec_max, zp] = cellfun (@__frequency_range__, sys_cell, {wbounds}, "uniformoutput", false); 54 55 if (strcmpi (wbounds, "std")) # plots with explicit frequencies 56 57 if (nargin == 4) # w = {wmin, wmax} 58 dec_min = log10 (wmin); 59 dec_max = log10 (wmax); 60 else 61 dec_min = min (cell2mat (dec_min)); 62 dec_max = max (cell2mat (dec_max)); 63 endif 64 65 zp = horzcat (zp{:}); 66 67 ## include zeros and poles for nice peaks in plots 68 idx = find (zp > 10^dec_min & zp < 10^dec_max); 69 zp = zp(idx); 70 71 w = logspace (dec_min, dec_max, N); 72 w = unique ([w, zp]); # unique also sorts frequency vector 73 74 w = repmat ({w}, 1, len); # return cell of frequency vectors 75 76 elseif (strcmpi (wbounds, "ext")) # plots with implicit frequencies 77 78 if (nargin == 4) 79 dec_min = repmat ({log10(wmin)}, 1, len); 80 dec_max = repmat ({log10(wmax)}, 1, len); 81 endif 82 83 idx = cellfun (@(zp, dec_min, dec_max) find (zp > 10^dec_min & zp < 10^dec_max), ... 84 zp, dec_min, dec_max, "uniformoutput", false); 85 zp = cellfun (@(zp, idx) zp(idx), zp, idx, "uniformoutput", false); 86 87 w = cellfun (@logspace, dec_min, dec_max, {N}, "uniformoutput", false); 88 w = cellfun (@(w, zp) unique ([w, zp]), w, zp, "uniformoutput", false); 89 ## unique also sorts frequency vector 90 91 else 92 error ("frequency_vector: invalid argument 'wbounds'"); 93 endif 94 95 if (! isc) # for __sys2frd__ methods 96 w = w{1}; 97 endif 98 99endfunction 100 101 102function [dec_min, dec_max, zp] = __frequency_range__ (sys, wbounds = "std") 103 104 if (isa (sys, "frd")) 105 w = get (sys, "w"); 106 dec_min = log10 (w(1)); 107 dec_max = log10 (w(end)); 108 zp = []; 109 return; 110 endif 111 112 zer = zero (sys); 113 pol = pole (sys); 114 tsam = abs (get (sys, "tsam")); # tsam could be -1 115 discrete = ! isct (sys); # static gains (tsam = -2) are assumed continuous 116 117 ## make sure zer, pol are row vectors 118 pol = reshape (pol, 1, []); 119 zer = reshape (zer, 1, []); 120 121 ## check for natural frequencies away from omega = 0 122 if (discrete) 123 ## The 2nd conditions prevents log(0) in the next log command 124 iiz = find (abs(zer-1) > norm(zer)*eps && abs(zer) > norm(zer)*eps); 125 iip = find (abs(pol-1) > norm(pol)*eps && abs(pol) > norm(pol)*eps); 126 127 ## avoid dividing empty matrices, it would work but looks nasty 128 if (! isempty (iiz)) 129 czer = log (zer(iiz))/tsam; 130 else 131 czer = []; 132 endif 133 134 if (! isempty (iip)) 135 cpol = log (pol(iip))/tsam; 136 else 137 cpol = []; 138 endif 139 else 140 ## continuous 141 iip = find (abs(pol) > norm(pol)*eps); 142 iiz = find (abs(zer) > norm(zer)*eps); 143 144 if (! isempty (zer)) 145 czer = zer(iiz); 146 else 147 czer = []; 148 endif 149 if (! isempty (pol)) 150 cpol = pol(iip); 151 else 152 cpol = []; 153 endif 154 endif 155 156 if (isempty (iip) && isempty (iiz)) 157 ## no poles/zeros away from omega = 0; pick defaults 158 dec_min = 0; # -1 159 dec_max = 2; # 3 160 else 161 dec_min = floor (log10 (min (abs ([cpol, czer])))); 162 dec_max = ceil (log10 (max (abs ([cpol, czer])))); 163 endif 164 165 ## expand to show the entirety of the "interesting" portion of the plot 166 switch (wbounds) 167 case "std" # standard 168 if (dec_min == dec_max) 169 dec_min -= 2; 170 dec_max += 2; 171 else 172 dec_min--; 173 dec_max++; 174 endif 175 case "ext" # extended (for nyquist) 176 if (any (abs (pol) < sqrt (eps))) # look for integrators 177 dec_min -= 0.5; 178 dec_max += 2; 179 else 180 dec_min -= 2; 181 dec_max += 2; 182 endif 183 otherwise 184 error ("frequency_range: second argument invalid"); 185 endswitch 186 187 ## run discrete frequency all the way to pi 188 if (discrete) 189 dec_max = log10 (pi/tsam); 190 endif 191 192 ## include zeros and poles for nice peaks in plots 193 zp = [abs(zer), abs(pol)]; 194 195endfunction 196