快速傅里叶变换（FFT）是数字信号处理里最基础的算法之一，它能把时域离散信号转成频域表示。在C++中自己实现一套FFT，既不依赖第三方库，也能帮助理解复数运算与分治思想。下面我们从数学原理出发，给出一个基于Cooley-Tukey基2算法的完整C++实现。
一、FFT的数学原理
离散傅里叶变换（DFT）的定义为：对长度为N的序列x[n]，其频域X[k]等于从n等于0到N减1的x[n]乘以e的负二倍πi kn除以N次方的求和。直接按这个公式算，每一个k都要循环N次，总体复杂度是O(N²)。当N大到几千或上万时，计算量就无法接受。
Cooley-Tukey算法发现，如果N是2的整数次幂，可以把序列按下标奇偶拆成两半。利用旋转因子W_N的周期性，原DFT变成两个长度为N/2的DFT组合，再乘上少许旋转因子。不断对半拆分，就得到深度为log₂N的递归树，总复杂度降到O(N log N)。这就是FFT提速的本质。
二、C++中的复数与基础结构
C++标准库在<complex>头文件中提供了std::complex模板，能直接用复数加减乘除以及实部和虚部访问。我们用std::complex<double>来表示信号采样点和频谱值，避免手写实数与虚数运算出错。
为了让代码清晰，先定义类型别名，并准备一个生成旋转因子的函数。旋转因子就是单位圆上的等分点，用std::polar生成模为1、角度为指定弧度的复数。下面这段代码展示了基础准备：
#include <complex>
#include <vector>
#include <cmath>
// 定义复数类型
typedef std::complex<double> Complex;
// 定义复数向量类型
typedef std::vector<Complex> CVector;
// 生成长度为n的旋转因子表
CVector genTwiddle(int n) {
CVector w(n);
for (int i = 0; i < n; ++i) {
// 第i个旋转因子角度
double angle = -2.0 * M_PI * i / n;
w[i] = std::polar(1.0, angle);
}
return w;
}
三、位反转排列
递归版FFT在拆分时自然形成了奇偶下标顺序，但迭代版需要先把输入按“位反转”重排。比如长度为8时，下标二进制001应放到100的位置。这样后续蝴蝶操作才能按相邻配对进行。
位反转可以用循环移位实现：从0开始，每次把已有反转值左移一位，并根据原下标最低位决定最高位填0还是1。下面给出迭代FFT前必须的位反转置换代码：
// 对长度为n（2的幂）的数组做位反转重排
void bitReverse(CVector& a) {
int n = a.size();
for (int i = 1, j = 0; i < n; ++i) {
int bit = n >> 1;
for (; j & bit; bit >>= 1) {
j ^= bit;
}
j ^= bit;
if (i < j) {
std::swap(a[i], a[j]);
}
}
}
四、迭代版FFT实现
做完位反转后，外层循环控制当前处理的蝴蝶长度len，从2一直翻倍到N；内层循环遍历每个长度为len的块，用预先算好的旋转因子做复数乘加。这种迭代写法没有递归调用，栈空间固定，适合嵌入式或高频调用场景。
注意正向变换不需要除以N，而逆变换要把结果除以N并取旋转因子共轭。下面给出完整正向FFT函数：
// 迭代版FFT，输入长度必须为2的幂
void fft(CVector& a) {
int n = a.size();
if (n <= 1) return;
bitReverse(a);
// 预先生成旋转因子
CVector w = genTwiddle(n);
for (int len = 2; len <= n; len <<= 1) {
int half = len >> 1;
for (int i = 0; i < n; i += len) {
for (int k = 0; k < half; ++k) {
Complex u = a[i + k];
// 取对应旋转因子
Complex t = a[i + k + half] * w[n / len * k];
a[i + k] = u + t;
a[i + k + half] = u - t;
}
}
}
}
五、使用示例与验证
我们可以用一个含两个频率的正弦波来测试。构造长度为1024的采样序列，分别加入50Hz与120Hz分量，跑完FFT后观察频谱幅值最大处对应的下标，换算频率应与设定一致。
下面示例演示如何填数据并调用上面写好的函数，最后打印前几个频谱点的模：
#include <iostream>
int main() {
const int N = 1024;
CVector sig(N);
double fs = 1024.0; // 采样率
for (int i = 0; i < N; ++i) {
double t = i / fs;
sig[i] = 0.5 * std::sin(2 * M_PI * 50 * t)
+ 0.3 * std::sin(2 * M_PI * 120 * t);
}
fft(sig);
// 输出前8点幅值
for (int i = 0; i < 8; ++i) {
std::cout << i << ": " << std::abs(sig[i]) << "n";
}
return 0;
}
六、常见错误与优化建议
初学者常把旋转因子符号弄反，导致频谱镜像或相位错乱。正向FFT用负指数，逆变换才取共轭。另外若输入长度不是2的幂，位反转和len翻倍会越界，调用前务必用零填充到最近二次幂。
在性能敏感处，可把genTwiddle结果缓存起来复用，或对蝴蝶循环做SIMD向量化。若只需实数信号，还能用实数FFT把内存和运算再省一半。掌握这些细节后，你写的C++ FFT就能稳妥地用在生产环境的信号处理模块里。