在SciPy中创建低通滤波器-了解方法和单位 [英] Creating lowpass filter in SciPy - understanding methods and units

问题描述

我正在尝试使用python过滤嘈杂的心率信号.因为心率永远不应该超过每分钟220次跳动，所以我想过滤掉所有220 bpm以上的噪音.我将220/min转换为3.66666666赫兹，然后将其转换为rad/s，以获得23.0383461 rad/sec.

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.

要采集数据的芯片的采样频率为30Hz，所以我将其转换为rad/s以获得188.495559 rad/s.

The sampling frequency of the chip that takes data is 30Hz so I converted that to rad/s to get 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)

我对此感到非常困惑，因为我非常确定黄油功能会以rad/s为单位进行截止和采样频率，但是我似乎得到了一个奇怪的输出.真的是赫兹吗?

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?

第二，这两行的目的是什么

Secondly, what is the purpose of these two lines:

    nyq = 0.5 * fs
    normalCutoff = cutOff / nyq

我知道这与标准化有关，但我认为奈奎斯特是采样频率的2倍，而不是一半.为什么要使用奈奎斯特作为归一化器?

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?

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

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

得到的结果显然不会以23 rad/s的速度截止:

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

推荐答案

一些评论:

• 奈奎斯特频率是采样率的一半.
• 您正在使用定期采样的数据，因此需要数字滤波器而不是模拟滤波器.这意味着您不应在对butter的调用中使用analog=True，而应使用scipy.signal.freqz(而不是freqs)来生成频率响应.
• 这些简短实用程序功能的一个目标是让您保留所有以Hz表示的频率.您不必转换为弧度/秒.只要您以一致的单位表示频率，实用程序功能中的缩放就可以为您进行归一化.

• 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 scaling in the utility functions takes care of the normalization 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):
    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()

这篇关于在SciPy中创建低通滤波器-了解方法和单位的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！