Filtre de calcul (b, a, x, zi) en utilisant FFT

Question

Je voudrais essayer de calculer y=filter(b,a,x,zi) et dy[i]/dx[j] en utilisant des FFT plutôt que dans le domaine temporel pour une accélération possible dans une implémentation GPU.

Je ne suis pas sûr que ce soit possible, en particulier lorsque zi est différent de zéro. J'ai regardé à travers comment scipy.signal.lfilter dans Scipy et filter dans l'octave sont mis en œuvre. Ils sont tous deux faits directement dans le domaine temporel, avec scipy utilisant la forme directe 2 et la forme directe d'octave 1 (en regardant le code en DLD-FUNCTIONS/filter.cc). Je n'ai vu nulle part une implémentation FFT analogue à fftfilt pour les filtres FIR dans MATLAB (c'est-à-dire a = [1.]).

J'ai essayé de faire y = ifft(fft(b)/fft(a) * fft(x)) mais cela semble être une erreur conceptuelle. En outre, je ne suis pas sûr de savoir comment gérer le transitoire initial zi. Toute référence, pointeur vers une implémentation existante, serait appréciée.

Exemple de 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) 

Answer 1

J'ai oublié peu que je connaissais TFR mais vous pouvez jeter un oeil à sedit.py et frequency.py à http://jc.unternet.net/src/ et voir si quelque chose il serait utile.

Answer 2

Essayez scipy.signal.lfiltic(b, a, y, x=None) pour obtenir les conditions initiales.

texte Doc pour 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

Vous pouvez réduire l'erreur dans votre implémentation existante en remplaçant P = T + MN - 1 avec P = T + 2*MN - 1. Ceci est purement intuitive, mais il me semble que la division de bf et af nécessitera 2*MN termes, en raison de l'enroulement. C.S. Burrus a un aperçu assez concis de la façon de considérer le filtrage, que FIR ou IIR, d'une manière orientée par bloc, here. Je ne l'ai pas lu en détail, mais je pense qu'il vous donne les équations dont vous avez besoin pour mettre en œuvre le filtrage RII par convolution, y compris les états intermédiaires.