别再死记硬背公式了！用Python手搓一个FFT，从分治到蝶形运算一次搞懂

用Python手搓FFT从分治到蝶形运算的视觉化之旅当我在信号处理课上第一次接触快速傅里叶变换FFT时那些复杂的数学公式让我望而生畏。直到有一天我决定用Python从零实现这个算法才真正理解了它的精妙之处——这就像用乐高积木搭建一座城堡每个蝶形运算单元都是不可或缺的零件。1. 为什么需要理解FFT的实现傅里叶变换是数字信号处理的基石而FFT将这个O(n²)的计算优化到O(n log n)。但教科书上的数学推导往往让人迷失在公式海洋中。通过代码实现你会发现分治思想的直观体现将大问题分解为相同的小问题蝶形运算的优雅结构复数旋转因子的精妙应用位反转置换的计算机友好特性内存访问模式的优化import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation2. 从简单案例开始8点FFT的解剖让我们从一个长度为8的实数序列开始这足够小到可以手动计算又足够大来展示所有关键概念。输入序列示例x np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])2.1 分治策略的可视化FFT的核心是将DFT分解为奇偶两部分。对于N点FFT$$ X[k] E[k] e^{-j2πk/N}O[k] $$其中E[k]是偶数点DFTO[k]是奇数点DFT。这个过程会递归进行直到最基本的2点DFT。递归分组过程递归层级分组情况初始[0,1,2,3,4,5,6,7]第1次[[0,2,4,6], [1,3,5,7]]第2次[[[0,4], [2,6]], [[1,5], [3,7]]]终止条件单个2点DFTdef visualize_division(x): fig, ax plt.subplots(figsize(10,4)) markers [o,s,D,^] for i in range(int(np.log2(len(x)))): groups [x[j::2**i] for j in range(2**i)] for j, group in enumerate(groups): ax.plot(group, markermarkers[i], linestyle, labelfLevel {i} Group {j}) ax.set_title(FFT Recursive Grouping Process) ax.legend() plt.show()3. 蝶形运算FFT的原子操作蝶形运算(Butterfly Operation)是FFT的基本计算单元得名于其数据流图 resembling a butterfly.基本运算规则A A W·B B A - W·B其中W是旋转因子(Twiddle Factor)$W_N^k e^{-j2πk/N}$3.1 旋转因子的计算与可视化def twiddle_factor(k, N): return np.exp(-2j * np.pi * k / N) theta np.linspace(0, 2*np.pi, 8, endpointFalse) W np.exp(-1j * theta) plt.figure(figsize(6,6)) plt.scatter(np.real(W), np.imag(W)) plt.title(Twiddle Factors on Unit Circle) plt.grid(True) for i, w in enumerate(W): plt.annotate(fW^{i}, (np.real(w), np.imag(w)))3.2 蝶形运算的Python实现def butterfly(x, W): N len(x) X np.zeros(N, dtypenp.complex128) half N // 2 for k in range(half): X[k] x[k] W[k] * x[k half] X[k half] x[k] - W[k] * x[k half] return X4. 完整FFT实现递归与迭代对比4.1 递归实现最直观的表达def recursive_fft(x): N len(x) if N 1: return x even recursive_fft(x[::2]) odd recursive_fft(x[1::2]) W np.exp(-2j * np.pi * np.arange(N//2) / N) return np.concatenate([ even W * odd, even - W * odd ])递归调用树示例N8FFT8 ├── FFT4(even) │ ├── FFT2(even) │ └── FFT2(odd) └── FFT4(odd) ├── FFT2(even) └── FFT2(odd)4.2 迭代实现更高效的版本迭代实现通过位反转重排和分层计算避免了递归开销。def iterative_fft(x): N len(x) X np.array(x, dtypenp.complex128) # 位反转重排 j 0 for i in range(1, N): bit N 1 while j bit: j - bit bit 1 j bit if i j: X[i], X[j] X[j], X[i] # 分层计算 L 2 while L N: W_L np.exp(-2j * np.pi / L) for k in range(0, N, L): w 1 for j in range(L//2): u X[k j] v w * X[k j L//2] X[k j] u v X[k j L//2] u - v w * W_L L 1 return X性能对比表实现方式时间复杂度空间复杂度适用场景递归实现O(n log n)O(n log n)教学、理解算法迭代实现O(n log n)O(n)实际生产环境使用5. 逆FFT的实现IFFT与FFT几乎相同只有两处关键区别旋转因子取共轭$W_N^k e^{j2πk/N}$最后结果要除以Ndef recursive_ifft(X): N len(X) if N 1: return X even recursive_ifft(X[::2]) odd recursive_ifft(X[1::2]) W np.exp(2j * np.pi * np.arange(N//2) / N) # 关键变化 return np.concatenate([ even W * odd, even - W * odd ]) / 2 # 合并时除以2最终相当于除以N6. 验证与可视化让数学变得可见6.1 验证FFT/IFFT的正确性x np.random.random(8) X recursive_fft(x) x_recon recursive_ifft(X) print(原始信号:, x) print(重建信号:, np.real(x_recon)) # 虚部应为近似0 print(最大误差:, np.max(np.abs(x - np.real(x_recon))))6.2 蝶形运算的动画演示def animate_butterfly(x): fig, ax plt.subplots(figsize(10,6)) N len(x) stages int(np.log2(N)) # 初始化绘图元素 lines [] for i in range(N): line, ax.plot([], [], o-, labelfX[{i}]) lines.append(line) def init(): ax.set_xlim(-1, stages1) ax.set_ylim(-1.5, 1.5) return lines def update(frame): # 这里简化了实际动画逻辑 for i in range(N): lines[i].set_data([frame, frame1], [0, 0]) return lines ani FuncAnimation(fig, update, framesstages, init_funcinit, blitTrue) plt.title(FFT Butterfly Operations) plt.show() return ani7. 实际应用音频频谱分析示例让我们用自实现的FFT分析真实音频信号from scipy.io import wavfile rate, data wavfile.read(audio.wav) mono data.mean(axis1) if len(data.shape) 1 else data N 1024 # 分析窗口大小 freq np.fft.fftfreq(N, 1/rate)[:N//2] # 使用我们实现的FFT window np.hamming(N) X recursive_fft(mono[:N] * window)[:N//2] magnitude 20 * np.log10(np.abs(X)) plt.plot(freq, magnitude) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude (dB)) plt.title(Audio Spectrum (Using Our FFT))8. 性能优化与进阶话题虽然我们的实现易于理解但与专业库相比还有差距预计算旋转因子避免重复计算SIMD指令利用现代CPU的并行计算能力多线程处理分治法的天然并行性混合基数FFT不限于基2的情况# 预计算旋转因子的优化示例 def optimized_fft(x): N len(x) logN int(np.log2(N)) # 预计算所有可能用到的旋转因子 W_dict {} for s in range(1, logN1): m 1 s W_dict[m] np.exp(-2j * np.pi * np.arange(m//2) / m) # ...剩余部分与迭代实现类似...在实现FFT的过程中最让我惊讶的是如此简单的蝶形操作组合起来竟能完成如此强大的变换。当第一次看到自己实现的FFT正确分析出音频信号的频率成分时那种成就感是单纯理解数学公式无法比拟的。