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