从零手写FFT用Python代码理解快速傅里叶变换的魔法傅里叶变换是数字信号处理领域的基石而快速傅里叶变换(FFT)则是让这个数学工具真正走向工程实践的关键突破。本文将带你从零开始用Python代码实现FFT算法通过实践理解这个改变世界的算法。1. 为什么需要理解FFT在开始编码之前让我们先思考一个基本问题为什么FFT如此重要想象一下你正在开发一个音乐识别应用需要分析音频信号的频率成分或者你是一名工程师需要从传感器数据中提取特征。这些场景都离不开频谱分析而FFT正是实现这一目标的利器。传统离散傅里叶变换(DFT)的时间复杂度是O(N²)而FFT通过巧妙的分解将其降低到O(N log N)。当N较大时这种效率提升是惊人的数据点数量(N)DFT计算量(N²)FFT计算量(N log₂N)加速比25665,5362,04832x10241,048,57610,240102x409616,777,21649,152341x2. 从DFT到FFT分而治之的智慧2.1 朴素DFT的实现我们先从最基础的DFT实现开始这将成为我们理解FFT的起点import numpy as np import matplotlib.pyplot as plt def naive_dft(x): 朴素的离散傅里叶变换实现 N len(x) n np.arange(N) k n.reshape((N, 1)) M np.exp(-2j * np.pi * k * n / N) return np.dot(M, x)这个实现直接根据DFT的数学定义 $$ X[k] \sum_{n0}^{N-1} x[n] \cdot e^{-i 2\pi k n / N} $$2.2 DFT的性能瓶颈让我们测试一下这个朴素实现的性能# 生成测试信号 N 256 # 信号长度 t np.linspace(0, 1, N) freq 5 # 信号频率 x np.sin(2 * np.pi * freq * t) 0.5 * np.random.randn(N) # 计算DFT并计时 %timeit naive_dft(x)在我的笔记本上对于N256的信号朴素DFT大约需要10毫秒。虽然看起来不错但当N增加到4096时计算时间将增加到约1.7秒——这对于实时处理来说太慢了。3. FFT的递归实现3.1 分治思想的核心FFT的魔力在于它利用了DFT的对称性和周期性将问题分解为更小的子问题。关键观察点是将序列分为奇数项和偶数项分别计算奇数项和偶数项的DFT合并结果得到完整的DFT数学表达式为 $$ X[k] E[k] e^{-i2πk/N} \cdot O[k] $$ $$ X[kN/2] E[k] - e^{-i2πk/N} \cdot O[k] $$其中E[k]是偶数项的DFTO[k]是奇数项的DFT。3.2 Python递归实现def recursive_fft(x): 递归实现的FFT N len(x) if N 1: return x # 分离奇偶项 even recursive_fft(x[::2]) odd recursive_fft(x[1::2]) # 计算旋转因子 T [np.exp(-2j * np.pi * k / N) * odd[k] for k in range(N // 2)] # 合并结果 return [even[k] T[k] for k in range(N // 2)] \ [even[k] - T[k] for k in range(N // 2)]让我们验证这个实现的正确性# 比较FFT和DFT的结果 fft_result recursive_fft(x) dft_result naive_dft(x) print(最大差异:, np.max(np.abs(fft_result - dft_result))) # 应该非常小3.3 递归实现的局限性虽然递归实现直观易懂但它有一些缺点Python的递归深度限制递归调用的开销需要额外的内存空间对于较大的N我们更需要迭代实现。4. FFT的迭代优化4.1 位反转排列迭代FFT的第一步是将输入数据重新排列为位反转顺序。这相当于递归实现中不断将数据分为奇偶项的自然结果。def bit_reverse_copy(x): 生成位反转排列 N len(x) bit_reversed np.zeros(N, dtypecomplex) for k in range(N): # 计算k的位反转 rev_k int({:0{}b}.format(k, int(np.log2(N)))[::-1], 2) bit_reversed[rev_k] x[k] return bit_reversed4.2 迭代FFT实现def iterative_fft(x): 迭代优化的FFT实现 N len(x) if (N (N - 1)) ! 0: raise ValueError(长度必须是2的幂次) # 位反转排列 x bit_reverse_copy(x) # 迭代计算 s 1 # 当前处理的子问题大小 while s N // 2: for k in range(0, N, 2 * s): for j in range(s): # 蝴蝶操作 w np.exp(-2j * np.pi * j / (2 * s)) u x[k j] v x[k j s] * w x[k j] u v x[k j s] u - v s * 2 return x4.3 性能对比让我们比较三种实现的性能N 1024 x np.random.randn(N) %timeit naive_dft(x) # 朴素DFT %timeit recursive_fft(x) # 递归FFT %timeit iterative_fft(x) # 迭代FFT %timeit np.fft.fft(x) # NumPy的优化FFT在我的测试中对于N1024朴素DFT约170ms递归FFT约12ms迭代FFT约5msNumPy FFT约0.02ms可以看到我们的迭代实现比朴素DFT快约34倍而NumPy的实现又比我们的快约250倍——这展示了高度优化的工业级实现的威力。5. 逆FFT实现有了FFT实现逆FFT(IFFT)非常简单只需稍作修改def iterative_ifft(X): 迭代优化的IFFT实现 N len(X) if (N (N - 1)) ! 0: raise ValueError(长度必须是2的幂次) # 位反转排列 X bit_reverse_copy(X) # 迭代计算 s 1 while s N // 2: for k in range(0, N, 2 * s): for j in range(s): # 蝴蝶操作使用共轭旋转因子 w np.exp(2j * np.pi * j / (2 * s)) # 注意符号变化 u X[k j] v X[k j s] * w X[k j] u v X[k j s] u - v s * 2 # 最后要除以N return X / N验证IFFT的正确性# 生成随机信号 x np.random.randn(256) 1j * np.random.randn(256) # FFT然后IFFT X iterative_fft(x) x_recon iterative_ifft(X) # 检查重建误差 print(最大重建误差:, np.max(np.abs(x - x_recon))) # 应该非常小6. 实际应用频谱分析现在让我们用自己实现的FFT来分析真实信号的频谱# 生成含噪声的信号 N 1024 t np.linspace(0, 1, N) freqs [10, 20, 30] # 信号包含的频率成分 x sum(np.sin(2 * np.pi * f * t) for f in freqs) 0.5 * np.random.randn(N) # 计算频谱 X iterative_fft(x) f np.linspace(0, N, N) # 频率轴 # 绘制结果 plt.figure(figsize(12, 6)) plt.subplot(121) plt.plot(t, x) plt.title(时域信号) plt.xlabel(时间) plt.ylabel(幅度) plt.subplot(122) plt.plot(f[:N//2], np.abs(X[:N//2])) # 只显示正频率部分 plt.title(频域谱) plt.xlabel(频率) plt.ylabel(幅度) plt.tight_layout() plt.show()这段代码会显示信号的时域波形和频域谱你应该能在频谱上清晰地看到我们在信号中加入的10Hz、20Hz和30Hz成分。7. 性能优化技巧虽然我们的迭代实现已经比朴素DFT快很多但还有优化空间7.1 预先计算旋转因子def optimized_fft(x): N len(x) logN int(np.log2(N)) # 预先计算所有可能的旋转因子 W [np.exp(-2j * np.pi * k / N) for k in range(N // 2)] # 位反转排列 x bit_reverse_copy(x) # 迭代计算 s 1 while s N // 2: w_step N // (2 * s) for k in range(0, N, 2 * s): for j in range(s): w W[j * w_step] u x[k j] v x[k j s] * w x[k j] u v x[k j s] u - v s * 2 return x7.2 使用内存视图对于大型数组使用内存视图可以减少复制开销def memory_efficient_fft(x): N len(x) x x.copy() # 避免修改原数组 # 位反转排列 for k in range(N): rev_k int({:0{}b}.format(k, int(np.log2(N)))[::-1], 2) if k rev_k: x[k], x[rev_k] x[rev_k], x[k] # 迭代计算 s 1 while s N // 2: w_step N // (2 * s) for k in range(0, N, 2 * s): for j in range(s): w np.exp(-2j * np.pi * j / (2 * s)) u x[k j] v x[k j s] * w x[k j] u v x[k j s] u - v s * 2 return x8. 常见问题与调试技巧实现FFT时可能会遇到各种问题以下是一些常见问题及解决方法结果不正确检查旋转因子的符号FFT用负IFFT用正确保输入长度是2的幂次验证位反转排列是否正确数值不稳定对于大N累积的浮点误差可能显著考虑使用更高精度的浮点类型性能不如预期避免在循环中重复计算旋转因子尽量减少数组复制操作对于非常大的N考虑分块处理9. 扩展与进阶掌握了基本FFT后你可以进一步探索多维FFT将一维FFT扩展到图像处理等二维场景非2的幂次长度研究混合基FFT或Chirp-Z变换稀疏FFT针对稀疏信号的优化算法并行FFT利用GPU或多核CPU加速计算10. 完整代码示例以下是本文讨论的所有FFT实现的完整代码import numpy as np import matplotlib.pyplot as plt def naive_dft(x): N len(x) n np.arange(N) k n.reshape((N, 1)) M np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) def recursive_fft(x): N len(x) if N 1: return x even recursive_fft(x[::2]) odd recursive_fft(x[1::2]) T [np.exp(-2j * np.pi * k / N) * odd[k] for k in range(N // 2)] return [even[k] T[k] for k in range(N // 2)] \ [even[k] - T[k] for k in range(N // 2)] def bit_reverse_copy(x): N len(x) bit_reversed np.zeros(N, dtypecomplex) for k in range(N): rev_k int({:0{}b}.format(k, int(np.log2(N)))[::-1], 2) bit_reversed[rev_k] x[k] return bit_reversed def iterative_fft(x): N len(x) if (N (N - 1)) ! 0: raise ValueError(长度必须是2的幂次) x bit_reverse_copy(x) s 1 while s N // 2: for k in range(0, N, 2 * s): for j in range(s): w np.exp(-2j * np.pi * j / (2 * s)) u x[k j] v x[k j s] * w x[k j] u v x[k j s] u - v s * 2 return x def iterative_ifft(X): N len(X) if (N (N - 1)) ! 0: raise ValueError(长度必须是2的幂次) X bit_reverse_copy(X) s 1 while s N // 2: for k in range(0, N, 2 * s): for j in range(s): w np.exp(2j * np.pi * j / (2 * s)) u X[k j] v X[k j s] * w X[k j] u v X[k j s] u - v s * 2 return X / N # 测试代码 if __name__ __main__: # 生成测试信号 N 256 t np.linspace(0, 1, N) x np.sin(2 * np.pi * 10 * t) 0.5 * np.random.randn(N) # 计算并比较不同实现 X_naive naive_dft(x) X_rec recursive_fft(x) X_iter iterative_fft(x) print(朴素DFT与递归FFT最大差异:, np.max(np.abs(X_naive - X_rec))) print(朴素DFT与迭代FFT最大差异:, np.max(np.abs(X_naive - X_iter))) # 绘制频谱 plt.plot(np.abs(X_iter)[:N//2]) plt.title(FFT频谱) plt.xlabel(频率) plt.ylabel(幅度) plt.show()通过这个完整的实现过程我们不仅理解了FFT的数学原理还掌握了如何将其转化为高效的代码。这种从理论到实践的转换能力正是工程师和科学家的核心技能之一。
别再死记硬背公式了！用Python手搓一个FFT，彻底搞懂快速傅里叶变换（附完整代码）
从零手写FFT用Python代码理解快速傅里叶变换的魔法傅里叶变换是数字信号处理领域的基石而快速傅里叶变换(FFT)则是让这个数学工具真正走向工程实践的关键突破。本文将带你从零开始用Python代码实现FFT算法通过实践理解这个改变世界的算法。1. 为什么需要理解FFT在开始编码之前让我们先思考一个基本问题为什么FFT如此重要想象一下你正在开发一个音乐识别应用需要分析音频信号的频率成分或者你是一名工程师需要从传感器数据中提取特征。这些场景都离不开频谱分析而FFT正是实现这一目标的利器。传统离散傅里叶变换(DFT)的时间复杂度是O(N²)而FFT通过巧妙的分解将其降低到O(N log N)。当N较大时这种效率提升是惊人的数据点数量(N)DFT计算量(N²)FFT计算量(N log₂N)加速比25665,5362,04832x10241,048,57610,240102x409616,777,21649,152341x2. 从DFT到FFT分而治之的智慧2.1 朴素DFT的实现我们先从最基础的DFT实现开始这将成为我们理解FFT的起点import numpy as np import matplotlib.pyplot as plt def naive_dft(x): 朴素的离散傅里叶变换实现 N len(x) n np.arange(N) k n.reshape((N, 1)) M np.exp(-2j * np.pi * k * n / N) return np.dot(M, x)这个实现直接根据DFT的数学定义 $$ X[k] \sum_{n0}^{N-1} x[n] \cdot e^{-i 2\pi k n / N} $$2.2 DFT的性能瓶颈让我们测试一下这个朴素实现的性能# 生成测试信号 N 256 # 信号长度 t np.linspace(0, 1, N) freq 5 # 信号频率 x np.sin(2 * np.pi * freq * t) 0.5 * np.random.randn(N) # 计算DFT并计时 %timeit naive_dft(x)在我的笔记本上对于N256的信号朴素DFT大约需要10毫秒。虽然看起来不错但当N增加到4096时计算时间将增加到约1.7秒——这对于实时处理来说太慢了。3. FFT的递归实现3.1 分治思想的核心FFT的魔力在于它利用了DFT的对称性和周期性将问题分解为更小的子问题。关键观察点是将序列分为奇数项和偶数项分别计算奇数项和偶数项的DFT合并结果得到完整的DFT数学表达式为 $$ X[k] E[k] e^{-i2πk/N} \cdot O[k] $$ $$ X[kN/2] E[k] - e^{-i2πk/N} \cdot O[k] $$其中E[k]是偶数项的DFTO[k]是奇数项的DFT。3.2 Python递归实现def recursive_fft(x): 递归实现的FFT N len(x) if N 1: return x # 分离奇偶项 even recursive_fft(x[::2]) odd recursive_fft(x[1::2]) # 计算旋转因子 T [np.exp(-2j * np.pi * k / N) * odd[k] for k in range(N // 2)] # 合并结果 return [even[k] T[k] for k in range(N // 2)] \ [even[k] - T[k] for k in range(N // 2)]让我们验证这个实现的正确性# 比较FFT和DFT的结果 fft_result recursive_fft(x) dft_result naive_dft(x) print(最大差异:, np.max(np.abs(fft_result - dft_result))) # 应该非常小3.3 递归实现的局限性虽然递归实现直观易懂但它有一些缺点Python的递归深度限制递归调用的开销需要额外的内存空间对于较大的N我们更需要迭代实现。4. FFT的迭代优化4.1 位反转排列迭代FFT的第一步是将输入数据重新排列为位反转顺序。这相当于递归实现中不断将数据分为奇偶项的自然结果。def bit_reverse_copy(x): 生成位反转排列 N len(x) bit_reversed np.zeros(N, dtypecomplex) for k in range(N): # 计算k的位反转 rev_k int({:0{}b}.format(k, int(np.log2(N)))[::-1], 2) bit_reversed[rev_k] x[k] return bit_reversed4.2 迭代FFT实现def iterative_fft(x): 迭代优化的FFT实现 N len(x) if (N (N - 1)) ! 0: raise ValueError(长度必须是2的幂次) # 位反转排列 x bit_reverse_copy(x) # 迭代计算 s 1 # 当前处理的子问题大小 while s N // 2: for k in range(0, N, 2 * s): for j in range(s): # 蝴蝶操作 w np.exp(-2j * np.pi * j / (2 * s)) u x[k j] v x[k j s] * w x[k j] u v x[k j s] u - v s * 2 return x4.3 性能对比让我们比较三种实现的性能N 1024 x np.random.randn(N) %timeit naive_dft(x) # 朴素DFT %timeit recursive_fft(x) # 递归FFT %timeit iterative_fft(x) # 迭代FFT %timeit np.fft.fft(x) # NumPy的优化FFT在我的测试中对于N1024朴素DFT约170ms递归FFT约12ms迭代FFT约5msNumPy FFT约0.02ms可以看到我们的迭代实现比朴素DFT快约34倍而NumPy的实现又比我们的快约250倍——这展示了高度优化的工业级实现的威力。5. 逆FFT实现有了FFT实现逆FFT(IFFT)非常简单只需稍作修改def iterative_ifft(X): 迭代优化的IFFT实现 N len(X) if (N (N - 1)) ! 0: raise ValueError(长度必须是2的幂次) # 位反转排列 X bit_reverse_copy(X) # 迭代计算 s 1 while s N // 2: for k in range(0, N, 2 * s): for j in range(s): # 蝴蝶操作使用共轭旋转因子 w np.exp(2j * np.pi * j / (2 * s)) # 注意符号变化 u X[k j] v X[k j s] * w X[k j] u v X[k j s] u - v s * 2 # 最后要除以N return X / N验证IFFT的正确性# 生成随机信号 x np.random.randn(256) 1j * np.random.randn(256) # FFT然后IFFT X iterative_fft(x) x_recon iterative_ifft(X) # 检查重建误差 print(最大重建误差:, np.max(np.abs(x - x_recon))) # 应该非常小6. 实际应用频谱分析现在让我们用自己实现的FFT来分析真实信号的频谱# 生成含噪声的信号 N 1024 t np.linspace(0, 1, N) freqs [10, 20, 30] # 信号包含的频率成分 x sum(np.sin(2 * np.pi * f * t) for f in freqs) 0.5 * np.random.randn(N) # 计算频谱 X iterative_fft(x) f np.linspace(0, N, N) # 频率轴 # 绘制结果 plt.figure(figsize(12, 6)) plt.subplot(121) plt.plot(t, x) plt.title(时域信号) plt.xlabel(时间) plt.ylabel(幅度) plt.subplot(122) plt.plot(f[:N//2], np.abs(X[:N//2])) # 只显示正频率部分 plt.title(频域谱) plt.xlabel(频率) plt.ylabel(幅度) plt.tight_layout() plt.show()这段代码会显示信号的时域波形和频域谱你应该能在频谱上清晰地看到我们在信号中加入的10Hz、20Hz和30Hz成分。7. 性能优化技巧虽然我们的迭代实现已经比朴素DFT快很多但还有优化空间7.1 预先计算旋转因子def optimized_fft(x): N len(x) logN int(np.log2(N)) # 预先计算所有可能的旋转因子 W [np.exp(-2j * np.pi * k / N) for k in range(N // 2)] # 位反转排列 x bit_reverse_copy(x) # 迭代计算 s 1 while s N // 2: w_step N // (2 * s) for k in range(0, N, 2 * s): for j in range(s): w W[j * w_step] u x[k j] v x[k j s] * w x[k j] u v x[k j s] u - v s * 2 return x7.2 使用内存视图对于大型数组使用内存视图可以减少复制开销def memory_efficient_fft(x): N len(x) x x.copy() # 避免修改原数组 # 位反转排列 for k in range(N): rev_k int({:0{}b}.format(k, int(np.log2(N)))[::-1], 2) if k rev_k: x[k], x[rev_k] x[rev_k], x[k] # 迭代计算 s 1 while s N // 2: w_step N // (2 * s) for k in range(0, N, 2 * s): for j in range(s): w np.exp(-2j * np.pi * j / (2 * s)) u x[k j] v x[k j s] * w x[k j] u v x[k j s] u - v s * 2 return x8. 常见问题与调试技巧实现FFT时可能会遇到各种问题以下是一些常见问题及解决方法结果不正确检查旋转因子的符号FFT用负IFFT用正确保输入长度是2的幂次验证位反转排列是否正确数值不稳定对于大N累积的浮点误差可能显著考虑使用更高精度的浮点类型性能不如预期避免在循环中重复计算旋转因子尽量减少数组复制操作对于非常大的N考虑分块处理9. 扩展与进阶掌握了基本FFT后你可以进一步探索多维FFT将一维FFT扩展到图像处理等二维场景非2的幂次长度研究混合基FFT或Chirp-Z变换稀疏FFT针对稀疏信号的优化算法并行FFT利用GPU或多核CPU加速计算10. 完整代码示例以下是本文讨论的所有FFT实现的完整代码import numpy as np import matplotlib.pyplot as plt def naive_dft(x): N len(x) n np.arange(N) k n.reshape((N, 1)) M np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) def recursive_fft(x): N len(x) if N 1: return x even recursive_fft(x[::2]) odd recursive_fft(x[1::2]) T [np.exp(-2j * np.pi * k / N) * odd[k] for k in range(N // 2)] return [even[k] T[k] for k in range(N // 2)] \ [even[k] - T[k] for k in range(N // 2)] def bit_reverse_copy(x): N len(x) bit_reversed np.zeros(N, dtypecomplex) for k in range(N): rev_k int({:0{}b}.format(k, int(np.log2(N)))[::-1], 2) bit_reversed[rev_k] x[k] return bit_reversed def iterative_fft(x): N len(x) if (N (N - 1)) ! 0: raise ValueError(长度必须是2的幂次) x bit_reverse_copy(x) s 1 while s N // 2: for k in range(0, N, 2 * s): for j in range(s): w np.exp(-2j * np.pi * j / (2 * s)) u x[k j] v x[k j s] * w x[k j] u v x[k j s] u - v s * 2 return x def iterative_ifft(X): N len(X) if (N (N - 1)) ! 0: raise ValueError(长度必须是2的幂次) X bit_reverse_copy(X) s 1 while s N // 2: for k in range(0, N, 2 * s): for j in range(s): w np.exp(2j * np.pi * j / (2 * s)) u X[k j] v X[k j s] * w X[k j] u v X[k j s] u - v s * 2 return X / N # 测试代码 if __name__ __main__: # 生成测试信号 N 256 t np.linspace(0, 1, N) x np.sin(2 * np.pi * 10 * t) 0.5 * np.random.randn(N) # 计算并比较不同实现 X_naive naive_dft(x) X_rec recursive_fft(x) X_iter iterative_fft(x) print(朴素DFT与递归FFT最大差异:, np.max(np.abs(X_naive - X_rec))) print(朴素DFT与迭代FFT最大差异:, np.max(np.abs(X_naive - X_iter))) # 绘制频谱 plt.plot(np.abs(X_iter)[:N//2]) plt.title(FFT频谱) plt.xlabel(频率) plt.ylabel(幅度) plt.show()通过这个完整的实现过程我们不仅理解了FFT的数学原理还掌握了如何将其转化为高效的代码。这种从理论到实践的转换能力正是工程师和科学家的核心技能之一。
