This page provides code and explanation on how to implement a first order low pass filter of the form,

Five implementations are provided. There are different ways to implement it because there are multiple discretization methods. Example Python code is provided for all methods. Below is an example of smoothing a noisy signal using one of these methods:

Python code for the above is as follows:

import numpy as np import matplotlib.pyplot as plt dt = 0.01 # Time step [sec] t = np.arange(0, 2, dt) # Input signal (curve with noise) u = 20*np.sin(t/1.6*2*np.pi) + 0.2*t + 50 u = np.random.normal(u, 15) # Add noise with standard deviation 15 # Output signal y = np.empty(len(t)) y[0] = u[0] # Initialize first sample # Low pass filter tau = 0.24 # Time constant [sec] alpha = dt / (tau + dt) for k in range(1, len(t)): y[k] = (1 - alpha) * y[k-1] + alpha * u[k] plt.plot(t, u) plt.plot(t, y) plt.show()

(Note: the plot will differ on every run due to the use of random number generator)

Before diving into the code, it is necessary to first define parameters that characterize a first order low pass filter. Continuous-time transfer function for a first order low pass filter can be written as follows:

where *ω _{c}* is the angular cutoff frequency. This is the frequency
at which the gain is 1/sqrt(2) of the DC gain (-3dB or about half the power).
This parameter describes which frequencies are passed and which are reduced.

Angular cutoff frequency *ω _{c}* (in rad/sec) is related to the
cutoff frequency

Another way to write the transfer function is:

where *τ* is the time constant in seconds. This is the time it takes for the system
to reach 63.2% of the final value to a step input. This parameter describes
how fast the output responds to a change in input.

The time constant and angular cutoff freqneucy are related by,

Finally, the discretization time step will be denoted *Δt* (in seconds). The
sampling frequency is then given by

It is not possible to remove frequencies above half the sampling rate
(Nyquist freqneucy) so *Δt* must be small enough and *ω _{c}*
also small enough to satisfy:

If *y* is the output and *u* is the input, the code to filter is.

where

Continuous transfer function can be discretized with forward Euler
method by replacing *s* with:

Doing this to *H(s)* and simplifying, we get,

Let *U(z)* be the input and *Y(z)* be the output,

Take the inverse Z-transform and we arrive at the difference equation:

Another way to reach the same result is by using state space method.
Let *Y(s)* be the output and *U(S)* the input. Then

Then, using the fact that

or

It can be discretize by forward Euler method:

where

so

If *y* is the output and *u* is the input, the code to filter is,

where

Continuous transfer function can be discretized with backward Euler
method by replacing *s* with:

Doing this to *H(s)* and simplifying, we get,

Let *U(z)* be the input and *Y(z)* be the output,

Take the inverse Z-transform and we arrive at the difference equation:

Another way to reach the same result is by using state space method.
Let *Y(s)* be the output and *U(S)* the input. Then

Then, using the fact that

or

It can be discretize by backward Euler method:

where

so

If *y* is the output and *u* is the input, the code to filter is,

where

Continuous transfer function can be discretized with zero-order hold using the following formula:

Substituting *H(s)* and simplifying, we get,

Let *U(z)* be the input and *Y(z)* be the output,

Take the inverse Z-transform and we arrive at the difference equation:

Another way to reach the same result is by using state space method.
Let *Y(s)* be the output and *U(S)* the input. Then

Then, discretize by zero-order hold:

where

so

This is also called trapezoidal approximation or Tustin method.
If *y* is the output and *u* is the input, the code to filter is,

where

Continuous transfer function can be discretized with bilinear transform
by replacing *s* with:

Doing this to *H(s)* and simplifying, we get,

Let *U(z)* be the input and *Y(z)* be the output,

Take the inverse Z-transform and we arrive at the difference equation:

*Y(s)* be the output and *U(S)* the input. Then

Then, using the fact that

It can be discretize by:

where

so

If *y* is the output and *u* is the input, the code to filter is,

where

Continuous transfer function can be discretized with bilinear transform
by replacing *s* with:

The rest of the derivation is similar to that of previous section.

import numpy as np import matplotlib.pyplot as plt dt = 0.01 # Time step [sec] t = np.arange(0, 2, dt) # Time [sec] # Input signal (curve with noise) u = 20*np.sin(t/1.6*2*np.pi) + 20*t + 50 u = np.random.normal(u, 15) # Add noise to input signal # Time constant [sec] tau = 0.24 def forward_euler(): alpha = dt / tau y = np.empty(len(t)) y[0] = u[0] for k in range(1, len(t)): y[k] = (1 - alpha) * y[k-1] + alpha * u[k-1] return y def backward_euler(): alpha = dt/ (tau + dt) y = np.empty(len(t)) y[0] = u[0] for k in range(1, len(t)): y[k] = (1 - alpha) * y[k-1] + alpha * u[k] return y def zero_order_hold(): alpha = np.exp(-dt / tau) y = np.empty(len(t)) y[0] = u[0] for k in range(1, len(t)): y[k] = alpha * y[k-1] + (1-alpha) * u[k-1] return y def bilinear_transform(): alpha = (2*tau - dt)/(2*tau + dt) beta = dt / (2*tau + dt) y = np.empty(len(t)) y[0] = u[0] for k in range(1, len(t)): y[k] = alpha * y[k-1] + beta * (u[k] + u[k-1]) return y def bilinear_transform_with_prewarp(): c = np.tan(dt / (2 * tau)) alpha = (1 - c)/(1 + c) beta = c/(1 + c) y = np.empty(len(t)) y[0] = u[0] for k in range(1, len(t)): y[k] = alpha * y[k-1] + beta * (u[k] + u[k-1]) return y omega_c = 1/tau # Cutoff angular frequency [rad/sec] f_c = omega_c / (2 * np.pi) # Cutoff frequency [Hz] f_s = 1 / dt # Sampling frequency if f_c >= 0.5 * f_s: raise Exception('Cutoff frequency is above Nyquist frequency') y1 = forward_euler() y2 = backward_euler() y3 = zero_order_hold() y4 = bilinear_transform() y5 = bilinear_transform_with_prewarp() plt.plot(t, u) plt.plot(t, y1) plt.plot(t, y2) plt.plot(t, y3) plt.plot(t, y4) plt.plot(t, y5) plt.show()

Output:

(Note: the plot will differ on every run due to the use of random number generator)

Above code was tested on below Python versions:

- Python 3.11.1
- numpy 1.26.4
- matplotlib 3.7.1