How to apply a FIR filter: convolve, fftconvolve, convolve1d or lfilter?
The following plot shows the time required to apply a finite impulse response (FIR) filter of varying length to a signal of length 131072 using several different functions that are available in numpy and scipy. The details of how this figure was created are given below.
The Details
There are several functions in the numpy and scipy libraries that can be used to apply a FIR filter to a signal. From scipy.signal, lfilter() is designed to apply a discrete IIR filter to a signal, so by simply setting the array of denominator coefficients to [1.0], it can be used to apply a FIR filter. Applying a FIR filter is equivalent to a discrete convolution, so one can also use convolve() from numpy, convolve() or fftconvolve() from scipy.signal, or convolve1d() from scipy.ndimage. In this page, we demonstrate each of these functions, and we look at how the computational time varies when the data signal size is fixed and the FIR filter length is varied. We'll use a data signal length of 131072, which is 2**17. We assume that we have m channels of data, so our input signal is an m by n array.
We assume our FIR filter coefficients are in a one-dimensional array b. The function numpy.convolve only accepts one-dimensional arrays, so we'll have to use a python loop over our input array to perform the convolution over all the channels. One way to do that is
np.array([np.convolve(xi, b, mode='valid') for xi in x])
We use a list comprehension to loop over the rows of x, and pass the result to np.array to reassemble the filtered data into a two-dimensional array.
Both signal.convolve and signal.fftconvolve perform a two-dimensional convolution of two-dimensional arrays. To filter our m by n array with either of these functions, we shape our filter to be a two-dimensional array, with shape 1 by len(b). The python code looks like this:
y = convolve(x, b[np.newaxis, :], mode='valid')
where x is a numpy array with shape (m, n), and b is the one-dimensional array of FIR filter coefficients. b[np.newaxis, :] is the two dimensional view of b, with shape 1 by len(b). y is the filtered data; it includes only those terms for which the full convolution was computed, so it has shape (m, n - len(b) + 1).
ndimage.convolve1d() is designed to do a convolution of a 1d array along the given axis of another n-dimensional array. It does not have the option mode='valid', so to extract the valid part of the result, we slice the result of the function:
y = convolve1d(x, b)[:, (len(b)-1)//2 : -(len(b)//2)]
signal.lfilter is designed to filter one-dimensional data. It can take a two-dimensional array (or, in general, an n-dimensional array) and filter the data in any given axis. It can also be used for IIR filters, so in our case, we'll pass in [1.0] for the denominator coefficients. In python, this looks like:
y = lfilter(b, [1.0], x)
To obtain exactly the same array as computed by convolve or fftconvolve (i.e. to get the equivalent of the 'valid' mode), we must discard the beginning of the array computed by lfilter. We can do this by slicing the array immediately after the call to filter:
y = lfilter(b, [1.0], x)[:, len(b) - 1:]
The following script computes and plots the results of applying a FIR filter to a 2 by 131072 array of data, with a series of FIR filters of increasing length.
1 import time
2
3 import numpy as np
4 from numpy import convolve as np_convolve
5 from scipy.signal import convolve as sig_convolve, fftconvolve, lfilter, firwin
6 from scipy.ndimage import convolve1d
7 from pylab import grid, show, legend, loglog, xlabel, ylabel, figure
8
9 # Create the m by n data to be filtered.
10 m = 4
11 n = 2 ** 17
12 x = np.random.random(size=(m, n))
13
14 conv_time = []
15 npconv_time = []
16 fftconv_time = []
17 conv1d_time = []
18 lfilt_time = []
19
20 diff_list = []
21 diff2_list = []
22 diff3_list = []
23
24 ntaps_list = 2 ** np.arange(2, 13)
25
26 for ntaps in ntaps_list:
27 # Create a FIR filter.
28 b = firwin(ntaps, [0.05, 0.95], width=0.05, pass_zero=False)
29
30 if ntaps <= 2 ** 9:
31 # --- signal.convolve ---
32 # We know this is slower than the others when ntaps is
33 # large, so we only compute it for small values.
34 tstart = time.time()
35 conv_result = sig_convolve(x, b[np.newaxis, :], mode='valid')
36 conv_time.append(time.time() - tstart)
37
38 # --- numpy.convolve ---
39 tstart = time.time()
40 npconv_result = np.array([np_convolve(xi, b, mode='valid') for xi in x])
41 npconv_time.append(time.time() - tstart)
42
43 # --- signal.fftconvolve ---
44 tstart = time.time()
45 fftconv_result = fftconvolve(x, b[np.newaxis, :], mode='valid')
46 fftconv_time.append(time.time() - tstart)
47
48 # --- convolve1d ---
49 tstart = time.time()
50 # convolve1d doesn't have a 'valid' mode, so we expliclity slice out
51 # the valid part of the result.
52 conv1d_result = convolve1d(x, b)[:, (len(b)-1)//2 : -(len(b)//2)]
53 conv1d_time.append(time.time() - tstart)
54
55 # --- lfilter ---
56 tstart = time.time()
57 lfilt_result = lfilter(b, [1.0], x)[:, len(b) - 1:]
58 lfilt_time.append(time.time() - tstart)
59
60 diff = np.abs(fftconv_result - lfilt_result).max()
61 diff_list.append(diff)
62
63 diff2 = np.abs(conv1d_result - lfilt_result).max()
64 diff2_list.append(diff2)
65
66 diff3 = np.abs(npconv_result - lfilt_result).max()
67 diff3_list.append(diff3)
68
69 # Verify that np.convolve and lfilter gave the same results.
70 print "Did np.convolve and lfilter produce the same results?",
71 check = all(diff < 1e-13 for diff in diff3_list)
72 if check:
73 print "Yes."
74 else:
75 print "No! Something went wrong."
76
77 # Verify that fftconvolve and lfilter gave the same results.
78 print "Did fftconvolve and lfilter produce the same results?",
79 check = all(diff < 1e-13 for diff in diff_list)
80 if check:
81 print "Yes."
82 else:
83 print "No! Something went wrong."
84
85 # Verify that convolve1d and lfilter gave the same results.
86 print "Did convolve1d and lfilter produce the same results?",
87 check = all(diff2 < 1e-13 for diff2 in diff2_list)
88 if check:
89 print "Yes."
90 else:
91 print "No! Something went wrong."
92
93 figure(1, figsize=(8, 5.5))
94 loglog(ntaps_list, npconv_time, 'c-s', label='numpy.convolve')
95 loglog(ntaps_list, conv1d_time, 'k-p', label='ndimage.convolve1d')
96 loglog(ntaps_list, fftconv_time, 'g-*', markersize=8, label='signal.fftconvolve')
97 loglog(ntaps_list[:len(conv_time)], conv_time, 'm-d', label='signal.convolve')
98 loglog(ntaps_list, lfilt_time, 'b-o', label='signal.lfilter')
99 legend(loc='best', numpoints=1)
100 grid(True)
101 xlabel('Number of taps')
102 ylabel('Time to filter (seconds)')
103 show()
The plot shows that, depending on the number of taps, either scipy.ndimage.convolve1d, numpy.convolve or scipy.signal.fftconvolve is the fastest. The above script can be used to explore variations of these results.