成都傅立叶（使用傅立叶变换清理时间序列数据噪声）

傅立叶变换是一种从完全不同的角度查看数据的强大方法：从时域到频域。 但是这个强大的运算用它的数学方程看起来很可怕。

将时域波变换为频域的公式如下：

下图很好地说明了傅立叶变换：将一个复杂的波分解成许多规则的正弦波。

这是完整的动画，解释了将时域波数据转换为频域视图时会发生什么。

我们可以轻松地处理频域中的数据，例如：去除噪声波。 之后，我们可以使用这个逆方程将频域数据转换回时域波：

让我们暂时忽略 FT 方程的复杂性。 假设我们已经完全理解数学方程的含义，让我们使用傅立叶变换在 Python 中做一些实际工作。

理解任何事物的最好方法就是使用它，就像学习游泳的最好方法是到进入到泳池中。

将干净的数据与噪声混合

创建两个正弦波并将它们合并为一个正弦波，然后故意用 np.random.randn(len(t)) 生成的数据污染干净的波。

import numpy as np import matplotlib.pyplot as plt plt.rcParams['figure.figsize'] = [16,10] plt.rcParams.update({'font.size':18}) #Create a simple signal with two frequencies data_step = 0.001 t = np.arange(start=0,stop=1,step=data_step) f_clean = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t) f_noise = f_clean + 2.5*np.random.randn(len(t)) plt.plot(t,f_noise,color='c',Linewidth=1.5,label='Noisy') plt.plot(t,f_clean,color='k',Linewidth=2,label='Clean') plt.legend()

（将两个信号组合成第三个信号也称为卷积或信号卷积。这是另一个有趣的话题，现在在神经网络模型中被大量使用）

带有带噪音的波浪，黑色是我们想要的波浪，绿线是噪音。

如果我隐藏图表中的颜色，我们几乎无法将噪声从干净的数据中分离出来，但是 傅立叶变换在这里可以提供帮助。我们需要做的就是将数据转换到另一个角度，从时间视图（x 轴）到频率视图（x 轴将是波频率）。

从时域到频域的转换

这里可以使用 numpy.fft 或 scipy.fft（pytorch1.8以后也增加了torch.fft这里就不详细说了）。 我发现 scipy.fft 非常方便且功能齐全，所以在本文中使用 scipy.fft，但是如果想使用其他模块或者根据公式构建自己的一个也是没问题的（代码见最后）。

from scipy.fft import rfft,rfftfreq n = len(t) yf = rfft(f_noise) xf = rfftfreq(n,data_step) plt.plot(xf,np.abs(yf))

在代码中，我使用 rfft 而不是 fft。 r 意味着reduce（我认为）只计算正频率。 所有负镜像频率将被省略。 因为他的速度更快。 rfft 函数的 yf 结果是一个复数，形式类似于 a+bj。 np.abs() 函数将为复数计算 √(a² + b²)。

这是我们原始波的神奇的频域视图。x轴表示频率。

一些在时域看起来很复杂的东西现在被转换成非常简单的频域数据。这两个峰代表两个正弦波的频率。一种波是50Hz，另一种是120Hz。再回顾一下生成正弦波的代码。

f_clean = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t)

其他频率就是噪音，并且在下一个步骤中很容易去除。

去除噪声频率

在Numpy的帮助下，我们可以很容易地将这些频率数据设置为0，除了50Hz和120Hz。

yf_abs = np.abs(yf) indices = yf_abs>300 # filter out those value under 300 yf_clean = indices * yf # noise frequency will be set to 0 plt.plot(xf,np.abs(yf_clean))

现在，所有的噪音都被清除了。

回到时域数据

最后一步就是逆转换了，将频域数据转换回时域数据

from scipy.fft import irfft new_f_clean = irfft(yf_clean) plt.plot(t,new_f_clean) plt.ylim(-6,8)

结果表明，所有的噪声波都被去除了。

这种转变是如何进行的

回到变换方程:

原始时域信号由小写 x 表示。 x[n] 表示第 n 个位置（时间）的时域数据点。

假设有10个数据点。

x = np.random.random(10)

N 应该是 10，所以，n 的范围是 0 到 9，10 个数据点。 k代表频率#，它的范围是0到9，为什么？ 极端情况是每个数据点代表一个独立的正弦波。

在传统的编程语言中，它将需要两个 for 循环，一个循环用于 k，另一个用于 n。 在 Python 中（其实使用了numpy）可以进行矢量化的操作替代循环。

Python 对复数的原生支持非常棒。 让我们构建傅立叶变换函数。

import numpy as np from scipy.fft import fft def DFT_slow(x): x = np.asarray(x, dtype=float)# ensure the data type N = x.shape[0] # get the x array length n = np.arange(N) # 1d array k = n.reshape((N, 1)) # 2d array, 10 x 1, aka column array M = np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) # [a,b] . [c,d] = ac + bd, it is a sum x = np.random.random(1024) np.allclose(DFT_slow(x), fft(x))

与来自numpy或scipy的函数相比，这个函数相对较慢，但对于理解FFT函数的工作原理来说已经足够了。

进一步的思考

傅立叶变换的思想是如此的深刻。它提醒我世界可能不是你所看到的，你的生活可能有一个完全不同的新面貌，只能通过一种变换才能看到，比如傅立叶变换。

你不仅可以转换声音数据，还可以转换图像，视频，电磁波，甚至股票交易数据(Kondratiev波)。

傅立叶变换也可以用描述运动来解释。

大圈就是我们的国家或者这个时代。 我们的个体是微小的内圈。 没有驱动一切的大圈，我们能做的很少。

工业革命发生在英国而不是其他国家不仅仅是因为蒸汽机，还有很多其他原因。 - 为什么英国首先工业化？因为当时的大圈只在英国出现了。

我们的成功可能并不全是因为自己的优点，而主要是因为你站在了风口上、你周围的人、你合作的好公司等等。 没有那些推动你前进的大圈子，小圈子再怎么转也是微乎其微的。

对傅里叶变换了解得越多，就越会觉得约瑟夫·傅里叶在 1822 年提出了这个令人难以置信的方程是有史以来最伟大的数学发现之一。

附录：四种傅里叶变换

本文中提到的所有傅里叶变换都是指离散傅里叶变换：

一般情况下我们使用电脑并尝试使用傅立叶变换做一些事情时，只会使用 DFT——本文正在讨论的变换。 如果你决定陷入数学深渊中，这里有两本书推荐阅读：

Frequency Domain and Fourier Transforms. Paul W. Cuff’s textbook from Princeton.
Digital Signal Processing — by Steven W.Smith 的第8章 DSP

引用

如果你想更深入的了解，可以看看以下的参考

1. An Interactive Guide To The Fourier Transform
2. Fourier Transforms With scipy.fft: Python Signal Processing
3. Understanding the FFT Algorithm
4. Frequency Domain and Fourier Transforms
5. Denoising Data with FFT [Python]
6. The FFT Algorithm — Simple Step by Step

作者：Andrew Zhu