Matlab Continuous to Discrete Transformation with Arbitrary Discretization Pattern

Matlab has c2d function to do the continuous to discrete time transformation. The c2d function has zoh foh impulse tustin matched least-squares discretization methods which meet most applications. However, sometimes you may want to implement your own discretization method. This is what this article aims for.

The c2d function

The commonly used discretization method is ZOH and Tustin Approximation (Bilinear Aproximation). The method are name as zoh and tustin in c2d function. The syntax of c2d is

sysd = c2d(sys,Ts)
sysd = c2d(sys,Ts,method)
sysd = c2d(sys,Ts,opts)
[sysd,G] = c2d(sys,Ts,method)
[sysd,G] = c2d(sys,Ts,opts)

ZOH (Backward Differentiation)

z=\cfrac{1}{1-sT_s}

s=\cfrac{1-z^{-1}}{T_s}

Tustin Approximation

without frequency prewarping

The Tustin or bilinear approximation yields the best frequency-domain match between the continuous-time and discretized systems.

z=e^{sT_s}\approx\cfrac{1+sT_s/2}{1-sT_s/2}

s=\cfrac{2}{T_s}\cfrac{1-z^{-1}}{1+z^{-1}}

with frequency prewarping

This method ensures a match between the continuous- and discrete-time responses at the prewarp frequency.

s=\cfrac{\omega_{prewarp}}{\tan(\omega_{prewarp}T_s/2)}\cfrac{1-z^{-1}}{1+z^{-1}}

Arbitrary Discretization

We can use Matlab symbolic operation to do arbitrary discretization. Take a one-order low-pass filter (LPF) for example

H_{LPF}(s) = \cfrac{\omega_c}{s+\omega_c}

Assuming \omega_c=200\pi\ \mathrm{rad/s}, T_s=10^{-3}\ s, with c2d

omega_s = 200*pi;
T_s     = 1e-3;
Hs_LPF  = tf(omega_s, [1 omega_s]);
Hz_LPF  = c2d(Hs_LPF,T_s,'tustin');

The result is

Hz_LPF =

  0.2391 z + 0.2391
  -----------------
     z - 0.5219

Sample time: 0.001 seconds
Discrete-time transfer function.

Symbolic Operation

syms s z
Ss_LPF = omega_s / (s + omega_s);
Sz_LPF = subs(Ss_LPF,s,2/T_s * (1-z^-1) / (1+z^-1));
[Snum,Sden] = numden(Sz_LPF);

num = flip(eval(coeffs(Snum)));
den = flip(eval(coeffs(Sden)));

num = num./den(1);
den = den./den(1);

Hz_LPF = tf(num,den,T_s);

The Hz_LPF is the same as c2d.

Hz_LPF =

  0.2391 z + 0.2391
  -----------------
     z - 0.5219

Sample time: 0.001 seconds
Discrete-time transfer function.

Explanation

subs

Symbolic substitution. subs(s,old,new) returns a copy of s, replacing all occurrences of old with new, and then evaluates s.

numden

Extract numerator and denominator. [N,D] = numden(A) converts A to a rational form where the numerator and denominator are relatively prime polynomials with integer coefficients. The function returns the numerator and denominator of the rational form of an expression.

If A is a symbolic or a numeric matrix, then N is the symbolic matrix of numerators, and D is the symbolic matrix of denominators. Both N and D are matrices of the same size as A.

coeffs

Coefficients of polynomial. C = coeffs(p) returns coefficients of the polynomial p with respect to all variables determined in p by symvar. The coefficients are ordered from the lowest degree to the highest degree. e.g.,

syms x
c = coeffs(16*x^2 + 19*x + 11)
>> c = [ 11, 19, 16]

Thus, a flip function should be used to convert the result to tf(num,den) format.