filtre passe-bas et FFT pour les débutants avec Python

23

Je suis nouveau dans le traitement du signal et en particulier dans la FFT, donc je ne sais pas si je fais la bonne chose ici et je suis un peu confus avec le résultat.

J'ai une fonction réelle discrète (données de mesure) et je veux mettre en place un filtre passe-bas à ce sujet. L'outil de choix est Python avec le paquet numpy. Je suis cette procédure:

  • calculer le fft de ma fonction
  • couper les hautes fréquences
  • effectuer le fft inverse

Voici le code que j'utilise:

import numpy as np
sampling_length = 15.0*60.0 # measured every 15 minutes
Fs = 1.0/sampling_length
ls = range(len(data)) # data contains the function
freq = np.fft.fftfreq(len(data), d = sampling_length)
fft = np.fft.fft(data)
x = freq[:len(data)/2] 
for i in range(len(x)):
if x[i] > 0.005: # cut off all frequencies higher than 0.005
    fft[i] = 0.0
    fft[len(data)/2 + i] = 0.0
inverse = np.fft.ifft(fft)

Est-ce la bonne procédure? Le résultat inversecontient des valeurs complexes, ce qui m'embrouille.

Jusqu'à B
la source
1
Lorsque j'apprenais la FFT, j'ai trouvé cet article de blog très utile. glowingpython.blogspot.com/2011/08/…
David Poole

Réponses:

23

Il faut s'attendre à ce que le résultat soit complexe. Je veux souligner deux ou trois choses:

Vous appliquez un filtre de domaine de fréquence de mur de briques aux données, tentant de mettre à zéro toutes les sorties FFT qui correspondent à une fréquence supérieure à 0,005 Hz, puis de transformer en inverse pour obtenir à nouveau un signal de domaine temporel. Pour que le résultat soit réel, l'entrée de la FFT inverse doit être conjuguée symétrique . Cela signifie que pour une longueur - FFT,N

X[k]=X[N-k],k=1,2,,N2-1(Neven)

X[k]=X[N-k],k=1,2,,N2(No)
  • Notez que pour pair, X [ 0 ] et X [ NNX[0]ne sont pas égaux en général, mais ils sont tous les deux réels. PourNimpair,X[0]doit être réel.X[N2]NX[0]

Je vois que vous avez tenté de faire quelque chose comme ça dans votre code ci-dessus, mais ce n'est pas tout à fait correct. Si vous appliquez la condition ci-dessus sur le signal que vous passez à la FFT inverse, vous devriez obtenir un vrai signal.

sjenc(X)sjenc

sjenc

tracé de la fonction sinc

sjencsjenc

Il existe des moyens plus pratiques d'appliquer des filtres passe-bas, à la fois dans les domaines temporel et fréquentiel. Les filtres à réponse impulsionnelle finie et à réponse impulsionnelle infinie peuvent être appliqués directement à l'aide de leur représentation par équation aux différences . Ou, si votre filtre a une réponse impulsionnelle suffisamment longue, vous pouvez souvent obtenir des avantages en termes de performances en utilisant des techniques de convolution rapide basées sur la FFT (application du filtre en multipliant dans le domaine fréquentiel au lieu de convolution dans le domaine temporel), comme le chevauchement - enregistrer et superposer-ajouter des méthodes.

Jason R
la source
La fonction sinc est un filtrage idéal, cependant, non? C'est ce que tous les autres filtres visent, mais n'atteignent pas. C'est mauvais pour le traitement d'image car les images ne sont pas anti-crénelées en premier, donc cela produit une sonnerie qui a l'air horrible, mais pour les signaux audio ou autres qui ont été filtrés anti-crénelage avant l'échantillonnage, n'est-ce pas le meilleur filtre que vous pouvez obtenir?
endolith
1
Oui, mon résultat n'était pas conjugué symétrique. J'ai corrigé le code, maintenant tout fonctionne bien. Merci!
Jusqu'au
3
@endolith - un Sinc est un interpolateur idéal pour certains types d'interpolation, mais peut être loin d'être idéal comme filtre pour la plupart des types d'exigences de filtre courantes, telles que la planéité de la réponse de la bande passante, le rejet de la bande d'arrêt, etc.
hotpaw2
+1 pour la belle explication sur "pourquoi les gens
n'implémentent
Vous devez utiliser un sinc fenêtré. Si vous n'êtes pas limité dans le temps, c'est le filtre optimal, bien meilleur que Chebichev.