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.
źródło
Odpowiedzi:
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 pasmabtype="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:
źródło
x[0]
? Próbowałem podobnych rzeczy z filtrem dolnoprzepustowym Cheby1 i mam ten sam problem.scipy.signal.lfilter_zi
izi
argumentu dolfilter
. Aby uzyskać szczegółowe informacje, zobacz dokumentację dotyczącąlfilter_zi
. TL; DR? Po prostu zmieńy = lfilter(b, a, data)
nazi = 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.)scipy.signal.lfiter()
wrt jest przesunięcie fazowe o 180 stopni do oryginalnego sygnału isignal.filtfilt()
wyjścia, dlaczego tak jest? Czy powinienem użyćfiltfilt()
zamiast tego, jeśli czas jest dla mnie ważny?filtfilt()
. Moja odpowiedź tutaj obejmuje przykład stosowaniafiltfilt()
w celu uniknięcia opóźnienia wywołanego przez filtr.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)
źródło
sosfilt
zsosfiltfilt
.sos
wyjście, ponieważ to zawsze pozwoli uniknąć niestabilności. Jeśli nie potrzebujesz przetwarzania w czasie rzeczywistym, zawsze powinieneś używaćsosfiltfilt
.fs
argumentu do wielu funkcji w programiescipy.signal
. Odpowiedź jest od dawna oczekiwana na aktualizację.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ą.
źródło
gpass=0.0
podnosi błąd dzielenia przez zero, więc zmieniłem go na 0,1 i błąd ustał. Poza tym, doktorzybutter
powiedzą: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ę.ws
iwp
posiada dwa elementy każdy filtr realizuje tylko niskie lub górnoprzepustowy (przezbtype
argument), ale nie pasmowo-przepustowy)