如何用Scipy.signal.butter实现带通巴特沃斯滤波器
更新：
我在这个问题上找到了一个西西食谱！因此，对于任何感兴趣的人，请直接转到：内容信号处理巴特沃斯带
我有一个困难的时间来实现，最初似乎是一个简单的任务，实现巴特沃斯带通滤波器的一维数字序列(时间序列)。
我必须包含的参数是sample_rate、赫兹的截止频率和可能的阶次(其他参数，如衰减、自然频率等，对我来说更加模糊，因此任何“默认”值都可以)。
我现在拥有的是这个，它看起来像一个高通过滤器，但我不确定我是否做得对：
def butter_highpass(interval, sampling_rate, cutoff, order=5):
nyq = sampling_rate * 0.5
stopfreq = float(cutoff)
cornerfreq = 0.4 * stopfreq # (?)
ws = cornerfreq/nyq
wp = stopfreq/nyq
# for bandpass:
# wp = [0.2, 0.5], ws = [0.1, 0.6]
N, wn = scipy.signal.buttord(wp, ws, 3, 16) # (?)
# for hardcoded order:
# N = order
b, a = scipy.signal.butter(N, wn, btype='high') # should 'high' be here for bandpass?
sf = scipy.signal.lfilter(b, a, interval)
return sf
文档和示例既混乱又晦涩，但我想实现标记为“用于带通”的推荐中所显示的表单。注释中的问号显示了我只是在哪里复制粘贴一些示例，而不了解正在发生什么。
我不是电气工程或科学家，只是一个医疗设备设计师，需要对EMG信号进行一些相当简单的带通滤波。
回答 3
Stack Overflow用户
发布于 2012-09-02 06:41:24
您可以跳过buttord的使用，而只是为过滤器选择一个订单，看看它是否符合您的过滤标准。要生成带通滤波器的滤波器系数，给出滤波器的阶数、截止频率Wn=[lowcut, highcut]、采样率fs (以与截止频率相同的单位表示)和带型btype="band"。
下面是一个脚本，它为使用Butterworth带通过滤器定义了几个方便的函数。当作为脚本运行时，它会形成两个情节。一种是在相同采样率和截止频率下，在几个滤波器阶数下的频率响应。另一幅图演示了过滤器(带有order=6)对样本时间序列的影响。
from scipy.signal import butter, lfilter
def butter_bandpass(lowcut, highcut, fs, order=5):
return butter(order, [lowcut, highcut], fs=fs, btype='band')
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
b, a = butter_bandpass(lowcut, highcut, fs, order=order)
y = lfilter(b, a, data)
return y
if __name__ == "__main__":
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz
# Sample rate and desired cutoff frequencies (in Hz).
fs = 5000.0
lowcut = 500.0
highcut = 1250.0
# Plot the frequency response for a few different orders.
plt.figure(1)
plt.clf()
for order in [3, 6, 9]:
b, a = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = freqz(b, a, fs=fs, worN=2000)
plt.plot(w, abs(h), label="order = %d" % order)
plt.plot([0, 0.5 * fs], [np.sqrt(0.5), np.sqrt(0.5)],
'--', label='sqrt(0.5)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain')
plt.grid(True)
plt.legend(loc='best')
# Filter a noisy signal.
T = 0.05
nsamples = T * fs
t = np.arange(0, nsamples) / fs
a = 0.02
f0 = 600.0
x = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
x += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
x += a * np.cos(2 * np.pi * f0 * t + .11)
x += 0.03 * np.cos(2 * np.pi * 2000 * t)
plt.figure(2)
plt.clf()
plt.plot(t, x, label='Noisy signal')
y = butter_bandpass_filter(x, lowcut, highcut, fs, order=6)
plt.plot(t, y, label='Filtered signal (%g Hz)' % f0)
plt.xlabel('time (seconds)')
plt.hlines([-a, a], 0, T, linestyles='--')
plt.grid(True)
plt.axis('tight')
plt.legend(loc='upper left')
plt.show()下面是由这个脚本生成的情节：
Stack Overflow用户
发布于 2018-02-08 04:01:15
接受答案中的滤波器设计方法是正确的，但也存在缺陷。用b，a设计的SciPy带通滤波器是不稳定，并可能导致错误滤波器在较高过滤阶数上。
相反，使用滤波器设计的sos (二阶段)输出。
from scipy.signal import butter, sosfilt, sosfreqz
def butter_bandpass(lowcut, highcut, fs, order=5):
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
sos = butter(order, [low, high], analog=False, btype='band', output='sos')
return sos
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
sos = butter_bandpass(lowcut, highcut, fs, order=order)
y = sosfilt(sos, data)
return y此外，您还可以通过更改
b, a = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = freqz(b, a, worN=2000)至
sos = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = sosfreqz(sos, worN=2000)Stack Overflow用户
发布于 2012-08-23 15:26:15
对于带通滤波器，ws是包含上下角频率的元组。它们代表滤波器响应比通带小3 dB的数字频率。
wp是一个包含停止带数字频率的元组。它们代表最大衰减开始的位置。
gpass是dB中通带的最大衰减，而gstop是阻带中的注意点。
例如，您想要设计一个采样率为8000采样/秒的滤波器，其拐角频率分别为300和3100 Hz。奈奎斯特频率是采样率除以两个，或在本例中，4000赫兹。等效数字频率为1.0。两个拐角频率分别为300/4000和3100/4000。
现在假设你希望从拐角频率下降30 dB +/- 100赫兹。因此，你的阻带将从200和3200赫兹开始，从而产生200/4000和3200/4000的数字频率。
要创建您的过滤器，您可以将按钮调用为
fs = 8000.0
fso2 = fs/2
N,wn = scipy.signal.buttord(ws=[300/fso2,3100/fso2], wp=[200/fs02,3200/fs02],
gpass=0.0, gstop=30.0)滤波器的长度将取决于停止带的深度和响应曲线的陡度，响应曲线的陡度取决于角频率和阻带频率之间的差异。
https://stackoverflow.com/questions/12093594
复制相似问题
腾讯云开发者
