I have some signal that sampled each 1 ns (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 10 MHz. 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, first order filter have 20 db/decade slope (power roll off) after cutoff frequency, second order filter have 40 db/dec slope after cutoff frequency and so on. High performance of code is important.
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.
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)
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.
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.
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.
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.