在 SciPy 中创建低通滤波器-理解方法和单位

我试图用巨蟒过滤一个嘈杂的心率信号。因为心率不应该超过每分钟220次，所以我想过滤掉所有超过每分钟220次的噪音。我把220/分转换成3.66666666赫兹，然后把这个赫兹转换成拉德/秒，得到23.0383461拉德/秒。

采集数据的芯片的采样频率是30Hz，因此我将其转换为 rad/s，得到188.495559 rad/s。

在网上查了一些东西之后，我找到了一些用于带通滤波器的函数，我想把它做成一个低通滤波器。这是带通码的链接，所以我把它转换成这样:

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 函数以 rad/s 表示截止值和采样频率，但我似乎得到了一个奇怪的输出。真的是赫兹吗？

其次，这两句话的目的是什么:

    nyq = 0.5 * fs
normalCutoff = cutOff / nyq

我知道这是关于正常化的东西，但我认为奈奎斯特是2倍的采样频率，而不是一半。你为什么要把 Nyquist 当成正常人？

有人能解释一下如何使用这些函数创建过滤器吗？

我用以下方法绘制了滤镜图:

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拉德/秒:

result

213720 次浏览
小开
最佳答案

A few comments:

  • The Nyquist frequency is half the sampling rate.
  • You are working with regularly sampled data, so you want a digital filter, not an analog filter. This means you should not use analog=True in the call to butter, and you should use scipy.signal.freqz (not freqs) to generate the frequency response.
  • One goal of those short utility functions is to allow you to leave all your frequencies expressed in Hz. You shouldn't have to convert to rad/sec. As long as you express your frequencies with consistent units, the fs parameter of the SciPy functions will take care of the scaling for you.

Here's my modified version of your script, followed by the plot that it generates.

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt




def butter_lowpass(cutoff, fs, order=5):
return butter(order, cutoff, fs=fs, btype='low', analog=False)


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, fs=fs, worN=8000)
plt.subplot(2, 1, 1)
plt.plot(w, 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()

lowpass example


提交答案