简介这份资源以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次谐波频率是2550Hz60Hz工频下是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系列主频168MHz1024点浮点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频率分辨率Δf51次谐波所在谱线特点1286400Hz50Hz第51根最省资源每根谱线恰为一次谐波25612800Hz50Hz第102根多4倍频谱细节分辨率仍对准50Hz102412800Hz12.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采样由ADCDMA自动完成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优化选项或者改用定点FFTTHD计算结果与商用测试仪对不上只算到某次就截断或未计入间谐波统一到51次确认旁瓣泄漏和底噪已经压低还有一个容易被忽略的点当作FFT的数据长度不是2的幂时或者你在做循环缓冲的拼接时数据不连续会导致频谱出现剧烈的“刺状”噪声。解决方法是确保送入FFT的1024个点在时间上是严格连续的采样序列不能有缺帧或重复。用1024点、12800Hz采样做出来的频谱如果你用它直接计算功率谱密度可以用每根谱线的幅值平方除以频率分辨率在功率谱密度定义为每Hz功率时。不过对于谐波分析而言我们通常更关心幅值和功率谱而不是功率谱密度因为标准要求的谐波指标以“百分比”为单位。只有做噪声分析时才需要用到PSD的积分概念。最后再说一个实际调测的技巧FFT的谐波结果虽然能准确反映电网质量但网格线上的首个谱峰往往不是基波而是直流附近的低频干扰。我建议在看频谱前先用一个20Hz的高通滤波器预处理数据比如用简单的差分滤波器可以极大改善低频段的信噪比。这个改动成本极低但效果立竿见影。上次项目交付后甲方又拿来一台进口变频器的输出电流做同样分析我只改了下采样率和基波频率参数程序半小时内完成适配。这就是把FFT谐波分析做成通用模块的好处——参数可配置算法不需要跟着硬件换。本文还有配套的精品资源点击获取