|
3431
|
1 ## Copyright (C) 1998 Auburn University. All rights reserved. |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
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 |
|
|
7 ## Free Software Foundation; either version 2, or (at your option) any |
|
|
8 ## later version. |
|
|
9 ## |
|
|
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 |
|
|
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
13 ## for more details. |
|
|
14 ## |
|
|
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 |
|
|
17 ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. |
|
|
18 |
|
|
19 ## -*- texinfo -*- |
|
3502
|
20 ## @deftypefn {Function File} {[@var{dsys}, @var{fidx}] =} dmr2d (@var{sys}, @var{idx}, @var{sprefix}, @var{ts2}, @var{cuflg}) |
|
3431
|
21 ## convert a multirate digital system to a single rate digital system |
|
3502
|
22 ## states specified by @var{idx}, @var{sprefix} are sampled at @var{ts2}, all |
|
|
23 ## others are assumed sampled at @var{ts1} = @code{sysgettsam (@var{sys})}. |
|
3431
|
24 ## |
|
|
25 ## @strong{Inputs} |
|
|
26 ## @table @var |
|
|
27 ## @item sys |
|
|
28 ## discrete time system; |
|
|
29 ## @code{dmr2d} exits with an error if @var{sys} is not discrete |
|
|
30 ## @item idx |
|
3462
|
31 ## indices or names of states with sampling time |
|
|
32 ## @code{sysgettsam(@var{sys})} (may be empty); see @code{listidx} |
|
3431
|
33 ## @item sprefix |
|
|
34 ## list of string prefixes of states with sampling time |
|
3462
|
35 ## @code{sysgettsam(@var{sys})} (may be empty) |
|
3502
|
36 ## @item ts2 |
|
3431
|
37 ## sampling time of states not specified by @var{idx}, @var{sprefix} |
|
|
38 ## must be an integer multiple of @code{sysgettsam(@var{sys})} |
|
|
39 ## @item cuflg |
|
|
40 ## "constant u flag" if @var{cuflg} is nonzero then the system inputs are |
|
3502
|
41 ## assumed to be constant over the revised sampling interval @var{ts2}. |
|
3431
|
42 ## Otherwise, since the inputs can change during the interval |
|
3502
|
43 ## @var{t} in @math{[k ts2, (k+1) ts2]}, an additional set of inputs is |
|
3431
|
44 ## included in the revised B matrix so that these intersample inputs |
|
|
45 ## may be included in the single-rate system. |
|
|
46 ## default @var{cuflg} = 1. |
|
|
47 ## @end table |
|
|
48 ## |
|
|
49 ## @strong{Outputs} |
|
|
50 ## @table @var |
|
|
51 ## @item dsys |
|
3502
|
52 ## equivalent discrete time system with sampling time @var{ts2}. |
|
3431
|
53 ## |
|
3502
|
54 ## The sampling time of sys is updated to @var{ts2}. |
|
3431
|
55 ## |
|
|
56 ## if @var{cuflg}=0 then a set of additional inputs is added to |
|
|
57 ## the system with suffixes _d1, ..., _dn to indicate their |
|
3502
|
58 ## delay from the starting time k @var{ts2}, i.e. |
|
3431
|
59 ## u = [u_1; u_1_d1; ..., u_1_dn] where u_1_dk is the input |
|
3502
|
60 ## k*ts1 units of time after u_1 is sampled. (@var{ts1} is |
|
|
61 ## the original sampling time of the discrete time system and |
|
|
62 ## @var{ts2} = (n+1)*ts1) |
|
3431
|
63 ## |
|
|
64 ## @item fidx |
|
|
65 ## indices of "formerly fast" states specified by @var{idx} and @var{sprefix}; |
|
3502
|
66 ## these states are updated to the new (slower) sampling interval @var{ts2}. |
|
3431
|
67 ## @end table |
|
|
68 ## |
|
|
69 ## @strong{WARNING} Not thoroughly tested yet; especially when |
|
|
70 ## @var{cuflg} == 0. |
|
|
71 ## @end deftypefn |
|
|
72 |
|
|
73 ## Adapted from c2d by a.s.hodel@eng.auburn.edu |
|
|
74 |
|
|
75 function [dsys, fidx] = dmr2d (sys, idx, sprefix, Ts2, cuflg) |
|
|
76 |
|
|
77 ## parse input arguments |
|
|
78 if(nargin != 4 | nargout > 2) |
|
|
79 usage("[dsys,fidx] = dmr2d (sys, idx, sprefix, Ts2 {,cuflg})"); |
|
|
80 |
|
|
81 elseif (!is_struct(sys)) |
|
|
82 error("sys must be in system data structure form"); |
|
|
83 |
|
|
84 elseif(!is_digital(sys)) |
|
|
85 error("sys must be discrete-time; continuous time passed"); |
|
3462
|
86 |
|
|
87 endif |
|
|
88 |
|
|
89 if(is_signal_list(idx) | isstr(idx)) |
|
|
90 idx = sysidx(sys,"st",idx); |
|
3431
|
91 |
|
|
92 elseif (!(is_vector(idx) | isempty(idx))) |
|
|
93 error(["idx(",num2str(rows(idx)),"x",num2str(columns(idx)), ... |
|
|
94 ") must be a vector"]); |
|
|
95 |
|
|
96 elseif (any(idx <= 0)) |
|
|
97 idv = find(idx <= 0); |
|
|
98 ii = idv(1); |
|
|
99 error(["idx(",num2str(ii),")=",num2str(idx(ii)), ... |
|
|
100 "; entries of idx must be positive"]); |
|
|
101 |
|
|
102 elseif(!(is_signal_list(sprefix) | isempty(sprefix))) |
|
|
103 error("sprefix must be a signal list (see is_signal_list) or empty"); |
|
|
104 |
|
|
105 elseif(!is_sample(Ts2)) |
|
|
106 error(["Ts2=",num2str(Ts2),"; invalid sampling time"]); |
|
|
107 |
|
|
108 endif |
|
|
109 |
|
|
110 ## optional argument: cuflg |
|
|
111 if(nargin <= 4) |
|
|
112 cuflg = 1; # default: constant inputs over Ts2 sampling interv. |
|
|
113 elseif( !is_scalar(cuflg) ) |
|
|
114 error("cuflg must be a scalar") |
|
|
115 elseif( cuflg != 0 | cuflg != 1) |
|
|
116 error(["cuflg = ",num2str(cuflg),", should be 0 or 1"]); |
|
|
117 endif |
|
|
118 |
|
|
119 ## extract state space information |
|
|
120 [da,db,dc,dd,Ts1,nc,nz,stname,inname,outname,yd] = sys2ss(sys); |
|
|
121 |
|
|
122 ## compute number of steps |
|
|
123 if(Ts1 > Ts2) |
|
3462
|
124 error(["Current sampling time=",num2str(Ts1)," > Ts2=",num2str(Ts2)]); |
|
3431
|
125 endif |
|
|
126 nstp = floor(Ts2/Ts1+0.5); |
|
|
127 if(abs((Ts2 - Ts1*nstp)/Ts1) > 1e-12) |
|
|
128 warning(["dmr2d: Ts1=",num2str(Ts1),", Ts2=",num2str(Ts2), ... |
|
|
129 ", selecting nsteps=",num2str(nstp),"; mismatch"]); |
|
|
130 endif |
|
|
131 |
|
|
132 if(isempty(sprefix) & isempty(idx)) |
|
|
133 warning("both sprefix and idx are empty; returning dsys=sys"); |
|
|
134 fidx = []; |
|
|
135 dsys = sys; |
|
|
136 return |
|
|
137 elseif(isempty(sprefix)) |
|
|
138 fidx = idx; |
|
|
139 else |
|
|
140 fidx = reshape(idx,1,length(idx)); |
|
|
141 ## find states whose name begins with any strings in sprefix. |
|
|
142 ns = length(sprefix); |
|
|
143 for kk=1:ns |
|
|
144 spk = nth(sprefix,kk); # get next prefix and length |
|
|
145 spl = length(spk); |
|
|
146 |
|
|
147 ## check each state name |
|
|
148 for ii=1:nz |
|
|
149 sti = nth(stname,ii); # compare spk with this state name |
|
|
150 if(length(sti) >= spl) |
|
|
151 ## if the prefix matches and ii isn't already in the list, add ii |
|
|
152 if(strcmp(sti(1:spl),spk) & !any(fidx == ii) ) |
|
|
153 fidx = sort([fidx,ii]); |
|
|
154 endif |
|
|
155 endif |
|
|
156 endfor |
|
|
157 endfor |
|
|
158 endif |
|
|
159 |
|
|
160 if(nstp == 0) |
|
|
161 warning("dmr2d: nstp = 0; setting tsam and returning"); |
|
|
162 dsys = syschtsam(sys,Ts2); |
|
|
163 return |
|
|
164 elseif(nstp < 0) |
|
|
165 error(["nstp = ", num2str(nstp)," < 0; this shouldn't be!"]); |
|
|
166 endif |
|
|
167 |
|
|
168 ## permute system matrices |
|
|
169 pv = sysreorder(nz,fidx); |
|
|
170 pv = pv(nz:-1:1); # reverse order to put fast modes in leading block |
|
|
171 |
|
|
172 ## construct inverse permutation |
|
|
173 Inz = eye(nz); |
|
|
174 pvi = (Inz(pv,:)'*[1:nz]')'; |
|
|
175 |
|
|
176 ## permute A, B (stname permuted for debugging only) |
|
|
177 da = da(pv,pv); |
|
|
178 db = db(pv,:); |
|
3462
|
179 stname = stname(pv); |
|
3431
|
180 |
|
|
181 ## partition A, B: |
|
|
182 lfidx = length(fidx); |
|
|
183 bki = 1:lfidx; |
|
|
184 a11 = da(bki,bki); |
|
|
185 b1 = db(bki,:); |
|
|
186 |
|
|
187 if(lfidx < nz) |
|
|
188 lfidx1 = lfidx+1; |
|
|
189 bki2 = (lfidx1):nz; |
|
|
190 a12 = da(bki,bki2); |
|
|
191 b2 = db(bki2,:); |
|
|
192 else |
|
|
193 warning("dmr2d: converting entire A,B matrices to new sampling rate"); |
|
|
194 lfidx1 = -1; |
|
|
195 bki2 = []; |
|
|
196 endif |
|
|
197 |
|
|
198 ## begin system conversion: nstp steps |
|
|
199 |
|
|
200 ## compute abar_{n-1}*a12 and appropriate b matrix stuff |
|
|
201 a12b = a12; # running total of abar_{n-1}*a12 |
|
|
202 a12w = a12; # current a11^n*a12 (start with n = 0) |
|
|
203 if(cuflg) |
|
|
204 b1b = b1; |
|
|
205 b1w = b1; |
|
|
206 else |
|
|
207 ## cuflg == 0, need to keep track of intersample inputs too |
|
|
208 nzdx = find(max(abs(b1)) != 0); # FIXME: check tolerance relative to ||b1|| |
|
|
209 b1w = b1(nzdx); |
|
|
210 innamenz = inname(nzdx); |
|
|
211 b1b = b1; # initial b1 must match columns in b2 |
|
|
212 endif |
|
|
213 |
|
|
214 ## compute a11h = a11^nstp by squaring |
|
|
215 a11h = eye(size(a11)); |
|
|
216 p2 = 1; |
|
|
217 a11p2 = a11; #a11^p2 |
|
|
218 |
|
|
219 nstpw = nstp; # workspace for computing a11^nstp |
|
|
220 while(nstpw > 0.5) |
|
|
221 oddv = rem(nstpw,2); |
|
|
222 if(oddv) |
|
|
223 a11h = a11h*a11p2; |
|
|
224 endif |
|
|
225 nstpw = (nstpw-oddv)/2; |
|
|
226 if(nstpw > 0.5) |
|
|
227 a11p2 = a11p2*a11p2; # a11^(next power of 2) |
|
|
228 endif |
|
|
229 endwhile |
|
|
230 |
|
|
231 ## FIXME: this part should probably also use squaring, but |
|
|
232 ## that would require exponentially growing memory. What do do? |
|
|
233 for kk=2:nstp |
|
|
234 ## update a12 block to sum(a12 + ... + a11^(kk-1)*a12) |
|
|
235 a12w = a11*a12w; |
|
|
236 a12b = a12b + a12w; |
|
|
237 |
|
|
238 ## similar for b1 block (checking for cuflg first!) |
|
|
239 b1w = a11*b1w; |
|
|
240 if(cuflg) |
|
|
241 b1b = b1b + b1w; # update b1 block just like we did a12 |
|
|
242 else |
|
|
243 b1b = [b1b, b1w]; # append new inputs |
|
|
244 newin = strappend(innamenz,["_d",num2str(kk-1)]); |
|
|
245 inname = append(inname,newin); |
|
|
246 endif |
|
|
247 endfor |
|
|
248 |
|
|
249 ## reconstruct system and return |
|
|
250 da(bki,bki) = a11h; |
|
|
251 db(bki,1:columns(b1b)) = b1b; |
|
|
252 if(!isempty(bki2)) |
|
|
253 da(bki,bki2) = a12b; |
|
|
254 endif |
|
|
255 |
|
|
256 da = da(pvi,pvi); |
|
|
257 db = db(pvi,:); |
|
3462
|
258 stname = stname(pvi); |
|
3431
|
259 |
|
|
260 ## construct new system and return |
|
|
261 dsys = ss2sys(da,db,dc,dd,Ts2,0,nz,stname,inname,outname,find(yd == 1)); |
|
|
262 |
|
|
263 endfunction |
