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.

## 0 Comments