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

问题描述

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

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

在网上查了一些东西后，我发现了一些带通滤波器的函数，我想把它变成低通滤波器.是采样率的一半.

您正在处理定期采样的数据，因此您需要一个数字滤波器，而不是模拟滤波器.这意味着您不应该在对 butter 的调用中使用 analog=True，而应该使用 scipy.signal.freqz(而不是 freqs) 以生成频率响应.

这些简短的效用函数的一个目标是让您保留以赫兹表示的所有频率.您不必转换为 rad/sec.只要您用一致的单位表示频率，效用函数中的缩放就会为您处理归一化.

这是我修改后的脚本版本，然后是它生成的情节.

将 numpy 导入为 np从 scipy.signal 导入黄油、lfilter、freqz导入 matplotlib.pyplot 作为 pltdef butter_lowpass(cutoff, fs, order=5):nyq = 0.5 * fsnormal_cutoff = 截止/nyqb, a = butter(order, normal_cutoff, btype='low',analog=False)返回 b, adef butter_lowpass_filter(data, cutoff, fs, order=5):b, a = butter_lowpass(cutoff, fs, order=order)y = lfilter(b, a, 数据)返回 y# 过滤要求.订单 = 6fs = 30.0 # 采样率，Hzcutoff = 3.667 # 滤波器所需的截止频率，Hz# 获取滤波器系数，以便我们可以检查其频率响应.b, a = butter_lowpass(cutoff, fs, order)# 绘制频率响应.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(截止, 0.5*np.sqrt(2), 'ko')plt.axvline(cutoff, color='k')plt.xlim(0, 0.5*fs)plt.title("低通滤波器频率响应")plt.xlabel('频率 [Hz]')plt.grid()# 演示过滤器的使用.# 首先制作一些要过滤的数据.T = 5.0 # 秒n = int(T * fs) # 样本总数t = np.linspace(0, T, n, 端点=假)# 嘈杂"数据.我们想从中恢复 1.2 Hz 信号.数据 = 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)# 过滤数据，并绘制原始信号和过滤后的信号.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('时间 [秒]')plt.grid()plt.legend()plt.subplots_adjust(hspace=0.35)plt.show()

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.

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)

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

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

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

解决方案

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 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屋！