Mercurial > hg > octave-lyh

comparison scripts/control/system/c2d.m @ 7135:8aa770b6c5bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-11-08 18:54:10 by jwe]
author jwe
date Thu, 08 Nov 2007 18:54:10 +0000
parents a1dbe9d80eee
children aeeb646f6538
comparison
equal deleted inserted replaced
7134:ff0b965b65bc 7135:8aa770b6c5bf
98 ## Author: R. Bruce Tenison <btenison@eng.auburn.edu> 98 ## Author: R. Bruce Tenison <btenison@eng.auburn.edu>
99 ## Created: October 1993 99 ## Created: October 1993
100 ## Updated by John Ingram for system data structure August 1996 100 ## Updated by John Ingram for system data structure August 1996
101 101
102 function dsys = c2d (sys, opt, T) 102 function dsys = c2d (sys, opt, T)
103
103 ## parse input arguments 104 ## parse input arguments
104 if(nargin < 1 | nargin > 3) 105 if (nargin < 1 || nargin > 3)
105 print_usage (); 106 print_usage ();
106 elseif (!isstruct(sys)) 107 elseif (! isstruct (sys))
107 error("sys must be a system data structure"); 108 error ("sys must be a system data structure");
108 elseif (nargin == 1) 109 elseif (nargin == 1)
109 opt = "ex"; 110 opt = "ex";
110 elseif (nargin == 2 & !ischar(opt) ) 111 elseif (nargin == 2 && ! ischar (opt))
111 T = opt; 112 T = opt;
112 opt = "ex"; 113 opt = "ex";
113 endif 114 endif
114 115
115 if (! ischar (opt)) 116 if (! ischar (opt))
116 error ("expecting option as a string"); 117 error ("expecting option as a string");
117 endif 118 endif
118 119
119 ## check if sampling period T was passed. 120 ## check if sampling period T was passed.
120 Ts = sysgettsam(sys); 121 Ts = sysgettsam (sys);
121 if(!exist("T")) 122 if (! exist ("T"))
122 T = Ts; 123 T = Ts;
123 if(T == 0) 124 if (T == 0)
124 error("sys is purely continuous; no sampling period T provided"); 125 error ("sys is purely continuous; no sampling period T provided");
125 endif 126 endif
126 elseif (T != Ts & Ts > 0) 127 elseif (T != Ts && Ts > 0)
127 warning(["c2d: T=",num2str(T),", system tsam==",num2str(Ts), ... 128 warning ("c2d: T=%g, system tsam=%g: using T=", T, Ts, min (T, Ts));
128 ": using T=", num2str(min(T,Ts))]); 129 T = min (T, Ts);
129 T = min(T,Ts); 130 endif
130 endif 131
131 132 if (! is_sample (T))
132 if (!is_sample(T)) 133 error ("sampling period T must be a positive, real scalar");
133 error("sampling period T must be a positive, real scalar");
134 elseif (! (strcmp (opt, "ex") 134 elseif (! (strcmp (opt, "ex")
135 || strcmp (opt, "bi") 135 || strcmp (opt, "bi")
136 || strcmp (opt, "matched"))) 136 || strcmp (opt, "matched")))
137 error ("invalid option passed: %s", opt); 137 error ("invalid option passed: %s", opt);
138 endif 138 endif
139 139
140 sys = sysupdate(sys,"ss"); 140 sys = sysupdate (sys, "ss");
141 [n,nz,m,p] = sysdimensions(sys); 141 [n, nz, m, p] = sysdimensions (sys);
142 if(n == 0) 142 if (n == 0)
143 dsys = syssetsignals(sys,"yd",ones(1:p)); 143 dsys = syssetsignals (sys, "yd", ones(1:p));
144 elseif(strcmp(opt,"ex")); 144 elseif (strcmp (opt, "ex"));
145 [aa,bb,cc,dd] = sys2ss(sys); 145 [aa, bb, cc, dd] = sys2ss (sys);
146 crng= 1:n; 146 crng= 1:n;
147 drng = n+(1:nz); 147 drng = n+(1:nz);
148 148
149 ## partition state equations into continuous, imaginary subsystems 149 ## partition state equations into continuous, imaginary subsystems
150 Ac = aa(crng,crng); 150 Ac = aa(crng,crng);
151 Bc = bb(crng,:); 151 Bc = bb(crng,:);
152 if(nz == 0) 152 if (nz == 0)
153 Acd = Adc = Add = Bd = 0; 153 Acd = Adc = Add = Bd = 0;
154 else 154 else
155 Acd = aa(crng,drng); 155 Acd = aa(crng,drng);
156 Adc = aa(drng,crng); 156 Adc = aa(drng,crng);
157 Add = aa(drng,drng); 157 Add = aa(drng,drng);
159 Bc = [Bc, Acd]; ## append discrete states as inputs to cont system 159 Bc = [Bc, Acd]; ## append discrete states as inputs to cont system
160 endif 160 endif
161 161
162 ## convert state equations 162 ## convert state equations
163 mat = [Ac, Bc; zeros(m+nz,n+nz+m)]; 163 mat = [Ac, Bc; zeros(m+nz,n+nz+m)];
164 matexp = expm(mat * T); 164 matexp = expm (mat * T);
165 165
166 ## replace Ac 166 ## replace Ac
167 aa(crng,crng) = matexp(crng,crng); ## discretized homegenous diff eqn 167 aa(crng,crng) = matexp(crng,crng); ## discretized homegenous diff eqn
168 168
169 ## replace Bc 169 ## replace Bc
170 bb(crng,:) = matexp(crng,n+(1:m)); 170 bb(crng,:) = matexp(crng,n+(1:m));
171 171
172 ## replace Acd 172 ## replace Acd
173 if(nz) 173 if (nz)
174 aa(crng,drng) = matexp(crng,n+m+(1:nz)); 174 aa(crng,drng) = matexp(crng,n+m+(1:nz));
175 end 175 end
176 176
177 stnames = sysgetsignals(sys,"st"); ## continuous states renamed below 177 stnames = sysgetsignals (sys, "st"); ## continuous states renamed below
178 innames = sysgetsignals(sys,"in"); 178 innames = sysgetsignals (sys, "in");
179 outnames = sysgetsignals(sys,"out"); 179 outnames = sysgetsignals (sys, "out");
180 outlist = 1:p; 180 outlist = 1:p;
181 dsys = ss(aa,bb,cc,dd,T,0,n+nz,stnames,innames, ... 181 dsys = ss (aa, bb, cc, dd, T, 0, n+nz, stnames, innames,
182 outnames,outlist); 182 outnames, outlist);
183 ## rename states 183 ## rename states
184 for ii=1:n 184 for ii = 1:n
185 strval = sprintf("%s_d",sysgetsignals(dsys,"st",ii,1)); 185 strval = sprintf ("%s_d", sysgetsignals (dsys, "st", ii, 1));
186 dsys = syssetsignals(dsys,"st",strval,ii); 186 dsys = syssetsignals (dsys, "st", strval, ii);
187 endfor 187 endfor
188 188
189 elseif(strcmp(opt,"bi")) 189 elseif (strcmp (opt, "bi"))
190 if(is_digital(sys)) 190 if (is_digital (sys))
191 error("c2d: system is already digital") 191 error ("c2d: system is already digital")
192 else 192 else
193 ## convert with bilinear transform 193 ## convert with bilinear transform
194 [a,b,c,d,tsam,n,nz,stname,inname,outname,yd] = sys2ss(sys); 194 [a, b, c, d, tsam, n, nz, stname, inname, outname, yd] = sys2ss (sys);
195 IT = (2/T)*eye(size(a)); 195 IT = (2/T) * eye (size (a));
196 A = (IT+a)/(IT-a); 196 A = (IT+a)/(IT-a);
197 iab = (IT-a)\b; 197 iab = (IT-a)\b;
198 tk=2/sqrt(T); 198 tk = 2 / sqrt (T);
199 B = tk*iab; 199 B = tk*iab;
200 C = tk*(c/(IT-a)); 200 C = tk*(c/(IT-a));
201 D = d + (c*iab); 201 D = d + (c*iab);
202 stnamed = strappend(stname,"_d"); 202 stnamed = strappend (stname, "_d");
203 dsys = ss(A,B,C,D,T,0,rows(A),stnamed,inname,outname); 203 dsys = ss (A, B, C, D, T, 0, rows (A), stnamed, inname, outname);
204 endif 204 endif
205 elseif(strcmp(opt,"matched")) 205 elseif (strcmp (opt, "matched"))
206 if(is_digital(sys)) 206 if (is_digital (sys))
207 error("c2d: system is already digital"); 207 error ("c2d: system is already digital");
208 elseif((length(sys.inname) != 1) || (length(sys.outname) != 1)) 208 elseif (length (sys.inname) != 1 || length (sys.outname) != 1)
209 error("c2d: system in not single input, single output"); 209 error ("c2d: system in not single input, single output");
210 else 210 else
211 sys = sysupdate(sys,"zp"); 211 sys = sysupdate (sys, "zp");
212 p = exp(sys.pol*T); 212 p = exp (sys.pol*T);
213 z = exp(sys.zer*T); 213 z = exp (sys.zer*T);
214 infinite_zeros = max(size(sys.pol))-max(size(sys.zer))-1; 214 infinite_zeros = max (size (sys.pol)) - max (size (sys.zer)) - 1;
215 for i = 1:infinite_zeros 215 for i = 1:infinite_zeros
216 z = [z ; -1]; 216 z = [z ; -1];
217 endfor 217 endfor
218 ## Should the freaquency we adjust around always be 1? 218 ## Should the freaquency we adjust around always be 1?
219 [cmag,cphase,cw] = bode(sys,1); 219 [cmag, cphase, cw] = bode (sys, 1);
220 [dmag,dpahse,dw] = bode(zp(z,p,1,T),1); 220 [dmag, dpahse, dw] = bode (zp (z, p, 1, T), 1);
221 dsys = zp(z,p,cmag/dmag,T); 221 dsys = zp (z, p, cmag/dmag, T);
222 endif 222 endif
223 else 223 else
224 error ("invalid option = %s", opt); 224 error ("invalid option = %s", opt);
225 endif 225 endif
226 226
227 endfunction 227 endfunction