scipy Python中的带通滤波器

7z5jn7bk  于 2023-08-05  发布在  Python
关注(0)|答案(3)|浏览(155)

我试图在python中得到一个带通滤波器,它具有128点汉明窗口,截止频率为0.7- 4 Hz。我从图像中获取信号样本。(1个样本= 1个图像)。FPS经常变化。
如何在Python中做到这一点?我读到了:http://mpastell.com/2010/01/18/fir-with-scipy/但我觉得firwin相当混乱。如何使用这个变量fps来实现这一点?

y4ekin9u

y4ekin9u1#

可以使用函数scipy.signal.firwinscipy.signal.firwin2创建带通FIR滤波器。还可以使用scipy.signal.remez设计FIR滤波器
以下代码提供了一些用于创建带通FIR滤波器的方便 Package 器。它使用这些来创建与问题中要求的数字相对应的带通滤波器。这假设采样是均匀进行的。如果采样不均匀,则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')

字符串
下面是脚本生成的情节。当然,在本地运行脚本要有用得多,这样您就可以放大细节。


的数据

9w11ddsr

9w11ddsr2#

尝试过滤采样率不一致的数据是非常困难的(不可能?)。因此,您将需要执行以下操作:
1.以固定采样率创建新信号。固定采样率应该是最大采样率或更高。通过设置一个新的“网格”来表示新样本的位置,并从现有数据中插入它们的值。根据所需的精度,存在多种插值方法。线性插值可能不是一个坏的起点,但它取决于你在做什么。如果你不确定的话,可以问https://dsp.stackexchange.com/
1.一旦你做到了这一点,你就可以对你的信号应用标准的信号处理方法,因为样本是均匀放置的,比如你链接的帖子中描述的那些。
1.如有必要,可能需要再次插值以恢复原始采样位置。
如果你只想分析你的数据,你可能会对Lomb Periodigram感兴趣。而不是带通你的数据,然后分析,你会使用Lomb Periodigram,然后只看相关的频率,或权重的结果,无论你想要的。(另见数字食谱系列。第13.8章被称为“非均匀间隔数据的谱分析”,这似乎是一个比维基百科页面更友好的介绍)

ubby3x7f

ubby3x7f3#

另一选项将是(异步)采样率转换,以在滤波之前将数据转换为恒定采样率(如“网格化”)。当然,这只有在你知道采样率的情况下才有效,而且只有在你真的需要过滤数据(而不仅仅是估计频谱)的情况下才有用。
为此,例如插值的UnivariateSpline从scipy.插值可以应用于逐帧,这是相当快的。

相关问题