Jak zaimplementować pasmowy filtr Butterwortha za pomocą Scipy.sygnał.masło
Aktualizacja:
Znalazłem przepis Scipy na podstawie tego pytania! Więc, dla wszystkich zainteresowanych, przejdź od razu do:
Http://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
Ciężko mi osiągnąć to, co początkowo wydawało się prostym zadaniem implementacji filtra pasmowo-przepustowego Butterwortha dla tablicy numpy 1 - d (szeregów czasowych).
Parametry, które muszę uwzględnić, to sample_rate, częstotliwości odcięcia w hercach i ewentualnie kolejność (inne parametry, takie jak tłumienie, częstotliwość naturalna itp. są dla mnie bardziej niejasne, więc każda" domyślna " wartość by wystarczyła).
Mam teraz 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
Dokumenty i przykłady są mylące i niejasne, ale chciałbym zaimplementować formularz przedstawiony w pochwale oznaczony jako "dla bandpass". Znaki zapytania w komentarzach pokazują gdzie właśnie skopiowałem-wkleiłem jakiś przykład bez zrozumienie, co się dzieje.
Nie jestem inżynierem elektrykiem ani naukowcem, tylko projektantem sprzętu medycznego, który musi wykonać proste filtrowanie pasmowe na sygnałach EMG.Dzięki za pomoc!
3 answers
Możesz pominąć użycie buttord, a zamiast tego po prostu wybrać zamówienie na filtr i sprawdzić, czy spełnia on Twoje kryterium filtrowania. Aby wygenerować współczynniki filtra dla filtra pasmowego, podaj butter () kolejność filtrów, częstotliwości odcięcia Wn=[low, high] (wyrażone jako ułamek częstotliwości Nyquista, która jest połową częstotliwości próbkowania) i typ pasma btype="band".
Oto skrypt, który definiuje kilka wygodnych funkcji do pracy z filtrem pasma Butterwortha. Kiedy działa jako skrypt, składa się z dwóch wątków. Jeden pokazuje Pasmo przenoszenia w kilku rzędach filtrów dla tej samej częstotliwości próbkowania i częstotliwości odcięcia. Drugi wykres pokazuje 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, które są generowane przez ten skrypt:
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2012-09-02 06:41:24
Metoda projektowania filtrów w zaakceptowanej odpowiedzi jest poprawna, ale ma wadę. Filtry pasmowe SciPy zaprojektowane z b, a są niestabilne i mogą powodować błędne filtry w wyższe zamówienia filtrów.
Zamiast tego użyj wyjścia sos (second-order sections) 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 również 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)
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2018-03-27 00:23:41
Dla filtra pasmowego, WS jest krotką zawierającą dolną i górną częstotliwość narożników. Reprezentują one częstotliwość cyfrową, w której odpowiedź filtra jest o 3 dB mniejsza od pasma Pass.
Wp jest krotką zawierającą cyfrowe częstotliwości pasma stop. Reprezentują one miejsce, w którym zaczyna się maksymalne tłumienie.
Gpass jest maksymalną atenutacją w passband w dB, podczas gdy gstop jest attentuation w passbands.
Powiedzmy na przykład, że chciałeś zaprojektować filtr dla częstotliwości próbkowania 8000 próbek / s o częstotliwościach narożnych 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 to 1.0. Dwie częstotliwości narożne to 300/4000 i 3100/4000.
Teraz powiedzmy, że chciałeś, aby Stopery były niższe o 30 dB +/- 100 Hz od częstotliwości narożnych. Tak więc, twoje Stopery zaczynałyby się od 200 i 3200 Hz, co skutkowałoby częstotliwościami cyfrowymi 200/4000 i 3200/4000.
Aby utworzyć filtr, nazwałbyś buttord jako
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 stopu i stromości krzywej odpowiedzi, która jest określona przez różnicę między częstotliwością narożną i częstotliwością pasma stopu.
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2012-08-23 15:26:15