数据集的傅里叶平滑处理

Question

我正在按照这个链接来平滑我的数据集。这个方法是基于一个原理，就是去掉信号的傅里叶变换中的高阶项，从而得到一个平滑的函数。这是我代码的一部分：

N = len(y)

y = y.astype(float)               # fix issue, see below
yfft = fft(y, N)

yfft[31:] = 0.0                   # set higher harmonics to zero
y_smooth = fft(yfft, N)

ax.errorbar(phase, y, yerr = err, fmt='b.', capsize=0, elinewidth=1.0)
ax.plot(phase, y_smooth/30, color='black') #arbitrary normalization, see below

不过，有些地方并没有正常工作。你可以查看生成的图：

蓝色的点是我的数据，而黑色的线应该是平滑后的曲线。

首先，我需要按照这个讨论来转换我的数据数组y。

其次，我只是随便进行了归一化，以便将曲线与数据进行比较，因为我不知道为什么原始曲线的值比数据点高得多。

最重要的是，这条曲线看起来像是“镜像”了数据点，我不知道为什么会这样。如果能得到一些建议，尤其是针对第三点的建议，那就太好了。还有更一般来说，如何针对我这个特定数据集的形状来优化平滑处理。

Answer 1

我会对使用这种方法非常谨慎。因为当你把快速傅里叶变换（FFT）中的频率成分归零时，其实是在频率领域里构建一个“砖墙”滤波器。这会导致在时间领域中与一个叫做sinc的函数进行卷积，这样可能会扭曲你想处理的信息。你可以查一下“吉布斯现象”来了解更多信息。

其实，设计一个低通滤波器或者使用简单的N点移动平均（这本身就是一个低通滤波器）来实现平滑效果，可能会更好。

Answer 2

根据我了解的情况，你想通过以下步骤来构建一个低通滤波器：

1. 先转换到频率域。（也就是做傅里叶变换）
2. 去掉不需要的频率。
3. 再转换回时间域。（也就是做逆傅里叶变换）

看你的代码，你在做第3步的时候，其实是又做了一次傅里叶变换。你应该尝试真正做一个逆傅里叶变换，这样才能回到时间域：

y_smooth = ifft(yfft, N)

可以看看 scipy signal，那里有很多现成的滤波器。

（编辑：我很想看看结果，记得分享哦！）

Answer 3

你的问题可能是因为标准的快速傅里叶变换（FFT）会导致数据的偏移。你可以在这里了解更多信息。

你的数据是实数，所以你可以利用傅里叶变换中的对称性，使用一个特别的函数 np.fft.rfft。

import numpy as np

x = np.arange(40)
y = np.log(x + 1) * np.exp(-x/8.) * x**2 + np.random.random(40) * 15
rft = np.fft.rfft(y)
rft[5:] = 0   # Note, rft.shape = 21
y_smooth = np.fft.irfft(rft)

plt.plot(x, y, label='Original')
plt.plot(x, y_smooth, label='Smoothed')
plt.legend(loc=0)
plt.show()

如果你画出 rft 的绝对值图，你会发现频率超过 5 的部分几乎没有信息，这就是我选择这个阈值的原因（当然也有一点实验的成分）。

这里是结果：