Python：仅使用 numpy 的低通滤波器

Question

我需要在 Python 中实现低通滤波器，但我唯一可以使用的模块是 numpy （不是 scipy）。我尝试np.fft.fft()在信号上使用，然后将所有高于截止频率的频率设置为 0，然后使用np.fft.ifft(). 然而这不起作用，我根本不知道如何应用过滤器。

编辑： 更改后np.abs()结果np.real()几乎是正确的。但在频谱图中，幅度小于原始参考和滤波参考中的幅度（相差 6dB）。所以看起来好像并不完全正确。有什么想法可以解决这个问题吗？

我的低通函数应采用以下参数：

signal: audio signal to be filtered
cutoff_freq: cout off frequency in Hz above which to cut off frequencies
sampling_rate: sampling rate in samples/second

应返回filterd 信号。

我当前的功能

signal: audio signal to be filtered
cutoff_freq: cout off frequency in Hz above which to cut off frequencies
sampling_rate: sampling rate in samples/second

Answer 1

我看到@Cris Luengo 的评论已经将您的解决方案发展到了正确的方向。您现在缺少的最后一件事是您从以下位置获得的频谱np.fft.fft由前半部分的正频率分量和后半部分的“镜像”负频率分量组成。

如果您现在将超出的所有分量设置bandlimit_index为零，则您将消除这些负频率分量。这解释了信号幅度下降了 6dB，您消除了一半的信号功率（+，因为您已经注意到每个真实信号都必须具有共轭对称频谱）。函数np.fft.ifft文档（ifft 文档）很好地解释了预期的格式。它指出：

“输入的顺序应与 fft 返回的顺序相同，即”

• a[0] 应包含零频率项，
• a[1:n//2] 应包含正频率项，
• a[n//2 + 1:] 应包含负频率项，从最负频率开始按升序排列。

这本质上就是您必须保持的对称性。因此，为了保留这些分量，只需将之间的分量设置bandlimit_index + 1 -> (len(fsig) - bandlimit_index)为零。

    def low_pass_filter(adata: np.ndarray, bandlimit: int = 1000, sampling_rate: int = 44100) -> np.ndarray:
        # translate bandlimit from Hz to dataindex according to sampling rate and data size
        bandlimit_index = int(bandlimit * adata.size / sampling_rate)
    
        fsig = np.fft.fft(adata)
        
        for i in range(bandlimit_index + 1, len(fsig) - bandlimit_index ):
            fsig[i] = 0
            
        adata_filtered = np.fft.ifft(fsig)
    
        return np.real(adata_filtered)

或者如果你想变得更“Pythonic”，你也可以将组件设置为零，如下所示：

fsig[bandlimit_index+1 : -bandlimit_index] = 0

这是一个可以浏览的合作实验室： https://colab.research.google.com/drive/1RR_9EYlApDMg4jAS2HuJIpSqwg5RLzGW ?usp=sharing

Answer 2

真实信号的完整 FFT 包含两个对称的一半。每一半将包含另一半的反射（纯虚数）复共轭：奈奎斯特之后的所有内容都不是真实信息。如果您将超过某个频率的所有内容设置为零，则必须尊重对称性。

这是一个以 1kHz 采样的人为信号及其 FFT：

import numpy as np
from matplotlib import pyplot as plt

t = np.linspace(0, 10, 10000, endpoint=False)
y = np.sin(2 * np.pi * 2 * t) * np.exp(-0.5 * ((t - 5) / 1.5)**2)
f = np.fft.fftfreq(t.size, 0.001)
F = np.fft.fft(y)

plt.figure(constrained_layout=True)

plt.subplot(1, 2, 1)
plt.plot(t, y)
plt.title('Time Domain')
plt.xlabel('time (s)')

plt.subplot(1, 2, 2)
plt.plot(np.fft.fftshift(f), np.fft.fftshift(F.imag), label='imag')
plt.xlim([-3, 3])
plt.title('Frequency Domain')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Imginatry Component')

频率轴如下所示：

>>> f
array([ 0. ,  0.1,  0.2, ..., -0.3, -0.2, -0.1])

请注意，除了直流分量（bin 0）之外，该轴关于中点对称，最高（奈奎斯特）频率位于中间。这就是我调用fftshift绘制绘图的原因：它重新排列数组，从最小到最大。

您很可能不需要将输入限制为整数。小数band_limit是完全可以接受的。请记住，要将频率转换为索引，您需要乘以大小和采样频率（除以 ram，而不是除以：

def low_pass_filter(data, band_limit, sampling_rate):
    cutoff_index = int(band_limit * data.size / sampling_rate)
    F = np.fft.fft(data)
    F[cutoff_index + 1 : -cutoff_index] = 0
    return np.fft.ifft(F).real

您仍然需要返回实数部分，因为 FFT 在最低位中始终会存在一些虚数舍入误差。

下面是在 2Hz 以上截止的样本信号图：

y2 = low_pass_filter(y, 2, 1000)
f2 = np.fft.fftfreq(t.size, 0.001)
F2 = np.fft.fft(y2)

plt.figure(constrained_layout=True)

plt.subplot(1, 2, 1)
plt.plot(t, y2)
plt.title('Time Domain')
plt.xlabel('time (s)')

plt.subplot(1, 2, 2)
plt.plot(np.fft.fftshift(f2), np.fft.fftshift(F2.imag), label='imag')
plt.xlim([-3, 3])
plt.title('Frequency Domain')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Imginatry Component')

还记得我们说过纯实数（或纯复数）信号的 FFT 中只有一半包含非冗余信息吗？Numpy 尊重这一点，并提供np.fft.rfft, np.fft.irfft,np.fft.rfftfreq来处理实值信号。您可以使用它来编写过滤器的更简单版本，因为不再存在对称约束。

def low_pass_filter(data, band_limit, sampling_rate):
     cutoff_index = int(band_limit * data.size / sampling_rate)
     F = np.fft.rfft(data)
     F[cutoff_index + 1:] = 0
     return np.fft.irfft(F, n=data.size).real

唯一需要注意的是，我们必须显式传入nto irfft，否则无论输入大小如何，输出大小始终是偶数。