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