简介：这份资源以C语言实现FFT快速傅里叶变换，可用于电力系统、音频处理与通信领域的谐波分析，能够计算从基波到第51次谐波的含量，帮助评估非线性负载导致的波形失真。压缩包内共3个文件，包括C源码、配套头文件以及一份Word格式的FFT说明文档，整体仅62KB，轻量易读；源码体现DIT蝶形运算结构，文档则详细解释输入数据准备、参数设置、函数调用方式与结果解读，便于学习者对照源码理解每个环节。已有2827人学习下载。借助该工具，既能掌握FFT算法的C语言实现原理，也能直接用于谐波含量计算与信号频谱分析，同时可作为数字信号处理课程的实验参考，适合嵌入式开发者、电气工程师及高校学生阅读研究。 FFT谐波分析这活儿，玩嵌入式或者电力电子的人早晚都得碰。我前阵子做一台变频器输出端的电能质量评估，甲方上来就问“谐波到几次”“THD多少”，我手里的手持式电能分析仪又刚好不在现场，最后干脆用板子上的MCU直接跑FFT，把51次以内的谐波含量全部算了出来，实测结果和送检报告对得上。今天就把这套做法完整拆开讲一遍。
1. 需求分析与整体方案设计
1.1 为什么必须算到51次谐波
国标GB/T 14549-1993对公用电网谐波的规定是算到50次，IEC 61000-4-7标准也是建议测量到50次谐波。我在实际工程里习惯把上限做到51次，主要原因是想覆盖奇数次的统计便利性：50Hz工频下51次谐波频率是2550Hz，60Hz工频下是3060Hz，这个范围足够覆盖大多数变频器、UPS、开关电源产生的低频传导谐波。
更关键的是，现在的PWM整流器和变频器开关频率普遍在2kHz到20kHz之间，开关频率附近的边带谐波往往落在第40次到第60次频段。如果只算到50次，正好把这些边带漏掉一半，THD（总谐波畸变率）的评估就会有明显偏差。51次不是一个理论推导出来的数，是工程实践中摸出来的保险值。
对于电力电子工程师来说，谐波分析的目标通常就三类：
- 评估并网电流/电压是否符合标准——需要各次谐波的含有率，核心是幅值精度
- 排查设备之间干扰——关心谐波频率分布，峰值位置找对了就行
- 计算THD和功率因数校正——需要从基波到高次谐波的全谱积分
这三种需求决定了FFT实现时的取舍方向：是重精度还是重速度，是算幅值谱还是算功率谱。我在这次项目里三种都要覆盖，所以直接用浮点FFT，避免定点实现带来的动态范围问题。
1.2 FFT方案选型：库函数、IP核还是手写
做FFT谐波分析，代码不是难点，真正的难点在于选对实现途径。目前主流的做法有三条路：
第一，MCU专用DSP库。比如ARM的CMSIS-DSP库，或者TI的DSP库，直接调arm_cfft_f32这类接口，底层已经是优化过的汇编或高度优化的C代码。这是嵌入式项目里性价比最高的方案，我是用了STM32F4系列，主频168MHz，1024点浮点FFT单次跑下来不到1ms，完全满足实时显示的需求。
第二，FPGA的FFT IP核。Xilinx Vivado里的FFT IP核，如果采样点数固定，计算延迟可以压到微秒级，通过AXI-Stream接口输入数据、输出频谱。这个方案适合需要极低延迟或超大点数（8192点以上）的场景，逻辑开发周期比MCU方案长很多，但性能上限完全不同。
第三，纯手写FFT。教学用可以，量产项目不建议。理由是：手写FFT的bug排查成本极高，而且性能很难超过芯片厂商或者MathWorks等专业机构多年的优化成果。我在大学的课程设计里写过一次基2时间抽取FFT，跑通没问题，放到产品里对比CMSIS库函数，同样的1024点，性能差了3倍以上，精度也有差距。
如果要在上位机或者Matlab里做离线分析，那就更简单了，fft()一行的活，但要注意Matlab的FFT输出规则和嵌入式的完全一致，核心是处理好归一化和频谱搬移。
2. 关键参数计算与原理解析
2.1 采样率与FFT点数的匹配计算
FFT谐波分析最容易被忽视、又最决定成败的环节，是采样参数的确定。这里有个标准的计算链条：
第一步：确定最高分析频率。基波50Hz，算到51次谐波，最高分析频率fmax = 51 × 50 = 2550Hz。
第二步：确定采样率。按奈奎斯特采样定理，采样率必须大于2倍最高频率，即fs > 5100Hz。实际工程中我习惯留20%到30%的裕量，同时兼顾ADC的时钟分频便利性，所以取了fs = 12800Hz。这个数值是128的整数倍，对后续整周期采样特别友好。
第三步：确定FFT点数。这里有一个决定频谱分辨率的公式：Δf = fs / N，其中N是FFT点数。频谱分辨率决定了你能不能分清相邻的两个频率分量。对于50Hz基波系统，理想的Δf正好等于50Hz（每个谐波恰好占一根谱线）或者50的整数分之一，这样各次谐波才能落到整数的谱线位置上，避免栅栏效应。
我做过的选型对比：
| FFT点数N | 采样率fs | 频率分辨率Δf | 51次谐波所在谱线 | 特点 |
|---|---|---|---|---|
| 128 | 6400Hz | 50Hz | 第51根 | 最省资源，每根谱线恰为一次谐波 |
| 256 | 12800Hz | 50Hz | 第102根 | 多4倍频谱细节，分辨率仍对准50Hz |
| 1024 | 12800Hz | 12.5Hz | 第204根 | 频谱细节丰富，能分辨基频附近的间谐波 |
我最终选了1024点、12800Hz采样率。原因很简单：工业现场的电网上不只有50Hz基波，还有大量的间谐波和噪声。如果分辨率是50Hz，那基波附近的49Hz、51Hz分量就会被直接抹平，看不出来。用12.5Hz分辨率，至少能看到这些异常分量的轮廓。代价是计算量增大了8倍，但前面说过，CMSIS-DSP跑1024点FFT不到1ms，这个代价完全可以接受。
2.2 频谱泄漏与窗函数选择
FFT的本质是假设被分析的信号是一个无限周期信号的一个周期切片。如果采样窗口内不是整数个基波周期，频谱就会向四周“泄漏”，原本应该在某一根谱线上的能量，散到了旁边很多根谱线上。
解决泄漏问题有两条路：一是做到整周期采样，二是加窗函数。
整周期采样的做法是让采样持续时间正好等于基波周期的整数倍。我上面的参数就是为此设计的：12800Hz采样率，1024点数据，采样持续时间为1024/12800 = 0.08秒，正好是50Hz周期的4倍。这种情况下，不加窗也不会泄漏。
但现实是电网频率会波动，49.8Hz、50.2Hz都是合法的。频率一旦偏离50Hz，整周期就破坏了，还是会泄漏。所以工程上必须加窗。窗函数的选型有个经典对比：
| 窗函数 | 主瓣宽度 | 旁瓣衰减 | 频率分辨率 | 适用场景 |
|---|---|---|---|---|
| 矩形窗 | 窄 | -13dB | 最高 | 整周期采样严格成立时 |
| 汉宁窗 | 中等 | -31dB | 中等 | 通用电力谐波分析，我最常用 |
| 布莱克曼窗 | 宽 | -58dB | 较低 | 强干扰下的弱信号检测 |
我测试过：在电网频率偏移0.5Hz的情况下，矩形窗测出的基波幅值误差可达2%以上，汉宁窗可以控制在0.5%以内。对谐波分析来说，这个差距已经足够决定THD是否超标了。所以我默认给所有数据加汉宁窗，除非用户明确要求矩形窗做频谱细看。
窗函数不是免费的：它会让主瓣变宽，导致相邻很近的分量区分度下降。好在谐波信号彼此间隔50Hz，汉宁窗的3dB主瓣宽度约为2个谱线间隔（在这个例子里是25Hz），不足以吞并相邻谐波。
3. 嵌入式FFT核心代码与谐波提取实现
3.1 基于CMSIS-DSP的FFT完整流程
我用的STM32F4，采样由ADC+DMA自动完成，ADC采样率和定时器触发精确到12800Hz。外设配置这里不展开，重点讲数据从DMA缓冲区到谐波结果输出的完整链路。
第一步，把ADC的原始数据填入FFT输入缓冲。CMSIS-DSP的arm_cfft_f32要求数据按实部、虚部交替排列，所以我们把采样值放在偶数位，奇数位填0：
#define FFT_SIZE 1024 #define SAMPLING_FS 12800.0f #define FUND_FREQ 50.0f #define HARM_MAX 51 float32_t fft_input[FFT_SIZE * 2]; float32_t fft_mag[FFT_SIZE / 2]; uint16_t adc_buffer[FFT_SIZE]; // DMA自动填充 for (uint16_t i = 0; i < FFT_SIZE; i++) { // 先做去直流：减去一个周期的平均值，避免直流分量淹没低频细节 fft_input[2 * i] = (float32_t)adc_buffer[i] - dc_offset; fft_input[2 * i + 1] = 0.0f; }这里有一个细节值得注意：dc_offset必须先算出来，而不是等于ADC量程中点。工业现场的信号调理电路可能引入直流偏置，这个偏置会在频谱的0Hz处形成一个大峰值，严重时会影响第1、2根谱线的读数。我是在进入FFT之前维护一个滑动平均来估计直流分量，然后做差。
第二步，加窗。以汉宁窗为例，窗系数需要提前算好存成表格，避免运行时重复计算正弦函数：
static float32_t hanning_win[FFT_SIZE]; void window_init(void) { for (uint16_t i = 0; i < FFT_SIZE; i++) { hanning_win[i] = 0.5f * (1.0f - arm_cos_f32(2.0f * PI * i / (FFT_SIZE - 1))); } } void apply_window(void) { for (uint16_t i = 0; i < FFT_SIZE; i++) { fft_input[2 * i] *= hanning_win[i]; // 虚部为0，不用管 } }第三步，调用FFT并计算幅值。CMSIS-DSP的FFT函数是原地操作，输入和输出共用同一个缓冲区：
arm_cfft_f32(&arm_cfft_sR_f32_len1024, fft_input, 0, 1); arm_cmplx_mag_f32(fft_input, fft_mag, FFT_SIZE / 2);arm_cfft_f32最后一个参数1表示做正变换，如果是逆变换则填0。arm_cmplx_mag_f32是求复数模长，也就是对每个谱线计算sqrt(re^2 + im^2)，但是注意它的输出长度是FFT_SIZE/2，因为我们只看0到奈奎斯特频率的正频率部分，负频率部分在实数输入的条件下是共轭对称的。
3.2 51次谐波分量的提取与THD计算
频谱算出来后，谐波提取反而是个细致活。由于我设置的采样参数，每个谐波峰值的谱线位置是可以直接计算出来的：
float harm_amp[HARM_MAX + 1]; // harm_amp[0]为直流，实际不使用 for (uint16_t k = 1; k <= HARM_MAX; k++) { // 换算谐波频率在频谱中的索引 float freq_k = k * FUND_FREQ; uint16_t idx = (uint16_t)(freq_k * FFT_SIZE / SAMPLING_FS + 0.5f); // 在该索引附近找局部最大值 uint16_t peak_idx = idx; float max_val = fft_mag[idx]; for (int16_t offset = -2; offset <= 2; offset++) { uint16_t ni = idx + offset; if (ni >= 1 && ni < FFT_SIZE / 2 && fft_mag[ni] > max_val) { max_val = fft_mag[ni]; peak_idx = ni; } } harm_amp[k] = max_val * 2.0f / FFT_SIZE; // 幅值还原 }这里面有两个容易出错的地方。第一个是2.0f / FFT_SIZE这个归一化系数。FFT输出的幅值与真实幅值之间不是直接的等号关系：如果不归一化，一个1000点的正弦信号经过1024点FFT后，峰值谱线的幅值大约是500，等于信号幅值的N/2倍。所以要除以N/2，也就是乘上2/N。注意这是针对单频正弦信号的峰值幅值，不是RMS值。如果要算有效值，再除以1.414。
第二个是你搜索峰值时不能只取idx这一根谱线。工程电网频率有波动，谐波峰不会严格落在理论谱线上，可能偏了1到2根。我上面代码里在idx ± 2范围内找局部最大，这个方法简单粗暴但非常有效。如果你不想牺牲精度，可以做频域插值，我放到第四章讲。
最后是THD计算。按IEC的定义，THD是各次谐波有效值与基波有效值的比值（这里用幅值代替有效值，因为同一信号下的比例关系不变）：
float thd = 0.0f; float sum_harm_sq = 0.0f; for (uint16_t k = 2; k <= HARM_MAX; k++) { sum_harm_sq += harm_amp[k] * harm_amp[k]; } thd = sqrtf(sum_harm_sq) / harm_amp[1] * 100.0f; // 单位为%算到51次这一点特别重要：如果我按传统的只算到25次，THD可能显示为8%，但算到51次后可能变成11%。这个2到3个百分点的差距，会直接影响设备是“合格”还是“超标”的判定。
4. 实测常见问题与排查经验
4.1 频率偏移下的栅栏效应和插值校正
嵌入式设备测电网，最头疼的就是电网频率波动。我调试时遇到过这么个情况：基波实际是49.8Hz，采样参数是按50Hz设计的，结果频谱图上基波峰值不在第128根谱线（理论位置）上，而是落到128和129之间的某个位置。这时直接取第128根谱线的幅值，会明显偏低——因为能量分散到了相邻谱线上。这就是栅栏效应，只通过离散谱线“看”信号，总有看不见的角落。
解决办法是双谱线插值，原理很简单：取峰值谱线两侧两根谱线的幅值，按权重估算真实峰的位置和幅度。经典的汉宁窗双谱线插值校正公式如下：
设峰值谱线索引为k0，相邻两根谱线的幅值分别为y1和y2（其中y1 >= y2），定义参数β = y2 / y1，那么频率偏移量δ可以通过查表或近似公式求得。汉宁窗下，δ与β之间有近似关系：
δ ≈ (2β - 1) / (β + 1)
真实幅值A = y1 × (2 + δ) / 1.5（汉宁窗的幅值恢复系数略有不同），真实频率f = (k0 + δ) × fs / N。
我实测下来，用这个校正后，在49.5Hz到50.5Hz的偏移范围内，幅值误差从超过2%降到0.3%以内。如果你的系统对精度要求较高，这步不能省。
4.2 我踩过的实战避坑清单
谐波分析这个功能，看着简单，实际跑起来问题一堆。我把踩过的坑整理成了表格，按出现频率排序：
| 现象 | 根因 | 解决办法 |
|---|---|---|
| 频谱第0根谱线异常高大，后面低频段模糊 | 输入信号含直流偏置 | 先做去直流处理，减去滑动平均 |
| 各次谐波峰值偏小，基波误差尤其明显 | FFT幅值归一化系数不对 | 检查系数是否为2/N，窗函数补偿系数是否引入 |
| 频谱出现宽大“裙边”，谐波峰连成一片 | 采样非整周期且未加窗 | 加汉宁窗，或在锁相环控制下进行同步采样 |
| 偶数次谐波也很大 | 波形正负半周不对称 | 检查变送器、运放供电和偏置，常见原因是信号调理电路单电源供电 |
| 51次以后仍有很多大分量 | 采样率不足或信号本身含有高频振荡 | 确认采样率大于2倍的目标频率，必要时提高采样率分析更高频段 |
| MCU资源占用过高，频谱刷新率不足 | 使用了大点数FFT场景却没用浮点加速 | 使用带FPU的MCU，打开编译器的FPU优化选项，或者改用定点FFT |
| THD计算结果与商用测试仪对不上 | 只算到某次就截断，或未计入间谐波 | 统一到51次，确认旁瓣泄漏和底噪已经压低 |
还有一个容易被忽略的点：当作FFT的数据长度不是2的幂时，或者你在做循环缓冲的拼接时，数据不连续会导致频谱出现剧烈的“刺状”噪声。解决方法是确保送入FFT的1024个点，在时间上是严格连续的采样序列，不能有缺帧或重复。
用1024点、12800Hz采样做出来的频谱，如果你用它直接计算功率谱密度，可以用每根谱线的幅值平方除以频率分辨率（在功率谱密度定义为每Hz功率时）。不过对于谐波分析而言，我们通常更关心幅值和功率谱，而不是功率谱密度，因为标准要求的谐波指标以“百分比”为单位。只有做噪声分析时才需要用到PSD的积分概念。
最后再说一个实际调测的技巧：FFT的谐波结果虽然能准确反映电网质量，但网格线上的首个谱峰往往不是基波，而是直流附近的低频干扰。我建议在看频谱前先用一个20Hz的高通滤波器预处理数据，比如用简单的差分滤波器，可以极大改善低频段的信噪比。这个改动成本极低，但效果立竿见影。
上次项目交付后，甲方又拿来一台进口变频器的输出电流做同样分析，我只改了下采样率和基波频率参数，程序半小时内完成适配。这就是把FFT谐波分析做成通用模块的好处——参数可配置，算法不需要跟着硬件换。
本文还有配套的精品资源，点击获取
