scipy.signal.firwin 低通滤波器的作用类似于高通滤波器

Posted 2023-03-10

【中文标题】scipy.signal.firwin 低通滤波器的作用类似于高通滤波器

【英文标题】：scipy.signal.firwin lowpass filter acts like highpass filter

【发布时间】：2021-06-12 12:41:29

【问题描述】：

在设计低/高/带通滤波器时，我遇到了它们无法按预期工作的问题（请参阅代码和输出）。 我想隔离三个频率峰值之一（具有最低频率的低通滤波器等）。但是，错误的峰被隔离了...有谁知道我的错误是什么？

    import numpy as np
    from scipy import signal
    import matplotlib.pyplot as plt

    def create_single_freq_state(N, w):
        """
        Creates a state that only contains the frequencies in w. Outputs array with N entries.
        """
        res = np.exp(2.0*np.pi*w[:,np.newaxis]*np.arange(N)/N*1j)
        return np.sum(res, axis=0)

    def band_pass_filter(f1, f2, numtaps, window="hamming"):

        fil = signal.firwin(numtaps, [f1, f2], pass_zero=False, fs=2*numtaps, window=window, scale=False)
        return fil

    def high_pass_filter(f, numtaps, window="hamming"):

        fil = signal.firwin(numtaps, f, pass_zero=False, fs=2*numtaps, window=window, scale=False)
        return fil

    def low_pass_filter(f, numtaps, window="hamming"):

        fil = signal.firwin(numtaps, f, pass_zero=True, fs=2*numtaps, window=window, scale=False)
        return fil

    def plot_freq_response(fil):

        fs = 2*len(fil)
        w, h = signal.freqz(fil, 1, worN=2000)
        plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="Hamming window")
        plt.show()

    def plot_func(f):

        x = np.arange(1, len(f)+1)  
        plt.plot(x, f, color ="red")  
        plt.show()

    #create state with several frequencies (those contained in the variable ws)
    N = 60
    ws = np.array([15, 30, 50])
    f = create_single_freq_state(N, ws)


    numtaps = 60-1
    freq1 = 19.9
    freq2 = 40.0

    #choose filter
    fil = low_pass_filter(freq1, numtaps)
    #fil = band_pass_filter(freq1, freq2, numtaps)
    #fil = high_pass_filter(freq2, numtaps)

    #plot signal in fourierspace
    plot_func(np.fft.fft(f))
    #plot frequency response
    plot_freq_response(fil)
    #plot result of the application of the filter
    plot_func(np.absolute(np.fft.fft(np.convolve(fil, f, mode="same"))))

输出：The output of my code. The wrong frequency is isolated.

【问题讨论】：

*1j 位于 create_single_freq_state 函数的分母中。这是故意的吗？另外，我建议不要在卷积中使用mode="same"，因为当信号和滤波器的长度大致相同时，边缘效应很重要。

哦，我以为这样写就意味着它在分子中...我将行更改为 res = np.exp(2.0*np.pi*1j*w[:,np. newaxis]*np.arange(N)/N)，但没有任何改变。关于卷积中的模式，如果我改变模式，结果只会变得更糟......

你是正确的1j 在分子中，对不起。无论如何，让f 成为真正唯一的信号的意义在于。然后，您可以使用np.fft.rfft，这样就无需考虑双面变换。还可以尝试使您的信号f 比您的过滤器fil 长得多，至少要长10 倍。这将消除开始和结束瞬态的影响，并使滤波器输出达到稳定状态。

【参考方案1】：

您的两个频率非常高，以至于它们显示为直流和低频（这称为混叠，是奈奎斯特采样定理的结果）。您正在处理一个复杂的输入信号，而且您的信号相对于您的滤波器来说很短，这让您很难看到这一点。

我将ws 更改为ws = np.array([7.5, 15, 25])，并将您对plot_func 的第一次调用更改为plot_func(np.abs(np.fft.fft(f)))（因为其中一个峰值完全在变换的复杂部分中）。现在输入和输出图如下所示：

【讨论】：

非常感谢您的回答！有没有一种方法可以实现峰值的隔离，而不是通过移动它们，而是通过操纵滤波器/奈奎斯特频率的大小（即我传递给 signal.firwin 的参数 fs）？样本数量也需要保持不变。这可能是一个非常愚蠢的问题，我感觉我还没有完全理解这个过程。

奈奎斯特频率将是采样率的一半，因此如果您想分析 0 到 20 kHz 的频率，例如，您必须以 > 40 kHz 的频率进行采样。您应该注意，在采样之前，可能会从奈奎斯特频率以上混叠到分析频带中的频率会被抗混叠滤波器拒绝。对于您的信号，采样设置是什么？您在测试中使用复杂信号是否有原因？

好的，我现在对所有事情的理解都好多了。关于您的评论：-)您知道如何在 python 中实现抗锯齿过滤器吗？ -)问题是我并没有真正测量。我正在自己构建信号。但是，由于问题的限制，我只能知道信号在时间 0, 1/N, 2/N, ..., (N-1)/N 的值。因此，我将永远无法过滤高于 N/2 的频率。这不是一个（大）问题。能够过滤较低频率就足够了。我希望使用抗锯齿过滤器可以做到这一点！

我没有意识到你正在构建信号。在这种情况下，您需要注意不要在构造中引入高于 Fs/2 的频率内容。在对连续数据进行采样以及更改已采样数据的采样率时，需要使用抗混叠滤波器。要在 python 中实现过滤器，您可以按照自己的方式设计过滤器，然后使用np.convolve 或scipy.signal.lfilter 应用过滤器。 lfilter 的好处是它返回过滤器状态，因此您可以继续过滤长信号。