3

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

我不确定这是否可能，尤其是在zi非零时。我查看了如何scipy.signal.lfilter在 scipy 和filteroctave 中实现。它们都直接在时域中完成，scipy 使用直接形式 2 和八度直接形式 1（通过查看 中的代码DLD-FUNCTIONS/filter.cc）。我在任何地方都没有看到类似于fftfiltMATLAB 中的 FIR 滤波器的 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)
4

3 回答 3

2

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

CS Burrus 有一篇关于如何以面向块的方式看待过滤（无论是 FIR 还是 IIR）的非常简洁的文章，请点击此处。我没有详细阅读它，但我认为它为您提供了通过卷积实现 IIR 滤波所需的方程式，包括中间状态。

于 2010-12-12T18:06:34.393 回答
1

我忘记了我对 FFT 知之甚少，但您可以查看http://jc.unternet.net/src/上的 sedit.py 和 frequency.py，看看是否有任何帮助。

于 2010-12-12T02:15:51.860 回答
1

尝试scipy.signal.lfiltic(b, a, y, x=None)获取初始条件。

文档文本lfiltic：

Given a linear filter (b,a) and initial conditions on the output y
and the input x, return the inital conditions on the state vector zi
which is used by lfilter to generate the output given the input.

If M=len(b)-1 and N=len(a)-1.  Then, the initial conditions are given
in the vectors x and y as

x = {x[-1],x[-2],...,x[-M]}
y = {y[-1],y[-2],...,y[-N]}

If x is not given, its inital conditions are assumed zero.
If either vector is too short, then zeros are added
  to achieve the proper length.

The output vector zi contains

zi = {z_0[-1], z_1[-1], ..., z_K-1[-1]}  where K=max(M,N).
于 2010-12-12T02:23:19.097 回答