Jak zaimplementować filtr pasmowy Butterwortha z Scipy.signal.butter

84

AKTUALIZACJA:

Znalazłem przepis na Scipy oparty na tym pytaniu! Zatem dla wszystkich zainteresowanych przejdź bezpośrednio do: Spis treści »Przetwarzanie sygnału» Pasmo Butterwortha


Mam trudności z osiągnięciem tego, co początkowo wydawało się prostym zadaniem, polegającym na zaimplementowaniu filtru pasmowo-przepustowego Butterwortha dla macierzy 1-D numpy (szeregi czasowe).

Parametry, które muszę uwzględnić, to sample_rate, częstotliwości odcięcia w HERTZ i ewentualnie kolejność (inne parametry, takie jak tłumienie, częstotliwość drgań własnych itp. Są dla mnie bardziej niejasne, więc wystarczyłaby każda wartość „domyślna”).

Teraz mam to, co wydaje się działać jako filtr górnoprzepustowy, ale nie jestem pewien, czy robię to dobrze:

def butter_highpass(interval, sampling_rate, cutoff, order=5):
    nyq = sampling_rate * 0.5

    stopfreq = float(cutoff)
    cornerfreq = 0.4 * stopfreq  # (?)

    ws = cornerfreq/nyq
    wp = stopfreq/nyq

    # for bandpass:
    # wp = [0.2, 0.5], ws = [0.1, 0.6]

    N, wn = scipy.signal.buttord(wp, ws, 3, 16)   # (?)

    # for hardcoded order:
    # N = order

    b, a = scipy.signal.butter(N, wn, btype='high')   # should 'high' be here for bandpass?
    sf = scipy.signal.lfilter(b, a, interval)
    return sf

wprowadź opis obrazu tutaj

Dokumentacja i przykłady są zagmatwane i niejasne, ale chciałbym zaimplementować formularz przedstawiony w poleceniu oznaczonym jako „for bandpass”. Znaki zapytania w komentarzach pokazują, gdzie właśnie skopiowałem i wkleiłem jakiś przykład bez zrozumienia, co się dzieje.

Nie jestem inżynierem elektrykiem ani naukowcem, tylko projektantem sprzętu medycznego, który musi wykonać dość proste filtrowanie pasmowoprzepustowe sygnałów EMG.

heltonbiker
źródło
Próbowałem czegoś na dsp.stackexchange, ale zbytnio koncentrują się (bardziej niż mogę znieść) na koncepcyjnych zagadnieniach inżynierii, a nie na używaniu funkcji scipy.
heltonbiker

Odpowiedzi:

119

Możesz pominąć użycie buttord, a zamiast tego po prostu wybrać kolejność dla filtra i sprawdzić, czy spełnia Twoje kryterium filtrowania. Aby wygenerować współczynniki filtru dla filtru pasmowoprzepustowego, należy podać butter () kolejność filtrów, częstotliwości odcięcia Wn=[low, high](wyrażone jako ułamek częstotliwości Nyquista, która stanowi połowę częstotliwości próbkowania) i typ pasma btype="band".

Oto skrypt, który definiuje kilka wygodnych funkcji do pracy z filtrem pasmowoprzepustowym Butterwortha. Uruchomiony jako skrypt tworzy dwa wykresy. Jeden przedstawia odpowiedź częstotliwościową przy kilku rzędach filtrów dla tej samej częstotliwości próbkowania i częstotliwości odcięcia. Drugi wykres przedstawia wpływ filtra (o kolejności = 6) na przykładowy szereg czasowy.

from scipy.signal import butter, lfilter


def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y


if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import freqz

    # Sample rate and desired cutoff frequencies (in Hz).
    fs = 5000.0
    lowcut = 500.0
    highcut = 1250.0

    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [3, 6, 9]:
        b, a = butter_bandpass(lowcut, highcut, fs, order=order)
        w, h = freqz(b, a, worN=2000)
        plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)

    plt.plot([0, 0.5 * fs], [np.sqrt(0.5), np.sqrt(0.5)],
             '--', label='sqrt(0.5)')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    # Filter a noisy signal.
    T = 0.05
    nsamples = T * fs
    t = np.linspace(0, T, nsamples, endpoint=False)
    a = 0.02
    f0 = 600.0
    x = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
    x += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
    x += a * np.cos(2 * np.pi * f0 * t + .11)
    x += 0.03 * np.cos(2 * np.pi * 2000 * t)
    plt.figure(2)
    plt.clf()
    plt.plot(t, x, label='Noisy signal')

    y = butter_bandpass_filter(x, lowcut, highcut, fs, order=6)
    plt.plot(t, y, label='Filtered signal (%g Hz)' % f0)
    plt.xlabel('time (seconds)')
    plt.hlines([-a, a], 0, T, linestyles='--')
    plt.grid(True)
    plt.axis('tight')
    plt.legend(loc='upper left')

    plt.show()

Oto wykresy generowane przez ten skrypt:

Pasmo przenoszenia dla kilku rzędów filtrów

wprowadź opis obrazu tutaj

Warren Weckesser
źródło
1
Czy wiesz, dlaczego filtrowane wyjście zawsze zaczyna się od wartości zero? Czy można dopasować go do rzeczywistej wartości wejściowej x[0]? Próbowałem podobnych rzeczy z filtrem dolnoprzepustowym Cheby1 i mam ten sam problem.
LWZ
2
@LWZ: Użyj funkcji scipy.signal.lfilter_zii ziargumentu do lfilter. Aby uzyskać szczegółowe informacje, zobacz dokumentację dotyczącą lfilter_zi. TL; DR? Po prostu zmień y = lfilter(b, a, data)na zi = lfilter_zi(b, a); y, zo = lfilter(b, a, data, zi=zi*data[0]). (Ale to może nie mieć znaczenia w przypadku filtra środkowoprzepustowego lub górnoprzepustowego.)
Warren Weckesser
1
Zauważyłem, że na wyjściu scipy.signal.lfiter()wrt jest przesunięcie fazowe o 180 stopni do oryginalnego sygnału i signal.filtfilt()wyjścia, dlaczego tak jest? Czy powinienem użyć filtfilt()zamiast tego, jeśli czas jest dla mnie ważny?
Jason
1
To jest opóźnienie fazowe filtra przy tej częstotliwości. Opóźnienie fazowe sinusoidy przez filtr Butterwortha zależy nieliniowo od częstotliwości. W przypadku opóźnienia fazy zerowej tak, możesz użyć filtfilt(). Moja odpowiedź tutaj obejmuje przykład stosowania filtfilt()w celu uniknięcia opóźnienia wywołanego przez filtr.
Warren Weckesser
1
Hej Jason, polecam zadawać pytania dotyczące teorii przetwarzania sygnałów na stronie dsp.stackexchange.com . Jeśli masz pytanie dotyczące napisanego przez Ciebie kodu, który nie działa zgodnie z oczekiwaniami, możesz rozpocząć nowe pytanie tutaj na temat stackoverflow.
Warren Weckesser
41

Metoda projektowania filtra w przyjętej odpowiedzi jest poprawna, ale ma wadę. Filtry pasmowoprzepustowe SciPy zaprojektowane z b, a są niestabilne i mogą powodować błędne filtry przy wyższych rzędach filtrów .

Zamiast tego użyj wyjścia sos (sekcje drugiego rzędu) projektu filtra.

from scipy.signal import butter, sosfilt, sosfreqz

def butter_bandpass(lowcut, highcut, fs, order=5):
        nyq = 0.5 * fs
        low = lowcut / nyq
        high = highcut / nyq
        sos = butter(order, [low, high], analog=False, btype='band', output='sos')
        return sos

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
        sos = butter_bandpass(lowcut, highcut, fs, order=order)
        y = sosfilt(sos, data)
        return y

Możesz także wykreślić pasmo przenoszenia, zmieniając

b, a = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = freqz(b, a, worN=2000)

do

sos = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = sosfreqz(sos, worN=2000)
user13107
źródło
+1, ponieważ w wielu przypadkach jest to teraz lepszy sposób. Podobnie jak w komentarzach do zaakceptowanej odpowiedzi, możliwe jest również wyeliminowanie opóźnienia fazy za pomocą filtrowania w przód-w tył. Wystarczy wymienić sosfiltz sosfiltfilt.
Mike
@Mike i user13107 Czy ten sam błąd wpływa również na górnoprzepustowe i dolnoprzepustowe filtry Butterwortha? Czy rozwiązanie jest takie samo?
dewarrn1
3
@ dewarrn1 Nie jest właściwie nazywanie tego "błędem"; algorytm jest poprawnie zaimplementowany, ale jest z natury niestabilny, więc jest to po prostu zły wybór algorytmu. Ale tak, wpływa na każdy filtr wyższego rzędu - nie tylko górnoprzepustowy lub dolnoprzepustowy i nie tylko filtry Butterwortha, ale także inne, takie jak Czebyszewa i tak dalej. W każdym razie, ogólnie rzecz biorąc, najlepiej jest zawsze wybierać soswyjście, ponieważ to zawsze pozwoli uniknąć niestabilności. Jeśli nie potrzebujesz przetwarzania w czasie rzeczywistym, zawsze powinieneś używać sosfiltfilt.
Mike
Przepraszam, że nie zauważyłem tej odpowiedzi dawno temu! @ user13107, tak, funkcja przenoszenia (lub „ba”) reprezentująca filtr liniowy ma poważne problemy numeryczne, gdy kolejność filtrów jest duża. Nawet filtry stosunkowo niskiego rzędu mogą mieć problemy, gdy pożądana szerokość pasma jest mała w porównaniu z częstotliwością próbkowania. Moja oryginalna odpowiedź została napisana przed dodaniem reprezentacji SOS do SciPy i przed dodaniem fsargumentu do wielu funkcji w programie scipy.signal. Odpowiedź jest od dawna oczekiwana na aktualizację.
Warren Weckesser
jakaś pomoc w tym? stackoverflow.com/q/60866193/5025009
seralouk
4

W przypadku filtru pasmowoprzepustowego ws jest krotką zawierającą częstotliwości dolnego i górnego rogu. Reprezentują one częstotliwość cyfrową, w której odpowiedź filtra jest o 3 dB mniejsza niż pasmo przepustowe.

wp to krotka zawierająca częstotliwości cyfrowe pasma zatrzymania. Reprezentują miejsce, w którym zaczyna się maksymalne tłumienie.

gpass to maksymalne tłumienie w pasmie przepustowym w dB, podczas gdy gstop to tłumienie w pasmach zatrzymania.

Załóżmy na przykład, że chciałeś zaprojektować filtr dla częstotliwości próbkowania 8000 próbek / s, z częstotliwościami narożnymi 300 i 3100 Hz. Częstotliwość Nyquista to częstotliwość próbkowania podzielona przez dwa, lub w tym przykładzie 4000 Hz. Równoważna częstotliwość cyfrowa wynosi 1,0. Dwie częstotliwości narożne wynoszą wówczas 300/4000 i 3100/4000.

Teraz powiedzmy, że chcesz, aby pasma zatrzymania spadły o 30 dB +/- 100 Hz od częstotliwości narożnych. W ten sposób twoje pasma zatrzymania zaczynałyby się od 200 i 3200 Hz, co daje częstotliwości cyfrowe 200/4000 i 3200/4000.

Aby utworzyć filtr, nazwałbyś buttord as

fs = 8000.0
fso2 = fs/2
N,wn = scipy.signal.buttord(ws=[300/fso2,3100/fso2], wp=[200/fs02,3200/fs02],
   gpass=0.0, gstop=30.0)

Długość powstałego filtra będzie zależna od głębokości pasm zaporowych i stromości krzywej odpowiedzi, która jest określona przez różnicę między częstotliwością narożną i częstotliwością zaporową.

sizzzzlerz
źródło
Próbowałem to wdrożyć, ale wciąż czegoś brakuje. Jedna rzecz jest taka, że gpass=0.0podnosi błąd dzielenia przez zero, więc zmieniłem go na 0,1 i błąd ustał. Poza tym, doktorzy butterpowiedzą: Passband and stopband edge frequencies, normalized from 0 to 1 (1 corresponds to pi radians / sample).mam wątpliwości, czy twoja odpowiedź dobrze wykonała obliczenia, więc nadal nad tym pracuję i wkrótce przekażę opinię.
heltonbiker
(Również, chociaż My wsi wpposiada dwa elementy każdy filtr realizuje tylko niskie lub górnoprzepustowy (przez btypeargument), ale nie pasmowo-przepustowy)
heltonbiker
1
Zgodnie z dokumentacją na docs.scipy.org/doc/scipy/reference/generated/… , buttord projektuje filtry dolnoprzepustowe , górne i pasmowe. Jeśli chodzi o gpass, przypuszczam, że buttord nie pozwala na tłumienie 0 dB w paśmie przepustowym. Następnie ustaw ją na jakąś niezerową wartość.
sizzzzlerz