Implementing 1st Order Low Pass Filter

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)

Transfer Function

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:

Forward Euler

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

where

Derivation 1

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:

Derivation 2

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

Backward Euler

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

where

Derivation 1

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:

Derivation 2

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

Zero Order Hold (ZOH)

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

where

Derivation 1

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:

Derivation 2

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

Bilinear Transform

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

Derivation 1

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:

Derivation 2

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

Bilinear Transform with Prewarp Correction

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

where

Derivation

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.

Sample Code

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)

Notes

Above code was tested on below Python versions: