Tag: Digital Filter

Band-pass Filter

Band-pass Filter

A band-pass filter attenuates both low and high frequency components of a signal and therefore, as the name suggests, allows for a band of frequencies to pass through. In this post, we will briefly go over an analog implementation of this type of filter by cascading high-pass and low-pass passive analog filters. I will also go through a digital implementation using Python, as seen on the post Joystick with Raspberry Pi.

The RC circuit diagram shows the high- and low-pass filters in a cascaded arrangement. The cutoff frequency (rad/s) for the high-pass is and for the low-pass. In terms of time constants (which will be helpful later with conciseness), we have respectively and .

For example, a band-pass filter with cutoff frequencies of approximately 0.1 and 2 Hz can be realized with, , , and . Or, using time constants, and .

That’s about it for the analog implementation of a first-order band-pass filter. Let’s see how we can derive the digital version of the same filter and have it ultimately implemented in Python code (at the bottom of this post if you want to skip the technical mumbo jumbo).

Continuous Filter

The first step is to come up with a mathematical representation of the band-pass filter in the analog (continuous) domain. More specifically, as already mentioned in the PID tuning post, let’s work with the Laplace transform (or s-domain), which conveniently turns differential equations into much easier to solve polynomial equations.

In the case of the RC circuit for the band-pass filter, we want to find the transfer function between the output voltageand the input voltage. In other words, between the raw and the filtered signals.

The relationships between voltage and current for the resistor and capacitor are given by:

and

In the s-domain, the derivative operatoris “replaced” with the algebraic s operator:

and

The filter can be separated into two filters (a high-pass and a low-pass), where distinct transfer functions betweenand, as well as betweenand , can be derived and then combined.

Solving the system’s first equation (below) for the current and substituting in the second equation gives:

And rearranging as a transfer function:

Solving the system’s first equation (below) for the current and substituting in the second equation gives:

And rearranging like the high-pass case:

Finally, solving the function for the low-pass filter forand substituting the result in the function for the high-pass filter, will lead to the transfer function of the band-pass filter below (also shown in terms of time constants and descending powers of s):

Discrete Filter

The next task is to convert the continuous system to a discrete one. Similarly to the continuous domain, where the derivative operatoris represented by the operator s, the discrete domain has the operator z. (or 1) corresponds to the current time step discrete value, corresponds to the previous value , to the value before the previous one, and so on.

As seen in the post about digital filtering, for a discrete system with sampling period, the derivative of a continuous variable can be approximated by the backward difference below:

And in terms of operators, we can write:

Substituting s in the continuous transfer function for the band-pass filter gives the equation in the z-domain below. Note also that, because we are talking about a filter for any type of signal, the voltages were replaced by more general input and output variables.

Carrying out the math so we can regroup in terms of descending powers of z, and multiplying both top and bottom by, leads to:

We can then rearrange the z transfer function above to resemble a difference equation:

Where:

At last, the equation above can be brought to a more useful difference equation form, by replacing the z operators by their corresponding discrete system values:

Python Code for Band-Pass Filter

The following Python program implements the digital filter with a fictitious signal comprised of three sinusoidal waves with frequencies 0.01, 0.5, and 25 Hz. The low and high cutoff frequencies of the filter are 0.1 and 2 Hz, therefore allowing only the 0.5 Hz component of the original signal to pass through.

By analyzing the code, it is straightforward to identify the filter coefficients determined in the previous section. The difference equation was implemented in the more general form

where the sets of coefficientsandare normalized by the coefficientof the output. This approach allows for any type of digital filter to be used in the code, as long as the normalized coefficients are defined.

# Importing modules and classes
import numpy as np
from utils import plot_line

# "Continuous" signal parameters
tstop = 20  # Signal duration (s)
Ts0 = 0.0002  # Time step (s)
fs0 = 1/Ts0  # Sampling frequency (Hz)
# Discrete signal parameters
tsample = 0.01  # Sampling period for code execution (s)
# Preallocating output arrays for plotting
tn = []
xn = []
yn = []

# Creating arbitrary signal with multiple sine functions
freq = [0.01, 0.5, 25]  # Sine frequencies (Hz)
ampl = [0.4, 1.0, 0.2]  # Sine amplitudes
t = np.arange(0, tstop+Ts0, Ts0)
xs = np.zeros(len(t))
for ai, fi in zip(ampl, freq):
    xs = xs + ai*np.sin(2*np.pi*fi*t)

# First order digital band-pass filter parameters
fc = np.array([0.1, 2])  # Band-pass cutoff frequencies (Hz)
tau = 1/(2*np.pi*fc)  # Filter time constants (s)
# Filter difference equation coefficients
a0 = tau[0]*tau[1]+(tau[0]+tau[1])*tsample+tsample**2
a1 = -(2*tau[0]*tau[1]+(tau[0]+tau[1])*tsample)
a2 = tau[0]*tau[1]
b0 = tau[0]*tsample
b1 = -tau[0]*tsample
# Defining normalized coefficients
a = np.array([1, a1/a0, a2/a0])
b = np.array([b0/a0, b1/a0])
# Initializing filter values
x = np.array([0.0]*(len(b)))  # x[n], x[n-1], x[n-2], ...
y = np.array([0.0]*(len(a)))  # y[n], y[n-1], y[n-2], ...

# Executing DAQ loop
tprev = 0
tcurr = 0
while tcurr <= tstop:
    # Doing I/O computations every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Simulating DAQ device signal acquisition
        x[0] = xs[int(tcurr/Ts0)]
        # Filtering signal
        y[0] = -np.sum(a[1::]*y[1::]) + np.sum(b*x)
        # Updating previous input values
        for i in range(len(b)-1, 0, -1):
                x[i] = x[i-1]
        # Updating previous output values
        for i in range(len(a)-1, 0, -1):
                y[i] = y[i-1]
    # Updating output arrays
    tn.append(tcurr)
    xn.append(x[0])
    yn.append(y[0])
    # Incrementing time step
    tprev = tcurr
    tcurr += Ts0

# Plotting results
plot_line(
    [t, tn], [xs, yn], yname=['X Input', 'Y Output'],
    legend=['Raw', 'Filtered'], figsize=(1300, 250)
)

The plot below shows the code output where the low frequency drift and the high frequency noise are almost fully removed from the signal, while the desired component of 0.5 Hz is preserved.

Digital Filtering

Digital Filtering

After a continuous signal goes through an analog-to-digital conversion, additional digital filtering can be applied to improve the signal quality. Whether the signal is being used in a real-time application or has been collected for a later analysis, implementing a digital filter via software is quite straightforward. We already saw them being casually used in some of my previous posts, such as MCP3008 with Raspberry Pi. Now, let’s go over these types of filters in more detail. For the examples, we will be using the signal processing functions from the the SciPy Python package.

One of the drawbacks of digital filtering (or analog filtering, for that matter) is the introduction of a phase delay in the filtered signal.

The plots on the left illustrate what an actual filtered signal and its phase lag would look like (top), in contrast with an ideal filtered signal with no lag (bottom). While the ideal situation cannot be achieved with real-time signal processing, it can be achieved when processing signals that are already sampled and stored as a discrete sequence of numerical values.

Since the sequence from start to finish is available, the idea is to filter it once (causing a phase shift) and then filter the resulting sequence backwards (causing a phase shift again). Since the filtering is done both forward and backwards, the phase shifts cancel each other, thus resulting in a zero-phase filtered signal!

To see how it’s done, check out the filtfilt function and the very good application examples on the SciPy documentation. You can also go for the code used to generate the plots above on my GitHub page. I should note that I use filtfilt all the time in my signal processing endeavors. There’s a lot of information that can be extracted from a signal after you clean it up a little bit. If the time alignment between events is critical, avoiding a phase shift is the way to go.

Before we get into the real-time application of digital filters, let’s talk briefly about how they’re implemented. I’ll focus on filters that can be described as the difference equation below:

whereis the input sequence (raw signal) andis the output sequence (filtered signal). Before you abandon this webpage, let me rewrite the equation so we can start bringing it to a more applicable form. Since the idea is to find the outputas a function of its previous values and the input, we can go with:

It’s common to normalize the coefficients so . Also, for simplicity, consider a first order filter which, in terms of difference equations, only depends on the values of the current and last time steps, or and . Thus:

So, to build a first-order low-pass digital filter, all we need to do is determine the coefficients , and. Luckily for us, the SciPy MATLAB-Style filter design functions return those coefficients, reducing our task to just the implementation of the filter using Python. Before we go over a couple of code examples, let’s examine a first-order filter that I use quite a bit.

Backward Difference Digital Filter

Going from the continuous domain to the discrete domain involves a transformation where we approximate a continuous variable by its discrete equivalent. I will start with a first-order low-pass filter, by analyzing its continuous behavior, more specifically the response to a unit step input.

The continuous system response to the constant input is given by , whereis the response time of the filter and is related to the filter cutoff frequencyby:

The response time of the filter is the time it takes for the output to reach approximately 63 % of the final value. In the case of the graph above, the response time is 0.1 seconds. If you’re wondering where that comes from, just calculate the system response using .

The differential equation that corresponds to the first order continuous system in question is:

And here is where our discrete approximation ofusing backward difference (with sampling period) comes into play:

By also transformingandinto their discrete counterpartsand, we can arrive to the difference equation below, which represents the first-order low-pass digital filter, where we used backward difference to approximate. Remember that and are the discrete current and previous time steps.

Finally, solving forgives us an equation very similar to the one we saw at the end of the previous section. Note that the coefficients andare a function of the sampling rate and the filter response time (in this particular case).

What I like about this filter implementation is that it’s fairly straightforward to see that the outputis a weighed average of its previous valueand the current input value. Furthermore, the smaller the response time (), the faster the filter is (higher cutoff frequency) and the more it follows the input value. Conversely, the slower the filter (), the more the output takes into account the previous value.

The Python code below shows how to implement this filter by placing it inside an execution loop and running a step input excitation through it. It will produce the plot shown at the beginning of this section. I invite you to experiment with different response times (cutoff frequencies) and sampling periods.

import numpy as np
import matplotlib.pyplot as plt

# Creating time array for "continuous" signal
tstop = 1  # Signal duration (s)
Ts0 = 0.001  # "Continuous" time step (s)
Ts = 0.02  # Sampling period (s)
t = np.arange(0, tstop+Ts0, Ts0)

# First order continuous system response to unit step input
tau = 0.1  # Response time (s)
y = 1 - np.exp(-t/tau)  # y(t)

# Preallocating signal arrays for digital filter
tf = []
yf = []

# Initializing previous and current values
xcurr = 1  # x[n] (step input)
yfprev = 0  # y[n-1]
yfcurr = 0  # y[n]

# Executing DAQ loop
tprev = 0
tcurr = 0
while tcurr <= tstop:
    # Doing filter computations every `Ts` seconds
    if (np.floor(tcurr/Ts) - np.floor(tprev/Ts)) == 1:
        yfcurr = tau/(tau+Ts)*yfprev + Ts/(tau+Ts)*xcurr
        yfprev = yfcurr
    # Updating output arrays
    tf.append(tcurr)
    yf.append(yfcurr)
    # Updating previous and current "continuous" time steps
    tprev = tcurr
    tcurr += Ts0

# Creating Matplotlib figure
fig = plt.figure(
    figsize=(6.3, 2.8),
    facecolor='#f8f8f8',
    tight_layout=True)
# Adding and configuring axes
ax = fig.add_subplot(
    xlim=(0, max(t)),
    xlabel='Time (s)',
    ylabel='Output ( - )',
    )
ax.grid(linestyle=':')
# Plotting signals
ax.plot(t, y, linewidth=1.5, label='Continuous', color='#1f77b4')
ax.plot(tf, yf, linewidth=1.5, label='Discrete', color='#ff7f0e')
ax.legend(loc='lower right')

SciPy Butterworth Digital Filter

The non-functional code snippet below shows how to import the SciPy signal processing module an use it to create a first-order digital Butterworth filter. Note that the implementation of the filter in the DAQ execution loop is very similar to what was done in the example code above. In this particular case however, you should be able to immediately identify the difference equation for the filter

where the arrays containing the coefficients , and, are the output of the SciPy function signal.butter for the corresponding cutoff frequency and discrete sampling period.

# Importing SciPy module
from scipy import signal

#
# Doing other initialization
#

# Discrete signal parameters
tsample = 0.01  # Sampling period for code execution (s)
fs = 1/tsample  # Sampling frequency (Hz)

# Finding a, b coefficients for digital butterworth filter
fc = 10  # Low-pass cutoff frequency (Hz)
b, a = signal.butter(1, fc, fs=fs)

# Preallocating variables
xcurr = 0  # x[n]
ycurr = 0  # y[n]
xprev = 0  # x[n-1]
yprev = 0  # y[n-1]

tprev = 0
tcurr = 0
# Executing DAQ loop
while tcurr <= tstop:
    # Doing I/O tasks every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Simulating DAQ device signal acquisition
        xcurr = mydaq.value
        # Filtering signal (digital butterworth filter)
        ycurr = -a[1]*yprev + b[0]*xcurr + b[1]*xprev
        yprev = ycurr
    # Updating previous values
    xprev = xcurr
    tprev = tcurr
    # Incrementing time step
    tcurr += Ts0

Finally, as I mentioned at the top, you can implement higher order or different types of filters. Take a look at the signal.iirnotch filter, which is a very useful design that removes specific frequencies from the original signal, such as the good old 60 Hz electric grid one.