应用 FIR 滤波器¶
| 日期 | 2016-04-10（最后修改），2011-06-06（创建） |
|---|
如何应用 FIR 滤波器：convolve、fftconvolve、convolve1d 或 lfilter？¶
下图显示了使用 numpy 和 scipy 中提供的几种不同函数，将不同长度的有限脉冲响应 (FIR) 滤波器应用于长度为 131072 的信号所需的时间。下面详细介绍了如何创建此图。
fig
详细信息¶
numpy 和 scipy 库中有多个函数可用于将 FIR 滤波器 应用于信号。从 scipy.signal 中，lfilter() 用于将离散 IIR 滤波器 应用于信号，因此只需将分母系数数组设置为 [1.0]，即可用于应用 FIR 滤波器。应用 FIR 滤波器等效于离散 卷积，因此也可以使用 numpy 中的 convolve()、scipy.signal 中的 convolve() 或 fftconvolve()，或 scipy.ndimage 中的 convolve1d()。在本页中，我们将演示这些函数中的每一个，并查看当数据信号大小固定而 FIR 滤波器长度变化时，计算时间如何变化。我们将使用 131072 的数据信号长度，即 2**17。我们假设我们有 m 个数据通道，因此我们的输入信号是一个 m x n 数组。对于一维计时，我们使用大小为 2*10^4 的源数组，并使用不同数量的滤波器系数。
我们假设我们的 FIR 滤波器系数在一个一维数组 b 中。numpy.convolve 函数只接受一维数组，因此我们必须使用 python 循环遍历我们的输入数组，以对所有通道执行卷积。一种方法是
y = np.array([np.convolve(xi, b, mode='valid') for xi in x])
我们使用列表推导来遍历x的行，并将结果传递给np.array，以将过滤后的数据重新组合成二维数组。需要注意的是，在单维情况下，scipy.signal.convolve可能会使用np.convolve。
signal.convolve和signal.fftconvolve都执行二维数组的二维卷积。为了用这两个函数中的任何一个过滤我们的m乘n数组，我们将我们的滤波器整形为一个二维数组，形状为1乘len(b)。python代码如下
y = convolve(x, b[np.newaxis, :], mode='valid')
其中x是一个形状为(m, n)的numpy数组，b是FIR滤波器系数的单维数组。b[np.newaxis, :]是b的二维视图，形状为1乘len(b)。y是过滤后的数据；它只包含计算了完整卷积的那些项，因此它的形状为(m, n - len(b) + 1)。
ndimage.convolve1d() 被设计为对另一个n维数组的给定轴上的1d数组进行卷积。它没有选项mode='valid'，因此要提取结果的有效部分，我们对函数的结果进行切片
y = convolve1d(x, b)[:, (len(b)-1)//2 : -(len(b)//2)]
signal.lfilter 被设计为过滤一维数据。它可以接受一个二维数组（或者更一般地说，一个n维数组）并在任何给定的轴上过滤数据。它也可以用于IIR滤波器，因此在我们的例子中，我们将为分母系数传递[1.0]。在python中，它看起来像这样
y = lfilter(b, [1.0], x)
为了获得与convolve或fftconvolve计算出的完全相同的数组（即获得'valid'模式的等效结果），我们必须丢弃lfilter计算出的数组的开头。我们可以在调用filter后立即对数组进行切片来做到这一点
y = lfilter(b, [1.0], x)[:, len(b) - 1:]以下脚本计算并绘制将FIR滤波器应用于2乘131072数据数组的结果，其中使用了一系列长度递增的FIR滤波器。
#!python
%matplotlib inline
from __future__ import print_function
import time
import numpy as np
from numpy import convolve as np_convolve
from scipy.signal import fftconvolve, lfilter, firwin
from scipy.signal import convolve as sig_convolve
from scipy.ndimage import convolve1d
import matplotlib.pyplot as plt
# Create the m by n data to be filtered.
m = 1
n = 2 ** 18
x = np.random.random(size=(m, n))
conv_time = []
npconv_time = []
fftconv_time = []
conv1d_time = []
lfilt_time = []
diff_list = []
diff2_list = []
diff3_list = []
ntaps_list = 2 ** np.arange(2, 14)
for ntaps in ntaps_list:
# Create a FIR filter.
b = firwin(ntaps, [0.05, 0.95], width=0.05, pass_zero=False)
# --- signal.convolve ---
tstart = time.time()
conv_result = sig_convolve(x, b[np.newaxis, :], mode='valid')
conv_time.append(time.time() - tstart)
# --- numpy.convolve ---
tstart = time.time()
npconv_result = np.array([np_convolve(xi, b, mode='valid') for xi in x])
npconv_time.append(time.time() - tstart)
# --- signal.fftconvolve ---
tstart = time.time()
fftconv_result = fftconvolve(x, b[np.newaxis, :], mode='valid')
fftconv_time.append(time.time() - tstart)
# --- convolve1d ---
tstart = time.time()
# convolve1d doesn't have a 'valid' mode, so we expliclity slice out
# the valid part of the result.
conv1d_result = convolve1d(x, b)[:, (len(b)-1)//2 : -(len(b)//2)]
conv1d_time.append(time.time() - tstart)
# --- lfilter ---
tstart = time.time()
lfilt_result = lfilter(b, [1.0], x)[:, len(b) - 1:]
lfilt_time.append(time.time() - tstart)
diff = np.abs(fftconv_result - lfilt_result).max()
diff_list.append(diff)
diff2 = np.abs(conv1d_result - lfilt_result).max()
diff2_list.append(diff2)
diff3 = np.abs(npconv_result - lfilt_result).max()
diff3_list.append(diff3)
# Verify that np.convolve and lfilter gave the same results.
print("Did np.convolve and lfilter produce the same results?",)
check = all(diff < 1e-13 for diff in diff3_list)
if check:
print( "Yes.")
else:
print( "No! Something went wrong.")
# Verify that fftconvolve and lfilter gave the same results.
print( "Did fftconvolve and lfilter produce the same results?")
check = all(diff < 1e-13 for diff in diff_list)
if check:
print( "Yes.")
else:
print( "No! Something went wrong.")
# Verify that convolve1d and lfilter gave the same results.
print( "Did convolve1d and lfilter produce the same results?",)
check = all(diff2 < 1e-13 for diff2 in diff2_list)
if check:
print( "Yes.")
else:
print( "No! Something went wrong.")
def timeit(fn, shape, lfilter=False, n_x=2e4, repeats=3):
x = np.random.rand(int(n_x))
y = np.random.rand(*shape)
args = [x, y] if not lfilter else [x, x, y]
times = []
for _ in range(int(repeats)):
start = time.time()
c = fn(*args)
times += [time.time() - start]
return min(times)
npconv_time2, conv_time2, conv1d_time2 = [], [], []
fftconv_time2, sig_conv_time2, lconv_time2 = [], [], []
Ns_1d = 2*np.logspace(0, 4, num=11, dtype=int)
for n in Ns_1d:
npconv_time2 += [timeit(np_convolve, shape=(n,))]
conv1d_time2 += [timeit(convolve1d, shape=(n,))]
fftconv_time2 += [timeit(fftconvolve, shape=(n,))]
sig_conv_time2 += [timeit(sig_convolve, shape=(n,))]
lconv_time2 += [timeit(lfilter, shape=(n,), lfilter=True)]
fig = plt.figure(1, figsize=(16, 5.5))
plt.subplot(1, 2, 1)
plt.loglog(ntaps_list, conv1d_time, 'k-p', label='ndimage.convolve1d')
plt.loglog(ntaps_list, lfilt_time, 'c-o', label='signal.lfilter')
plt.loglog(ntaps_list, fftconv_time, 'm-*', markersize=8, label='signal.fftconvolve')
plt.loglog(ntaps_list[:len(conv_time)], conv_time, 'g-d', label='signal.convolve')
plt.loglog(ntaps_list, npconv_time, 'b-s', label='numpy.convolve')
plt.legend(loc='best', numpoints=1)
plt.grid(True)
plt.xlabel('Number of taps')
plt.ylabel('Time to filter (seconds)')
plt.title('Multidimensional timing')
plt.subplot(1, 2, 2)
plt.loglog(Ns_1d, conv1d_time2, 'k-p', label='ndimage.convolve1d')
plt.loglog(Ns_1d, lconv_time2, 'c-o', label='signal.lfilter')
plt.loglog(Ns_1d, fftconv_time2, 'm-*', markersize=8, label='signal.fftconvolve')
plt.loglog(Ns_1d, sig_conv_time2, 'g-d', label='signal.convolve')
plt.loglog(Ns_1d, npconv_time2, 'b-s', label='np.convolve')
plt.grid()
plt.xlabel('Number of taps')
plt.ylabel('Time to filter (seconds)')
plt.title('One dimensional timing')
plt.legend(loc='best')
plt.show()
该图显示，根据抽头的数量，numpy.convolve或scipy.signal.fftconvolve/scipy.signal.convolve是最快的。上面的脚本可以用来探索这些结果的变化。
需要注意的是，scipy.signal.convolve会根据输入数组的大小/形状选择使用scipy.signal.fftconvolve、np.convolve或scipy.signal.convolve中自然找到的卷积实现。
章节作者：WarrenWeckesser，WarrenWeckesser，Scott Sievert