- 编注：这篇译文由@unblock 和@jingliang 共同完成。
- 再次推荐：《如果看了此文你还不懂傅里叶变换，那就过来掐死我吧【完整版】》
快速傅里叶变换（Fast Fourier Transform）是信号处理与数据分析领域里最重要的算法之一。没有正规计算机科学课程背景的我，使用这个算法多年，但这周我却突然想起自己从没思考过为什么FFT能如此快速地计算离散傅里叶变换。我打开一本老旧的算法书，欣赏了JW Cooley 和 John Tukey 在1965年的文章中，以看似简单的计算技巧来讲解这个东西。
本文的目标是，深入Cooley-Tukey FFT 算法，解释作为其根源的“对称性”，并以一些直观的python代码将其理论转变为实际。我希望这次研究能使数据科学家（例如我），对这个算法的背景原理有更全面的认识。
FFT（快速傅里叶变换）本身就是离散傅里叶变换（Discrete Fourier Transform）的快速算法，使算法复杂度由原本的O(N^2) 变为 O(NlogN)，离散傅里叶变换DFT，如同更为人熟悉的连续傅里叶变换，有如下的正、逆定义形式：
xn 到 Xk 的转化就是空域到频域的转换，这个转换有助于研究信号的功率谱，和使某些问题的计算更有效率。事实上，你还可以查看一下我们即将推出的天文学和统计学的图书的第十章（这里有一些图示和python代码）。作为一个例子，你可以查看下我的文章《用python求解薛定谔方程》，是如何利用FFT将原本复杂的微分方程简化。
正因为FFT在那么多领域里如此有用，python提供了很多标准工具和封装来计算它。NumPy 和 SciPy 都有经过充分测试的封装好的FFT库，分别位于子模块 numpy.fft 和 scipy.fftpack 。我所知的最快的FFT是在 FFTW包中 ，而你也可以在python的pyFFTW 包中使用它。
虽然说了这么远，但还是暂时先将这些库放一边，考虑一下怎样使用原始的python从头开始计算FFT。
计算离散傅里叶变换
简单起见，我们只关心正变换，因为逆变换也只是以很相似的方式就能做到。看一下上面的DFT表达式，它只是一个直观的线性运算：向量x的矩阵乘法，
矩阵M可以表示为
这么想的话，我们可以简单地利用矩阵乘法计算DFT：
|
1 2 3 4 5 6 7 8 9 |
import numpy as np def DFT_slow(x): """Compute the discrete Fourier Transform of the 1D array x""" x = np.asarray(x, dtype=float) N = x.shape[0] n = np.arange(N) k = n.reshape((N, 1)) M = np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) |
对比numpy的内置FFT函数，我们来对结果进行仔细检查，
|
1 2 |
x = np.random.random(1024) np.allclose(DFT_slow(x), np.fft.fft(x)) |
输出：
|
1 |
True |
现在为了验证我们的算法有多慢，对比下两者的执行时间
|
1 2 3 |
%timeit DFT_slow(x) %timeit np.fft.fft(x) |
输出：
|
1 2 |
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是怎么提速完事的呢？答案就在于他利用了对称性。
离散傅里叶变换中的对称性
算法设计者所掌握的最重要手段之一，就是利用问题的对称性。如果你能清晰地展示问题的某一部分与另一部分相关，那么你就只需计算子结果一次，从而节省了计算成本。
Cooley 和 Tukey 正是使用这种方法导出FFT。 首先我们来看下
的值。根据上面的表达式，推出：
对于所有的整数n，exp[2π i n]=1。
最后一行展示了DFT很好的对称性：
简单地拓展一下：
同理对于所有整数 i 。正如下面即将看到的，这个对称性能被利用于更快地计算DFT。
DFT 到 FFT：
利用对称性 Cooley 和 Tukey 证明了，DFT的计算可以分为两部分。从DFT的定义得：
我们将单个DFT分成了看起来相似两个更小的DFT。一个奇，一个偶。目前为止，我们还没有节省计算开销，每一部分都包含(N/2)∗N的计算量，总的来说，就是N^2 。
技巧就是对每一部分利用对称性。因为 k 的范围是 0≤k<N ， 而 n 的范围是 0≤n<M≡N/2 ， 从上面的对称性特点来看，我们只需对每个子问题作一半的计算量。O(N^2) 变成了 O(M^2) 。
但我们不是到这步就停下来，只要我们的小傅里叶变换是偶倍数，就可以再作分治，直到分解出来的子问题小到无法通过分治提高效率，接近极限时，这个递归是 O(n logn) 级的。
这个递归算法能在python里快速实现，当子问题被分解到合适大小时，再用回原本那种“慢方法”。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
def FFT(x): """A recursive implementation of the 1D Cooley-Tukey FFT""" x = np.asarray(x, dtype=float) N = x.shape[0] if N % 2 > 0: raise ValueError("size of x must be a power of 2") elif N <= 32: # this cutoff should be optimized return DFT_slow(x) else:
快速傅里叶变换（Fast Fourier Transform）是信号处理与数据分析领域里最重要的算法之一。没有正规计算机科学课程背景的我，使用这个算法多年，但这周我却突然想起自己从没思考过为什么FFT能如此快速地计算离散傅里叶变换。我打开一本老旧的算法书，欣赏了JW Cooley 和 John Tukey 在1965年的文章中，以看似简单的计算技巧来讲解这个东西。 本文的目标是，深入Cooley-Tukey FFT 算法，解释作为其根源的“对称性”，并以一些直观的python代码将其理论转变为实际。我希望这次研究能使数据科学家（例如我），对这个算法的背景原理有更全面的认识。 FFT（快速傅里叶变换）本身就是离散傅里叶变换（Discrete Fourier Transform）的快速算法，使算法复杂度由原本的O(N^2) 变为 O(NlogN)，离散傅里叶变换DFT，如同更为人熟悉的连续傅里叶变换，有如下的正、逆定义形式：
xn 到 Xk 的转化就是空域到频域的转换，这个转换有助于研究信号的功率谱，和使某些问题的计算更有效率。事实上，你还可以查看一下我们即将推出的天文学和统计学的图书的第十章（这里有一些图示和python代码）。作为一个例子，你可以查看下我的文章《用python求解薛定谔方程》，是如何利用FFT将原本复杂的微分方程简化。 正因为FFT在那么多领域里如此有用，python提供了很多标准工具和封装来计算它。NumPy 和 SciPy 都有经过充分测试的封装好的FFT库，分别位于子模块 numpy.fft 和 scipy.fftpack 。我所知的最快的FFT是在 FFTW包中 ，而你也可以在python的pyFFTW 包中使用它。 虽然说了这么远，但还是暂时先将这些库放一边，考虑一下怎样使用原始的python从头开始计算FFT。 计算离散傅里叶变换简单起见，我们只关心正变换，因为逆变换也只是以很相似的方式就能做到。看一下上面的DFT表达式，它只是一个直观的线性运算：向量x的矩阵乘法，
矩阵M可以表示为
这么想的话，我们可以简单地利用矩阵乘法计算DFT：
对比numpy的内置FFT函数，我们来对结果进行仔细检查，
输出：
现在为了验证我们的算法有多慢，对比下两者的执行时间
输出：
使用这种简化的实现方法，正如所料，我们慢了一千多倍。但问题不是这么简单。对于长度为N的输入矢量，FFT是O(N logN)级的，而我们的慢算法是O(N^2)级的。这就意味着，FFT用50毫秒能干完的活，对于我们的慢算法来说，要差不多20小时！ 那么FFT是怎么提速完事的呢？答案就在于他利用了对称性。 离散傅里叶变换中的对称性算法设计者所掌握的最重要手段之一，就是利用问题的对称性。如果你能清晰地展示问题的某一部分与另一部分相关，那么你就只需计算子结果一次，从而节省了计算成本。 Cooley 和 Tukey 正是使用这种方法导出FFT。 首先我们来看下
对于所有的整数n，exp[2π i n]=1。 最后一行展示了DFT很好的对称性：
简单地拓展一下：
同理对于所有整数 i 。正如下面即将看到的，这个对称性能被利用于更快地计算DFT。 DFT 到 FFT： 利用对称性 Cooley 和 Tukey 证明了，DFT的计算可以分为两部分。从DFT的定义得：
我们将单个DFT分成了看起来相似两个更小的DFT。一个奇，一个偶。目前为止，我们还没有节省计算开销，每一部分都包含(N/2)∗N的计算量，总的来说，就是N^2 。 技巧就是对每一部分利用对称性。因为 k 的范围是 0≤k<N ， 而 n 的范围是 0≤n<M≡N/2 ， 从上面的对称性特点来看，我们只需对每个子问题作一半的计算量。O(N^2) 变成了 O(M^2) 。 但我们不是到这步就停下来，只要我们的小傅里叶变换是偶倍数，就可以再作分治，直到分解出来的子问题小到无法通过分治提高效率，接近极限时，这个递归是 O(n logn) 级的。 这个递归算法能在python里快速实现，当子问题被分解到合适大小时，再用回原本那种“慢方法”。
|
