python - 使用 FFT 计算滤波器 (b,a,x,zi)

Question

我想尝试使用 FFT 而不是计算 y=filter(b,a,x,zi) 和 dy[i]/dx[j]在时域中可能加速 GPU 实现。

我不确定这是否可能，尤其是当 zi 不为零时。我查看了 scipy 中的 scipy.signal.lfilter 和 octave 中的 filter 是如何实现的。它们都是直接在时域中完成的，scipy 使用直接形式 2 和 Octave 直接形式 1(通过查看 DLD-FUNCTIONS/filter.cc 中的代码)。我还没有在任何地方看到 FFT 实现类似于 MATLAB 中 FIR 滤波器的 fftfilt(即 a = [1.])。

我试过 y = ifft(fft(b)/fft(a) * fft(x)) 但这在概念上似乎是错误的。另外，我不确定如何处理初始 transient zi。任何引用，指向现有实现的指针，将不胜感激。

示例代码，

import numpy as np
import scipy.signal as sg
import matplotlib.pyplot as plt

# create an IRR lowpass filter
N = 5
b, a = sg.butter(N, .4)
MN = max(len(a), len(b))

# create a random signal to be filtered
T = 100
P = T + MN - 1
x = np.random.randn(T)
zi = np.zeros(MN-1)

# time domain filter
ylf, zo = sg.lfilter(b, a, x, zi=zi)

# frequency domain filter
af = sg.fft(a, P)
bf = sg.fft(b, P)
xf = sg.fft(x, P)
yfft = np.real(sg.ifft(bf/af * xf))[:T]

# error
print np.linalg.norm(yfft - ylf)

# plot, note error is larger at beginning and with larger N
plt.figure(1)
plt.clf()
plt.plot(ylf)
plt.plot(yfft)

Answer 1

您可以通过将 P = T + MN - 1 替换为 P = T + 2*MN - 1 来减少现有实现中的错误。这纯粹是直观的，但在我看来，由于环绕，bf 和 af 的划分将需要 2*MN 项。

C.S. Burrus 有一篇关于如何以面向 block 的方式考虑过滤(无论是 FIR 还是 IIR)的相当简洁的文章，here .我没有详细阅读它，但我认为它为您提供了通过卷积实现 IIR 滤波所需的方程式，包括中间状态。