FFT谐波分析：51次谐波为何算不准与C实现避坑

简介：面向电力测量与数字信号处理初学者的 FFT 谐波分析 C 语言实现包，围绕基波到第 51 次谐波含量的计算展开，可用于评估非线性负载引入的波形失真与谐波污染程度，也适合作为教学与研究中的算法实例。压缩包共 3 个文件，约 62KB，含 1 个 doc 说明文档、1 个 c 源文件与 1 个 h 头文件：头文件声明接口，源文件给出基于 DIT 结构的蝶形运算实现，文档则交代输入数据的准备方式、参数设置、函数调用方法与结果解读。目前已有 2833 人学习下载。读者可从代码中看清倒序重排、旋转因子与逐级蝶形组合的完整流程，理解如何把长序列的离散傅里叶变换拆成小规模运算逐层合并，并借助文档把算法接入自己的采样数据，完成各次谐波幅值的统计与结果校验，省去从公式推导到可编译代码之间的摸索成本。

1. 算不准 51 次谐波，九成问题不在蝶形运算

很多人第一次做 FFT 谐波分析，代码是从别处拷来的、蝶形运算背得滚瓜烂熟，接上 ADC 一跑，基波对得上，可 51 次谐波的读数在 0 和某个随机值之间跳。第一反应是怀疑 FFT 写错了，于是反复检查旋转因子的正负号、位反转顺序、归一化系数，折腾一周还是那个数。

真正的原因往往藏在采样这一侧。51 次谐波在基波 50 Hz 时对应 2550 Hz，这个频率离折叠点有多远、谱线能不能正好落上去、截断处相位连不连续，才是决定这条谱线是锐利还是糊成一团的关键。蝶形运算只要写对一次就永远是对的，采样参数配错却是每换一个工况就换一种错法。

搞懂这套东西的价值在于：你以后看到任何一条谐波读数异常，能立刻判断是算法问题还是参数问题，而不是把整条链路重写一遍。下面从蝶形运算的原理讲起，一路讲到 51 次谐波那条谱线是怎么被精确定位出来的。

2. 从 DFT 到蝶形运算：O(N²) 到 O(N log N) 省下的到底是什么

2.1 旋转因子的对称性：把重复乘法折叠掉

N 点 DFT 的定义是 X[k] = Σ x[n]·W_N^{nk}，其中 W_N = e^{-j2π/N} 叫旋转因子。直接按定义算，每个 k 都要遍历 N 个 n，一共 N² 次复数乘法。N 取 1024 时是 104 万次复数乘法，在带 FPU 的核上约几十毫秒，做实时谐波分析根本来不及。

FFT 的突破口是 W_N 的两个性质。周期性：W_N^{k+N} = W_N^k，旋转因子只有 N 个不同的值。对称性：W_N^{k+N/2} = -W_N^k，隔半个圆周的两点只差一个负号，乘法可以省掉。

基 2 分解就是把下标按奇偶拆开：X[k] = E[k] + W_N^k·O[k]，其中 E[k] 是偶数序号点的 N/2 点 DFT，O[k] 是奇数序号点的 N/2 点 DFT，k 取 0 到 N/2-1。同时 X[k+N/2] = E[k] - W_N^k·O[k]。两个 N/2 点 DFT 的中间结果被复用两次，一次加法一次减法，这就是一个蝶形。

递归下去，N 点 DFT 拆成两个 N/2 点，再拆成四个 N/4 点，一直到 2 点。总层数是 log2N，每层有 N/2 个蝶形，每个蝶形 1 次复数乘法、2 次复数加法。复数乘法降到 (N/2)·log2N。N=1024 时是 5120 次，比 104 万次少了约 205 倍。这个数量级差距就是 FFT 能进嵌入式实时系统的根本原因。

注意：所谓「省下的是乘法」不准确。加法次数从 N² 降到了 N·log2N，乘法降得更狠，但真正决定定点实现能不能跑的是加法的溢出路径，不是乘法器数量。

2.2 基 2 蝶形的三种写法与原地计算的地址规律

同一个蝶形有三种常见组织形式，选哪种直接决定内存占用和缓存行为。

第一种是递归版，函数自己调自己，代码最短，但每层都要分配临时数组，栈开销大，嵌入式里基本不用。第二种是迭代版、输入自然序输出倒位序，读写地址规律差，需要额外索引表。第三种是迭代版、输入倒位序输出自然序，全部原地完成，只需要一个 N 点复数数组，是最常用的形态。

原地计算能成立，靠的是蝶形的一个地址性质：第 stage 层（stage 从 1 开始）的跨度 m = 2^stage，同一层内参与蝶形的两个点下标相差 half = m/2，且这两个点只在本层内被读写一次，之后不再被同层其他蝶形触碰。所以可以在原数组上直接覆盖，不需要影子缓冲。

旋转因子的复用也是在这一层榨出来的。第 stage 层里，跨度为 m 的两点之间的旋转因子是 W_m^k，而这张表在整张 N 点旋转因子表里的步长是 N/m。也就是说全篇只需要一张 N/2 项的旋转因子表，各层按不同步长去取，不必为每层单独生成。这一条直接把表的内存从 N·log2N 项降到 N/2 项。

2.3 位反转：输入序列为什么要先倒着排

奇偶拆分是递归进行的，第一层按最低位分，第二层按次低位分，最后一层按最高位分。如果输入按自然序存，每一层的访问都会跳跃，无法原地。把输入按照下标的二进制位反转重排之后，每一层的奇偶分组就变成了连续的两段，地址连续、可以原地。

位反转的朴素做法是开一个 N 项索引表，查表置换。省内存的做法是用一个游标 j 增量维护，逐位进位，复杂度 O(N)，不需要额外数组。

/* 位反转置换：把自然序输入重排成倒位序，为原地蝶形做准备 */ static void bit_reverse(complex_t *x, int n) { int i, j, bit; for (i = 1, j = 0; i < n; i++) { bit = n >> 1; /* j 是 i 的位反转累加器：从最高位开始逐位尝试置 1 */ for (; j & bit; bit >>= 1) j ^= bit; /* 该位已为 1，进位复位 */ j ^= bit; /* 把第一个为 0 的位翻成 1 */ if (i < j) { /* 只在 i<j 时交换，避免换回去 */ complex_t t = x[i]; x[i] = x[j]; x[j] = t; } } }

bit从 n/2 开始，对应下标的最高位；内层循环负责处理进位，遇到已经是 1 的位就清零并右移一位，遇到 0 就置 1 并停下。i < j这个判断是必需的，否则每一对会被交换两次，等于没换。

蝶形主体如下，输入已经是位反转序，输出直接是自然序：

/* 原地基 2 DIT-FFT：n 必须是 2 的幂，输入已按位反转排好 */ static void fft_radix2(complex_t *x, const complex_t *twiddle, int n) { int stage, m, half, k, i; for (stage = 1, m = 2; m <= n; stage++, m <<= 1) { half = m >> 1; for (k = 0; k < half; k++) { /* 本级第 k 个旋转因子在整张表里的步长是 n/m */ complex_t w = twiddle[k * (n / m)]; for (i = k; i < n; i += m) { complex_t t = cmul(w, x[i + half]); x[i + half] = csub(x[i], t); /* 蝶形下半支：E - W·O */ x[i] = cadd(x[i], t); /* 蝶形上半支：E + W·O */ } } } }

三层循环的分工：外层stage控制本级的蝶形跨度m；中层k遍历本级的 half 个旋转因子；内层i以m为步长扫过整个数组，做原地读写。w在中层循环里只算一次，内层全程复用，这是把复数乘法从 N·log2N 压到 (N/2)·log2N 的关键动作。

3. 51 次谐波的分辨率账本：采样率、点数与栅栏效应

3.1 Δf = fs/N 决定了谱线落在哪里

FFT 输出的是离散谱线，第 k 根谱线对应的频率是 k·fs/N，相邻谱线的间隔 Δf = fs/N 就是频率分辨率。这个公式是整套谐波分析的地基，51 次谐波能不能被准确读出，先看 2550 Hz 这个频率能不能落在某根谱线上。

以 fs = 6400 Hz、N = 1024 为例，Δf = 6.25 Hz。50 Hz 落在第 8 根谱线（50/6.25 = 8），51 次谐波 2550 Hz 落在第 408 根谱线（2550/6.25 = 408），都是整数。因为 6.25 恰好整除 50，所有整数次谐波都正好落在谱线上，不存在栅栏效应。

再看折叠频率。fs/2 = 3200 Hz，对应第 512 根谱线，也就是第 64 次谐波。理论上能分辨到 63 次，51 次留了 13 次的余量。

换一组参数就露馅了。fs = 10000 Hz、N = 1024，Δf = 9.7656 Hz，50 Hz 落在第 5.12 根谱线上，能量被劈到第 5 和第 6 两根谱线，各分走一部分，幅值读数直接掉三成。这就是栅栏效应：谱线之间的信号既不在谱线上、又不完全在谱线外，读数完全取决于泄漏的分布。

参数选取的顺序应该是反的：先定最高要分析的谐波次数，得到所需最高频率；再让 fs 至少是这个频率的 2 倍并留出滤波过渡带；最后让 Δf 整除基波频率，解出 N。51 次这个需求，2550 Hz 的两倍是 5100 Hz，取 6400 Hz 留了 1300 Hz 余量，N = 6400/6.25 = 1024，正好是 2 的幂次，可以直接用基 2 FFT。这一串数字不是巧合，是反推出来的。

3.2 电网频率不是 50.00 Hz：整周期采样的现实裂缝

上面的推导有个隐含前提——基波正好是 50.00 Hz。实际电网频率在 49.5 到 50.5 Hz 之间缓慢漂移，偏差 1% 就意味着 1024 点采样窗内不是整数个工频周期，截断处相位不连续，频谱泄漏立刻出现。

泄漏的表现很有欺骗性：基波幅值基本不变，但 51 次谐波这种高次谐波本来幅值就小，被泄漏过来的基波旁瓣一盖，读数能从真实值浮上来好几倍。

工程上有两条路。第一条是硬件同步采样：用锁相环跟踪基波过零点，把 ADC 的采样时钟按倍频锁定到实际工频上，采样窗永远等于整数个周期。代价是需要额外的过零检测电路和 PLL，成本上去了。第二条是软件重采样：先粗测基波相位，用三次样条或线性插值把原始采样点重映射到理想时间栅格上，再做 FFT。代价是插值误差和额外算力，好处是纯软件、成本几乎为零。

对 51 次谐波这种高次分量，插值误差会被放大——高次谐波对应的高频成分对时间抖动最敏感。所以软件重采样的插值阶数不能太低，线性插值在高次谐波上引入的幅值误差可能到几个百分点。

3.3 窗函数与幅值修正：汉宁窗的 0.5 从哪来

当整周期采样做不到时，加窗是最后一道防线。矩形窗等于不加窗，主瓣最窄（2Δf），但旁瓣只衰减 13 dB，远处的高次谐波很容易被基波的旁瓣淹没。汉宁窗主瓣展宽到 4Δf，旁瓣压到 -31 dB，代价是频率分辨率下降。

加窗之后幅值必须做补偿，否则读数会偏小。补偿系数叫相干增益，定义为窗函数在时域的平均值：矩形窗是 1，汉宁窗是 0.5。原因很直白：汉宁窗把信号两端压到了零，等效于把信号整体乘了一个平均值为 0.5 的包络，谱线幅值自然掉一半。

所以单边幅度谱的正确算法是：

公式里的系数 2 是单边谱的折算（负频率那一半的能量并到正频率），除 N 是 DFT 的归一化。汉宁窗那格写成 4/N，就是 2/(N×0.5)。

注意：加窗会改变主瓣宽度，如果窗的主瓣宽到把相邻两次谐波盖在一起（比如高次谐波间隔只有 50 Hz，而窗主瓣展宽到 4Δf = 25 Hz 就开始重叠），读数会互相污染。51 次谐波距离 50 次只有 50 Hz，这个间距在 6.25 Hz 的分辨率下是 8 根谱线，安全，但若把分辨率降到 25 Hz，就必须重新评估。

4. C 实现：旋转因子表、定点取舍与谐波提取

4.1 旋转因子表的预计算与精度衰减

旋转因子表是全篇唯一需要调用三角函数的操。N = 1024 时表长 512 项，每项两个 float，共 4 KB。表在启动时算一次，之后 FFT 全程查表。

/* 预计算 N/2 项旋转因子表，w[i] = exp(-j*2*pi*i/N) */ static void build_twiddle(complex_t *w, int n) { const float pi = 3.14159265358979f; int i; for (i = 0; i < n / 2; i++) { float angle = -2.0f * pi * (float)i / (float)n; w[i].r = cosf(angle); w[i].i = sinf(angle); } }

角度前面的负号不能丢，丢掉之后算出来的等效于逆变换，相位全反。用 float 还是 double 取决于目标平台：float 有 24 位有效位，log2(1024) = 10 级运算累积的舍入误差大概在 10⁻⁵ 量级，对谐波分析要求的分辨率（通常 -60 dB，即 10⁻³）绰绰有余。double 只是徒增内存和算力，除非平台本身没有单精度硬件。

定点实现更麻烦。Q15 把满量程映射到 ±32767，蝶形里的加法 x[i] + t 随时可能超过 32767，溢出后符号位翻掉，频谱上直接出现鬼影。常见的解法是块浮点：每一级算完之后统计本级的最大绝对值，整体右移一位，用一个计数器记录总共移了几位，最后统一左移回去。这样动态范围只损失几位，比固定 Q15 安全得多，代价是每级多一次求最大值的遍历。

4.2 实数输入的 N/2 点技巧

谐波分析的输入永远是实数，虚部为 0。直接做 N 点复数 FFT 等于白扔一半算力。标准做法是把 N 个实数打包成 N/2 个复数（偶数点作实部、奇数点作虚部），做 N/2 点复数 FFT，再用后处理公式把谱拆开。

后处理的核心是两个对称关系：Z[k] 和 Z[N/2-k] 的共轭组合出正频率部分的谱。这一步省下来的是实打实的一半乘法和一半内存，N = 1024 时等效降到 512 点复数 FFT，在 100 MHz 级别的 MCU 上从毫秒级掉到亚毫秒级。

代价是代码复杂度上升，而且后处理公式写错的话症状特别隐蔽——谱线位置对、幅值差一点、相位全乱，非常难查。如果算力富裕、开发周期紧，老老实实跑 N 点复数 FFT 是可以接受的取舍。

4.3 从复数谱到 51 次谐波含量

FFT 输出之后，谐波提取本身只是一次索引加一次归一化。基波落在第 base_bin = 50/Δf 根谱线，第 k 次谐波落在第 k·base_bin 根谱线——因为所有谐波频率都是基波的整数倍，谱线序号自然也是整数倍。

/* 从单边幅度谱中提取各次谐波幅值并计算 THD */ float harmonic_analysis(const complex_t *spec, int n, int base_bin, int max_order, float win_gain, float *hbuf) { int k; float h1, sum_sq = 0.0f; /* 单边幅值 = 2*|X[k]|/(N*相干增益)，直流分量不乘 2 */ h1 = 2.0f * hypotf(spec[base_bin].r, spec[base_bin].i) / ((float)n * win_gain); hbuf[0] = h1; for (k = 2; k <= max_order; k++) { int bin = base_bin * k; if (bin >= n / 2) break; /* 越过折叠频率直接截断 */ hbuf[k - 1] = 2.0f * hypotf(spec[bin].r, spec[bin].i) / ((float)n * win_gain); sum_sq += hbuf[k - 1] * hbuf[k - 1]; } /* THD 按基波归一：各次谐波有效值的平方和开根再除基波 */ return sqrtf(sum_sq) / h1; }

参数需要对齐：n是 FFT 点数，base_bin由基波频率除以分辨率得到，max_order取 51，win_gain矩形窗传 1.0、汉宁窗传 0.5。bin >= n/2这条截断必须写，否则当基波频率较高、采样率又不够时，k·base_bin 会越界读取数组后面的内存，读出来的值看着像谐波，实际是别的变量。

hbuf里存的是各次谐波的幅值，索引 0 是基波，索引 50 对应第 51 次。如果下游还要算各次谐波含量百分比，直接拿 hbuf[k]/h1 即可，不必再进一次 FFT。

5. 避坑：谐波分析里最容易翻车的五个地方

5.1 幅值与频率类问题

现象一：所有谐波幅值都只有真实值的 1/N。原因：忘了做 DFT 归一化，直接把 |X[k]| 当幅值用了。N 点 FFT 的 |X[k]| 是 N 倍量级，不除 N 读数全错。如果是单边谱还要再乘 2（直流和折叠频率点除外）。解决：把归一化系数抽成一个宏AMP_SCALE(n, win_gain)，全工程只此一处定义，别在每个调用点手写。

现象二：基波读数正常，51 次谐波忽大忽小，随工况跳变。原因：三种可能——采样率不够导致 2550 Hz 被混叠（检查 fs/2 是否大于 2550 且抗混叠滤波器在 2550 Hz 处有足够衰减）；非整周期截断导致泄漏（看基波谱线有没有能量扩散到相邻谱线）；加窗后没做相干增益补偿（汉宁窗读数会系统性偏小一半）。解决：先把输入改成理想双音信号（50 Hz 和 2550 Hz 各占一半幅值）做离线验证，理想输入下读数对不上，就一定是参数或归一化问题，与算法无关。

现象三：直流分量附近的谱线被抬高，基波幅值虚高。原因：原始信号里带直流偏置（ADC 零点漂移、运放失调），直流在谱上是一根很高的谱线，它的旁瓣会泄露到基波位置，基波幅值被抬升，THD 算出来反而偏小。解决：进 FFT 前先减掉均值，或者加一个截止频率远低于 50 Hz 的高通滤波器。减均值最省事，一行代码。

5.2 实现与数据类问题

现象四：频谱左右镜像，或者谱线位置整体错位。原因：位反转写错，或者蝶形里旋转因子的步长算错。位反转的典型错误是漏了i < j判断导致交换两次；步长错误的典型表现是低频成分跑到高频去。解决：用一个已知频率的单音信号（比如只有 50 Hz）做测试，正确结果应该只有第 8 根谱线上有值，其余全接近 0。如果出现两根对称的谱线，就是位反转或共轭方向错了。

现象五：定点实现跑一段时间后某个频段突然冒出大幅值噪声。原因：蝶形加法溢出。Q15 的蝶形 x[i] + t 可能超过 32767，溢出翻转后变成一个大负数，经过后面几级放大，在频谱上表现为随机频段的鬼影。解决：改用块浮点，每级统计最大值做整体移位；或者把输入预缩放到满量程的四分之一，用动态范围换安全裕度。后者的代价是高次谐波这种小信号的分辨率下降，如果 51 次谐波幅值只有基波的千分之一，预缩放后可能已经淹没在量化噪声里，只能上块浮点。

6. 设计取舍：51 这个数字背后的余量账

为什么是 51 次而不是 50 次，也不是 63 次？这个数字不是标准凑出来的，是采样链路上几个约束互相挤压之后的落脚点。

先看上限。fs = 6400 Hz 时折叠频率 3200 Hz 对应第 64 次谐波，理论上能分析的极限是 63 次。但折叠频率附近恰恰是抗混叠滤波器最难做的地方：滤波器从通带到阻带需要一个过渡带，如果要求 63 次（3150 Hz）处仍然平坦、64 次以上必须衰减到 -60 dB，过渡带只剩 50 Hz，意味着滤波器阶数要拉到很高，成本和延迟都受不了。把上限退到 51 次（2550 Hz），过渡带从 2550 到 3200 有 650 Hz 的宽度，一个中等阶数的模拟低通就能做到，这是最直接的工程余量。

再看分辨率。N = 1024 时窗长 0.16 秒，等于 8 个工频周期，Δf = 6.25 Hz。这个窗长在实时性上是有讲究的：电网谐波监测通常要求 200 ms 级别的刷新率，0.16 秒的窗加上计算时间刚好卡在这个量级。如果把 N 提到 2048 换取更高的频率分辨率，窗长翻倍到 0.32 秒，刷新率掉一半，对动态谐波（比如电弧炉、变频器负载）的跟踪就跟不上了。

然后是点数约束。Δf 必须整除 50 Hz，同时 N 必须是 2 的幂——这两个条件同时成立时 N 只能取 2 的幂且 fs = 50·N/2^m 的形式。fs = 6400、N = 1024 恰好满足，fs = 3200、N = 512 也满足但折叠频率只有 1600 Hz，51 次谐波 2550 Hz 直接超了。所以 51 次这个需求和 1024 点这个规模是被同一个采样率锁死的，动一个就得全盘重算。

最后是算法侧。1024 点基 2 FFT 是 5120 次复数乘法，定点在无 FPU 的核上大约几十万周期，100 MHz 主频下几毫秒，实时性没压力。如果为了 51 次去上混合基或者任意长度 FFT（比如为了凑 1280 点整周期采样），代码复杂度会翻好几倍，收益却只是省掉一个重采样环节，不划算。

我最早做这块的时候吃过一次亏：为了让 51 次谐波读数更「稳」，把 N 从 1024 加到 4096，结果刷新率掉到 5 Hz，现场一台设备的谐波是阶跃式变化的，FFT 窗口横跨了跳变点，读数反而更乱。那次之后我形成了一个习惯——每定一组采样参数，先拿理想双音信号跑一遍离线验证，确认谱线位置、幅值归一化、窗增益补偿三件事同时正确，再接到真实信号上。如果理想输入下读数都对不上，那一定是参数问题，别去改蝶形运算。

这套从蝶形到谱线到归一化的完整链路，配上 C 实现和参数表，基本覆盖了从 50 Hz 基波一路算到 51 次谐波的全部关键点。资源里把旋转因子表、窗函数、谐波提取这些环节都拆成了独立模块，改采样率或者改谐波次数只需要动一处配置，适合拿去对着源码逐段验证。希望帮到你。

本文还有配套的精品资源，点击获取