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

我不确定是否有可能,尤其是当zi为非零时.我研究了如何实现scipy中的scipy.signal.lfilter和八度中的filter.它们都在时域中直接完成,使用直接格式2和八度直接格式1 scipy(通过在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)

解决方法:

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

CS Burrus有一篇非常简洁的文章,介绍了如何以面向块的方式考虑滤波,无论是FIR还是IIR,为here.我没有详细阅读它,但是我认为它为您提供了实现IIR滤波所需的方程式卷积,包括中间状态.