3

I would like to try to compute y=filter(b,a,x,zi) and dy[i]/dx[j] using FFTs rather than in the time domain for possible speedup in a GPU implementation.

I am not sure it's possible, particularly when zi is non-zero. I looked through how scipy.signal.lfilter in scipy and filter in octave are implemented. They are both done directly in the time domain, with scipy using direct form 2 and octave direct form 1 (from looking through code in DLD-FUNCTIONS/filter.cc). I haven't seen anywhere an FFT implementation analogous to fftfilt for FIR filters in MATLAB (i.e. 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 に答える