d2cm.m

经典通信系统仿真书籍《通信系统与 MATLAB (Proakis)》源代码

function [acout,bc,cc,dc] = d2cm(a,b,c,d,Ts,method,w)
%D2CM	Conversion of discrete LTI systems to continuous-time.
%	[Ac,Bc,Cc,Dc] = D2CM(A,B,C,D,Ts,'method') converts the discrete-
%	time state-space system to continuous-time using 'method':
%	  'zoh'	        Convert to continuous time assuming a zero order
%	                hold on the inputs.
%	  'tustin'      Convert to continuous time using the bilinear 
%	                (Tustin) approximation to the derivative.
%	  'prewarp'     Convert to continuous time using the bilinear 
%	                (Tustin) approximation with frequency prewarping.
%	                Specify the critical frequency with an additional
%	                argument, i.e. D2CM(A,B,C,D,Ts,'prewarp',Wc)
%	  'matched'     Convert the SISO system to continuous time using
%	                the matched pole-zero method.
%
%	[NUMc,DENc] = D2CM(NUM,DEN,Ts,'method') converts the discrete-time
%	polynomial transfer function G(z) = NUM(z)/DEN(z) to continuous
%	time, G(s) = NUMc(s)/DENc(s), using 'method'.
%
%	Note:  'foh' is no longer available.
%
%	See also: D2C, and C2DM.

%	Clay M. Thompson  7-19-90
%	Copyright (c) 1986-93 by the MathWorks, Inc.

error(nargchk(3,7,nargin));

tol = 1.e-4; % Tolerance for 'foh' integrator and zero checking

tf = 0;
% --- Determine which syntax is being used ---
if (nargin==3),		% Transfer function without method, assume 'zoh'
  [num,den] = tfchk(a,b);
  Ts = c;
  method = 'zoh';
  [a,b,c,d] = tf2ss(num,den);
  tf = 1;

elseif (nargin==4),	% Transfer function with method.
  [num,den] = tfchk(a,b);
  Ts = c;
  method = d;
  [a,b,c,d] = tf2ss(num,den);
  tf = 1;

elseif (nargin==5),
  if isstr(d),		% Transfer function with method and prewarp const.
    [num,den] = tfchk(a,b);
    w = Ts;
    Ts = c;
    method = d;
    [a,b,c,d] = tf2ss(num,den); 
    tf = 1;
  else			% State space system without method, assume 'zoh'
    error(abcdchk(a,b,c,d));
    method = 'zoh';
  end

else			% State space system with method.
  error(abcdchk(a,b,c,d));

end
[nx,na] = size(a);
[nb,nu] = size(b);

% --- Determine conversion method ---
if method(1)=='z',	% Zero order hold approximation.
  [ac,bc] = d2c(a,b,Ts);
  cc = c; dc = d;

elseif method(1)=='f',	% First order hold approximation.
  error('Conversion to continuous time with ''foh'' is not available.')

elseif method(1)=='t',	% Tustin approximation.
  I = eye(nx);
  r = 2/Ts;
  P = inv(I + a);
  ac = r*(a-I)*P;
  bc = 2*P*b;
  cc = r*c*P;
  dc = d - c*P*b;

elseif method(1)=='p',	% Tustin approximation with frequency prewarping.
  if ~((nargin==5)|(nargin==7)),
    error('The critical frequency must be specified when using ''prewarp''.');
  end
  T = 2*tan(w*Ts/2)/w;		% Prewarp
  I = eye(nx);
  r = 2/T;
  P = inv(I + a);
  ac = r*(a-I)*P;
  bc = 2*P*b;
  cc = r*c*P;
  dc = d - c*P*b;
  
elseif method(1)=='m',	% Matched pole-zero approximation.
  [ny,nu] = size(d);
  if (ny>1)|(nu>1),
    error('System must be SISO for matched pole-zero method.');
  end
  if tf,
    z = roots(num); p = roots(den);
  else
    [z,p] = ss2zp(a,b,c,d,1);
  end
  z=esort(z); p = esort(p);
  pc = log(p)./Ts;
  zc = zeros(length(z),1);
  zinf = (abs(z+1)<sqrt(eps));	% Fuzzy compare (z==-1)
  if ~isempty(zc),
    zc(~zinf) = log(z(~zinf))./Ts;
    zc(zinf) = inf*ones(length(z(zinf)),1);
  end
  [ac,bc,cc,dc] = zp2ss(zc,pc,1);

  % Match D.C. gain or gain at s=1 for singular systems.
  if any(abs(p-1)<sqrt(eps)), % Match gain at s = 1.
    w = exp(sqrt(-1)*Ts);
    if tf,
      kd = abs(polyval(num,w))/abs(polyval(den,w));
    else
      kd = abs(c/(eye(nx)*w-a)*b + d);
    end
    kc = cc/(eye(nx)-ac)*bc + dc;
  else
    if tf,
      kd = sum(num')'/sum(den);
    else
      kd = c/(eye(nx)-a)*b + d;
    end
    kc = -cc/ac*bc + dc;
  end
  km = sqrt(abs(kd/kc));
  sm = sign(kd/kc);
  bc = bc.*km;
  cc = cc.*km.*sm;
  dc = dc.*km.*km.*sm;

else
  error('Conversion method is unknown.');

end

if nargout==0,		% Compare Bode or Singular value plots
  [ny,nc] = size(c);
  if (ny==1)&(nu==1),	% Plot Bode plots
    [magd,phased,wd] = dbode(a,b,c,d,Ts,1);
    [magc,phasec,wc] = bode(ac,bc,cc,dc,1);
    semilogx(wd,20*log10(magd),'-',wc,20*log10(magc),'--')
    xlabel('Frequency (rad/sec)'), ylabel('Gain dB')
    title('D2CM comparison plot')
  else
    [svd,wd] = dsigma(a,b,c,d,Ts);
    [svc,wc] = sigma(ac,bc,cc,dc);
    semilogx(wd,20*log10(svd),'-',wc,20*log10(svc),'--')
    xlabel('Frequency (rad/sec)'), ylabel('Singular Values dB')
    title('D2CM comparison plot')
  end
  return
end

if tf,		% Convert to TF form for output
  [ac,bc] = ss2tf(ac,bc,cc,dc,1);
end

acout = ac;