Jak napisać filtr dolnoprzepustowy dla próbkowanego sygnału w Pythonie?

16

Mam pewien sygnał, który próbkował co 1 ns (1e-9 sekund) i mam, powiedzmy, 1e4 punkty. Muszę odfiltrować wysokie częstotliwości z tego sygnału. Powiedzmy, że muszę filtrować częstotliwości wyższe niż 10 MHz. Chcę, aby dla częstotliwości niższych niż odcięty sygnał częstotliwości był przekazywany bez zmian. Oznacza to, że wzmocnienie filtra wyniesie 1 dla częstotliwości niższych niż częstotliwość graniczna. Chciałbym móc określić kolejność filtrów. Mam na myśli, że filtr pierwszego rzędu ma nachylenie 20 db / dekadę (wyłączanie zasilania) po częstotliwości odcięcia, filtr drugiego rzędu ma nachylenie 40 db / dec po częstotliwości odcięcia i tak dalej. Ważna jest wysoka wydajność kodu.

Alex
źródło

Odpowiedzi:

19

Pasmo przenoszenia dla filtra zaprojektowanego przy użyciu funkcji masła wynosi:

Odpowiedź filtra Butterwortha

Ale nie ma powodu, aby ograniczać filtr do stałej konstrukcji filtra monotonicznego. Jeśli potrzebujesz wyższego tłumienia w paśmie zatrzymania i bardziej stromym paśmie przejścia, istnieją inne opcje. Aby uzyskać więcej informacji na temat określania filtra za pomocą iirdesing, zobacz to . Jak pokazują wykresy odpowiedzi częstotliwości dla projektu masła, częstotliwość graniczna (punkt -3 dB) jest daleka od celu. Można to złagodzić poprzez próbkowanie w dół przed filtrowaniem (funkcje projektowe będą miały trudny czas przy tak wąskim filtrze, 2% szerokości pasma). Przyjrzyjmy się filtrowaniu oryginalnej częstotliwości próbkowania z określonym punktem odcięcia.

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') 

Oryginalne filtry częstotliwości próbkowania

Jak wspomniano, ponieważ próbujemy przefiltrować tak mały procent szerokości pasma, filtr nie będzie miał ostrego odcięcia. W tym przypadku filtr dolnoprzepustowy możemy zmniejszyć przepustowość, aby uzyskać lepiej wyglądający filtr. Funkcja ponownego próbkowania python / scipy.signal może być wykorzystana do zmniejszenia przepustowości.

Uwaga: funkcja ponownego próbkowania wykona filtrowanie, aby zapobiec aliasingowi. Filtrowanie wstępne może być również wykonane (w celu zmniejszenia aliasingu), w tym przypadku możemy po prostu przeskalować ponownie o 100 i wykonać , ale zadaje się pytanie o tworzenie filtrów. W tym przykładzie zmniejszymy próbkę o 25 i utworzymy nowy filtr

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

Jeśli zaktualizujemy parametry projektu dla filtra FIR, nowa odpowiedź to.

# 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) 

Odpowiedź filtra z obniżonym próbkowaniem

Filtr działający na danych o zmniejszonej próbce ma lepszą odpowiedź. Kolejną zaletą stosowania filtra FIR jest to, że będziesz miał liniową odpowiedź fazową.

Christopher Felton
źródło
1
Dziękuję Ci. Jak tworzysz wykres widma sygnału?
Alex
Dziękuję bardzo za doskonałą odpowiedź! Zastanawiam się, czy mógłbyś wyjaśnić, jak zastosować filtr FIR na podstawie współczynników obliczonych za pomocą Remeza? Mam problem ze zrozumieniem, czego filtfiltchce aparametr.
ali_m
Po uzyskaniu współczynników z projektu filtra ( b dla FIR b i a dla IIR) możesz użyć kilku różnych funkcji do wykonania filtrowania: lfilter , zwoje , filtry . Zazwyczaj wszystkie te funkcje działają podobnie: y = filtfilt (b, a, x) Jeśli masz filtr FIR po prostu ustaw a = 1 , x to sygnał wejściowy, b to współczynniki FIR. Ten post może również pomóc.
Christopher Felton
5

czy to działa?

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)

Masz rację, dokumentacja nie jest zbyt kompletna. Wygląda butterto na opakowanie iirfilter, które jest lepiej udokumentowane :

N: int Kolejność filtra. Wn: array_like Sekwencja skalarna lub długość-2, dająca częstotliwości krytyczne.

Większość tych rzeczy jest jednak sklonowanych z Matlaba, więc możesz też przejrzeć ich dokumentację :

znormalizowana częstotliwość odcięcia Wn musi być liczbą od 0 do 1, gdzie 1 odpowiada częstotliwości Nyquista, π radianów na próbkę.

Aktualizacja:

Dodałem dokumentację dla tych funkcji. :) Github ułatwia.

endolit
źródło
1

Nie jestem pewien, jaka jest twoja aplikacja, ale możesz sprawdzić Gnuradio: http://gnuradio.org/doc/doxygen/classgr__firdes.html

Bloki przetwarzania sygnałów są napisane w C ++ (chociaż wykresy przepływu Gnuradio są w języku Python), ale powiedziałeś, że wysoka wydajność jest ważna.

wrapperapps
źródło
1

Mam dobre wyniki z tym filtrem FIR. Zauważa, że ​​stosuje filtr dwa razy, przesuwając „do przodu” i „do tyłu”, aby zrekompensować przesunięcie sygnału ( filtfiltfunkcja nie działała, nie wiem dlaczego):

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

Świetnym źródłem do filtrowania projektu i użycia, skąd wziąłem ten kod i skąd można wziąć przykłady filtrów pasmowoprzepustowych i hi-passowych, jest TO .

heltonbiker
źródło
Nie sądzę, że jest wiele korzyści z filtrowania do przodu i do tyłu filtru FIR. Filtr IIR może korzystać z przewijania do przodu / do tyłu (filtfilt), ponieważ można uzyskać fazę liniową z nieliniowego filtra fazy przez filtrowanie do tyłu.
Christopher Felton,
2
@ChristopherFelton Właśnie odwracam, aby zsynchronizować sygnał elektromiograficzny RAW z jego wygładzoną wersją. Wiem, że mogę po prostu przesunąć sygnał, ale dwukrotne filtrowanie kończy się mniejszym problemem. Warto zauważyć, że drugie przejście prawie nie zmienia już przefiltrowanego pierwszego przejścia ... Dziękujemy za odnotowanie!
heltonbiker
Ach tak. Aby usunąć opóźnienie (opóźnienie grupowe), warto zwrócić uwagę.
Christopher Felton,
1

Nie mam uprawnień do komentowania ...

@endolith: Używam tego samego, co ty, z wyjątkiem scipy.signal.filtfilt (B, A, x), gdzie x jest wektorem wejściowym do filtrowania - np. numpy.random.normal (size = (N)) . filtfilt wykonuje przekazywanie sygnału do przodu i do tyłu. Ze względu na kompletność (większość jest taka sama jak @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, jak sugeruje również @heltonbiker, wymaga tablic współczynników, które moim zdaniem. W przypadku, gdy konieczne jest przeprowadzenie filtrowania pasmowoprzepustowego w złożonym paśmie podstawowym, wymagana jest bardziej zaangażowana konfiguracja, ale tutaj nie wydaje się to stanowić problemu.

Lars1
źródło