home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / D2CM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.7 KB  |  168 lines
  1. function [acout,bc,cc,dc] = d2cm(a,b,c,d,Ts,method,w)
  2. %D2CM    Conversion of discrete LTI systems to continuous-time.
  3. %    [Ac,Bc,Cc,Dc] = D2CM(A,B,C,D,Ts,'method') converts the discrete-
  4. %    time state-space system to continuous-time using 'method':
  5. %      'zoh'            Convert to continuous time assuming a zero order
  6. %                    hold on the inputs.
  7. %      'tustin'      Convert to continuous time using the bilinear 
  8. %                    (Tustin) approximation to the derivative.
  9. %      'prewarp'     Convert to continuous time using the bilinear 
  10. %                    (Tustin) approximation with frequency prewarping.
  11. %                    Specify the critical frequency with an additional
  12. %                    argument, i.e. D2CM(A,B,C,D,Ts,'prewarp',Wc)
  13. %      'matched'     Convert the SISO system to continuous time using
  14. %                    the matched pole-zero method.
  15. %
  16. %    [NUMc,DENc] = D2CM(NUM,DEN,Ts,'method') converts the discrete-time
  17. %    polynomial transfer function G(z) = NUM(z)/DEN(z) to continuous
  18. %    time, G(s) = NUMc(s)/DENc(s), using 'method'.
  19. %
  20. %    Note:  'foh' is no longer available.
  21. %
  22. %    See also: D2C, and C2DM.
  23.  
  24. %    Clay M. Thompson  7-19-90
  25. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  26.  
  27. error(nargchk(3,7,nargin));
  28.  
  29. tol = 1.e-4; % Tolerance for 'foh' integrator and zero checking
  30.  
  31. tf = 0;
  32. % --- Determine which syntax is being used ---
  33. if (nargin==3),        % Transfer function without method, assume 'zoh'
  34.   [num,den] = tfchk(a,b);
  35.   Ts = c;
  36.   method = 'zoh';
  37.   [a,b,c,d] = tf2ss(num,den);
  38.   tf = 1;
  39.  
  40. elseif (nargin==4),    % Transfer function with method.
  41.   [num,den] = tfchk(a,b);
  42.   Ts = c;
  43.   method = d;
  44.   [a,b,c,d] = tf2ss(num,den);
  45.   tf = 1;
  46.  
  47. elseif (nargin==5),
  48.   if isstr(d),        % Transfer function with method and prewarp const.
  49.     [num,den] = tfchk(a,b);
  50.     w = Ts;
  51.     Ts = c;
  52.     method = d;
  53.     [a,b,c,d] = tf2ss(num,den); 
  54.     tf = 1;
  55.   else            % State space system without method, assume 'zoh'
  56.     error(abcdchk(a,b,c,d));
  57.     method = 'zoh';
  58.   end
  59.  
  60. else            % State space system with method.
  61.   error(abcdchk(a,b,c,d));
  62.  
  63. end
  64. [nx,na] = size(a);
  65. [nb,nu] = size(b);
  66.  
  67. % --- Determine conversion method ---
  68. if method(1)=='z',    % Zero order hold approximation.
  69.   [ac,bc] = d2c(a,b,Ts);
  70.   cc = c; dc = d;
  71.  
  72. elseif method(1)=='f',    % First order hold approximation.
  73.   error('Conversion to continuous time with ''foh'' is not available.')
  74.  
  75. elseif method(1)=='t',    % Tustin approximation.
  76.   I = eye(nx);
  77.   r = 2/Ts;
  78.   P = inv(I + a);
  79.   ac = r*(a-I)*P;
  80.   bc = 2*P*b;
  81.   cc = r*c*P;
  82.   dc = d - c*P*b;
  83.  
  84. elseif method(1)=='p',    % Tustin approximation with frequency prewarping.
  85.   if ~((nargin==5)|(nargin==7)),
  86.     error('The critical frequency must be specified when using ''prewarp''.');
  87.   end
  88.   T = 2*tan(w*Ts/2)/w;        % Prewarp
  89.   I = eye(nx);
  90.   r = 2/T;
  91.   P = inv(I + a);
  92.   ac = r*(a-I)*P;
  93.   bc = 2*P*b;
  94.   cc = r*c*P;
  95.   dc = d - c*P*b;
  96.   
  97. elseif method(1)=='m',    % Matched pole-zero approximation.
  98.   [ny,nu] = size(d);
  99.   if (ny>1)|(nu>1),
  100.     error('System must be SISO for matched pole-zero method.');
  101.   end
  102.   if tf,
  103.     z = roots(num); p = roots(den);
  104.   else
  105.     [z,p] = ss2zp(a,b,c,d,1);
  106.   end
  107.   z=esort(z); p = esort(p);
  108.   pc = log(p)./Ts;
  109.   zc = zeros(length(z),1);
  110.   zinf = (abs(z+1)<sqrt(eps));    % Fuzzy compare (z==-1)
  111.   if ~isempty(zc),
  112.     zc(~zinf) = log(z(~zinf))./Ts;
  113.     zc(zinf) = inf*ones(length(z(zinf)),1);
  114.   end
  115.   [ac,bc,cc,dc] = zp2ss(zc,pc,1);
  116.  
  117.   % Match D.C. gain or gain at s=1 for singular systems.
  118.   if any(abs(p-1)<sqrt(eps)), % Match gain at s = 1.
  119.     w = exp(sqrt(-1)*Ts);
  120.     if tf,
  121.       kd = abs(polyval(num,w))/abs(polyval(den,w));
  122.     else
  123.       kd = abs(c/(eye(nx)*w-a)*b + d);
  124.     end
  125.     kc = cc/(eye(nx)-ac)*bc + dc;
  126.   else
  127.     if tf,
  128.       kd = sum(num')'/sum(den);
  129.     else
  130.       kd = c/(eye(nx)-a)*b + d;
  131.     end
  132.     kc = -cc/ac*bc + dc;
  133.   end
  134.   km = sqrt(abs(kd/kc));
  135.   sm = sign(kd/kc);
  136.   bc = bc.*km;
  137.   cc = cc.*km.*sm;
  138.   dc = dc.*km.*km.*sm;
  139.  
  140. else
  141.   error('Conversion method is unknown.');
  142.  
  143. end
  144.  
  145. if nargout==0,        % Compare Bode or Singular value plots
  146.   [ny,nc] = size(c);
  147.   if (ny==1)&(nu==1),    % Plot Bode plots
  148.     [magd,phased,wd] = dbode(a,b,c,d,Ts,1);
  149.     [magc,phasec,wc] = bode(ac,bc,cc,dc,1);
  150.     semilogx(wd,20*log10(magd),'-',wc,20*log10(magc),'--')
  151.     xlabel('Frequency (rad/sec)'), ylabel('Gain dB')
  152.     title('D2CM comparison plot')
  153.   else
  154.     [svd,wd] = dsigma(a,b,c,d,Ts);
  155.     [svc,wc] = sigma(ac,bc,cc,dc);
  156.     semilogx(wd,20*log10(svd),'-',wc,20*log10(svc),'--')
  157.     xlabel('Frequency (rad/sec)'), ylabel('Singular Values dB')
  158.     title('D2CM comparison plot')
  159.   end
  160.   return
  161. end
  162.  
  163. if tf,        % Convert to TF form for output
  164.   [ac,bc] = ss2tf(ac,bc,cc,dc,1);
  165. end
  166.  
  167. acout = ac;
  168.