信息处理 - 适用于 Python 初学者的低通滤波器和 FFT - 吾爱随笔录

适用于 Python 初学者的低通滤波器和 FFT

信息处理 fft 低通滤波器 Python

2022-01-09 23:20:06

我是信号处理的新手，尤其是 FFT，因此我不确定我在这里做的是否正确，我对结果有点困惑。

我有一个离散的实函数（测量数据），想在上面设置一个低通滤波器。选择的工具是带有 numpy 包的 Python。我遵循以下程序：

计算我的函数的 fft
切断高频
执行逆 fft

这是我正在使用的代码：

import numpy as np
sampling_length = 15.0*60.0 # measured every 15 minutes
Fs = 1.0/sampling_length
ls = range(len(data)) # data contains the function
freq = np.fft.fftfreq(len(data), d = sampling_length)
fft = np.fft.fft(data)
x = freq[:len(data)/2] 
for i in range(len(x)):
if x[i] > 0.005: # cut off all frequencies higher than 0.005
    fft[i] = 0.0
    fft[len(data)/2 + i] = 0.0
inverse = np.fft.ifft(fft)

这是正确的程序吗？结果inverse包含复杂的值，这让我感到困惑。

1个回答

结果是复杂的事实是可以预料的。我想指出几点：

您正在对数据应用砖墙频域滤波器，尝试将所有对应于大于 0.005 Hz 频率的 FFT 输出归零，然后进行逆变换以再次获得时域信号。为了使结果为实数，逆 FFT 的输入必须是共轭对称的。这意味着对于一个长度- $N$ 快速傅里叶变换，

X [k] = X^{*} [N - k], k = 1, 2, \dots, \frac{N}{2} - 1 (N e v e n)

$X[k] = X^*[N-k], k = 1, 2, \ldots , \frac{N}{2} - 1\;\;\;\;\;\;\;(N\;\; even)$

X [k] = X^{*} [N - k], k = 1, 2, \dots, ⌊ \frac{N}{2} ⌋ (N o d d)

$X[k] = X^*[N-k], k = 1, 2, \ldots , \lfloor \frac{N}{2} \rfloor\;\;\;\;\;\;\;\;(N\;\; odd)$

请注意，对于 $N$ 甚至， $X[0]$ 和 $X[\frac{N}{2}]$ 通常不相等，但它们都是真实的。对于奇数 $N$ , $X[0]$ 必须是真实的。

我看到你试图在上面的代码中做这样的事情，但它并不完全正确。如果您对传递给反向 FFT 的信号强制执行上述条件，那么您应该得到一个真实信号。

我的第二点更多是哲学上的：你正在做的事情会起作用，因为它会抑制你不想要的频域内容。然而，这通常不是在实践中实现低通滤波器的方式。正如我之前提到的，您所做的实际上是应用具有砖墙（即完美矩形）幅度响应的滤波器。这种滤波器的脉冲响应具有 $sinc(x)$ 形状。由于频域中的乘法等效于时域中的卷积（在使用 DFT 的情况下为循环），因此该操作等效于将时域信号与 $sinc$ 功能。

为什么这是个问题？回想一下 $sinc$ 函数在时域中看起来像（下图厚颜无耻地从维基百科借来的）：

sinc 函数图

这 $sinc$ 函数在时域上有非常广泛的支持；当您及时远离其主瓣时，它会非常缓慢地衰减。对于许多应用程序来说，这不是一个理想的属性。当您将信号与 $sinc$ ，缓慢衰减的旁瓣的影响通常在滤波输出信号的时域形式中很明显。这种效果通常被称为振铃。如果您知道自己在做什么，那么在某些情况下这种类型的过滤可能是合适的，但在一般情况下，这不是您想要的。

在时域和频域都有更实用的应用低通滤波器的方法。有限脉冲响应和无限脉冲响应滤波器可以使用它们的差分方程表示直接应用。或者，如果您的滤波器具有足够长的脉冲响应，您通常可以使用基于 FFT 的快速卷积技术（通过在频域中相乘而不是在时域中的卷积来应用滤波器）来获得性能优势，例如重叠-保存和重叠添加方法。

其它你可能感兴趣的问题

上一篇什么是 Walsh Hadamard 变换？如何在图像处理中使用它？下一篇如何从一组图像中重建 3D 结构？