Laagdoorlaatfilter maken in SciPy – methoden en eenheden begrijpen

Ik probeer een luidruchtig hartslagsignaal te filteren met Python. Omdat de hartslag nooit boven de 220 slagen per minuut mag komen, wil ik alle ruis boven 220 bpm eruit filteren. Ik heb 220/minuut omgezet in 3.66666666 Hertz en die Hertz vervolgens omgezet in rad/s om 23.0383461 rad/sec te krijgen.

De bemonsteringsfrequentie van de chip die gegevens opneemt is 30 Hz, dus ik heb dat geconverteerd naar rad/s om 188.495559 rad/s te krijgen.

Na wat dingen online te hebben gezocht, vond ik een aantal functies voor een banddoorlaatfilter die ik in een laagdoorlaatfilter wilde maken. Hier is de link naar de bandpass-code, dus ik heb het omgezet naar dit:

from scipy.signal import butter, lfilter
from scipy.signal import freqs
def butter_lowpass(cutOff, fs, order=5):
    nyq = 0.5 * fs
    normalCutoff = cutOff / nyq
    b, a = butter(order, normalCutoff, btype='low', analog = True)
    return b, a
def butter_lowpass_filter(data, cutOff, fs, order=4):
    b, a = butter_lowpass(cutOff, fs, order=order)
    y = lfilter(b, a, data)
    return y
cutOff = 23.1 #cutoff frequency in rad/s
fs = 188.495559 #sampling frequency in rad/s
order = 20 #order of filter
#print sticker_data.ps1_dxdt2
y = butter_lowpass_filter(data, cutOff, fs, order)
plt.plot(y)

Ik ben echter erg in de war omdat ik er vrij zeker van ben dat de butter-functie de cutoff en bemonsteringsfrequentie in rad/s opneemt, maar ik krijg een rare output. Is het eigenlijk in Hz?

Ten tweede, wat is het doel van deze twee regels:

   nyq = 0.5 * fs
    normalCutoff = cutOff / nyq

Ik weet dat het iets met normalisatie te maken heeft, maar ik dacht dat de nyquist 2 keer de bemonsteringsfrequentie was, niet de helft. En waarom gebruik je de nyquist als normalizer?

Kan iemand meer uitleggen over het maken van filters met deze functies?

Ik heb het filter geplot met:

w, h = signal.freqs(b, a)
plt.plot(w, 20 * np.log10(abs(h)))
plt.xscale('log')
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(100, color='green') # cutoff frequency
plt.show()

en kreeg dit dat duidelijk niet afsnijdt bij 23 rad/s:

Antwoord 1, autoriteit 100%

Een paar opmerkingen:

• De Nyquist-frequentieis de helft van de bemonsteringsfrequentie.
• U werkt met regelmatig gesamplede gegevens, dus u wilt een digitaal filter, geen analoog filter. Dit betekent dat u analog=Trueniet moet gebruiken in de aanroep van butter, en dat u scipy.signal.freqzmoet gebruiken (niet freqs) om de frequentierespons te genereren.
• Eén doel van die korte hulpprogramma’s is dat je al je frequenties uitgedrukt in Hz kunt laten staan. U hoeft niet te converteren naar rad/sec. Zolang u uw frequenties uitdrukt met consistente eenheden, zorgt de schaling in de hulpprogramma’s voor de normalisatie voor u.

Hier is mijn aangepaste versie van je script, gevolgd door de plot die het genereert.

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a
def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y
# Filter requirements.
order = 6
fs = 30.0       # sample rate, Hz
cutoff = 3.667  # desired cutoff frequency of the filter, Hz
# Get the filter coefficients so we can check its frequency response.
b, a = butter_lowpass(cutoff, fs, order)
# Plot the frequency response.
w, h = freqz(b, a, worN=8000)
plt.subplot(2, 1, 1)
plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')
plt.plot(cutoff, 0.5*np.sqrt(2), 'ko')
plt.axvline(cutoff, color='k')
plt.xlim(0, 0.5*fs)
plt.title("Lowpass Filter Frequency Response")
plt.xlabel('Frequency [Hz]')
plt.grid()
# Demonstrate the use of the filter.
# First make some data to be filtered.
T = 5.0         # seconds
n = int(T * fs) # total number of samples
t = np.linspace(0, T, n, endpoint=False)
# "Noisy" data.  We want to recover the 1.2 Hz signal from this.
data = np.sin(1.2*2*np.pi*t) + 1.5*np.cos(9*2*np.pi*t) + 0.5*np.sin(12.0*2*np.pi*t)
# Filter the data, and plot both the original and filtered signals.
y = butter_lowpass_filter(data, cutoff, fs, order)
plt.subplot(2, 1, 2)
plt.plot(t, data, 'b-', label='data')
plt.plot(t, y, 'g-', linewidth=2, label='filtered data')
plt.xlabel('Time [sec]')
plt.grid()
plt.legend()
plt.subplots_adjust(hspace=0.35)
plt.show()