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 fc (in Hz) by
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:
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
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: