Digital Low-Pass Filter

In this section, we are going to do the following activities:

  • discretize a continuous-time low-pass filter by using bilinear transform (trapezoidal or Tustin) method

  • implement the discretized low-pass filter into MATLAB Simulink

  • compare the results from the implemented discrete low-pass filter to the shipped discrete low-pass filter in MATLAB Simulink

Required Imports

from IPython.core.display import HTML
from sympy import *
from mathprint import *

Ts, tau = symbols('T_s tau', positive=True)
s = symbols('s', complex=True)
z = symbols('z')
wc = symbols('omega_c', positive=True)

x, y = symbols('x y')
x0, x1, x2, x3 = symbols('x_k x_{k-1} x_{k-2} x_{k-3}')
y0, y1, y2, y3 = symbols('y_k y_{k-1} y_{k-2} y_{k-3}')

First-Order Low-Pass Filter

Transfer function of a first order lowpass filter (time-constant filter):

\[ H(s) = \frac{1} {\tau s+1} \]

where \(\tau\) is the filter time constant (in seconds).

Discretization with Bilinear Transformation

Next, we transform \(s\) into \(z\) by applying the following substitution.

\[ \frac{1}{s} \longleftarrow \frac{T_s}{2} \frac{z+1}{z-1} \]
H = 1 / (tau*s+1)
mprint('H=',latex(H))

H = H.subs(1/s, Ts/2 * (z+1)/(z-1))
mprint('H=',latex(H))
\[\displaystyle H=\frac{1}{s \tau + 1}\]
\[\displaystyle H=\frac{1}{1 + \frac{2 \tau \left(z - 1\right)}{T_{s} \left(z + 1\right)}}\]

Let us define \(x\) as the input to the filter and \(y\) as the output (filtered input).

eq = Eq(y, H * x)
mprint(latex(eq))

eq = simplify(eq)
mprint(latex(eq))

eq = Eq(numer(eq.rhs), expand(eq.lhs * denom(eq.rhs)))
mprint(latex(eq))

eq =expand(Eq(numer(eq.rhs)/z/Ts, eq.lhs * denom(eq.rhs)/z/Ts))
mprint(latex(eq))
\[\displaystyle y = \frac{x}{1 + \frac{2 \tau \left(z - 1\right)}{T_{s} \left(z + 1\right)}}\]
\[\displaystyle y = \frac{T_{s} x \left(z + 1\right)}{T_{s} \left(z + 1\right) + 2 \tau \left(z - 1\right)}\]
\[\displaystyle T_{s} x \left(z + 1\right) = T_{s} y z + T_{s} y + 2 \tau y z - 2 \tau y\]
\[\displaystyle y + \frac{y}{z} + \frac{2 \tau y}{T_{s}} - \frac{2 \tau y}{T_{s} z} = x + \frac{x}{z}\]

Apply the following substitutions:

  • \(y\) becomes \(y_{k}\)

  • \(y/z\) becomes \(y_{k-1}\)

  • \(x\) becomes \(x_{k}\)

  • \(x/z\) becomes \(x_{k-1}\)

eq = eq.subs(x/z**3, x3).subs(x/z**2, x2).subs(x/z, x1).subs(x, x0).subs(y/z**3, y3).subs(y/z**2, y2).subs(y/z, y1).subs(y, y0)
mprint("\\small ", latex(eq))
\[\displaystyle \small y_{k} + y_{k-1} + \frac{2 \tau y_{k}}{T_{s}} - \frac{2 \tau y_{k-1}}{T_{s}} = x_{k} + x_{k-1}\]

Finally, by grouping the variables, we obtain:

eq = Eq(collect(eq.rhs, [x0, x1, x2, x3]), collect(eq.lhs, [y0, y1, y2, y3]) )
mprintb("\\small ", latex(eq))
\[\displaystyle \boxed{\small x_{k} + x_{k-1} = y_{k} \left(1 + \frac{2 \tau}{T_{s}}\right) + y_{k-1} \left(1 - \frac{2 \tau}{T_{s}}\right)}\]

Second Order Low-Pass Butterworth Filter

Transfer function of a second order lowpass filter (Butterworth):

\[ H(s) = \frac{\omega_{c}^{2}}{\omega_{c}^{2} + \sqrt{2} \omega_{c} s + s^{2}} \]

where \(\omega_c\) is the filter cut-off frequency (in Hz).

Discretization with Bilinear Transformation

Similiar to the previous section, here we also transform \(s\) into \(z\) by applying the following substitution.

\[ \frac{1}{s} \longleftarrow \frac{T_s}{2} \frac{z+1}{z-1} \]
H = wc**2 / (wc**2+sqrt(2)*s*wc+s**2)
mprint('H=',latex(H))

H = H.subs(1/s, Ts/2 * (z+1)/(z-1))
mprint('H=',latex(H))
\[\displaystyle H=\frac{\omega_{c}^{2}}{\omega_{c}^{2} + \sqrt{2} \omega_{c} s + s^{2}}\]
\[\displaystyle H=\frac{\omega_{c}^{2}}{\omega_{c}^{2} + \frac{2 \sqrt{2} \omega_{c} \left(z - 1\right)}{T_{s} \left(z + 1\right)} + \frac{4 \left(z - 1\right)^{2}}{T_{s}^{2} \left(z + 1\right)^{2}}}\]

Let us define \(x\) as the input to the filter and \(y\) as the output (filtered input).

eq = Eq(y, H * x)
mprint(latex(eq))

eq = simplify(eq)
mprint(latex(eq))

eq = Eq(numer(eq.rhs), expand(eq.lhs * denom(eq.rhs)))
mprint(latex(eq))

eq =expand(Eq(numer(eq.rhs)/z**2/Ts**2/wc**2, eq.lhs * denom(eq.rhs)/z**2/Ts**2/wc**2))
mprint(latex(eq))
\[\displaystyle y = \frac{\omega_{c}^{2} x}{\omega_{c}^{2} + \frac{2 \sqrt{2} \omega_{c} \left(z - 1\right)}{T_{s} \left(z + 1\right)} + \frac{4 \left(z - 1\right)^{2}}{T_{s}^{2} \left(z + 1\right)^{2}}}\]
\[\displaystyle y = \frac{T_{s}^{2} \omega_{c}^{2} x \left(z + 1\right)^{2}}{T_{s}^{2} \omega_{c}^{2} \left(z + 1\right)^{2} + 2 \sqrt{2} T_{s} \omega_{c} \left(z - 1\right) \left(z + 1\right) + 4 \left(z - 1\right)^{2}}\]
\[\displaystyle T_{s}^{2} \omega_{c}^{2} x \left(z + 1\right)^{2} = T_{s}^{2} \omega_{c}^{2} y z^{2} + 2 T_{s}^{2} \omega_{c}^{2} y z + T_{s}^{2} \omega_{c}^{2} y + 2 \sqrt{2} T_{s} \omega_{c} y z^{2} - 2 \sqrt{2} T_{s} \omega_{c} y + 4 y z^{2} - 8 y z + 4 y\]
\[\displaystyle y + \frac{2 y}{z} + \frac{y}{z^{2}} + \frac{2 \sqrt{2} y}{T_{s} \omega_{c}} - \frac{2 \sqrt{2} y}{T_{s} \omega_{c} z^{2}} + \frac{4 y}{T_{s}^{2} \omega_{c}^{2}} - \frac{8 y}{T_{s}^{2} \omega_{c}^{2} z} + \frac{4 y}{T_{s}^{2} \omega_{c}^{2} z^{2}} = x + \frac{2 x}{z} + \frac{x}{z^{2}}\]

Apply the following substitutions:

  • \(y\) becomes \(y_{k}\)

  • \(y/z\) becomes \(y_{k-1}\)

  • \(y/z^2\) becomes \(y_{k-2}\)

  • \(x\) becomes \(x_{k}\)

  • \(x/z\) becomes \(x_{k-1}\)

  • \(x/z^2\) becomes \(x_{k-2}\)

eq = eq.subs(x/z**3, x3).subs(x/z**2, x2).subs(x/z, x1).subs(x, x0).subs(y/z**3, y3).subs(y/z**2, y2).subs(y/z, y1).subs(y, y0)
mprint(latex(eq))
\[\displaystyle y_{k} + 2 y_{k-1} + y_{k-2} + \frac{2 \sqrt{2} y_{k}}{T_{s} \omega_{c}} - \frac{2 \sqrt{2} y_{k-2}}{T_{s} \omega_{c}} + \frac{4 y_{k}}{T_{s}^{2} \omega_{c}^{2}} - \frac{8 y_{k-1}}{T_{s}^{2} \omega_{c}^{2}} + \frac{4 y_{k-2}}{T_{s}^{2} \omega_{c}^{2}} = x_{k} + 2 x_{k-1} + x_{k-2}\]

Finally, by grouping the variables, we obtain:

eq = Eq(collect(eq.rhs, [x0, x1, x2, x3]), collect(eq.lhs, [y0, y1, y2, y3]) )
mprintb(latex(eq))
\[\displaystyle \boxed{x_{k} + 2 x_{k-1} + x_{k-2} = y_{k} \left(1 + \frac{2 \sqrt{2}}{T_{s} \omega_{c}} + \frac{4}{T_{s}^{2} \omega_{c}^{2}}\right) + y_{k-1} \left(2 - \frac{8}{T_{s}^{2} \omega_{c}^{2}}\right) + y_{k-2} \left(1 - \frac{2 \sqrt{2}}{T_{s} \omega_{c}} + \frac{4}{T_{s}^{2} \omega_{c}^{2}}\right)}\]