Mercurial > octave-nkf

comparison scripts/control/base/__freqresp__.m @ 5568:e9cde940b271

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-12-08 02:28:22 by jwe]
author jwe
date Thu, 08 Dec 2005 02:28:22 +0000
parents 4c8a2e4e0717
children 34f96dd5441b
comparison
equal deleted inserted replaced
5567:80e629357483 5568:e9cde940b271
45 45
46 function [ff, w] = __freqresp__ (sys, USEW, w); 46 function [ff, w] = __freqresp__ (sys, USEW, w);
47 47
48 ## SYS_INTERNAL accesses members of system data structure 48 ## SYS_INTERNAL accesses members of system data structure
49 49
50 save_warn_empty_list_elements = warn_empty_list_elements; 50 ## Check Args
51 unwind_protect 51 if ((nargin < 2) || (nargin > 4))
52 warn_empty_list_elements = 0; 52 usage ("[ff, w] = __freqresp__ (sys, USEW, w)");
53 elseif (USEW & (nargin < 3) )
54 error ("USEW = 1 but w was not passed.");
55 elseif (USEW & isempty(w))
56 warning("USEW = 1 but w is empty; setting USEW=0");
57 USEW = 0;
58 endif
53 59
54 ## Check Args 60 DIGITAL = is_digital(sys);
55 if ((nargin < 2) || (nargin > 4)) 61
56 usage ("[ff, w] = __freqresp__ (sys, USEW, w)"); 62 ## compute default w if needed
57 elseif (USEW & (nargin < 3) ) 63 if(!USEW)
58 error ("USEW = 1 but w was not passed."); 64 if(is_siso(sys))
59 elseif (USEW & isempty(w)) 65 sys = sysupdate(sys,"zp");
60 warning("USEW = 1 but w is empty; setting USEW=0"); 66 [zer,pol] = sys2zp(sys);
61 USEW = 0; 67 else
68 zer = tzero(sys);
69 pol = eig(sys2ss(sys));
62 endif 70 endif
63 71
64 DIGITAL = is_digital(sys); 72 ## get default frequency range
73 [wmin,wmax] = bode_bounds(zer,pol,DIGITAL,sysgettsam(sys));
74 w = logspace(wmin,wmax,50);
75 else
76 w = reshape(w,1,length(w)); # make sure it's a row vector
77 endif
65 78
66 ## compute default w if needed 79 ## now get complex values of s or z
67 if(!USEW) 80 if(DIGITAL)
68 if(is_siso(sys)) 81 jw = exp(i*w*sysgettsam(sys));
69 sys = sysupdate(sys,"zp"); 82 else
70 [zer,pol] = sys2zp(sys); 83 jw = i*w;
71 else 84 endif
72 zer = tzero(sys);
73 pol = eig(sys2ss(sys));
74 endif
75 85
76 ## get default frequency range 86 [nn,nz,mm,pp] = sysdimensions(sys);
77 [wmin,wmax] = bode_bounds(zer,pol,DIGITAL,sysgettsam(sys));
78 w = logspace(wmin,wmax,50);
79 else
80 w = reshape(w,1,length(w)); # make sure it's a row vector
81 endif
82 87
83 ## now get complex values of s or z 88 ## now compute the frequency response - divide by zero yields a warning
84 if(DIGITAL) 89 if (strcmp(sysgettype(sys),"zp"))
85 jw = exp(i*w*sysgettsam(sys)); 90 ## zero-pole form (preferred)
86 else 91 [zer,pol,sysk] = sys2zp(sys);
87 jw = i*w; 92 ff = ones(size(jw));
88 endif 93 l1 = min(length(zer)*(1-isempty(zer)),length(pol)*(1-isempty(pol)));
94 for ii=1:l1
95 ff = ff .* (jw - zer(ii)) ./ (jw - pol(ii));
96 endfor
89 97
90 [nn,nz,mm,pp] = sysdimensions(sys); 98 ## require proper transfer function, so now just get poles.
99 for ii=(l1+1):length(pol)
100 ff = ff ./ (jw - pol(ii));
101 endfor
102 ff = ff*sysk;
91 103
92 ## now compute the frequency response - divide by zero yields a warning 104 elseif (strcmp(sysgettype(sys),"tf"))
93 if (strcmp(sysgettype(sys),"zp")) 105 ## transfer function form
94 ## zero-pole form (preferred) 106 [num,den] = sys2tf(sys);
95 [zer,pol,sysk] = sys2zp(sys); 107 ff = polyval(num,jw)./polyval(den,jw);
96 ff = ones(size(jw)); 108 elseif (mm==pp)
97 l1 = min(length(zer)*(1-isempty(zer)),length(pol)*(1-isempty(pol))); 109 ## The system is square; do state-space form bode plot
98 for ii=1:l1 110 [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys);
99 ff = ff .* (jw - zer(ii)) ./ (jw - pol(ii)); 111 n = sysn + sysnz;
100 endfor 112 for ii=1:length(jw);
113 ff(ii) = det(sysc*((jw(ii).*eye(n)-sysa)\sysb)+sysd);
114 endfor;
115 else
116 ## Must be state space... bode
117 [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys);
118 n = sysn + sysnz;
119 for ii=1:length(jw);
120 ff(ii) = norm(sysc*((jw(ii)*eye(n)-sysa)\sysb)+sysd);
121 endfor
101 122
102 ## require proper transfer function, so now just get poles. 123 endif
103 for ii=(l1+1):length(pol)
104 ff = ff ./ (jw - pol(ii));
105 endfor
106 ff = ff*sysk;
107 124
108 elseif (strcmp(sysgettype(sys),"tf")) 125 w = reshape(w,1,length(w));
109 ## transfer function form 126 ff = reshape(ff,1,length(ff));
110 [num,den] = sys2tf(sys);
111 ff = polyval(num,jw)./polyval(den,jw);
112 elseif (mm==pp)
113 ## The system is square; do state-space form bode plot
114 [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys);
115 n = sysn + sysnz;
116 for ii=1:length(jw);
117 ff(ii) = det(sysc*((jw(ii).*eye(n)-sysa)\sysb)+sysd);
118 endfor;
119 else
120 ## Must be state space... bode
121 [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys);
122 n = sysn + sysnz;
123 for ii=1:length(jw);
124 ff(ii) = norm(sysc*((jw(ii)*eye(n)-sysa)\sysb)+sysd);
125 endfor
126
127 endif
128
129 w = reshape(w,1,length(w));
130 ff = reshape(ff,1,length(ff));
131
132 unwind_protect_cleanup
133 warn_empty_list_elements = save_warn_empty_list_elements;
134 end_unwind_protect
135 127
136 endfunction 128 endfunction
137 129