Python中的带通滤波器

Question

我正在尝试在Python中使用128点汉明窗口获得0.7-4Hz的带通滤波器。我从图像中获取信号的样本(1个样本=1个图像)。fps经常会改变。

如何在Python中完成这个任务？我阅读了这篇文章：http://mpastell.com/2010/01/18/fir-with-scipy/，但我觉得firwin相当困惑。如何考虑这个可变的fps？

- Olivier_s_j

---

你的采样率是多少？如果不知道，就无法指定数字截止频率。 - sizzzzlerz

既然它等于每秒帧数，那就是这样。它不断变化...从20fps到60fps不等...因此采样频率也会改变。 - Olivier_s_j

你是说你尝试过滤的数据没有以均匀的速率采样吗？如果是这种情况，你将无法使用FIR滤波器或任何其他数字滤波器。它们需要一个均匀、一致的采样率。 - sizzzzlerz

大约是63帧每秒，所以是63赫兹。 - Olivier_s_j

Answer 1

你可以使用函数scipy.signal.firwin或scipy.signal.firwin2来创建带通FIR滤波器。你也可以使用scipy.signal.remez来设计一个FIR滤波器。

以下代码提供了一些方便的包装器，用于创建带通FIR滤波器。它使用这些包装器来创建与问题中所请求的数字相对应的带通滤波器。这假设采样是均匀的。如果采样不均匀，则不适合使用FIR滤波器。

from scipy.signal import firwin, remez, kaiser_atten, kaiser_beta

# Several flavors of bandpass FIR filters.

def bandpass_firwin(ntaps, lowcut, highcut, fs, window='hamming'):
    taps = firwin(ntaps, [lowcut, highcut], fs=fs, pass_zero=False,
                  window=window, scale=False)
    return taps

def bandpass_kaiser(ntaps, lowcut, highcut, fs, width):
    atten = kaiser_atten(ntaps, width/(0.5*fs))
    beta = kaiser_beta(atten)
    taps = firwin(ntaps, [lowcut, highcut], fs=fs, pass_zero=False,
                  window=('kaiser', beta), scale=False)
    return taps

def bandpass_remez(ntaps, lowcut, highcut, fs, width):
    delta = 0.5 * width
    edges = [0, lowcut - delta, lowcut + delta,
             highcut - delta, highcut + delta, 0.5*fs]
    taps = remez(ntaps, edges, [0, 1, 0], fs=fs)
    return taps


if __name__ == "__main__":
    import matplotlib.pyplot as plt
    from scipy.signal import freqz

    # Sample rate and desired cutoff frequencies (in Hz).
    fs = 63.0
    lowcut = 0.7
    highcut = 4.0

    ntaps = 128
    taps_hamming = bandpass_firwin(ntaps, lowcut, highcut, fs=fs)
    taps_kaiser16 = bandpass_kaiser(ntaps, lowcut, highcut, fs=fs, width=1.6)
    taps_kaiser10 = bandpass_kaiser(ntaps, lowcut, highcut, fs=fs, width=1.0)
    remez_width = 1.0
    taps_remez = bandpass_remez(ntaps, lowcut, highcut, fs=fs,
                                width=remez_width)

    # Plot the frequency responses of the filters.
    plt.figure(1, figsize=(12, 9))
    plt.clf()

    # First plot the desired ideal response as a green(ish) rectangle.
    rect = plt.Rectangle((lowcut, 0), highcut - lowcut, 1.0,
                         facecolor="#60ff60", alpha=0.2)
    plt.gca().add_patch(rect)

    # Plot the frequency response of each filter.
    w, h = freqz(taps_hamming, 1, worN=2000, fs=fs)
    plt.plot(w, abs(h), label="Hamming window")

    w, h = freqz(taps_kaiser16, 1, worN=2000, fs=fs)
    plt.plot(w, abs(h), label="Kaiser window, width=1.6")

    w, h = freqz(taps_kaiser10, 1, worN=2000, fs=fs)
    plt.plot(w, abs(h), label="Kaiser window, width=1.0")

    w, h = freqz(taps_remez, 1, worN=2000, fs=fs)
    plt.plot(w, abs(h), '--',
             label="Remez algorithm, width=%.1f" % remez_width)

    plt.xlim(0, 8.0)
    plt.ylim(0, 1.1)
    plt.grid(True)
    plt.legend(shadow=True, framealpha=1)
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.title('Frequency response of several FIR filters, %d taps' % ntaps)

    plt.show()
    # plt.savefig('plot.png')

这是脚本生成的情节。当然，最好在本地运行脚本，这样你可以放大细节。

- Warren Weckesser

Answer 2

尝试使用不一致的采样率来过滤数据非常困难（不可能？）。因此，您需要执行以下操作：

1. 创建一个具有固定采样率的新信号。固定采样率应该是最大采样率或更高。通过设置新的“网格”来表示新样本应该放置在哪里，并从现有数据中插值其值来完成此操作。根据需要准确度而存在各种插值方法。线性插值可能是不错的起点，但这取决于您正在做什么。如果不确定，请在https://dsp.stackexchange.com/上咨询。

2. 完成此操作后，您可以对信号应用标准信号处理方法，因为样本是均匀分布的，例如您链接的文章中所描述的方法。

3. 如有必要，您可能需要再次插值以恢复原始的采样位置。

如果您只想分析数据，可能会对Lomb Periodigram感兴趣。不需要进行带通滤波和分析，使用Lomb Periodigram即可，然后只查看相关频率，或按照需要加权结果。(另请参见数值计算系列。第13.8章名为“不均匀间隔数据的谱分析”，似乎比维基百科页面更友好)

- Bjorn Roche

---

FYI：SciPy有Lomb-Scargle周期图的实现，scipy.signal.lombscargle。 - Warren Weckesser

Answer 3

另一个选择是（异步）采样率转换，将数据转换为恒定的采样率（如“网格化”），然后再进行滤波。当然，这仅在您知道采样率并且只有在您确实需要过滤数据（而不仅仅是估计频谱）时才有用。

为此，例如可以使用scipy.interpolate中的InterpolatedUnivariateSpline逐帧应用，速度相当快。

- Chipmuenk