Erstellen eines Tiefpassfilters in SciPy - Methoden und Einheiten verstehen

78

Ich versuche, ein verrauschtes Herzfrequenzsignal mit Python zu filtern. Da die Herzfrequenz niemals über 220 Schlägen pro Minute liegen sollte, möchte ich alle Geräusche über 220 Schlägen pro Minute herausfiltern. Ich wandelte 220 / Minute in 3,66666666 Hertz um und wandelte dieses Hertz dann in rad / s um, um 23,0383461 rad / s zu erhalten.

Die Abtastfrequenz des Chips, der Daten aufnimmt, beträgt 30 Hz, daher habe ich diese in rad / s konvertiert, um 188,495559 rad / s zu erhalten.

Nachdem ich einige Dinge online nachgeschlagen hatte, fand ich einige Funktionen für einen Bandpassfilter, den ich in einen Tiefpass verwandeln wollte. Hier ist der Link zum Bandpasscode , also habe ich ihn folgendermaßen konvertiert:

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)

Ich bin sehr verwirrt darüber, weil ich ziemlich sicher bin, dass die Butterfunktion die Grenz- und Abtastfrequenz in rad / s berücksichtigt, aber ich scheine eine seltsame Ausgabe zu bekommen. Ist es tatsächlich in Hz?

Zweitens, was ist der Zweck dieser beiden Zeilen:

    nyq = 0.5 * fs
    normalCutoff = cutOff / nyq

Ich weiß, dass es um Normalisierung geht, aber ich dachte, der Nyquist war das Zweifache der Abtastfrequenz, nicht die Hälfte. Und warum benutzt du den Nyquist als Normalisierer?

Kann jemand mehr darüber erklären, wie man mit diesen Funktionen Filter erstellt?

Ich habe den Filter gezeichnet mit:

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()

und habe dies bekommen, was bei 23 rad / s eindeutig nicht abschneidet:

Ergebnis

user3123955
quelle

Antworten:

155

Ein paar Kommentare:

  • Die Nyquist-Frequenz beträgt die Hälfte der Abtastrate.
  • Sie arbeiten mit regelmäßig abgetasteten Daten, möchten also einen digitalen Filter, keinen analogen Filter. Dies bedeutet, dass Sie analog=Trueim Aufruf von nicht verwenden sollten butterund Sie sollten scipy.signal.freqz(nicht freqs) verwenden, um den Frequenzgang zu generieren.
  • Ein Ziel dieser kurzen Dienstprogrammfunktionen ist es, Ihnen zu ermöglichen, alle Ihre Frequenzen in Hz auszudrücken. Sie sollten nicht in rad / sec konvertieren müssen. Solange Sie Ihre Frequenzen mit konsistenten Einheiten ausdrücken, übernimmt die Skalierung in den Dienstprogrammfunktionen die Normalisierung für Sie.

Hier ist meine modifizierte Version Ihres Skripts, gefolgt von der Handlung, die es generiert.

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()

Tiefpass Beispiel

Warren Weckesser
quelle
4
Ja, ich bin mir sicher. Betrachten Sie den Wortlaut "Sie müssen mit der doppelten Bandbreite abtasten". Das heißt, die Abtastrate muss doppelt so groß sein wie die Bandbreite des Signals.
Warren Weckesser
3
Mit anderen Worten: Sie tasten mit 30 Hz. Dies bedeutet, dass die Bandbreite Ihres Signals nicht mehr als 15 Hz betragen sollte. 15 Hz ist die Nyquist-Frequenz.
Warren Weckesser
1
Dies wiederbeleben : Ich habe in einem anderen Thread einen Vorschlag zur Verwendung von filtfilt gesehen, der eine Rückwärts- / Vorwärtsfilterung anstelle von lfilter durchführt. Wie ist Ihre Meinung dazu?
Bar
1
@ Bar: Diese Kommentare sind nicht der richtige Ort, um Ihre Frage zu beantworten. Sie könnten eine neue Stackoverflow-Frage erstellen, die Moderatoren könnten sie jedoch schließen, da es sich nicht um eine Programmierfrage handelt . Der beste Ort dafür ist wahrscheinlich dsp.stackexchange.com
Warren Weckesser
1
@samjewell Es ist eine Zufallszahl, mit der ich immer zu verwenden scheine freqz. Ich mag einen glatten Plot mit einer übermäßigen Auflösung, die ich ein wenig vergrößern kann, ohne den Plot neu generieren zu müssen, und 8000 ist groß genug, um dies in den meisten Fällen zu erreichen.
Warren Weckesser