Matlab's bandpass filter equivalent in python

Question

In matlab there is a function called bandpass that I often use. The doc of the function can be found here: https://ch.mathworks.com/help/signal/ref/bandpass.html

I am looking for a way to apply a bandpass filter in Python and get the same or almost the same output filtered signal.

My signal can be downloaded from here: https://gofile.io/?c=JBGVsH

Matlab code:

load('mysignal.mat')
y = bandpass(x, [0.015,0.15], 1/0.7);
plot(x);hold on; plot(y)

Python code:

import matplotlib.pyplot as plt
import scipy.io
from scipy.signal import butter, lfilter

x = scipy.io.loadmat("mysignal.mat")['x']

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a


def butter_bandpass_filter(data, lowcut, highcut, fs, order=6):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

y = butter_bandpass_filter(x, 0.015, 0.15, 1/0.7, order=6)

plt.plot(x);plt.plot(y);plt.show()

I need to find a way in python to apply similar filtering as in the Matlab example code block.

Gilstrap

Your butter_bandpass function is not using fs, for one... – Purificator

That said, isn't the huge DC offset of your original signal a problem for Matlab? – Purificator

In matlab the desired output is achieved as shown in the first plot. For the butter_bandpass i indeed pass fs as input argument. – Gilstrap

You don't use fs within butter_bandpass(). – Purificator

fixed but same results :/ – Gilstrap

Well, your code is creating a fifth-order filter, and Matlab's bandpass, quote, "uses a minimum-order filter with a stopband attenuation of 60 dB"... you'd probably have to generate such a filter instead. Since you do have Matlab, you could probably also peek at the bandpass() function's source? – Purificator

Answer 1

My favorite solution is here on Creating lowpass filter in SciPy - understanding methods and units which I changed to a band-pass example:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
""" 
bandpass example, from
https://mcmap.net/q/194580/-creating-lowpass-filter-in-scipy-understanding-methods-and-units
pl, 22.05.2022
"""

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt


def butter_bandpass(lowcut, highcut, fs, order=5):
    b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=6):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

# Filter requirements.
order = 6
fs = 60.0       # sample rate, Hz
#cutoff = 3.667  # desired cutoff frequency of the filter, Hz
cutoffLo = 2.0
cutoffHi = 4.0

# Demonstrate the use of the filter.
# First make some data to be filtered.
T = 5.0         # seconds
n = int(T * fs) # total number of samples
t = np.linspace(0, T, n, endpoint=False)
# "Noisy" data.  We want to recover the 1.2 Hz signal from this.
fl=1.2      # Hz
f0=3.0      # Hz
fh=7.0      # Hz
al=2.5
a0=1.0
ah=3.5
w=2*np.pi   # Omega

data = f0 * np.sin(w*f0*t) \
     + fl * np.sin(w*fl*t) \
     + fh * np.sin(w*fh*t) \
     + 8.0

# Get the filter coefficients so we can check its frequency response.
b, a = butter_bandpass(cutoffLo, cutoffHi, fs, order)

# Plot the frequency response.
w, h = freqz(b, a, fs=fs, worN=8000)
plt.subplot(2, 1, 1)
plt.plot(w, np.abs(h), 'b')
plt.plot(cutoffLo, 0.5*np.sqrt(2), 'ko')
plt.plot(cutoffHi, 0.5*np.sqrt(2), 'ko')
plt.axvline(cutoffLo, color='k')
plt.axvline(cutoffHi, color='k')
plt.axvline(fl, color='r',linestyle=":")
plt.axvline(f0, color='r',linestyle="-.")
plt.axvline(fh, color='r',linestyle=":")
plt.xlim(0, 2.0*cutoffHi)
plt.title("Bandpass Filter Frequency Response")
plt.xlabel('Frequency [Hz]')
plt.grid()


# Filter the data, and plot both the original and filtered signals.
y = butter_bandpass_filter(data, cutoffLo, cutoffHi, fs, order)

plt.subplot(2, 1, 2)
plt.plot(t, data, 'b-', label='data')
plt.plot(t, y, 'g-', linewidth=2, label='filtered data')
plt.xlabel('Time [sec]')
plt.grid()
plt.legend()

plt.subplots_adjust(hspace=0.35)
plt.show()

Lientery