大家好，我们又见面了，我是你们的朋友全战军。
快速傅立叶变换（Fast）是信号处理和数据分析领域最重要的算法之一。 我打开一本旧的算法书，很欣赏 JW 和 John Tukey 1965 年的文章，用看似简单的计算技术解释了这个东西。
本文的目标是深入研究 Tukey FFT 算法，从根本上解释“对称性”，并通过一些直观的代码将其理论转化为实践。 希望本次研究能够让大家更全面地了解该算法的背景原理。
FFT（快速傅里叶变换）本身是离散傅里叶变换（）的一种快速算法，它将算法复杂度从原来的O(N^2)变为O(NlogN)。 离散傅里叶变换 DFT 就像一种较新的算法。 熟悉的连续傅里叶变换有以下正向和逆向定义：
快速傅里叶变换（FFT）算法【详细解释】【简单易懂】
xn到Xk的转换就是从空间域到频域的转换。 这种转换有助于研究信号的功率谱，并使某些问题的计算更加高效。 事实上，您还可以查看我们即将出版的天文学和统计学书籍的第 10 章（这里有一些图表和代码）。 作为例子，你可以查看我的文章《求解薛定谔方程》，其中展示了如何使用FFT来简化原本复杂的微分方程。
由于 FFT 在许多领域都很有用，因此提供了许多标准工具和软件包来计算它。 NumPy 和 SciPy 都有经过充分测试的封装 FFT 库，位于子模块 numpy.fft 和 scipy 中。 分别。 我所知道的最快的 FFT 是在 FFTW 包中，您也可以在 FFTW 包中使用它。
到目前为止，让我们暂时将这些库放在一边，考虑如何使用原始库从头开始计算 FFT。
计算离散傅里叶变换
为了简单起见，我们只关心正向变换，因为反向变换只能以非常相似的方式完成。 看上面的DFT表达式，只是一个直观的线性运算：向量x的矩阵乘法，
快速傅里叶变换（FFT）算法【详细解释】【简单易懂】
矩阵M可以表示为
快速傅里叶变换（FFT）算法【详细解释】【简单易懂】
这样想，我们可以使用矩阵乘法简单地计算 DFT：
1 import numpy as np
2 def DFT_slow(x):
3 """Compute the discrete Fourier Transform of the 1D array x"""
4 x = np.asarray(x, dtype=float)
5 N = x.shape[0]
6 n = np.arange(N)
7 k = n.reshape((N, 1))
8 M = np.exp(-2j * np.pi * k * n / N)
9 return np.dot(M, x)复制
对比numpy内置的FFT函数，我们仔细检查一下结果
x = np.random.random(1024)
np.allclose(DFT_slow(x), np.fft.fft(x))复制
输出：
True复制
现在要验证我们的算法有多慢，请比较两者的执行时间
%timeit DFT_slow(x)
%timeit np.fft.fft(x)复制
输出：
10 loops, best of 3: 75.4 ms per loop
10000 loops, best of 3: 25.5 µs per loop复制
正如预期的那样，使用这种简化的实现，我们的速度慢了一千倍以上。 但问题并没有那么简单。 对于长度为 N 的输入向量，FFT 为 O(N logN)，而我们的慢速算法为 O(N^2)。 这意味着 FFT 可以在 50 毫秒内完成的事情对于我们的慢速算法来说将花费近 20 个小时！ 那么FFT如何加速这个过程呢？ 答案在于他对对称性的运用。
离散傅立叶变换中的对称性
算法设计者可以使用的最重要的工具之一是利用问题的对称性。 如果你能清楚地表明问题的一部分与另一部分相关，那么你只需要计算一次子结果，从而节省计算成本。
Tukey正是使用这种方法来推导FFT。首先让我们看一下
快速傅里叶变换（FFT）算法【详细解释】【简单易懂】
价值。 根据上面的表达式，可以推导出：
快速傅里叶变换（FFT）算法【详细解释】【简单易懂】
对于所有整数 n，exp[2π in]=1。
最后一行显示了 DFT 的良好对称性：
快速傅里叶变换（FFT）算法【详细解释】【简单易懂】
简单地扩展它：
快速傅里叶变换（FFT）算法【详细解释】【简单易懂】
对于所有整数 i 也是如此。 如下所示，利用这种对称性可以更快地计算 DFT。
DFT 到 FFT：
利用对称性和Tukey证明，DFT的计算可以分为两部分。 由DFT的定义可知：
快速傅里叶变换（FFT）算法【详细解释】【简单易懂】
我们将单个 DFT 拆分为两个看起来相似的较小 DFT。 一个是奇数，一个是偶数。 到目前为止，我们还没有节省计算开销，每个部分包含(N/2)*N次计算，总共N^2。
诀窍是对每个部分使用对称性。 因为k的范围是0≤k
但我们不会止步于此。 只要我们的小傅里叶变换是偶数倍，我们就可以继续分而治之，直到分解出来的子问题太小而无法通过分而治之来提高效率。 当接近极限时，这个递归是O(n logn)级别的。
这种递归算法在这里可以很快实现，当子问题分解到合适的大小时，可以再次使用原来的“慢速方法”。
1 def FFT(x):
2 """A recursive implementation of the 1D Cooley-Tukey FFT"""
3 x = np.asarray(x, dtype=float)
4 N = x.shape[0]
5
6 if N % 2 > 0:
7 raise ValueError("size of x must be a power of 2")
8 elif N <= 32: # this cutoff should be optimized
9 return DFT_slow(x)
10 else:
11 X_even = FFT(x[::2])
12 X_odd = FFT(x[1::2])
13 factor = np.exp(-2j * np.pi * np.arange(N) / N)
14 return np.concatenate([X_even + factor[:N / 2] * X_odd,
15 X_even + factor[N / 2:] * X_odd])复制
现在让我们快速检查一下结果是否正确：
x = np.random.random(1024)
np.allclose(FFT(x), np.fft.fft(x))复制
输出：
True复制
然后与“慢速方法”的运行时间进行比较：
%timeit DFT_slow(x)
%timeit FFT(x)
%timeit np.fft.fft(x)复制
输出：
10 loops, best of 3: 77.6 ms per loop
100 loops, best of 3: 4.07 ms per loop
10000 loops, best of 3: 24.7 µs per loop复制
当前的算法比以前的算法快一个数量级。 此外，我们的递归算法是渐进的 O(n logn)。 我们实现了 FFT。
需要注意的是，我们还没有实现 numpy 的内置 FFT 算法，这是预期的。 numpy 的 fft 背后的算法是用 实现的，并且已经调整了很多年。 此外，我们的 NumPy 解决方案涉及堆栈递归和许多临时数组的分配，这显着增加了计算时间。
如果您仍然想加快速度，一个好方法是在使用/NumPy 时尽可能对重复计算进行向量化。 我们可以通过消除计算过程中的递归并使 FFT 更加高效来实现这一点。
向量化 NumPy
请注意，在上面的递归 FFT 实现中，在最低递归级别，我们执行 N/32 矩阵向量乘积。 我们的算法将受益于将这些矩阵向量乘积转换为一次性计算的矩阵矩阵乘积。 在每个递归级别，重复计算也可以向量化。 因为NumPy非常擅长这种运算，我们可以利用这一点来实现向量化FFT
1 def FFT_vectorized(x):
2 """A vectorized, non-recursive version of the Cooley-Tukey FFT"""
3 x = np.asarray(x, dtype=float)
4 N = x.shape[0]
5
6 if np.log2(N) % 1 > 0:
7 raise ValueError("size of x must be a power of 2")
8
9 # N_min here is equivalent to the stopping condition above,
10 # and should be a power of 2
11 N_min = min(N, 32)
12
13 # Perform an O[N^2] DFT on all length-N_min sub-problems at once
14 n = np.arange(N_min)
15 k = n[:, None]
16 M = np.exp(-2j * np.pi * n * k / N_min)
17 X = np.dot(M, x.reshape((N_min, -1)))
18
19 # build-up each level of the recursive calculation all at once
20 while X.shape[0] < N:
21 X_even = X[:, :X.shape[1] / 2]
22 X_odd = X[:, X.shape[1] / 2:]
23 factor = np.exp(-1j * np.pi * np.arange(X.shape[0])
24 / X.shape[0])[:, None]
25 X = np.vstack([X_even + factor * X_odd,
26 X_even - factor * X_odd])
27
28 return X.ravel()复制
x = np.random.random(1024)
np.allclose(FFT_vectorized(x), np.fft.fft(x))复制
输出：
True复制
因为我们算法的效率有了很大的提升，所以我们来做一个更大的测试（不包括在内）
x = np.random.random(1024 * 16)
%timeit FFT(x)
%timeit FFT_vectorized(x)
%timeit np.fft.fft(x)复制
输出：
10 loops, best of 3: 72.8 ms per loop
100 loops, best of 3: 4.11 ms per loop
1000 loops, best of 3: 505 µs per loop复制
我们的实施已进入新的水平。 在这里，我们仅使用几十行 + NumPy 代码，对大约 10 个因子进行基准测试。 尽管没有相应的计算来证明，但该版本远远优于源代码，您可以在此处浏览。
那么如何获得最后一点加速度呢？ 也许这只是一本详细的日志，花费了大量的时间来确保任何子计算都可以被重用。 我们这里的numpy版本涉及额外的内存分配和复制，可以通过一些低级语言（例如）轻松控制和最小化。 此外，Tukey算法也可以分为两个以上的部分（正如我们在这里使用的-Tukey FFT base 2算法），并且还可以应用其他更先进的FFT算法，包括基于基本卷积的算法。 关于不同的方法（例如算法和 Rader 算法）。 结合上述思想和方法的扩展，即使数组大小不满足2的幂，也可以快速执行FFT。
