pass - signal processing python

Creando filtros de paso bajo en SciPy-entendiendo métodos y unidades

Question

Estoy tratando de filtrar una señal de ritmo cardíaco ruidosa con python. Debido a que la frecuencia cardíaca nunca debe ser de aproximadamente 220 latidos por minuto, quiero filtrar todo el ruido por encima de los 220 bpm. Convertí 220 / minuto en 3.66666666 Hertz y luego convertí ese Hertz en rad / s para obtener 23.0383461 rad / seg.

La frecuencia de muestreo del chip que toma los datos es de 30 Hz, de modo que la convertí a rad / s para obtener 188.495559 rad / s.

Después de buscar algunas cosas en línea, encontré algunas funciones para un filtro de paso de banda que quería convertir en un paso bajo. Aquí está el enlace del código de paso de banda , así que lo convertí en este:

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)

Sin embargo, estoy muy confundido por esto porque estoy bastante seguro de que la función de manteca admite el corte y la frecuencia de muestreo en rad / s pero parece que obtengo una salida extraña. ¿Está realmente en Hz?

En segundo lugar, cuál es el propósito de estas dos líneas:

    nyq = 0.5 * fs
    normalCutoff = cutOff / nyq

Sé que es algo acerca de la normalización, pero pensé que el nyquist era 2 veces la frecuencia de muestreo, no la mitad. ¿Y por qué estás usando el nyquist como un normalizador?

¿Se puede explicar más sobre cómo crear filtros con estas funciones?

Traje el filtro usando

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

y obtuve esto que claramente no se corta a 23 rad / s:

Answer 1

Algunos comentarios:

• La frecuencia de Nyquist es la mitad de la frecuencia de muestreo.
• Está trabajando con datos muestreados regularmente, por lo que desea un filtro digital, no un filtro analógico. Esto significa que no debe usar analog=True en la llamada a butter , y debe usar scipy.signal.freqz (no freqs ) para generar la respuesta de frecuencia.
• Uno de los objetivos de esas funciones cortas de utilidad es permitirle dejar todas sus frecuencias expresadas en Hz. No deberías tener que convertir a rad / seg. Mientras expreses tus frecuencias con unidades consistentes, la escala en las funciones de utilidad se encarga de la normalización por ti.

Aquí está mi versión modificada de su script, seguida de la trama que genera.

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