数据集的傅里叶平滑答案

【问题标题】：Fourier smoothing of data set数据集的傅里叶平滑
【发布时间】：2014-05-29 11:04:03
【问题描述】：

我正在关注this link 对我的数据集进行平滑处理。该技术基于去除信号傅里叶变换的高阶项的原理，从而获得平滑函数。这是我的代码的一部分：

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

但是有些东西不能正常工作。实际上，您可以检查结果图：蓝点是我的数据，黑线应该是平滑曲线。

首先，我必须通过关注this discussion 来转换我的数据数组y。

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

最重要的是，曲线对于数据点来说就像“镜面反射”，我不知道为什么会这样。如果能提供一些建议，尤其是对第三点的建议，以及更一般地说，如何针对我的特定数据集形状使用这种技术优化平滑，那就太好了。

【问题讨论】：

标签： python scipy fft smoothing

【解决方案1】：

您的问题可能是由于标准 FFT 所做的转换。你可以阅读它here.

您的数据是真实的，因此您可以利用 FT 中的对称性并使用特殊功能 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 的频率几乎没有信息，所以这就是我选择这个阈值的原因（也有点玩弄）。

结果如下：

【讨论】：

【解决方案2】：

据我所知，您希望通过执行以下操作来构建低通滤波器：

移至频域。（傅里叶变换）
删除不需要的频率。
回到时域。（逆傅立叶变换）

查看您的代码，而不是执行 3) 您只是在执行另一个傅立叶变换。相反，请尝试进行实际的傅立叶逆变换以返回时域：

y_smooth = ifft(yfft, N)

查看scipy signal 以查看一堆已经可用的过滤器。

（编辑：我很想看看结果，请分享！）

【讨论】：

在数学上，FT 是它自己的逆，取模一些取决于定义的负号，我认为这就是 OP 使用它的原因。当然，魔鬼在细节中，看起来像使用的那个，加上频移，意味着$FT[FT[f](x)](k) = f(-x)$。有关图像输出，请参阅我的答案。

【解决方案3】：

我会非常谨慎地使用这种技术。通过将 FFT 的频率分量归零，您可以有效地在频域中构建砖墙滤波器。这将导致与时域中的 sinc 卷积，并可能扭曲您要处理的信息。更多信息请查阅“吉布斯现象”。

您最好设计一个低通滤波器或使用简单的 N 点移动平均线（它本身就是一个 LPF）来完成平滑。

【讨论】：