Я пытаюсь отфильтровать шумовой сигнал сердечного ритма с помощью Python. Поскольку частота сердечных сокращений никогда не должна превышать 220 ударов в минуту, я хочу отфильтровать весь шум выше 220 ударов в минуту. Я конвертировал 220/мин в 3,66666666 герц, а затем конвертировал этот герц в рад/с, чтобы получить 23,0383461 рад/с.
Частота дискретизации чипа, принимающего данные, составляет 30 Гц, поэтому я преобразовал его в рад/с, чтобы получить 188,495559 рад/с.
Посмотрев некоторые вещи в Интернете, я нашел некоторые функции для полосового фильтра, которые я хотел превратить в низкочастотный. Вот ссылка на код пропускной способности, поэтому я преобразовал ее в следующее:
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)
Это меня очень смущает, потому что я почти уверен, что функция "butter" принимает частоту среза и частоты дискретизации в рад/с, но, похоже, получаю странный вывод. Это на самом деле в Гц?
Во-вторых, какова цель этих двух строк:
nyq = 0.5 * fs
normalCutoff = cutOff / nyq
Я знаю кое-что о нормализации, но я подумал, что nyquist был в 2 раза больше, чем выборка, а не наполовину. И почему вы используете никвист в качестве нормализатора?
Может кто-нибудь объяснить больше о том, как создавать фильтры с этими функциями?
Я построил фильтр, используя:
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()
и получил это, которое явно не обрезается при 23 рад/с: