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