我正在尝试使用Python过滤嘈杂的心率信号。因为心率不应超过每分钟220次,所以我想过滤掉所有高于220 bpm的噪音。我将每分钟220次转换为3.66666666赫兹,然后将该赫兹转换为弧度/秒,得到23.0383461弧度/秒。
采集数据的芯片的采样频率为30Hz,因此我将其转换为弧度/秒,得到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
我知道这与归一化有关,但我以为奈奎斯特频率是采样频率的两倍,而不是一半。为什么要使用奈奎斯特频率作为归一化器呢?
有人能解释一下如何使用这些函数创建滤波器吗?
我使用以下方法绘制了滤波器:
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 rad/s处截断:
freqz
一起使用的随机数。我喜欢 平滑 的图形,具有足够的多余分辨率,以便我可以稍微放大一点而不必重新生成图形,8000 在大多数情况下都足够实现这一点。 - Warren Weckesser