做信号处理的工程师十有八九都跟FFT打过交道。无论是分析振动数据、排查设备故障还是处理电力系统的谐波问题FFT几乎是通用语言。我最初接触FFT是在一次电机电流信号分析中用MATLAB写了一段几十行的程序把时域波形转换成频谱立刻看到了隐藏在高噪声背后的特征频率。那次经历让我意识到FFT不是教科书上的数学公式而是实实在在能帮你找到问题根源的工具。这篇博文要聊的正是一个基于MATLAB的FFT分析与滤波程序——它能对数据信号做频谱分析把波形中的谐波分量拆出来逐个看清同时利用数字滤波技术按需处理频段。不管你是刚入门信号处理的在校学生还是需要在现场快速分析数据的一线工程师这套思路和代码都能直接上手。我会把FFT分析、谐波识别、滤波实现这三块的核心原理和实操细节都掰开揉碎讲清楚附上可以直接复制的完整代码。1. 为什么信号处理离不开FFT与滤波的组合拳1.1 时域波形看不懂的痛先回忆一个常见场景你采集到一段加速度传感器的数据时域波形是杂乱无章的振荡振幅忽大忽小完全看不出规律。如果你只盯着时域图去看除了能看出来信号在变化其余什么都得不到。这是所有时域信号的通病——信息被混叠在一起特征频率被噪声淹没。FFT解决的核心问题是把信号从时域视角切换到频域视角。做个简单的类比你听一首交响乐时域信号是所有人同时演奏的混响而频谱分析就像是把录音拆分成各个乐器的分轨能单独看清小提琴在拉什么音、鼓在敲什么节奏。FFT就是干这个的。它能把一个时域波形分解成一系列正弦波的叠加每个频率成分都有对应的幅度和相位。1.2 FFT与滤波的分工关系FFT和滤波是一对黄金搭档但它们的定位完全不同。FFT解决的是看清有什么——把信号中的频率成分展示出来滤波解决的是按需取舍——把不想要的频率成分削弱或滤除。实际项目中这两个功能常常配合使用先用FFT分析信号的频谱结构确定哪些频率是有效信号、哪些是噪声干扰再根据分析结果设计滤波器把干扰成分滤掉滤完之后再做一次FFT验证滤波效果。我经常把这个流程称为先诊断、后治疗、再复查。举一个实际例子。在处理逆变器输出电流信号时基波是50Hz但还会有大量的5次、7次、11次谐波以及开关频率附近的噪声。通过FFT分析能快速确定各次谐波的幅值占比然后设计一个低通滤波器把开关噪声滤除或者设计陷波器针对特定次谐波做处理。整个过程环环相扣缺一不可。2. 频谱分析前的数据准备工作2.1 采样率、采样点数与频率分辨率的关系我见过太多人拿到数据就急着调用fft函数结果频谱图出来后一头雾水——频率轴不对、峰值位置不准、幅度看起来也不合理。这些问题十有八九是采样参数没算清楚就上手了。在按下运行按钮之前有几个参数必须弄清楚。第一个是采样率fs。采样定理告诉我们要无混叠地分析信号采样率必须大于最高分析频率的两倍。比如信号里可能有1000Hz的分量采样率至少要在2000Hz以上。实际工程中通常留出余量取最高频率的5到10倍。第二个是采样点数N。频率分辨率df fs / N也就是频谱上相邻两条谱线之间的间隔。这个值决定了你区分两个相近频率的能力——如果两个频率相差不到一个df它们在频谱上是分不开的会融成一个峰。比如采样率是1000Hz采样点数是1000个那么频率分辨率就是1Hz这意味着只能区分相差1Hz以上的频率成分。想要分辨得更细只能增加采样点数没有捷径。第三个是FFT变换的点数NFFT。MATLAB的fft函数允许你指定计算点数通常取等于采样点数N也可以做补零处理获得更密的频谱显示。但这里有个容易误解的地方——补零不会提升物理分辨率它只是让频谱曲线看起来更平滑峰值位置的估计精度并不会因为补零而提升。2.2 数据预处理去直流与去除异常点很多AD采集卡出来的信号都带直流偏置不做处理就做FFT分析会在0Hz处产生一个巨大的分量导致频谱图被压扁其他频率成分的细节全部看不清。处理办法很简单——减掉信号的均值。原始信号减去均值相当于把直流分量去掉。除直流偏置外数据中常混有野值点。有些是传感器瞬间超量程产生的削顶有些是通信干扰导致的数据跳变。这些异常点如果在FFT分析前不处理会表现为频谱中宽频段的背景抬升导致小信号特征被淹没。野值点的检测可以用滑动窗口内中值滤波的思路——超出窗口内中值3倍标准差的点标记为异常用相邻正常值的插值替换。预处理这一步我建议做成一个独立函数因为它在滤波、FFT分析、特征提取等后续所有环节中都会用到。一次写好后反复调用能省不少事。2.3 窗函数的选择与频谱泄露的抑制频谱泄露这个词几乎所有做过FFT的人都会遇到。它指的是某个频率的能量漏到了旁边的频率上去造成频谱上的旁瓣。产生的原因很简单FFT要求处理的信号是周期信号且截取长度内恰好包含整数个周期。实际数据不满足这个条件截断会导致信号在边界处不连续从数学上看就是给信号乘了一个矩形窗而矩形窗的频谱带有旁瓣于是泄漏就发生了。对付频谱泄露最直接的办法是加窗。加窗的本质是对截取范围内的信号做加权让边界处的幅度逐渐衰减到零从而减小截断带来的突变。常用的窗函数有汉宁窗、海明窗、布莱克曼窗等它们的主要区别在于主瓣宽度与旁瓣衰减的权衡。汉宁窗是我最常用的选择它在主瓣宽度和旁瓣抑制之间取得了一个不错的平衡。如果你需要精确测量幅度而且信号中含有比较强的频率成分可以考虑用平顶窗它能提高幅度测量的精度。如果分析的是瞬态信号使用矩形窗反而更合适。选窗函数没有绝对正确只有适不适合当前场景讲清楚这一点比记住任何一张对照表都重要。3. FFT频谱分析的核心代码与谐波识别3.1 基础FFT分析完整代码下面这段代码是我在多个项目中反复使用的模板它可以接受任意一维信号自动完成去直流、加窗、FFT计算、单边谱转换和峰值标注。为了让代码具备通用性这里把采样率、窗函数类型都作为参数传入。function [f, amp, phase] fft_analysis_utils(x, fs, winType) % FFT频谱分析工具函数 % 输入 % x - 输入时域信号一维数组 % fs - 采样率单位Hz % winType - 窗函数类型rect,hann,hamming,blackman % 输出 % f - 频率轴 % amp - 单边幅值谱 % phase - 相位谱弧度 % 去直流 x x(:) - mean(x); N length(x); % 加窗 switch winType case rect win ones(N,1); case hann win hann(N); case hamming win hamming(N); case blackman win blackman(N); otherwise win hann(N); end xw x .* win; % FFT计算 NFFT N; X fft(xw, NFFT); % 单边谱 f (0:NFFT/2-1) * fs / NFFT; amp abs(X(1:NFFT/2)) * 2 / N; phase angle(X(1:NFFT/2)); % 加窗幅度修正 % 加窗会使信号总能量衰减需要乘以修正系数恢复真实幅值 % 矩形窗修正系数为1汉宁窗为2海明窗约为1.853 switch winType case rect amp amp; case hann amp amp * 2; case hamming amp amp * 1.853; case blackman amp amp * 2.26; end end注意幅值修正这一步。很多人直接对原始信号做FFT得到幅值谱这没错但一旦加了窗如果不做幅度修正得到的峰值幅度会比真实值偏小。汉宁窗修正系数为2海明窗约为1.853布莱克曼窗约为2.26。这个修正系数很多人不知道很容易导致幅度测量偏差。验证这段代码的效果可以构造一个标准信号来测试。假设信号包含50Hz、幅值1和120Hz、幅值0.5的两个正弦分量加上白噪声设置采样率1000Hz采样点数2000。运行分析后50Hz处峰值应接近1考虑噪声干扰后有一定误差120Hz处应接近0.5。3.2 谐波分量的自动识别与参数提取FFT频谱图上能看到很多峰值但光用眼睛看是不够的尤其是当谐波数量多、噪声水平高时人工识别效率太低。更可靠的做法是写一个峰值搜索算法自动找出频谱中的显著分量并提取每个分量的频率、幅值和相位。峰值识别的基本思路是设置一个幅度阈值只保留幅度高于阈值的谱线再对超过阈值的区域做局部最大值搜索剔除旁瓣引起的虚假峰值。这里需要注意的一点是两个邻近的频率成分如果距离小于频率分辨率df它们会混叠成一个宽峰峰值搜索算法会把它们识别为一个分量。这种物理极限无法靠算法弥补只能通过增加采样点数来提高分辨率。提取频率时可以采用抛物线插值提高精度。FFT输出的离散谱线中真实峰值可能落在两条谱线之间单凭最大谱线的位置只能给出分辨率级别的估计。用峰值点及其相邻两点做抛物线拟合可以估计出更精确的峰值频率与幅值。这个技巧在谐波分析中非常实用尤其是测量电力系统中的间谐波和次谐波时精度提升明显。相位信息的提取是一个容易被忽略的细节。FFT结果中每个谱线都带有相位角但这里有一个隐含条件——如果信号加过窗相位会被窗函数的相位特性影响。矩形窗是零相位窗不存在这个问题但汉宁窗等会改变相位读数。对于只关心频率和幅值的场景可以不用在意这一点但如果需要分析各次谐波之间的相位关系建议用矩形窗或者做相位校正。3.3 用FFT结果计算谐波失真率THD谐波失真率THD是衡量信号质量的重要指标常用于电力系统、音频设备、逆变器输出质量的评估。THD的定义是全部谐波能量与基波能量之比的平方根用百分比表示。有了FFT分析结果THD的计算就很简单。假定基波频率为f0那么2倍、3倍……N倍频率处的幅值平方和除以基波幅值平方再开方乘100%就是THD。代码实现如下% 从频谱中提取THD function thd_value thd_from_fft(f, amp, f0) % 找到基波频率对应谱线索引 [~, idx0] min(abs(f - f0)); % 取基波幅值 amp_fund amp(idx0); % 谐波阶次上限 maxHarmonic floor(max(f) / f0); % 计算谐波能量 harmonic_power 0; for k 2:maxHarmonic [~, idxk] min(abs(f - k*f0)); % 搜索pk附近一定带宽内的最大幅值避免噪声干扰 half_bw 3; % 3根谱线范围 if idxk - half_bw 1 idxk half_bw length(amp) seg amp(idxk-half_bw:idxkhalf_bw); harmonic_power harmonic_power max(seg)^2; end end thd_value sqrt(harmonic_power) / amp_fund * 100; end这段代码有一个关键点搜索谐波幅值时不能用固定的一个谱线索引而应该在预估频率附近做一个带宽搜索。因为实际信号频率往往不是严格的整数倍关系且存在频率漂移固定的谱线索引很容易取错。在谐波频率附近3到5根谱线范围内搜索最大值能显著提升THD计算结果的稳定性。4. 滤波器的设计与实现完成频谱分析后接下来的任务是滤波。滤波的目的很明确把有用的频段保留把噪声和谐波干扰抑掉。下面从频域滤波和时域滤波两个角度展开。4.1 频域直接滤波最直观的滤波方式是在频域做乘法。思路是对信号做FFT把不需要的频段幅值置零再用逆FFT变回时域。这种方法的优点是概念清晰、实现简单但要注意它有几个严重的副作用。第一个副作用是吉布斯效应。对频谱做硬截断相当于在频域乘一个矩形窗对应到时域是信号与sinc函数的卷积这会在信号边界处产生振铃现象。信号长度越短振铃越明显。处理办法是不要用硬截止而是在频域使用过渡带让边界处的幅值平滑地从1变到0。第二个副作用是相位问题。如果直接在频域把负频率部分全部置零再逆变换回时域获得的是解析信号实部才是原始滤波结果。在MATLAB中处理这个环节要特别仔细正负频率要对称操作否则逆变换出来的信号会有虚部残留。频域直接滤波虽然看起来不如时域滤波器常用但在一些特定场景下非常好用。比如你只想去掉某个非常窄的频带如50Hz工频干扰用IIR陷波器设计就很麻烦而频域直接把这个频带清零干净利落。我的做法是把它封装成一个函数在需要做频谱定制时优先使用。4.2 FIR与IIR滤波器的设计如果需要在实际系统中实时应用滤波器还是应该使用标准的FIR或IIR滤波器。MATLAB中设计这两类滤波器都有现成工具但需要对两者特性有清楚认识。IIR滤波器的优势是阶数低、计算量小相同过渡带下可以用更少的阶数达到同样的衰减效果。但它有一个固有缺陷——非线性相位。这意味着不同频率成分通过滤波器后产生的延迟不一致信号波形会发生相位失真。如果后续处理需要保留波形特征比如峰值检测这可能会是个问题。FIR滤波器则恰恰相反它有严格的线性相位特性不会造成相位失真代价是阶数高很多。在相同设计指标下FIR滤波器的阶数往往是IIR滤波器的5到10倍。对于离线分析来说这个代价完全可以接受对于实时系统要评估一下DSP或MCU的计算能力是否够用。我个人的倾向是离线分析优先用FIR实时系统频段间隔充裕时用IIR追求相位精度时用FIR。MATLAB中设计一个巴特沃斯低通滤波器的基本流程如下% IIR滤波器设计 fc 100; % 截止频率100Hz fs 1000; % 采样率1000Hz order 4; % 阶数 [b, a] butter(order, fc/(fs/2), low); % 零相位滤波 filtered_signal filtfilt(b, a, raw_signal);注意我用了filtfilt而不是filter。filtfilt是零相位滤波也就是对信号正向滤波一遍再反向滤波一遍两次滤波的相位失真相互抵消。对于离线数据分析这个函数几乎是必备的它能避免滤波导致的时间偏移。如果是实时系统无法使用filtfilt只能使用filter此时滤波延时需要计入整体系统延时预算。设计截止频率时还有一个容易踩的坑截止频率f与MATLAB中的归一化频率不是一回事。butter函数的归一化频率是fc除以奈奎斯特频率fs/2。比如采样率1000Hz、截止频率100Hz归一化频率是100/5000.2。弄错这个单位关系会导致滤波器实际截止频率完全不符合预期。4.3 中值滤波、滑动平均与S-G滤波除了经典频率响应型滤波器还有一类时域平滑滤波方法在信号处理中也极其常用。这些方法的核心思想不是基于频率而是基于数值统计它们对脉冲噪声和高频毛刺有很好的抑制效果。中值滤波的作用是把某个采样点替换为窗口内所有点的中值。这个操作对脉冲噪声非常有效几个孤立的异常尖峰经中值滤波后会被剔除。窗口宽度一般取奇数常见的取值是3、5、7。窗口越大平滑效果越强但响应越迟钝细节保留越少。在处理加速度信号或压力信号时我会先用中值滤波去除传感器尖峰再做FFT分析这样频谱图会更干净。滑动平均滤波的数学形式是窗口内所有数据的算术平均值它本质上是FIR滤波器的一个特例。实现很简单MATLAB中有movmean函数一行搞定。滑动平均对高斯白噪声的抑制效果不错但对脉冲噪声和野值不敏感而且会有比较明显的波形钝化效应。窗口宽度建议通过实际试验确定从5开始逐步增大直到噪声水平可接受为止。S-G滤波Savitzky-Golay滤波是一种更高级的平滑方法它在滑动窗口内做多项式最小二乘拟合用拟合值替换窗口中心点的值。S-G滤波器比滑动平均的亮点在于它能更好地保留信号的峰值和谷值不会把波形整体削平。它有一个多项式阶数参数和窗口宽度参数需要调。经验是多项式阶数不宜超过窗口宽度的一半否则容易过拟合。对于峰值检测需求S-G滤波几乎是首选。下面这段代码对比了这三种滤波方法的效果% 构造带噪声的测试信号 t 0:0.001:1; x sin(2*pi*50*t) 0.2*sin(2*pi*150*t) 0.1*randn(size(t)); x(100:105) 5 * ones(1,6); % 加入脉冲噪声 % 中值滤波 med_filt medfilt1(x, 5); % 滑动平均 mean_filt movmean(x, 5); % S-G滤波 sg_filt sgolayfilt(x, 3, 11);实际操作时可以同时画出三种滤波结果的对比图。我做完对比之后常说一句话滤波方法没有绝对的好坏只有合不合适。中值滤波适合去脉冲滑动平均适合去随机噪声S-G滤波适合在去噪声的同时保留波形细节。4.4 滤波效果的评估滤波效果不能靠看起来平滑了来判断量化评估是必须的。我常用的评估手段包括时域和频域两个维度。时域维度看的是滤波前后的波形对比重点观察有没有相位偏移、幅值衰减、边缘效应。用一个简单的办法验证对一段纯正弦信号做滤波计算滤波前后的幅值比和相位差。如果幅值比不符合预期需要检查滤波器设计参数如果存在明显相位差但用filtfilt处理了却仍然存在那就得检查滤波器是否设计错误。频域维度是更本质的评估方式——把滤波前和滤波后的信号分别做FFT对比频谱图。滤波后噪声频段的幅值应该明显下降有效频段的幅值应基本不变。除此之外波形能量守恒也是一个可以检查的指标滤波前后信号的总能量时域平方和和频域能量幅值平方和应遵循Parseval定理的关系。5. 常见问题与排查技巧实录5.1 FFT频谱异常与频谱泄露问题症状频谱图上本应是一个尖锐峰的位置却出现了一个宽宽的山包两侧还有很多裙摆一样的旁瓣。原因几乎都是信号截断造成的频谱泄露。解决方法就是加窗具体窗口选择因人而异但汉宁窗是最稳妥的起点。还有一种情况是信号频率落在两条谱线之间这时主瓣会明显展宽峰值幅度也会下降可以通过前面提到过的抛物线插值来恢复峰值频率和幅度的精度。5.2 谐波识别不准或THD计算偏差症状某些高次谐波分量在频谱图上看得到但自动识别算法找不到或者THD计算结果明显比理论值偏大。排查方向有两个一是频率分辨率不够谐波峰值和旁瓣或者相邻谐波之间没有分开阈值设置过高会丢掉小谐波二是噪声底抬高了整个频谱导致谐波峰被淹没。处理方式可以选择更长的采样时长来提升分辨率或者对多次FFT做平均处理来压低噪声底。我做电力信号分析时经常用10次以上的平均来保证高次谐波识别的稳定性。5.3 滤波结果出现边界剧烈抖动症状用filtfilt滤波后信号的开头和结尾有几段明显的大幅振荡像振铃。这是滤波器初始条件处理不当的典型表现。filtfilt函数理论上会估计并补偿初始条件但信号开头如果含有很大的直流偏置或异常跳变滤波器状态估算会失准。解决思路是给信号做预处理先去除均值掐掉开头一小段暂态再进滤波器。如果你用的是filter而不是filtfilt边界振铃会更严重必须在滤波前对初始状态做合理的赋值或者丢弃滤波输出信号前几百个点。5.4 实时处理时滤波延时不可接受症状在线监测系统里滤波后的信号相比原始信号有明显滞后导致控制逻辑判断迟钝。任何因果滤波器即只能利用当前和过去数据都会引入滞后这是物理规律。应对方式是重新审视需求如果对延时敏感就应该放弃高阶滤波器改用一个低阶滤波器并接受其过渡带较宽的代价或者使用零相位方案比如引入固定长度的延迟缓存配合后处理方式逼近零相位效果。延时和滤波效果必须做折中不存在两全其美的选择这个观念要尽早建立。5.5 MATLAB仿真与真实数据差异巨大的原因症状在MATLAB中用仿真信号测试的滤波和FFT分析效果非常好一换到真实采集数据就完全不对了。这通常不是因为算法错了而是因为真实数据的特性远比仿真复杂真实信号往往是非平稳的频率会随时间漂移传感器本身带有非线性失真环境噪声不是理想白噪声而是有色噪声。应对方法很简单在设计算法时不要只测仿真信号要拿到一段真实数据作为测试样本把算法的参数窗长、阶数、截止频率在真实数据上验证一遍。我一直坚持算法必须用真实数据验收的原则这个习惯帮我避开了无数次方案被推翻的尴尬。6. 项目整体中我在实践中的体会做FFT分析和滤波中间有无数次我以为是信号问题的排查最终还是回到采样参数、窗函数选择、边界处理这些细枝末节上。把每一个环节的关键参数写清楚、把每一步的手法记录下来比写好一整段代码更重要。基于这个项目的积累我有几点操作层面的建议作为一个实践者而非教程作者分享给大家。FFT分析前花5分钟把采样率和分辨率的关系算清楚比盲目调窗函数节省的时间多得多。花1分钟把信号均值减掉频谱会变得干净很多。滤波参数的选取不要照搬任何推荐值一定要结合自己信号的频带范围现场试验。预备好一个小工具箱把FFT分析、谐波提取、THD计算、各种滤波方法都封装成独立函数项目到了后期你会发现这套工具在几乎所有信号类任务里都能复用。最后再分享一个我在项目中后期追加的技巧把FFT分析的结果与滤波的过程做成一条龙的处理管线用同一组参数自动完成分析-滤波-再分析的闭环验证。这样每次调整滤波器参数后我可以立即看到频谱上的变化而不是通过肉眼观察时域波形来推测效果。在实际工程中这种可见的反馈比任何理论计算都更有说服力。