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

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.

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

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

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) 

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ą.

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.

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.

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 .

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.