Discretization - Continuous-Time to Discrete-Time System Conversion with MATLAB

[@WonYYang]

The related files can be downloaded at kr.mathworks.com/matlabcentral/fileexchange/176363-discretization-ct-to-dt-system-conversion-with-matlab.
%ct2dt_1.m
clear, clf
% Analytical solution using the Laplace traasform
syms s a1 y0
ya=ilaplace((1/s+y0)/(s+a1))
ya_=@(t)(1-exp(-a1*t))/a1+y0*exp(-a1*t);
a1=1; y0=0;
% Solve the differential equation
t0=0; tf=4; tspan=[t0 tf];
u=@(t)t>=0; % Unit step function
dy=@(t,y)-a1*y+u(t); % Differential eq
% Numerical solution using ode45()
[tt,y]=ode45(dy,tspan,y0);
% Plot the numerical/analytical solutions
subplot(141)
plot(tt,y, tt,eval(ya_(tt)),'r:')
% Solve the difference equation
Ts=[0.8 0.2 0.1];
ud=@(n)n>=0;
yd(1)=y0;
for i=1:numel(Ts)
T=Ts(i); Nt=ceil((tf-t0)/T);
nn=[0:Nt];
% Solve the difference equation
for n=1:Nt
yd(n+1)=(1-a1*T)*yd(n)+T*ud(n); % Difference equation
end
subplot(141+i)
% Plot the solutions of differential/difference eqs
plot(tt,y), hold on, stem(nn*T,yd,'.')
end
%ct2dt_2.m
clear, clf
%syms s z
a1=200; a2=5e4;
%a1=2; a2=2;
B=a2; A=[1 a1 a2]; % Numerator/Denominator of an LPF transfer ftn
%B=[a1 0]; A=[1 a1 a2]; %BPF
%B=[1 0 0]; A=[1 a1 a2]; % HPF
Gs=tf(B,A) % Transfer function of a CT system
% Discretization
T=1e-3; % Sampling interval
wp=1/T; %sqrt(a2);
Gz_zoh=c2d(Gs,T,'zoh')
Gz_blt=c2d(Gs,T,'tustin')
%opt=c2dOptions('Method','tustin','PrewarpFrequency',wp);
%Gz_bltp=c2d(Gs,T,opt)
% Step responses of the DT systems
subplot(121)
step(Gs,'k', Gz_zoh,'g', Gz_blt,'b') %, Gz_bltp,'r')
legend('Gs','Gz-zoh','Gz-blt') %,'Gz-bltp')
% Frequency responses of the CT system
W=0:0.001:pi; w=W/T; % DT/CT frequency ranges
Gw=freqs(B,A,w);
GwmdB=20*log10(abs(Gw));
%GwmdB=abs(Gw);
% Frequency responses of the DT systems
[Bz_zoh,Az_zoh]=tfdata(Gz_zoh,'v');
Gw_zoh=freqz(Bz_zoh,Az_zoh,W);
GwmdB_zoh=20*log10(abs(Gw_zoh));
%GwmdB_zoh=abs(Gw_zoh);
[Bz_blt,Az_blt]=tfdata(Gz_blt,'v');
Gw_blt=freqz(Bz_blt,Az_blt,W);
GwmdB_blt=20*log10(abs(Gw_blt));
%GwmdB_blt=abs(Gw_blt);
%[Bz_bltp,Az_bltp]=tfdata(Gz_bltp,'v');
%Gw_bltp=freqz(Bz_bltp,Az_bltp,W);
%GwmdB_bltp=20*log10(abs(Gw_bltp));
%GwmdB_bltp=abs(Gw_bltp);
subplot(122)
plot(w,GwmdB,'k', w,GwmdB_zoh,'g', w,GwmdB_blt,'b') %, w,GwmdB_bltp,'r')
legend('GwmdB','GwmdB-zoh','GwmdB-blt') %,'GwmdB-bltp')
xlabel('CT radian frequency[rad/s]')
title('Frequency responses')
%hold on, plot([wp wp],[-100 0],'m:')
xlim([0 pi/T]); ylim([-50 5]) % for LPF
%xlim([0 3*wp]); ylim([-10 5]) % for HPF