使用FFTs计算滤波器（b，a，x，zi）

Question

我想尝试使用fft而不是时域来计算y=filter(b,a,x,zi)和{}，以便在GPU实现中加速。在

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

我试着做y = ifft(fft(b) / fft(a) * fft(x))，但这似乎在概念上是错误的。另外，我不知道如何处理初始瞬态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和{}的划分将需要2*MN项，因为它是环绕的。在

C.S.Burrus有一篇非常简洁的文章，介绍了如何以面向块的方式来看待滤波，不管是FIR还是IIR。我还没有详细读过，但我想它给出了通过卷积实现IIR滤波所需的方程，包括中间状态。在