当前位置：首页 > news >正文

别再死记硬背公式了！用Python手搓一个FFT，从蝴蝶操作到代码实现（附完整源码）

news 2026/6/13 15:57:25

从蝴蝶操作到完整实现：用Python手搓FFT的实战指南

为什么我们需要自己实现FFT？

在信号处理领域，快速傅里叶变换（FFT）就像是一把瑞士军刀，它能够将时域信号转换到频域，让我们看到信号背后的频率成分。虽然NumPy等库提供了现成的fft函数，但真正理解FFT的最佳方式就是亲手实现它。

我记得第一次在项目中需要使用FFT时，只是简单地调用了np.fft.fft()，结果出来了却不知道如何解释那些复数输出的实际意义。直到我亲手实现了FFT算法，才真正理解了时域和频域之间的转换关系。

FFT的核心：蝴蝶操作解析

蝴蝶操作（Butterfly Operation）是FFT算法中最关键的部分，得名于其数据流图看起来像蝴蝶的形状。这个操作本质上是一个简单的复数运算，它将两个复数通过旋转因子联系起来。

def butterfly(a, b, factor): """基础蝴蝶操作实现""" add = a + factor * b sub = a - factor * b return add, sub

理解这个操作的关键在于旋转因子（twiddle factor），它实际上是一个单位圆上的复数：

def get_twiddle_factor(k, N): """计算旋转因子 e^(-2πjk/N)""" angle = -2 * math.pi * k / N return complex(math.cos(angle), math.sin(angle))

递归实现：最直观的FFT算法

递归实现是最符合FFT分治思想的实现方式。它的基本步骤是：

将输入序列分成奇数位和偶数位两部分
对两部分分别递归调用FFT
将结果通过蝴蝶操作组合起来

def fft_recursive(x): N = len(x) if N <= 1: return x # 分治：奇偶分组 even = fft_recursive(x[0::2]) odd = fft_recursive(x[1::2]) # 合并结果 result = [0] * N for k in range(N//2): factor = get_twiddle_factor(k, N) result[k] = even[k] + factor * odd[k] result[k + N//2] = even[k] - factor * odd[k] return result

虽然递归实现直观易懂，但在实际应用中，我们更常用迭代实现，因为它的效率更高，且避免了递归调用的开销。

迭代实现：高效的FFT算法

迭代实现FFT需要先对输入序列进行位反转排列（Bit-reversal permutation），然后自底向上地进行蝴蝶操作。

位反转排列

def bit_reverse(n, bits): """计算n的位反转值""" reversed_n = 0 for i in range(bits): reversed_n = (reversed_n << 1) | (n & 1) n >>= 1 return reversed_n def fft_iterative(x): N = len(x) bits = N.bit_length() - 1 # 位反转排列 result = [x[bit_reverse(i, bits)] for i in range(N)] # 迭代计算FFT size = 2 while size <= N: half_size = size // 2 step = N // size for i in range(0, N, size): for j in range(half_size): factor = get_twiddle_factor(j * step, N) a = result[i + j] b = result[i + j + half_size] * factor result[i + j] = a + b result[i + j + half_size] = a - b size *= 2 return result

验证我们的实现：与NumPy对比

实现完成后，我们需要验证它的正确性。最直接的方式就是与NumPy的实现结果进行对比。

import numpy as np # 生成测试信号 t = np.linspace(0, 1, 256, endpoint=False) signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 10 * t) # 计算FFT fft_numpy = np.fft.fft(signal) fft_custom = fft_iterative(signal.tolist()) # 比较结果 print("最大差异:", np.max(np.abs(fft_numpy - fft_custom)))

性能优化技巧

虽然我们的实现已经可以工作，但在处理大数据量时可能会遇到性能问题。以下是几个优化方向：

预计算旋转因子：避免重复计算相同的旋转因子
使用内存视图：减少列表操作的开销
并行计算：利用多核CPU加速计算

def optimized_fft(x): N = len(x) bits = N.bit_length() - 1 # 预计算所有旋转因子 factors = [get_twiddle_factor(k, N) for k in range(N//2)] # 位反转排列 result = [x[bit_reverse(i, bits)] for i in range(N)] # 优化后的迭代计算 size = 2 while size <= N: half_size = size // 2 step = N // size for i in range(0, N, size): for j in range(half_size): factor = factors[j * step] a = result[i + j] b = result[i + j + half_size] * factor result[i + j] = a + b result[i + j + half_size] = a - b size *= 2 return result

实际应用案例：音频频谱分析

让我们用一个实际的例子来展示FFT的应用。我们将分析一段音频信号的频谱。

import wave import struct # 读取音频文件 with wave.open('audio.wav', 'rb') as wav: frames = wav.readframes(wav.getnframes()) samples = struct.unpack('{n}h'.format(n=wav.getnframes()*wav.getnchannels()), frames) sample_rate = wav.getframerate() # 计算FFT fft_result = optimized_fft(samples[:1024]) # 取前1024个样本 freqs = np.fft.fftfreq(1024, 1/sample_rate) magnitude = np.abs(fft_result) # 绘制频谱图 import matplotlib.pyplot as plt plt.plot(freqs[:512], magnitude[:512]) # 只显示正频率部分 plt.xlabel('Frequency (Hz)') plt.ylabel('Magnitude') plt.title('Audio Spectrum Analysis') plt.show()

常见问题与调试技巧

在实现FFT的过程中，可能会遇到以下问题：

结果不正确：检查位反转排列是否正确，旋转因子的计算是否有误
性能低下：确保避免了重复计算，使用更高效的数据结构
数值不稳定：对于大N值，浮点误差可能会累积，考虑使用更高精度的数据类型

提示：当FFT结果看起来不对时，可以从小的输入（如N=2或4）开始逐步调试，手动计算预期结果并与程序输出对比。

完整实现代码

以下是完整的FFT实现，包括递归和迭代两种方式，以及相应的逆变换：

import math import cmath from typing import List, Union def fft_recursive(x: List[complex]) -> List[complex]: """递归实现的FFT""" N = len(x) if N <= 1: return x even = fft_recursive(x[0::2]) odd = fft_recursive(x[1::2]) result = [0] * N for k in range(N//2): factor = cmath.exp(-2j * math.pi * k / N) result[k] = even[k] + factor * odd[k] result[k + N//2] = even[k] - factor * odd[k] return result def ifft_recursive(x: List[complex]) -> List[complex]: """递归实现的IFFT""" N = len(x) if N <= 1: return x even = ifft_recursive(x[0::2]) odd = ifft_recursive(x[1::2]) result = [0] * N for k in range(N//2): factor = cmath.exp(2j * math.pi * k / N) result[k] = even[k] + factor * odd[k] result[k + N//2] = even[k] - factor * odd[k] return [val / N for val in result] def fft_iterative(x: List[Union[float, complex]]) -> List[complex]: """迭代实现的FFT""" N = len(x) if N & (N - 1) != 0: raise ValueError("输入长度必须是2的幂") # 位反转排列 logN = int(math.log2(N)) result = [x[int('{:0{width}b}'.format(i, width=logN)[::-1], 2)] for i in range(N)] # 迭代计算 size = 2 while size <= N: half_size = size // 2 step = N // size for i in range(0, N, size): for j in range(half_size): factor = cmath.exp(-2j * math.pi * j * step / N) a = result[i + j] b = result[i + j + half_size] * factor result[i + j] = a + b result[i + j + half_size] = a - b size *= 2 return result def ifft_iterative(x: List[complex]) -> List[complex]: """迭代实现的IFFT""" N = len(x) if N & (N - 1) != 0: raise ValueError("输入长度必须是2的幂") # 位反转排列 logN = int(math.log2(N)) result = [x[int('{:0{width}b}'.format(i, width=logN)[::-1], 2)] for i in range(N)] # 迭代计算 size = 2 while size <= N: half_size = size // 2 step = N // size for i in range(0, N, size): for j in range(half_size): factor = cmath.exp(2j * math.pi * j * step / N) a = result[i + j] b = result[i + j + half_size] * factor result[i + j] = a + b result[i + j + half_size] = a - b size *= 2 return [val / N for val in result]