<p>我想尝试使用fft而不是时域来计算<code>y=filter(b,a,x,zi)</code>和{<cd2>},以便在GPU实现中加速。在</p>
<p>我不确定这是否可能,尤其是当<code>zi</code>非零时。我研究了scipy中的<code>scipy.signal.lfilter</code>和倍频程中的<code>filter</code>是如何实现的。它们都是直接在时域中完成的,scipy使用直接形式2和倍频程直接形式1(通过查看<code>DLD-FUNCTIONS/filter.cc</code>中的代码)。我在MATLAB中还没有见过类似于FIR滤波器的<code>fftfilt</code>的FFT实现(即a=[1.])。在</p>
<p>我试着做<code>y = ifft(fft(b) / fft(a) * fft(x))</code>,但这似乎在概念上是错误的。另外,我不知道如何处理初始瞬态<code>zi</code>。任何引用,现有实现的指针,将不胜感激。在</p>
<p>示例代码</p>
<pre><code>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)
</code></pre>