如何实现带通巴特沃斯滤波器

更新:

我发现了一个基于这个问题的 Scipy 食谱! 所以，对任何感兴趣的人，直接去: 内容“信号处理”巴特沃斯带通

我很难完成最初看起来很简单的任务，那就是实现一维数字数组(时间序列)的 Butterworth 带通滤波器。

我必须包含的参数是 sample _ rate、 HERTZ 的截止频率和可能的顺序(其他参数，比如衰减、自然频率等对我来说更加模糊，所以任何“默认”值都可以)。

我现在拥有的是这个，它看起来像是一个高通滤波器，但我不确定我做得对不对:

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

enter image description here

文档和示例令人困惑且晦涩难懂，但我想实现标记为“ for bandpass”的推荐中提供的表单。注释中的问号表示我只是复制粘贴了一些示例，而没有理解正在发生什么。

我不是电气工程师或科学家，只是一个医疗设备设计师需要执行一些相当简单的带通滤波肌电信号。

202756 次浏览
小开

For a bandpass filter, ws is a tuple containing the lower and upper corner frequencies. These represent the digital frequency where the filter response is 3 dB less than the passband.

wp is a tuple containing the stop band digital frequencies. They represent the location where the maximum attenuation begins.

gpass is the maximum attenutation in the passband in dB while gstop is the attentuation in the stopbands.

Say, for example, you wanted to design a filter for a sampling rate of 8000 samples/sec having corner frequencies of 300 and 3100 Hz. The Nyquist frequency is the sample rate divided by two, or in this example, 4000 Hz. The equivalent digital frequency is 1.0. The two corner frequencies are then 300/4000 and 3100/4000.

Now lets say you wanted the stopbands to be down 30 dB +/- 100 Hz from the corner frequencies. Thus, your stopbands would start at 200 and 3200 Hz resulting in the digital frequencies of 200/4000 and 3200/4000.

To create your filter, you'd call buttord as

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)

The length of the resulting filter will be dependent upon the depth of the stop bands and the steepness of the response curve which is determined by the difference between the corner frequency and stopband frequency.

小开
最佳答案

You could skip the use of buttord, and instead just pick an order for the filter and see if it meets your filtering criterion. To generate the filter coefficients for a bandpass filter, give butter() the filter order, the cutoff frequencies Wn=[lowcut, highcut], the sampling rate fs (expressed in the same units as the cutoff frequencies) and the band type btype="band".

Here's a script that defines a couple convenience functions for working with a Butterworth bandpass filter. When run as a script, it makes two plots. One shows the frequency response at several filter orders for the same sampling rate and cutoff frequencies. The other plot demonstrates the effect of the filter (with order=6) on a sample time series.

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

Here are the plots that are generated by this script:

Frequency response for several filter orders

enter image description here

小开

The filter design method in accepted answer is correct, but it has a flaw. SciPy bandpass filters designed with b, a are unstable and may result in erroneous filters at higher filter orders.

Instead, use sos (second-order sections) output of filter design.

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

Also, you can plot frequency response by changing 

b, a = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = freqz(b, a, worN=2000)

to

sos = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = sosfreqz(sos, worN=2000)

提交答案