快速傅里叶变换（FFT）是离散傅里叶变换（DFT）的高效实现算法，核心思想是通过分治策略减少重复计算，将原本O(n²)复杂度的DFT运算降低到O(n log n)。在信号处理、频谱分析等场景中应用十分广泛，其实现的核心是蝶形运算单元和频域转换的位逆序重排逻辑。
FFT核心原理与蝶形运算逻辑
FFT基于DFT的对称性和周期性拆分计算，最常见的基2 FFT要求输入数据长度为2的整数次幂。蝶形运算是FFT的基本计算单元，每一次蝶形运算会处理两个输入数据，输出两个结果，核心公式如下：
设旋转因子为W_n^k = cos(2πk/n) - i sin(2πk/n)，两个输入为a、b，则蝶形运算输出为：
- 输出1：a + W_n^k * b
- 输出2：a - W_n^k * b
蝶形运算的层级与数据长度相关，长度为N=2^m的序列需要进行m级蝶形运算，每级的蝶形分组间隔和旋转因子索引都会发生变化。
位逆序重排
FFT输入数据需要先进行位逆序重排，才能保证蝶形运算后输出是自然顺序的频域结果。位逆序指的是将数据的索引二进制位反转，比如长度为8的序列，索引3（二进制011）反转后是6（二进制110）。
C++实现FFT完整代码
下面的代码实现了基2的按时间抽取FFT，包含复数定义、位逆序重排、蝶形运算和最终的FFT函数，能够将时域复数序列转换为频域序列。
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
// 定义复数结构体
struct Complex {
double real; // 实部
double imag; // 虚部
Complex(double r = 0, double i = 0) : real(r), imag(i) {}
};
// 复数乘法
Complex mul(const Complex& a, const Complex& b) {
return Complex(
a.real * b.real - a.imag * b.imag,
a.real * b.imag + a.imag * b.real
);
}
// 复数加法
Complex add(const Complex& a, const Complex& b) {
return Complex(a.real + b.real, a.imag + b.imag);
}
// 复数减法
Complex sub(const Complex& a, const Complex& b) {
return Complex(a.real - b.real, a.imag - b.imag);
}
// 计算旋转因子 W_N^k
Complex getTwiddle(int n, int k) {
double theta = 2 * M_PI * k / n;
return Complex(cos(theta), -sin(theta));
}
// 位逆序重排函数
void bitReverse(vector<Complex>& data) {
int n = data.size();
int rev = 0;
for (int i = 0; i < n; ++i) {
if (i < rev) {
swap(data[i], data[rev]);
}
// 计算下一个逆序索引
int mask = n >> 1;
while (rev & mask) {
rev ^= mask;
mask >>= 1;
}
rev |= mask;
}
}
// FFT主函数，输入为时域复数序列，输出为频域复数序列
void fft(vector<Complex>& data) {
int n = data.size();
// 先进行位逆序重排
bitReverse(data);
// 迭代每一级蝶形运算，共log2(n)级
for (int len = 2; len <= n; len <<= 1) {
int halfLen = len >> 1;
// 当前级的旋转因子步长
for (int i = 0; i < n; i += len) {
for (int j = 0; j < halfLen; ++j) {
Complex w = getTwiddle(len, j);
Complex u = data[i + j];
Complex v = mul(w, data[i + j + halfLen]);
// 蝶形运算
data[i + j] = add(u, v);
data[i + j + halfLen] = sub(u, v);
}
}
}
}
int main() {
// 测试数据：长度为8的时域序列
vector<Complex> testData = {
Complex(1, 0), Complex(2, 0), Complex(3, 0), Complex(4, 0),
Complex(5, 0), Complex(6, 0), Complex(7, 0), Complex(8, 0)
};
cout << "时域输入数据：" << endl;
for (auto& c : testData) {
cout << "(" << c.real << ", " << c.imag << ") ";
}
cout << endl;
// 执行FFT
fft(testData);
cout << "频域输出结果：" << endl;
for (auto& c : testData) {
cout << "(" << c.real << ", " << c.imag << ") ";
}
cout << endl;
return 0;
}
代码逻辑解析
上述代码的核心流程分为三步：
- 首先通过
bitReverse函数对输入序列进行位逆序重排，保证后续蝶形运算的顺序正确。 - 外层循环控制蝶形运算的级数，从长度为2的分组开始，每次分组长度翻倍，直到等于输入序列长度。
- 内层循环处理每一组内的蝶形运算，根据当前分组长度和偏移量计算旋转因子，完成两个输入数据的加减运算，得到新的输出结果。
运行代码后，输出的结果就是输入时域序列对应的频域复数表示，每个复数点的模长对应对应频率分量的幅值，相位对应频率分量的相位信息。
注意事项
- 输入数据长度必须是2的整数次幂，否则需要先进行补零操作到最近的2的整数次幂长度。
- 旋转因子的计算使用了三角函数，实际工程中可以通过查表的方式预存旋转因子，进一步提升运算效率。
- 如果需要实现逆FFT（从频域转回时域），只需要将旋转因子的虚部取反，最后结果除以序列长度即可。