Filtrage de l'espace de Fourier

Question

J'ai une vraie série de vecteurs x de longueur T et un filtre h de longueur t < < T. h est un filtre dans l'espace de Fourier, réel et symétrique. C'est environ 1/f.

Je voudrais filtrer x avec h pour obtenir y. Supposons que t == T et les FFT de longueur T puissent entrer en mémoire (aucune des deux ne sont vraies). Pour obtenir mes filtrés x en python, je le ferais:

import numpy as np 
from scipy.signal import fft, ifft 

y = np.real(np.ifft(np.fft(x) * h))) 

Étant donné que les conditions ne sont pas, j'ai essayé le hack:

1. Sélectionnez une taille de remplissage P < t/2, sélectionnez une taille de bloc B telle que B + 2P est une bonne taille FFT
2. Échelle h par interpolation spline pour être de taille B + 2P> t (h_scaled)
3. y = []; Loop:
   • Take bloc de longueur B + 2P de x (appelé X_B)
   • Effectuer y_b = IFFT (fft (X_B) * h_scaled)
   • goutte rembourrage P de part et d'autre de y_b et concaténer avec y
   • Sélectionnez ensuite X_B avec chevauchement dernier par P

Est-ce une bonne stratégie? Comment sélectionner le padding P dans le bon sens? Quelle est la bonne façon de faire cela? Je ne connais pas beaucoup le traitement du signal. C'est une bonne chance d'apprendre. J'utilise CUFFT pour accélérer les choses, donc ce serait génial si la plupart des opérations sont des FFT. Le problème réel est 3D. En outre, je ne suis pas préoccupé par les artefacts d'un filtre acausal.

Merci, Paul.

Answer 1

Vous êtes sur la bonne voie. La technique est appelée overlap-save processing. Est-ce que t est suffisamment court pour que les FFT de cette longueur soient en mémoire? Si oui, vous pouvez choisir votre taille de bloc B tel que B > 2*min(length(x),length(h)) et fait pour une transformation rapide. Ensuite, lorsque vous traitez, vous laissez tomber la première moitié de y_b, plutôt que de tomber des deux extrémités. Pour voir pourquoi vous déposez la première moitié, n'oubliez pas que la multiplication spectrale est la même que la circonvolution circulaire dans le domaine temporel. Convolving avec le zéro-rembourré h crée étranges transitoires glitchy dans la première moitié du résultat, mais par la seconde moitié tous les transitoires ont disparu parce que le point d'enroulement circulaire x est aligné avec la partie zéro de h. Il y a une bonne explication de ceci, avec des images, en "Theory and Application of Digital Signal Processing" by Lawrence Rabiner and Bernard Gold.

Il est important que votre filtre de domaine temporel s'amincit à 0 au moins à une extrémité de sorte que vous n'obteniez pas d'artefacts de sonnerie. Vous mentionnez que h est réel dans le domaine fréquentiel, ce qui implique qu'il a tous 0 phase. Habituellement, un tel signal ne sera continu que de manière circulaire, et lorsqu'il est utilisé comme un filtre, il crée une distorsion tout au long de la bande de fréquence. Un moyen facile de créer un filtre raisonnable est de le concevoir dans le domaine fréquentiel avec 0 phase, transformation inverse et rotation.Par exemple:

def OneOverF(N): 
    import numpy as np 
    N2 = N/2; #N has to be even! 
    x = np.hstack((1, np.arange(1, N2+1), np.arange(N2-1, 0, -1))) 
    hf = 1/(2*np.pi*x/N2) 
    ht = np.real(np.fft.ifft(hf)) # discard tiny imag part from numerical error 
    htrot = np.roll(ht, N2) 
    htwin = htrot * np.hamming(N) 
    return ht, htrot, htwin 

(Je suis assez nouveau à Python, s'il vous plaît laissez-moi savoir s'il y a une meilleure façon de coder cela).

Si vous comparez les réponses en fréquence de ht, htrot et htwin, vous voyez ce qui suit (axe x est la fréquence normalisée jusqu'à pi):

ht, en haut, a beaucoup d'ondulation . Ceci est dû à la discontinuité au bord. htrot, au milieu, c'est mieux, mais a encore ondulation. htwin est agréable et lisse, au détriment de l'aplatissement à une fréquence légèrement plus élevée. Notez que vous pouvez étendre la longueur de la section linéaire en utilisant une valeur plus grande pour N.

J'ai écrit sur le problème de discontinuité, et j'ai également écrit un exemple Matlab/Octave dans another SO question si vous voulez voir plus de détails.