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

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

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.

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:

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)

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