FFTを使用したフィルタ（b、a、x、zi）の計算

Question

y=filter(b,a,x,zi)およびdy[i]/dx[j]を、GPUの実装で可能なスピードアップのために時間領域ではなくFFTを使用して計算したいと考えています。

特に、ziがゼロ以外の場合は、それが可能かどうかはわかりません。私はscipeyでscipy.signal.lfilter、オクターブでfilterがどのように実装されているか調べました。それらは、ダイレクトフォーム2とオクターブ直接フォーム1（DLD-FUNCTIONS/filter.ccのコードを見て）を使用してscipyで、時間領域で直接行われます。私は、MATLABのFIRフィルター（つまりa = [1.]）の場合はfftfiltに類似したFFT実装を見たことがありません。

私は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

私は少し私はFFTのを知っていたが、あなたはhttp://jc.unternet.net/src/でsedit.pyとfrequency.pyを見て、何が役立つかどうかを見ることができるものを忘れてしまいました。

Answer 2

初期条件を得るために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). 

Answer 3

あなたはP = T + 2*MN - 1でP = T + MN - 1を交換することにより、既存の実装の誤差を低減することができます。これは純粋に直感的ですが、周りの回り込みのためbfとafの分割には2*MNという用語が必要です。

C.S. Burrusは、ブロック指向の方法で、FIRまたはIIRのいずれかのフィルタリングをどのように考慮するかのかなり簡潔な書き方を持っています。here。私は詳細を読んでいませんが、中間状態を含む畳み込みによってIIRフィルタリングを実装するのに必要な式が得られたと思います。