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

Creating lowpass filter in SciPy - understanding methods and units

提问by user3123955

I am trying to filter a noisy heart rate signal with python. Because heart rates should never be above about 220 beats per minute, I want to filter out all noise above 220 bpm. I converted 220/minute into 3.66666666 Hertz and then converted that Hertz to rad/s to get 23.0383461 rad/sec.

我正在尝试用 python 过滤嘈杂的心率信号。因为心率不应超过每分钟 220 次左右，所以我想过滤掉所有高于 220 bpm 的噪音。我将 220/分钟转换为 3.66666666 赫兹，然后将该赫兹转换为弧度/秒以获得 23.0383461 弧度/秒。

The sampling frequency of the chip that takes data is 30Hz so I converted that to rad/s to get 188.495559 rad/s.

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

After looking up some stuff online I found some functions for a bandpass filter that I wanted to make into a lowpass. Here is the link the bandpass code, so I converted it to be this:

在网上查了一些东西后，我发现了一些带通滤波器的功能，我想把它做成低通。这是带通代码的链接，所以我将其转换为：

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)

I am very confused by this though because I am pretty sure the butter function takes in the cutoff and sampling frequency in rad/s but I seem to be getting a weird output. Is it actually in Hz?

我对此感到非常困惑，因为我很确定黄油函数以 rad/s 为单位接收截止频率和采样频率，但我似乎得到了一个奇怪的输出。它实际上是赫兹吗？

Secondly, what is the purpose of these two lines:

其次，这两行的目的是什么：

    nyq = 0.5 * fs
    normalCutoff = cutOff / nyq

I know it's something about normalization but I thought the nyquist was 2 times the sampling requency, not one half. And why are you using the nyquist as a normalizer?

我知道这与标准化有关，但我认为 nyquist 是采样频率的 2 倍，而不是一半。为什么你使用 nyquist 作为规范化器？

Can someone explain more about how to create filters with these functions?

有人可以解释更多有关如何使用这些功能创建过滤器的信息吗？

I plotted the filter using:

我使用以下方法绘制了过滤器：

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

and got this which clearly does not cut-off at 23 rad/s:

并得到这个显然不会在 23 rad/s 下截止的：

采纳答案by Warren Weckesser

A few comments:

几点意见：

• The Nyquist frequencyis 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=Truein 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 scaling in the utility functions takes care of the normalization for you.

• 的奈奎斯特频率是采样率的一半。
• 您正在处理定期采样的数据，因此您需要数字滤波器，而不是模拟滤波器。这意味着您不应analog=True在对 的调用中使用butter，而应使用scipy.signal.freqz(not freqs) 来生成频率响应。
• 这些简短的效用函数的一个目标是允许您保留所有频率以 Hz 表示。您不必转换为 rad/sec。只要您用一致的单位表示频率，效用函数中的缩放就会为您处理归一化。

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