Python：仅使用numpy实现低通滤波器

Question

我需要在Python中实现一个低通滤波器，但我只能使用numpy模块（不是scipy）。我尝试在信号上使用np.fft.fft()，然后将所有高于截止频率的频率设置为0，然后使用np.fft.ifft()。然而这并没有起作用，我也不确定如何应用滤波器。

编辑后，更改np.abs()为np.real()，结果几乎正确了。但在频谱图中，振幅比原始和过滤参考值要小（差6dB）。所以看起来它还没有完全正确。有什么方法可以解决这个问题吗？

我的Lowpass函数应该接受以下参数：

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

过滤后的信号应该被返回。

我的当前函数。

    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)):
            fsig[i] = 0
            
        adata_filtered = np.fft.ifft(fsig)
    
        return np.real(adata_filtered)

- TheBaum

---

请提供一个最小可复现示例。 - a_guest

“this didn't work” - 具体是什么问题？ - DYZ

如果逆傅里叶变换返回复数（其中虚部与0显著不同），则表示您的滤波器存在问题。在您的情况下，它没有在频率域中保持共轭对称性。如果正确执行，您的逆傅里叶变换将返回一个复杂数组，但虚部应该近似为零。取其实部而非幅值即可。 - Cris Luengo

你不需要上传文件。你应该绘制步骤并首先向自己解释流程在哪里出了问题。 - Mad Physicist

@CrisLuengo。当然，如果您不对零进行对称设置并尊重额外的零，则会出现这种情况。 - Mad Physicist

Answer 1

我看到@Cris Luengo的评论已经把你的解决方案引向了正确的方向。现在你所缺少的仅是使用np.fft.fft获得的频谱由第一半部分的正频率成分和第二半部分的“镜像”负频率成分组成。

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

“输入应按与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)

或者如果你想更符合Python风格，也可以像这样将组件设置为零：

fsig[bandlimit_index+1 : -bandlimit_index] = 0

这是一个可供参考的Colab： https://colab.research.google.com/drive/1RR_9EYlApDMg4jAS2HuJIpSqwg5RLzGW?usp=sharing

- lwohlhart

---

非常感谢，您的解释真的很有帮助。现在我的结果也符合参考文献，并且Colab演示真的有助于我理解 :) - TheBaum

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

唯一的注意事项是，我们必须明确传递 n 给 irfft，否则输出大小将始终为偶数，无论输入大小如何。

- Mad Physicist

---

很有帮助！我之前不知道还有这个方便的 np.fft.rfft :-) - lwohlhart