How to write lowpass filter for sampled signal in Python?

Question

Sorry for question that may be trivial, but I am new in signal processing, and I didn't find python code with good explanation. So I will be glad to get explanation for dummy :)

I have some signal that sampled each 1 nsec (1e-9 sec) and have, let say, 1e4 points. I need to filter high frequencies from this signal. Let say I need to filter frequencies higher than 10MHz. I want that for frequencies lower than cutoff frequency signal will be passed unchanged. It means gain of the filter will be 1 for frequencies lower than cutoff frequency. I would like to be able to specify filter order. I mean, fitsr order filter have 20 db/decade slope (power rolloff) after cutoff frequency, second order filter have 40 db/dec slope after cutoff frequency and so on. High perfomance of code is importent.

Thank you very much for help. Alex.

Answer 1

The frequency response for the filter designed using the butter function is:

But there is no reason to limit the filter to a constant monotonic filter design. If you desire a higher attenuation in the stopband and steeper transition band, other options exist. For more information on specifying a filter using iirdesing see this. As shown by the frequency response plots for the butter design the cutoff frequency (-3dB point) is far from the goal. This can be alleviated by down-sampling before filtering (the design functions will have a difficult time with such a narrow filter, 2% of the bandwidth). Lets look at filtering the original sample rate with the cutoff specified.

import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

from scipy.signal import fir_filter_design as ffd
from scipy.signal import filter_design as ifd

# setup some of the required parameters
Fs = 1e9           # sample-rate defined in the question, down-sampled

# remez (fir) design arguements
Fpass = 10e6       # passband edge
Fstop = 11.1e6     # stopband edge, transition band 100kHz
Wp = Fpass/(Fs)    # pass normalized frequency
Ws = Fstop/(Fs)    # stop normalized frequency

# iirdesign agruements
Wip = (Fpass)/(Fs/2)
Wis = (Fstop+1e6)/(Fs/2)
Rp = 1             # passband ripple
As = 42            # stopband attenuation

# Create a FIR filter, the remez function takes a list of 
# "bands" and the amplitude for each band.
taps = 4096
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000) 

# The iirdesign takes passband, stopband, passband ripple, 
# and stop attenuation.
bc, ac = ifd.iirdesign(Wip, Wis, Rp, As, ftype='ellip')  
bb, ab = ifd.iirdesign(Wip, Wis, Rp, As, ftype='cheby2') 

As mentioned, because we are trying to filter such a small percent of the bandwidth the filter will not have a sharp cutoff. In this case, lowpass filter, we can reduce the bandwidth to get a better looking filter. The python/scipy.signal resample function can be used to reduce the bandwidth.

Note the resample function will perform filtering to prevent aliasing. Prefiltering can also be perfomed (to reduce aliasing) and in this case we could simply resample by 100 and be done, but the question asked about creating filters. For this example we will downsample by 25 and create a new filter

R = 25;            # how much to down sample by
Fsr = Fs/25.       # down-sampled sample rate
xs = signal.resample(x, len(x)/25.)

If we update the design parameters for the FIR filter the new response is.

# Down sampled version, create new filter and plot spectrum
R = 25.             # how much to down sample by
Fsr = Fs/R          # down-sampled sample rate
Fstop = 11.1e6      # modified stopband
Wp = Fpass/(Fsr)    # pass normalized frequency
Ws = Fstop/(Fsr)    # stop normalized frequency
taps = 256
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000) 

The filter operating on the downsampled data has a better response. Another benefit of using a FIR filter is that you will have linear phase response.

Answer 2

Does this work?

from __future__ import division
from scipy.signal import butter, lfilter

fs = 1E9 # 1 ns -> 1 GHz
cutoff = 10E6 # 10 MHz
B, A = butter(1, cutoff / (fs / 2), btype='low') # 1st order Butterworth low-pass
filtered_signal = lfilter(B, A, signal, axis=0)

You're right, though, the documentation is not very complete. It looks like butter is a wrapper for iirfilter, which is better documented:

N : int The order of the filter. Wn : array_like A scalar or length-2 sequence giving the critical frequencies.

Most of this stuff is cloned from matlab, though, so you can look at their documentation too:

the normalized cutoff frequency Wn must be a number between 0 and 1, where 1 corresponds to the Nyquist frequency, π radians per sample.

Update:

I added documentation for these functions. :) Github makes it easy.

Answer 3

Not sure what your application is, but you may want to check out Gnuradio: http://gnuradio.org/doc/doxygen/classgr__firdes.html

The signal processing blocks are written in C++ (although the Gnuradio flow graphs are in Python), but you did say high performance is important.

Answer 4

I'm having good results with this FIR filter. Notices it applies the filter twice, going "forward" and "reverse", so as to compensate for signal offset (filtfilt function didn't work, don't know why):

def firfilt(interval, freq, sampling_rate):
    nfreq = freq/(0.5*sampling_rate)
    taps =  sampling_rate + 1
    a = 1
    b = scipy.signal.firwin(taps, cutoff=nfreq)
    firstpass = scipy.signal.lfilter(b, a, interval)
    secondpass = scipy.signal.lfilter(b, a, firstpass[::-1])[::-1]
    return secondpass

A great resource to filter design and use, from where I took this code, and from where band-pass and hi-pass filter examples can be taken, is THIS.

Answer 5

I do not have comment rights ...

@endolith: I use the same as you except using the scipy.signal.filtfilt(B, A, x) where x is the input vector to be filtered - e.g. numpy.random.normal(size=(N)). filtfilt makes a forward and reverse pass of the signal. For the sake of completeness (most being the same as @endolith):

import numpy as np
import scipy.signal as sps

input = np.random.normal(size=(N)) # Random signal as example
bz, az = sps.butter(FiltOrder, Bandwidth/(SamplingFreq/2)) # Gives you lowpass Butterworth as default
output = sps.filtfilt(bz, az, input) # Makes forward/reverse filtering (linear phase filter)

filtfilt as also suggested by @heltonbiker requires arrays of coefficients I believe. In case you need to perform bandpass filtering at complex baseband a more involved config is needed but this does not appear to be a problem here.