Comment écrire un filtre passe-bas pour le signal échantillonné en Python?

16

J'ai un signal qui échantillonne chaque 1 ns (1e-9 sec) et j'ai, disons, 1e4 points. J'ai besoin de filtrer les hautes fréquences de ce signal. Disons que je dois filtrer les fréquences supérieures à 10 MHz. Je veux que pour les fréquences inférieures à la fréquence de coupure, le signal soit transmis sans changement. Cela signifie que le gain du filtre sera de 1 pour les fréquences inférieures à la fréquence de coupure. Je voudrais pouvoir spécifier l'ordre des filtres. Je veux dire, le filtre du premier ordre a une pente de 20 dB / décade (coupure de courant) après la fréquence de coupure, le filtre du deuxième ordre a une pente de 40 dB / déc après la fréquence de coupure et ainsi de suite. La haute performance du code est importante.

Alex
la source

Réponses:

19

La réponse en fréquence du filtre conçu à l'aide de la fonction beurre est:

Réponse de Butterworth Filter

Mais il n'y a aucune raison de limiter le filtre à une conception de filtre monotone constante. Si vous souhaitez une atténuation plus élevée dans la bande d'arrêt et une bande de transition plus raide, d'autres options existent. Pour plus d'informations sur la spécification d'un filtre à l'aide d' iirdesing, consultez ceci . Comme le montrent les graphiques de réponse en fréquence pour la conception du beurre , la fréquence de coupure (point -3 dB) est loin de l'objectif. Cela peut être atténué par un sous-échantillonnage avant le filtrage (les fonctions de conception auront du mal avec un filtre aussi étroit, 2% de la bande passante). Permet de filtrer la fréquence d'échantillonnage d'origine avec la coupure spécifiée.

import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

from scipy.signal import fir_filter_design as ffd
from scipy.signal import filter_design as ifd

# setup some of the required parameters
Fs = 1e9           # sample-rate defined in the question, down-sampled

# remez (fir) design arguements
Fpass = 10e6       # passband edge
Fstop = 11.1e6     # stopband edge, transition band 100kHz
Wp = Fpass/(Fs)    # pass normalized frequency
Ws = Fstop/(Fs)    # stop normalized frequency

# iirdesign agruements
Wip = (Fpass)/(Fs/2)
Wis = (Fstop+1e6)/(Fs/2)
Rp = 1             # passband ripple
As = 42            # stopband attenuation

# Create a FIR filter, the remez function takes a list of 
# "bands" and the amplitude for each band.
taps = 4096
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000) 

# The iirdesign takes passband, stopband, passband ripple, 
# and stop attenuation.
bc, ac = ifd.iirdesign(Wip, Wis, Rp, As, ftype='ellip')  
bb, ab = ifd.iirdesign(Wip, Wis, Rp, As, ftype='cheby2') 

Filtres de taux d'échantillonnage d'origine

Comme mentionné, parce que nous essayons de filtrer un si petit pourcentage de la bande passante, le filtre n'aura pas de coupure nette. Dans ce cas, filtre passe-bas, nous pouvons réduire la bande passante pour obtenir un meilleur filtre. La fonction de rééchantillonnage python / scipy.signal peut être utilisée pour réduire la bande passante.

Notez que la fonction de rééchantillonnage effectuera un filtrage pour empêcher l'aliasing. Le préfiltrage peut également être effectué (pour réduire l'aliasing) et dans ce cas, nous pourrions simplement rééchantillonner par 100 et être fait , mais la question posée sur la création de filtres. Pour cet exemple, nous allons sous-échantillonner de 25 et créer un nouveau filtre

R = 25;            # how much to down sample by
Fsr = Fs/25.       # down-sampled sample rate
xs = signal.resample(x, len(x)/25.)

Si nous mettons à jour les paramètres de conception du filtre FIR, la nouvelle réponse est.

# Down sampled version, create new filter and plot spectrum
R = 25.             # how much to down sample by
Fsr = Fs/R          # down-sampled sample rate
Fstop = 11.1e6      # modified stopband
Wp = Fpass/(Fsr)    # pass normalized frequency
Ws = Fstop/(Fsr)    # stop normalized frequency
taps = 256
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000) 

Réponse du filtre sous-échantillonné

Le filtre opérant sur les données sous-échantillonnées a une meilleure réponse. Un autre avantage de l'utilisation d'un filtre FIR est que vous aurez une réponse de phase linéaire.

Christopher Felton
la source
1
Je vous remercie. Comment créez-vous un graphique du spectre du signal?
Alex
Merci beaucoup pour une excellente réponse! Je me demande si vous pourriez éventuellement expliquer comment appliquer un filtre FIR basé sur les coefficients calculés à l'aide de Remez? J'ai du mal à comprendre ce filtfiltque l' on veut pour le aparamètre.
ali_m
Une fois que vous avez les coefficients d'une conception de filtre ( b pour FIR b et a pour IIR), vous pouvez utiliser quelques fonctions différentes pour effectuer le filtrage: lfilter , convolve , filtfilt . Typiquement, toutes ces fonctions fonctionnent de manière similaire: y = filtfilt (b, a, x) Si vous avez un filtre FIR, définissez simplement a = 1 , x est le signal d'entrée, b est les coefficients FIR. Ce message pourrait également aider.
Christopher Felton
5

Est-ce que ça marche?

from __future__ import division
from scipy.signal import butter, lfilter

fs = 1E9 # 1 ns -> 1 GHz
cutoff = 10E6 # 10 MHz
B, A = butter(1, cutoff / (fs / 2), btype='low') # 1st order Butterworth low-pass
filtered_signal = lfilter(B, A, signal, axis=0)

Vous avez raison, cependant, la documentation n'est pas très complète. Il ressemble à butterun wrapper pour iirfilter, qui est mieux documenté :

N: int Ordre du filtre. Wn: array_like Une séquence scalaire ou de longueur 2 donnant les fréquences critiques.

La plupart de ces éléments sont clonés depuis matlab, vous pouvez donc également consulter leur documentation :

la fréquence de coupure normalisée Wn doit être un nombre compris entre 0 et 1, où 1 correspond à la fréquence de Nyquist, π radians par échantillon.

Mise à jour:

J'ai ajouté de la documentation pour ces fonctions. :) Github vous facilite la tâche.

endolith
la source
1

Je ne sais pas quelle est votre application, mais vous voudrez peut-être consulter Gnuradio: http://gnuradio.org/doc/doxygen/classgr__firdes.html

Les blocs de traitement du signal sont écrits en C ++ (bien que les graphiques de flux Gnuradio soient en Python), mais vous avez dit que les hautes performances sont importantes.

wrapperapps
la source
1

J'ai de bons résultats avec ce filtre FIR. Remarque qu'il applique le filtre deux fois, en "avant" et "en arrière", afin de compenser le décalage du signal (la filtfiltfonction ne fonctionnait pas, je ne sais pas pourquoi):

def firfilt(interval, freq, sampling_rate):
    nfreq = freq/(0.5*sampling_rate)
    taps =  sampling_rate + 1
    a = 1
    b = scipy.signal.firwin(taps, cutoff=nfreq)
    firstpass = scipy.signal.lfilter(b, a, interval)
    secondpass = scipy.signal.lfilter(b, a, firstpass[::-1])[::-1]
    return secondpass

CECI est une excellente ressource pour filtrer la conception et l'utilisation, d'où j'ai pris ce code, et d'où des exemples de filtres passe-bande et passe-haut .

heltonbiker
la source
Je ne pense pas qu'il y ait beaucoup d'avantages à filtrer en avant et en arrière un filtre FIR. Un filtre IIR peut bénéficier de l'avant / arrière (filtfilt) car vous pouvez obtenir une phase linéaire à partir d'un filtre de phase non linéaire par filtrage inverse.
Christopher Felton
2
@ChristopherFelton Je viens de renverser afin de synchroniser un signal électromyographique RAW avec la version lissée de lui-même. Je sais que je pourrais simplement décaler le signal, mais le filtrage deux fois finit par être moins difficile. Il est à noter que le deuxième passage ne change presque pas le premier passage déjà filtré ... Merci de noter!
heltonbiker
Ahh, oui. Pour supprimer le retard (délai de groupe), bon point.
Christopher Felton
1

Je n'ai pas de droit de commentaire ...

@endolith: j'utilise la même chose que vous sauf en utilisant scipy.signal.filtfilt (B, A, x)x est le vecteur d'entrée à filtrer - par exemple numpy.random.normal (taille = (N)) . filtfilt effectue une passe avant et arrière du signal. Par souci d'exhaustivité (la plupart étant identiques à @endolith):

import numpy as np
import scipy.signal as sps

input = np.random.normal(size=(N)) # Random signal as example
bz, az = sps.butter(FiltOrder, Bandwidth/(SamplingFreq/2)) # Gives you lowpass Butterworth as default
output = sps.filtfilt(bz, az, input) # Makes forward/reverse filtering (linear phase filter)

filtfilt comme l'a également suggéré @heltonbiker nécessite des tableaux de coefficients, je crois. Si vous devez effectuer un filtrage passe-bande sur une bande de base complexe, une configuration plus complexe est nécessaire, mais cela ne semble pas être un problème ici.

Lars1
la source