很多做信号处理的人应该都有过这种经历手里有一段实测数据频谱图上看得到几个明显的峰值但就是想不出来怎么把对应的分量干脆利落地拆出来。用固定带通滤波吧你得先知道频率在哪用EMD吧模态混叠能让你调一整天用VMD吧又在纠结模态数K到底设几。我之前分析一段振动信号时就卡在这个问题上后来折腾出一套基于傅里叶分析的3级自适应信号分解方法先通过FFT谱峰定位自动找出当前最强分量用频域掩码把它剥出去剩下的残差继续做同样的事最多迭代三级弱分量也能稳定分离。整个过程在MATLAB里实现起来并不复杂但参数细节和边界条件相当多。这个思路很适合处理少数窄带分量叠加背景噪声的信号比如轴承故障诊断、电力谐波分析、生物电信号分离。下面把我的实现过程和踩坑经验完整写出来供做信号分析、状态监测和MATLAB编程的同学参考。1. 为什么用剥洋葱思路傅里叶谱定位逐级剥离的优势1.1 传统傅里叶分析的局限傅里叶变换的最大价值是能直观告诉我们信号里有哪些频率成分但它的输出是一个全局平均结果。举个例子你用FFT看到频谱图上在3Hz、27Hz、83Hz三个位置有尖峰这只能说明信号里存在这几个周期性分量却拿不出每个分量在时域里长什么样、幅值怎么变化、包络有没有调制。很多后续分析工作恰恰需要时域波形比如计算瞬时相位、做包络解调、提取故障特征频率这时候傅里叶变换就帮不上忙了。另一个麻烦是傅里叶分析用固定基函数拟合信号分辨率受到数据长度的限制。频率间隔小于fs/N的两个分量在频谱上会混成一个大鼓包你既分不清有几个峰也没法准确定位峰的位置。这也是后面第4章要聊的栅栏效应问题的根源。1.2 EMD/VMD等自适应分解的痛点既然固定基不行大家自然想到自适应分解。经验模态分解EMD不用预设基函数能根据信号本身的时间尺度逐层筛出内禀模态函数听起来很美好用起来却有不少坑。模态混叠是大问题不同频率的成分可能被分到同一个模态里或者一个成分被拆到好几个模态里。端点效应也很头疼数据短的时候首尾会出现大幅度摆动。VMD整体上比EMD稳健它是把信号分解成若干个窄带模态每个模态中心频率和带宽通过优化自动确定。但VMD有个绕不开的门槛你得预先设定模态数K和惩罚参数α。K设小了大分量拆不干净K设大了会把一个分量劈成两半。问题在于很多工程信号里到底有几个可分离的窄带分量事先根本不知道往往是看了频谱图才大概知道。这就形成了一个尴尬的循环VMD要求你先知道答案再分解而FFT能给你答案却拿不出分量本身。1.3 3级自适应的基本判断我当时的想法很朴素既然FFT谱峰代表了当前信号里最强势的周期性成分那就每次只提取这一个最强成分剥掉之后再看剩余残差的频谱直到残差里没有明显的窄带尖峰为止。这个过程就像剥洋葱一层一层来。为什么设计成3级而不是任意级或者干脆预设一个固定数量因为我发现工程信号中可辨识的窄带分量通常在2到4个左右第一级把最强的基波或主频成分拿走第二级分离次强的倍频或故障特征频率第三级处理弱小的残余分量。三级之后如果残差里还有大面积能量说明剩余部分已经不是稀疏窄带结构继续分解只会拆出噪声。所以这里的3级是一个可配置的上限而不是必须拆满三个。实际级数由停止准则决定这一点和VMD预设模态数有本质区别。2. 方法总体框架从粗到细的三级处理流程2.1 每一级的具体任务3级自适应分解的每一级做的工作其实是同一件事只是参数尺度越来越精细。第一级是全局频谱粗定位。对原始信号做FFT搜索幅值谱峰值跳过直流估计主峰周围的谱峰宽度构建频域掩码把最强分量c1提取出来得到残差r1 x - c1。第二级是残差信号的精细分解。对r1做FFT但峰值搜索前先做补零FFT进行频率细化因为残差里的分量弱栅栏效应造成的频率估计偏差会更明显。同时带宽策略要收缩以更窄的频率范围提取次强分量c2残差r2 r1 - c2。第三级是弱分量确认与收敛判断。对r2继续做细化FFT用更窄的带宽提取c3然后检查残差能量占比以及c3的谱峰显著性判断分解是否可以停止。2.2 为什么用频域掩码而不是滤波器传统做法是用带通滤波器把目标频带滤出来但我选择直接在频域构造掩码原因是频域掩码天然是零相位的。先用FFT把信号变换到频域乘上一个只在目标频带附近为1的掩码再用IFFT变回时域整个过程没有IIR/FIR滤波器那种群延迟分量和残差从第一级开始就能在时间轴上严格对齐。这对后面做重构验证和瞬时特征提取非常关键。可能有人会说filtfilt也能做到近似零相位。filtfilt确实可以但它对边界做的是反射延拓在短数据上会引起额外的边界失真。频域掩码的边界效应主要由窗函数旁瓣决定物理上更可控。当然矩形的硬截断掩码会在时域引入Gibbs振铃这个问题的解决办法我放在第6章详细讲。2.3 整体调用结构整个流程抽象出来是这样对每一级 k 1,2,3: 1. 对当前输入信号 u第一级ux之后ur_{k-1}做FFT 2. 搜索幅值谱主峰跳过直流必要时跳过已被提取的频带 3. 估计主峰半功率宽度乘以安全系数得到提取带宽 4. 构建带过渡带的频域掩码单边谱 → 共轭对称扩展 5. c_k ifft(fft(u) .* mask)r_k u - c_k 6. 检查停止准则不满足则进入下一级这里每一级输入的是上一级的残差而不是每次都用原始信号减去所有分量的累和。两者在数学上等价但逐级更新在数值上更干净因为残差里减掉的是当级重构得到的分量可以避免把前一级的微小误差累积传递。3. 第一级全局FFT峰值检测与宽频带分离3.1 数据准备和频谱预处理我用一个经典的仿真信号来演示三个正弦分量分别位于3Hz、27Hz、83Hz幅值分别为2.0、0.7、0.3叠加标准差为0.05的白噪声。这样设计有两个目的一是模拟强弱分量共存的真实场景二是有已知真值可以验证分解结果。%% 第1级全局频谱峰值定位 fs 500; % 采样率(Hz) N 2048; % 采样点数 t (0:N-1) / fs; % 时间列向量 rng(42); f1 3; A1 2.0; f2 27; A2 0.7; f3 83; A3 0.3; x A1*sin(2*pi*f1*t 0.5) ... A2*sin(2*pi*f2*t 1.0) ... A3*sin(2*pi*f3*t 2.2) 0.05*randn(N,1); X fft(x); df fs / N; % 0.2441 Hz freq (0:N-1) * df; half floor(N/2) 1; % 1025 amp abs(X(1:half)); amp(1) 0; % 去掉直流 [~, k1] max(amp); % 谱峰索引 f_peak1 freq(k1); fprintf(第1级主峰频率: %.3f Hz\n, f_peak1);这里有几个预处理细节。第一amp(1) 0把直流分量的谱线强制置零否则峰值搜索很容易搜到0Hz附近。第二如果信号有明显的趋势项最好在FFT前做一次去均值甚至多项式去趋势不然0Hz附近的低频泄漏会污染附近的谱结构。第三如果主峰很可能是极低频建议直接把频率小于0.5Hz的谱线全部置零后再搜索避免把传感器零漂当成有效分量。3.2 半功率宽度估计与带宽选择找到谱峰后需要估计这个峰在频谱上占据的宽度。我采用的方案是半功率宽度法先取峰值的1/sqrt(2)作为阈值然后从峰值位置向左右两边搜索找到频谱幅值第一次低于阈值的位置这两个位置之间的频率间隔就是-3dB带宽。peak_amp amp(k1); th peak_amp / sqrt(2); li k1; ri k1; while li 2 amp(li-1) th li li - 1; end while ri half amp(ri1) th ri ri 1; end bw3dB freq(ri) - freq(li); bw_extract1 bw3dB * 1.8; % 安全系数为什么带宽要乘以1.8而不是直接用半功率宽度原因是半功率宽度只覆盖了主瓣的中心区域如果按这个宽度直接截取相当于给提取的分量加了一个过窄的频域窗时域包络会被拉长展宽重构出来和真实分量对不上。乘上1.5到2.0的系数把主瓣大部分能量包进来提取的分量在时域上才更接近原始分量。第一级我倾向用1.8因为最强分量的谱峰一般很干净带宽取宽一点不太担心混入别的分量前提是其他分量离得足够远。3.3 频域掩码构造带宽确定后构造一个带余弦过渡带的频域掩码。直接在通带内赋1、通带外赋0是矩形掩码等效于在时域乘以sinc函数重构波形的首尾会有明显的振铃。因此我在掩码边缘加了一段余弦渐变让掩码幅值从0平滑过渡到1。mask1 zeros(half, 1); [~, p_lo] min(abs(freq(1:half) - (f_peak1 - bw_extract1/2))); [~, p_hi] min(abs(freq(1:half) - (f_peak1 bw_extract1/2))); p_lo max(2, p_lo); p_hi min(half, p_hi); mask1(p_lo:p_hi) 1; % 过渡带宽度取通带宽度的20%或至少4条谱线 trans max(4, round((p_hi - p_lo) * 0.2)); if p_lo - trans 1 ramp 0.5 - 0.5*cos(linspace(0, pi, trans2)); mask1(p_lo-trans:p_lo-1) ramp(2:end-1); end if p_hi trans half ramp 0.5 0.5*cos(linspace(0, pi, trans2)); mask1(p_hi1:p_hitrans) ramp(2:end-1); end % 共轭对称扩展 mask1_full [mask1; mask1(end-1:-1:2)]; Xf1 X .* mask1_full; c1 real(ifft(Xf1)); r1 x - c1;这里有一个必须注意的细节实数信号的FFT频谱是共轭对称的正频率部分的掩码确定后负频率部分必须按镜像扩展。mask1_full [mask1; mask1(end-1:-1:2)]就是做这件事把正频部分的掩码反转复制到负频部分。如果不做这一步IFFT结果是复数虚部不为零提取出的时域信号就不对了。3.4 第一级的验证运行完这段代码后可以验证c1是否与真实的3Hz分量对应。首先看c1的频谱是否只在3Hz附近有能量其次计算c1与x的相关系数再看r1的频谱里3Hz的谱峰是否被去掉、27Hz和83Hz的谱峰是否完好保留。我实测的结果是第一级提取的c1幅值估计约1.98和真实的2.0非常接近。r1的频谱里3Hz附近的谱峰基本消失但27Hz和83Hz的谱线几乎不受影响。这就是频域掩码的好处掩码只改变目标频带的频谱成分其他频带保持原样不会像时域滤波器那样产生通带外失真。4. 第二级残差信号的Zoom-FFT频率细化与窄带提取4.1 栅栏效应为什么第二级要先细化频率再选掩码第一级提取的是最强分量峰值很突出栅栏效应引起的频率偏差可能不算严重。但到了第二级残差里的分量幅度较弱谱峰可能只高出噪声几个dB此时峰值索引对应的频率很可能偏离真实频率达半个频率分辨率。具体来说N2048、fs500时频率分辨率约0.244Hz。如果27Hz的分量真实频率落在两条谱线中间粗搜索给出的峰值频率可能偏了0.1Hz以上。用这个有偏差的频率作为带通中心弱分量会被切掉一部分能量因为中心位置和真实谱峰位置不重合掩码会把有效信号推到过渡带边缘甚至截断。这一点在第一级影响不大因为强分量峰宽很大偏差0.1Hz可以忽略但第二级和第三级处理弱分量时不细化频率会明显降低重构精度。4.2 补零FFT与抛物线插值最直观的频率细化手段是补零FFT。对残差信号做8倍补零FFT相当于在原来的频谱上插入更密的采样网格峰值位置的读数就精细多了。%% 第2级残差频率细化 N_zpad 8 * N; % 8倍补零 Xz fft(r1, N_zpad); fz (0:N_zpad-1) * fs / N_zpad; half_z N_zpad/2 1; amp_z abs(Xz(1:half_z)); amp_z(1) 0; [~, kz] max(amp_z); f_coarse2 fz(kz); % 三点抛物线插值 if kz 1 kz half_z y0 amp_z(kz-1); y1 amp_z(kz); y2 amp_z(kz1); delta 0.5 * (y0 - y2) / (y0 - 2*y1 y2); f_peak2 f_coarse2 delta * (fs / N_zpad); else f_peak2 f_coarse2; end fprintf(第2级主峰频率: %.4f Hz\n, f_peak2);补零FFT之后再做三点抛物线插值可以把峰值位置估计到亚网格精度。我用过不同数据测试插值后的频率估计偏差通常只有原始分辨率的一两成相对纯补零又能进一步改进。这里必须澄清一个常见误区补零FFT并没有提高真实频率分辨率。真实分辨率由观测时长决定N个点、采样率fs时能分辨的最小频率间隔就是fs/N。补零只是把DTFT采样得更密让原来卡在两条谱线中间的峰值不再被两边拉低而不是把靠得很近的两个峰分开。如果两个分量频率只差0.1Hz而你的数据长度对应的分辨率是0.244Hz那无论补多少零都分不开它们频谱上只会是一个鼓包。这种情况只能加长观测时间或者用参数化方法比如子空间法去估计。4.3 第二级带宽的自适应收缩残差中次强分量的谱峰通常比第一级更瘦因为主分量被剥离后没有强分量的旁瓣压着它峰形更接近窗函数的主瓣。因此第二级的带宽策略要做两件事一是继续用半功率宽度自适应估计二是把安全系数从1.8降到1.4左右。% 基于细化频谱估计半功率宽度 li2 kz; ri2 kz; th2 amp_z(kz) / sqrt(2); while li2 2 amp_z(li2-1) th2 li2 li2 - 1; end while ri2 half_z amp_z(ri21) th2 ri2 ri2 1; end bw3dB2 fz(ri2) - fz(li2); % 安全系数1.4并保留下限 bw_extract2 max(bw3dB2 * 1.4, 4*fs/N_zpad);下限4*fs/N_zpad的意思是至少覆盖4条原始频率分辨率对应的带宽也就是约1Hz避免因为谱峰太瘦导致提取带过窄、把真实分量的一部分能量切掉。这个下限值不能设得太大否则会把旁边的噪声频带包进来。4.4 提取与验证第二级的掩码构造和第一级相同只是把中心频率换成细化后的f_peak2带宽换成bw_extract2。提取后得到c2残差r2 r1 - c2。在实际运行中第二级提取的分量c2和真实27Hz正弦分量的相关系数通常能达到0.99以上。我专门做过对比如果直接用粗峰值频率不做补零细化构造掩码相关系数可能掉到0.86左右细化之后能到0.995。差距主要体现在弱分量上所以这一步不是可有可无的优化而是弱分量能否干净拆出的关键。5. 第三级弱分量提取与停止准则判断5.1 第三级提取的流程第三级处理的是r2此刻残差里只剩下83Hz弱分量和噪声。流程和第二级几乎一样但有两个参数进一步收缩补零倍数从8倍提到16倍频率网格更细带宽安全系数从1.4降到1.3下限压到3条原始谱线对应的频宽。%% 第3级弱分量提取 Nz3 16 * N; % 16倍补零 Xz3 fft(r2, Nz3); fz3 (0:Nz3-1) * fs / Nz3; half3 Nz3/2 1; amp3 abs(Xz3(1:half3)); amp3(1) 0; [~, k3] max(amp3); if k3 1 k3 half3 y0 amp3(k3-1); y1 amp3(k3); y2 amp3(k31); delta 0.5 * (y0 - y2) / (y0 - 2*y1 y2); f_peak3 fz3(k3) delta * (fs / Nz3); else f_peak3 fz3(k3); end th3 amp3(k3) / sqrt(2); li3 k3; ri3 k3; while li3 2 amp3(li3-1) th3 li3 li3 - 1; end while ri3 half3 amp3(ri31) th3 ri3 ri3 1; end bw_extract3 max((fz3(ri3) - fz3(li3)) * 1.3, 3*fs/Nz3);5.2 谱峰显著性判断防止把噪声当分量到第三级残差里可能已经主要是噪声了。如果不管三七二十一都按最大峰值去提取很可能把某一根噪声谱线当成分量拆出来得到的c3其实是噪声的窄带实现没有物理意义。所以第三级必须加显著性检验。我的做法是用残差频谱的中位数来估计噪声底再计算峰值相对噪声底的倍数noise_floor median(amp3(2:half3)); peak_to_floor amp3(k3) / noise_floor; if peak_to_floor 8 % 通过显著性检验继续提取c3 else fprintf(第3级谱峰不显著(峰值/噪声底%.1f 8)停止分解\n, peak_to_floor); end为什么用中位数而不是均值因为白噪声FFT的谱线幅值分布是右偏的少数强谱线会拉高均值中位数更稳健。8倍这个阈值的统计背景是对白噪声做N2048点FFT最大的谱线幅值与中位数的比通常不超过五六倍。所以8倍可以认为不是纯噪声造成的峰值。数据越长噪声谱越平稳这个比值阈值可以适当放宽到7数据短时要取高一些比如10。5.3 重构验证完成三级分解后一定要做一次完整重构验证x_recon c1 c2 c3 r3; recon_error max(abs(x_recon - x)); fprintf(重构最大误差: %.3e\n, recon_error); for k 1:3 ck eval([c num2str(k)]); est_amp sqrt(2) * rms(ck); fprintf(分量%d 估计幅值: %.3f, 能量占比: %.2f%%\n, ... k, est_amp, sum(ck.^2)/sum(x.^2)*100); end重构误差理论上应该在10的负12次方量级因为x_recon c1 c2 c3 r3在数值上恒等于x。真正需要关注的是各分量的幅值估计和能量占比。在我的仿真里三个分量的幅值估计应该分别接近2.000、0.700、0.300残差能量占比基本就是噪声能量占比大约在0.1%左右。如果c3的估计幅值明显大于0.3说明前两级带宽选择偏宽、把别人的能量也包进来了需要回看第一级和第二级的掩码宽度。如果c3幅值偏小多半是第三级带宽太窄把分量边缘切掉了。5.4 如果第三级之后还有明显谱峰怎么办这不是方法失效而是信号本身不是三分量窄带叠加的形态。比如遇到五六个谐波分量共存或者强噪声中的连续谱结构三级剥完残差谱仍然有尖峰。我有两个建议。一是把级数上限从3调到5继续按同样流程剥。我在处理转轴振动信号时遇到过一次七个分量的情况调到5级后依然稳定。二是检查是不是前两级把某个分量拆到了相邻频带导致残差互调。如果级数调到很高还在出峰那基本说明信号包含宽带成分或者分量之间频率间隔小于傅里叶分辨率这时候用VMD或EMD会更合适。这套逐级剥离的框架天生适合分量数量少、窄带特征明确的信号不适合当万能工具用。6. 参数调试与实测经验跑数据时踩过的坑6.1 峰值搜索被直流或泄漏污染我最早跑实测数据时第一级峰值搜索经常搜到0Hz附近去提取出来的最强分量是一条趋势项。原因就是传感器零漂和直流偏移没去掉。FFT的0Hz谱线幅值极大即使置零了旁边的低频泄漏也远高于真实信号分量。处理办法是先去均值必要时做一个多项式趋势去除搜索前把低于0.5Hz的谱线全部置零再考虑加窗比如Hamming窗压制频谱泄漏。但加窗会改变主瓣宽度影响半功率宽度的估计所以我在仿真演示里没有加窗实际处理时要根据信号特性权衡。6.2 掩码过渡带宽度对重构波形的影响矩形掩码导致的Gibbs振铃是最容易忽略的问题。我用一个简单例子说明正弦信号在频域就一根谱线矩形掩码刚好包住它逆变换回去几乎无失真但如果信号频率不在整数谱线上矩形截断会把频谱边缘切在不该切的位置时域波形首尾就出现抖动。解决方法是把掩码边缘做成余弦过渡带让频谱幅值从通带到阻带平滑渐变。过渡带太窄一两条谱线效果不明显我一般取通带宽度的15%到20%最少四到六条谱线。代价是掩码选择性略降过渡带内可能有少许泄漏但和振铃相比这点泄漏影响小得多。6.3 补零FFT的误区插值不等于提高分辨率这一点我在第4章提过但还是要再强调一次。补零FFT只是让频谱采样更密峰值位置读数更精细并不能分辨原本混叠的两个频率分量。如果你的数据长度对应分辨率是0.244Hz而两个分量只差0.2Hz补再多的零也分不开它们。正确做法是增加观测时长。如果你只能处理已有的一段数据那可以考虑用MUSIC或ESPRIT这类子空间方法去估计频率它们不受傅里叶分辨率的限制代价是对模型阶数敏感、计算量大。我把补零FFT定位成细读而不是分辨就不会踩这个坑。6.4 关于停止阈值别把噪声当信号峰值/噪声底比值取8、残差能量阈值取2%到5%这些数值是我在不同实测信号里试出来的比较稳妥的初始值。一个更稳的补充判据是把提取出的分量c_k的包络画出来如果包络变化剧烈、没有稳定的周期结构大概率是拆出了噪声。在自动化处理流程里可以计算包络的调制指数包络标准差除以包络均值超过一定阈值就判定为噪声分量不保留。不同数据长度的阈值调整经验参数第1级第2级第3级备注频率细化倍数1原始FFT816补零倍数半功率宽度乘数1.81.41.3带宽安全系数最小带宽下限6*df4*df3*dfdf为原始频率分辨率过渡带宽度20%通带宽或≥5条谱线同左同左抑制Gibbs振铃峰值/噪声底阈值不检查≥10≥8防噪声误判残差能量停止阈值--≤2%~5%噪声水平决定6.5 边界效应短数据要注意首尾失真频域掩码虽然零相位但在数据很短时掩码截断引起的边界效应依然存在。一个实用技巧分解完成后把各分量首尾各1%长度的采样点用淡入淡出窗口处理再重新计算残差防止边界失真传播到下一级。如果数据长度少于几百个点建议直接在边缘丢弃1%到2%的样本再做分析因为边界信息本来就不可靠留着只会误导后续的瞬时频率计算。7. 适用范围与扩展方向7.1 适合的信号类型这套方法最舒服的场景是窄带分量数量少2到5个、频率间隔大于三四倍傅里叶分辨率、各分量幅度差距可以很大比如差20dB以上、噪声是宽带平滑背景。典型例子包括齿轮箱振动信号轴频、啮合频率及其边带的分离、电力谐波分析基波加3次、5次谐波的提取、含呼吸调制的心电信号分离、水声信号中的单频成分提取等。在这些场景下它比固定带通滤波器省去了频率先验比EMD更不容易模态混叠比VMD少了反复调K的麻烦。7.2 不适合的情况如果信号本身是宽频瞬态冲击比如故障轴承初期那种周期性冲击逐级剥离窄带成分的意义就不大。因为冲击在频谱上是连续宽带结构用包络谱或者直接做时域冲击提取更合适。如果信号是强非平稳调频信号频率随时间快速变化单次FFT的主频概念本身就会失效应该用短时傅里叶或时频分析来指导分解。7.3 可扩展的几个方向第一提取出的c1、c2、c3可以继续做Hilbert变换直接得到各分量的瞬时幅值和瞬时频率构造完整的时频特征图这在机械故障诊断和生理信号分析里都很有用。第二带宽收缩系数和过渡带宽度可以通过优化算法自动搜索比如用遗传算法找使残差能量最小化的参数组合。但这类算法计算量大一般只在批量离线分析时值得做在线实时处理还是用固定经验值稳妥。第三配合MUSIC或ESPRIT估计频率可以突破傅里叶分辨率的限制处理频率非常接近的工程信号。我自己的体会是先把这套三级傅里叶剥离流程作为第一遍快速侦察筛查大多数信号如果发现确实存在频谱上挤在一起的峰再针对性地用子空间方法去精细估计这样多数问题已经能解决一大半。这套方法最初只是我处理振动信号时临时搭的方案B后来在好几个项目里反而成了首选因为它流程简单、参数透明、每一步的频谱图都能直接拿来做人工检查。如果让我总结一条最重要的经验就是永远不要只看最终分解结果要把每一级的残差频谱和提取分量的时域波形打印出来看一眼大部分问题都藏在中间级里。希望这篇实现记录能帮你少走几步弯路如果你在自己的数据上试出了不同的系数规律欢迎一起交流。
基于FFT谱峰定位的3级自适应信号分解方法详解
很多做信号处理的人应该都有过这种经历手里有一段实测数据频谱图上看得到几个明显的峰值但就是想不出来怎么把对应的分量干脆利落地拆出来。用固定带通滤波吧你得先知道频率在哪用EMD吧模态混叠能让你调一整天用VMD吧又在纠结模态数K到底设几。我之前分析一段振动信号时就卡在这个问题上后来折腾出一套基于傅里叶分析的3级自适应信号分解方法先通过FFT谱峰定位自动找出当前最强分量用频域掩码把它剥出去剩下的残差继续做同样的事最多迭代三级弱分量也能稳定分离。整个过程在MATLAB里实现起来并不复杂但参数细节和边界条件相当多。这个思路很适合处理少数窄带分量叠加背景噪声的信号比如轴承故障诊断、电力谐波分析、生物电信号分离。下面把我的实现过程和踩坑经验完整写出来供做信号分析、状态监测和MATLAB编程的同学参考。1. 为什么用剥洋葱思路傅里叶谱定位逐级剥离的优势1.1 传统傅里叶分析的局限傅里叶变换的最大价值是能直观告诉我们信号里有哪些频率成分但它的输出是一个全局平均结果。举个例子你用FFT看到频谱图上在3Hz、27Hz、83Hz三个位置有尖峰这只能说明信号里存在这几个周期性分量却拿不出每个分量在时域里长什么样、幅值怎么变化、包络有没有调制。很多后续分析工作恰恰需要时域波形比如计算瞬时相位、做包络解调、提取故障特征频率这时候傅里叶变换就帮不上忙了。另一个麻烦是傅里叶分析用固定基函数拟合信号分辨率受到数据长度的限制。频率间隔小于fs/N的两个分量在频谱上会混成一个大鼓包你既分不清有几个峰也没法准确定位峰的位置。这也是后面第4章要聊的栅栏效应问题的根源。1.2 EMD/VMD等自适应分解的痛点既然固定基不行大家自然想到自适应分解。经验模态分解EMD不用预设基函数能根据信号本身的时间尺度逐层筛出内禀模态函数听起来很美好用起来却有不少坑。模态混叠是大问题不同频率的成分可能被分到同一个模态里或者一个成分被拆到好几个模态里。端点效应也很头疼数据短的时候首尾会出现大幅度摆动。VMD整体上比EMD稳健它是把信号分解成若干个窄带模态每个模态中心频率和带宽通过优化自动确定。但VMD有个绕不开的门槛你得预先设定模态数K和惩罚参数α。K设小了大分量拆不干净K设大了会把一个分量劈成两半。问题在于很多工程信号里到底有几个可分离的窄带分量事先根本不知道往往是看了频谱图才大概知道。这就形成了一个尴尬的循环VMD要求你先知道答案再分解而FFT能给你答案却拿不出分量本身。1.3 3级自适应的基本判断我当时的想法很朴素既然FFT谱峰代表了当前信号里最强势的周期性成分那就每次只提取这一个最强成分剥掉之后再看剩余残差的频谱直到残差里没有明显的窄带尖峰为止。这个过程就像剥洋葱一层一层来。为什么设计成3级而不是任意级或者干脆预设一个固定数量因为我发现工程信号中可辨识的窄带分量通常在2到4个左右第一级把最强的基波或主频成分拿走第二级分离次强的倍频或故障特征频率第三级处理弱小的残余分量。三级之后如果残差里还有大面积能量说明剩余部分已经不是稀疏窄带结构继续分解只会拆出噪声。所以这里的3级是一个可配置的上限而不是必须拆满三个。实际级数由停止准则决定这一点和VMD预设模态数有本质区别。2. 方法总体框架从粗到细的三级处理流程2.1 每一级的具体任务3级自适应分解的每一级做的工作其实是同一件事只是参数尺度越来越精细。第一级是全局频谱粗定位。对原始信号做FFT搜索幅值谱峰值跳过直流估计主峰周围的谱峰宽度构建频域掩码把最强分量c1提取出来得到残差r1 x - c1。第二级是残差信号的精细分解。对r1做FFT但峰值搜索前先做补零FFT进行频率细化因为残差里的分量弱栅栏效应造成的频率估计偏差会更明显。同时带宽策略要收缩以更窄的频率范围提取次强分量c2残差r2 r1 - c2。第三级是弱分量确认与收敛判断。对r2继续做细化FFT用更窄的带宽提取c3然后检查残差能量占比以及c3的谱峰显著性判断分解是否可以停止。2.2 为什么用频域掩码而不是滤波器传统做法是用带通滤波器把目标频带滤出来但我选择直接在频域构造掩码原因是频域掩码天然是零相位的。先用FFT把信号变换到频域乘上一个只在目标频带附近为1的掩码再用IFFT变回时域整个过程没有IIR/FIR滤波器那种群延迟分量和残差从第一级开始就能在时间轴上严格对齐。这对后面做重构验证和瞬时特征提取非常关键。可能有人会说filtfilt也能做到近似零相位。filtfilt确实可以但它对边界做的是反射延拓在短数据上会引起额外的边界失真。频域掩码的边界效应主要由窗函数旁瓣决定物理上更可控。当然矩形的硬截断掩码会在时域引入Gibbs振铃这个问题的解决办法我放在第6章详细讲。2.3 整体调用结构整个流程抽象出来是这样对每一级 k 1,2,3: 1. 对当前输入信号 u第一级ux之后ur_{k-1}做FFT 2. 搜索幅值谱主峰跳过直流必要时跳过已被提取的频带 3. 估计主峰半功率宽度乘以安全系数得到提取带宽 4. 构建带过渡带的频域掩码单边谱 → 共轭对称扩展 5. c_k ifft(fft(u) .* mask)r_k u - c_k 6. 检查停止准则不满足则进入下一级这里每一级输入的是上一级的残差而不是每次都用原始信号减去所有分量的累和。两者在数学上等价但逐级更新在数值上更干净因为残差里减掉的是当级重构得到的分量可以避免把前一级的微小误差累积传递。3. 第一级全局FFT峰值检测与宽频带分离3.1 数据准备和频谱预处理我用一个经典的仿真信号来演示三个正弦分量分别位于3Hz、27Hz、83Hz幅值分别为2.0、0.7、0.3叠加标准差为0.05的白噪声。这样设计有两个目的一是模拟强弱分量共存的真实场景二是有已知真值可以验证分解结果。%% 第1级全局频谱峰值定位 fs 500; % 采样率(Hz) N 2048; % 采样点数 t (0:N-1) / fs; % 时间列向量 rng(42); f1 3; A1 2.0; f2 27; A2 0.7; f3 83; A3 0.3; x A1*sin(2*pi*f1*t 0.5) ... A2*sin(2*pi*f2*t 1.0) ... A3*sin(2*pi*f3*t 2.2) 0.05*randn(N,1); X fft(x); df fs / N; % 0.2441 Hz freq (0:N-1) * df; half floor(N/2) 1; % 1025 amp abs(X(1:half)); amp(1) 0; % 去掉直流 [~, k1] max(amp); % 谱峰索引 f_peak1 freq(k1); fprintf(第1级主峰频率: %.3f Hz\n, f_peak1);这里有几个预处理细节。第一amp(1) 0把直流分量的谱线强制置零否则峰值搜索很容易搜到0Hz附近。第二如果信号有明显的趋势项最好在FFT前做一次去均值甚至多项式去趋势不然0Hz附近的低频泄漏会污染附近的谱结构。第三如果主峰很可能是极低频建议直接把频率小于0.5Hz的谱线全部置零后再搜索避免把传感器零漂当成有效分量。3.2 半功率宽度估计与带宽选择找到谱峰后需要估计这个峰在频谱上占据的宽度。我采用的方案是半功率宽度法先取峰值的1/sqrt(2)作为阈值然后从峰值位置向左右两边搜索找到频谱幅值第一次低于阈值的位置这两个位置之间的频率间隔就是-3dB带宽。peak_amp amp(k1); th peak_amp / sqrt(2); li k1; ri k1; while li 2 amp(li-1) th li li - 1; end while ri half amp(ri1) th ri ri 1; end bw3dB freq(ri) - freq(li); bw_extract1 bw3dB * 1.8; % 安全系数为什么带宽要乘以1.8而不是直接用半功率宽度原因是半功率宽度只覆盖了主瓣的中心区域如果按这个宽度直接截取相当于给提取的分量加了一个过窄的频域窗时域包络会被拉长展宽重构出来和真实分量对不上。乘上1.5到2.0的系数把主瓣大部分能量包进来提取的分量在时域上才更接近原始分量。第一级我倾向用1.8因为最强分量的谱峰一般很干净带宽取宽一点不太担心混入别的分量前提是其他分量离得足够远。3.3 频域掩码构造带宽确定后构造一个带余弦过渡带的频域掩码。直接在通带内赋1、通带外赋0是矩形掩码等效于在时域乘以sinc函数重构波形的首尾会有明显的振铃。因此我在掩码边缘加了一段余弦渐变让掩码幅值从0平滑过渡到1。mask1 zeros(half, 1); [~, p_lo] min(abs(freq(1:half) - (f_peak1 - bw_extract1/2))); [~, p_hi] min(abs(freq(1:half) - (f_peak1 bw_extract1/2))); p_lo max(2, p_lo); p_hi min(half, p_hi); mask1(p_lo:p_hi) 1; % 过渡带宽度取通带宽度的20%或至少4条谱线 trans max(4, round((p_hi - p_lo) * 0.2)); if p_lo - trans 1 ramp 0.5 - 0.5*cos(linspace(0, pi, trans2)); mask1(p_lo-trans:p_lo-1) ramp(2:end-1); end if p_hi trans half ramp 0.5 0.5*cos(linspace(0, pi, trans2)); mask1(p_hi1:p_hitrans) ramp(2:end-1); end % 共轭对称扩展 mask1_full [mask1; mask1(end-1:-1:2)]; Xf1 X .* mask1_full; c1 real(ifft(Xf1)); r1 x - c1;这里有一个必须注意的细节实数信号的FFT频谱是共轭对称的正频率部分的掩码确定后负频率部分必须按镜像扩展。mask1_full [mask1; mask1(end-1:-1:2)]就是做这件事把正频部分的掩码反转复制到负频部分。如果不做这一步IFFT结果是复数虚部不为零提取出的时域信号就不对了。3.4 第一级的验证运行完这段代码后可以验证c1是否与真实的3Hz分量对应。首先看c1的频谱是否只在3Hz附近有能量其次计算c1与x的相关系数再看r1的频谱里3Hz的谱峰是否被去掉、27Hz和83Hz的谱峰是否完好保留。我实测的结果是第一级提取的c1幅值估计约1.98和真实的2.0非常接近。r1的频谱里3Hz附近的谱峰基本消失但27Hz和83Hz的谱线几乎不受影响。这就是频域掩码的好处掩码只改变目标频带的频谱成分其他频带保持原样不会像时域滤波器那样产生通带外失真。4. 第二级残差信号的Zoom-FFT频率细化与窄带提取4.1 栅栏效应为什么第二级要先细化频率再选掩码第一级提取的是最强分量峰值很突出栅栏效应引起的频率偏差可能不算严重。但到了第二级残差里的分量幅度较弱谱峰可能只高出噪声几个dB此时峰值索引对应的频率很可能偏离真实频率达半个频率分辨率。具体来说N2048、fs500时频率分辨率约0.244Hz。如果27Hz的分量真实频率落在两条谱线中间粗搜索给出的峰值频率可能偏了0.1Hz以上。用这个有偏差的频率作为带通中心弱分量会被切掉一部分能量因为中心位置和真实谱峰位置不重合掩码会把有效信号推到过渡带边缘甚至截断。这一点在第一级影响不大因为强分量峰宽很大偏差0.1Hz可以忽略但第二级和第三级处理弱分量时不细化频率会明显降低重构精度。4.2 补零FFT与抛物线插值最直观的频率细化手段是补零FFT。对残差信号做8倍补零FFT相当于在原来的频谱上插入更密的采样网格峰值位置的读数就精细多了。%% 第2级残差频率细化 N_zpad 8 * N; % 8倍补零 Xz fft(r1, N_zpad); fz (0:N_zpad-1) * fs / N_zpad; half_z N_zpad/2 1; amp_z abs(Xz(1:half_z)); amp_z(1) 0; [~, kz] max(amp_z); f_coarse2 fz(kz); % 三点抛物线插值 if kz 1 kz half_z y0 amp_z(kz-1); y1 amp_z(kz); y2 amp_z(kz1); delta 0.5 * (y0 - y2) / (y0 - 2*y1 y2); f_peak2 f_coarse2 delta * (fs / N_zpad); else f_peak2 f_coarse2; end fprintf(第2级主峰频率: %.4f Hz\n, f_peak2);补零FFT之后再做三点抛物线插值可以把峰值位置估计到亚网格精度。我用过不同数据测试插值后的频率估计偏差通常只有原始分辨率的一两成相对纯补零又能进一步改进。这里必须澄清一个常见误区补零FFT并没有提高真实频率分辨率。真实分辨率由观测时长决定N个点、采样率fs时能分辨的最小频率间隔就是fs/N。补零只是把DTFT采样得更密让原来卡在两条谱线中间的峰值不再被两边拉低而不是把靠得很近的两个峰分开。如果两个分量频率只差0.1Hz而你的数据长度对应的分辨率是0.244Hz那无论补多少零都分不开它们频谱上只会是一个鼓包。这种情况只能加长观测时间或者用参数化方法比如子空间法去估计。4.3 第二级带宽的自适应收缩残差中次强分量的谱峰通常比第一级更瘦因为主分量被剥离后没有强分量的旁瓣压着它峰形更接近窗函数的主瓣。因此第二级的带宽策略要做两件事一是继续用半功率宽度自适应估计二是把安全系数从1.8降到1.4左右。% 基于细化频谱估计半功率宽度 li2 kz; ri2 kz; th2 amp_z(kz) / sqrt(2); while li2 2 amp_z(li2-1) th2 li2 li2 - 1; end while ri2 half_z amp_z(ri21) th2 ri2 ri2 1; end bw3dB2 fz(ri2) - fz(li2); % 安全系数1.4并保留下限 bw_extract2 max(bw3dB2 * 1.4, 4*fs/N_zpad);下限4*fs/N_zpad的意思是至少覆盖4条原始频率分辨率对应的带宽也就是约1Hz避免因为谱峰太瘦导致提取带过窄、把真实分量的一部分能量切掉。这个下限值不能设得太大否则会把旁边的噪声频带包进来。4.4 提取与验证第二级的掩码构造和第一级相同只是把中心频率换成细化后的f_peak2带宽换成bw_extract2。提取后得到c2残差r2 r1 - c2。在实际运行中第二级提取的分量c2和真实27Hz正弦分量的相关系数通常能达到0.99以上。我专门做过对比如果直接用粗峰值频率不做补零细化构造掩码相关系数可能掉到0.86左右细化之后能到0.995。差距主要体现在弱分量上所以这一步不是可有可无的优化而是弱分量能否干净拆出的关键。5. 第三级弱分量提取与停止准则判断5.1 第三级提取的流程第三级处理的是r2此刻残差里只剩下83Hz弱分量和噪声。流程和第二级几乎一样但有两个参数进一步收缩补零倍数从8倍提到16倍频率网格更细带宽安全系数从1.4降到1.3下限压到3条原始谱线对应的频宽。%% 第3级弱分量提取 Nz3 16 * N; % 16倍补零 Xz3 fft(r2, Nz3); fz3 (0:Nz3-1) * fs / Nz3; half3 Nz3/2 1; amp3 abs(Xz3(1:half3)); amp3(1) 0; [~, k3] max(amp3); if k3 1 k3 half3 y0 amp3(k3-1); y1 amp3(k3); y2 amp3(k31); delta 0.5 * (y0 - y2) / (y0 - 2*y1 y2); f_peak3 fz3(k3) delta * (fs / Nz3); else f_peak3 fz3(k3); end th3 amp3(k3) / sqrt(2); li3 k3; ri3 k3; while li3 2 amp3(li3-1) th3 li3 li3 - 1; end while ri3 half3 amp3(ri31) th3 ri3 ri3 1; end bw_extract3 max((fz3(ri3) - fz3(li3)) * 1.3, 3*fs/Nz3);5.2 谱峰显著性判断防止把噪声当分量到第三级残差里可能已经主要是噪声了。如果不管三七二十一都按最大峰值去提取很可能把某一根噪声谱线当成分量拆出来得到的c3其实是噪声的窄带实现没有物理意义。所以第三级必须加显著性检验。我的做法是用残差频谱的中位数来估计噪声底再计算峰值相对噪声底的倍数noise_floor median(amp3(2:half3)); peak_to_floor amp3(k3) / noise_floor; if peak_to_floor 8 % 通过显著性检验继续提取c3 else fprintf(第3级谱峰不显著(峰值/噪声底%.1f 8)停止分解\n, peak_to_floor); end为什么用中位数而不是均值因为白噪声FFT的谱线幅值分布是右偏的少数强谱线会拉高均值中位数更稳健。8倍这个阈值的统计背景是对白噪声做N2048点FFT最大的谱线幅值与中位数的比通常不超过五六倍。所以8倍可以认为不是纯噪声造成的峰值。数据越长噪声谱越平稳这个比值阈值可以适当放宽到7数据短时要取高一些比如10。5.3 重构验证完成三级分解后一定要做一次完整重构验证x_recon c1 c2 c3 r3; recon_error max(abs(x_recon - x)); fprintf(重构最大误差: %.3e\n, recon_error); for k 1:3 ck eval([c num2str(k)]); est_amp sqrt(2) * rms(ck); fprintf(分量%d 估计幅值: %.3f, 能量占比: %.2f%%\n, ... k, est_amp, sum(ck.^2)/sum(x.^2)*100); end重构误差理论上应该在10的负12次方量级因为x_recon c1 c2 c3 r3在数值上恒等于x。真正需要关注的是各分量的幅值估计和能量占比。在我的仿真里三个分量的幅值估计应该分别接近2.000、0.700、0.300残差能量占比基本就是噪声能量占比大约在0.1%左右。如果c3的估计幅值明显大于0.3说明前两级带宽选择偏宽、把别人的能量也包进来了需要回看第一级和第二级的掩码宽度。如果c3幅值偏小多半是第三级带宽太窄把分量边缘切掉了。5.4 如果第三级之后还有明显谱峰怎么办这不是方法失效而是信号本身不是三分量窄带叠加的形态。比如遇到五六个谐波分量共存或者强噪声中的连续谱结构三级剥完残差谱仍然有尖峰。我有两个建议。一是把级数上限从3调到5继续按同样流程剥。我在处理转轴振动信号时遇到过一次七个分量的情况调到5级后依然稳定。二是检查是不是前两级把某个分量拆到了相邻频带导致残差互调。如果级数调到很高还在出峰那基本说明信号包含宽带成分或者分量之间频率间隔小于傅里叶分辨率这时候用VMD或EMD会更合适。这套逐级剥离的框架天生适合分量数量少、窄带特征明确的信号不适合当万能工具用。6. 参数调试与实测经验跑数据时踩过的坑6.1 峰值搜索被直流或泄漏污染我最早跑实测数据时第一级峰值搜索经常搜到0Hz附近去提取出来的最强分量是一条趋势项。原因就是传感器零漂和直流偏移没去掉。FFT的0Hz谱线幅值极大即使置零了旁边的低频泄漏也远高于真实信号分量。处理办法是先去均值必要时做一个多项式趋势去除搜索前把低于0.5Hz的谱线全部置零再考虑加窗比如Hamming窗压制频谱泄漏。但加窗会改变主瓣宽度影响半功率宽度的估计所以我在仿真演示里没有加窗实际处理时要根据信号特性权衡。6.2 掩码过渡带宽度对重构波形的影响矩形掩码导致的Gibbs振铃是最容易忽略的问题。我用一个简单例子说明正弦信号在频域就一根谱线矩形掩码刚好包住它逆变换回去几乎无失真但如果信号频率不在整数谱线上矩形截断会把频谱边缘切在不该切的位置时域波形首尾就出现抖动。解决方法是把掩码边缘做成余弦过渡带让频谱幅值从通带到阻带平滑渐变。过渡带太窄一两条谱线效果不明显我一般取通带宽度的15%到20%最少四到六条谱线。代价是掩码选择性略降过渡带内可能有少许泄漏但和振铃相比这点泄漏影响小得多。6.3 补零FFT的误区插值不等于提高分辨率这一点我在第4章提过但还是要再强调一次。补零FFT只是让频谱采样更密峰值位置读数更精细并不能分辨原本混叠的两个频率分量。如果你的数据长度对应分辨率是0.244Hz而两个分量只差0.2Hz补再多的零也分不开它们。正确做法是增加观测时长。如果你只能处理已有的一段数据那可以考虑用MUSIC或ESPRIT这类子空间方法去估计频率它们不受傅里叶分辨率的限制代价是对模型阶数敏感、计算量大。我把补零FFT定位成细读而不是分辨就不会踩这个坑。6.4 关于停止阈值别把噪声当信号峰值/噪声底比值取8、残差能量阈值取2%到5%这些数值是我在不同实测信号里试出来的比较稳妥的初始值。一个更稳的补充判据是把提取出的分量c_k的包络画出来如果包络变化剧烈、没有稳定的周期结构大概率是拆出了噪声。在自动化处理流程里可以计算包络的调制指数包络标准差除以包络均值超过一定阈值就判定为噪声分量不保留。不同数据长度的阈值调整经验参数第1级第2级第3级备注频率细化倍数1原始FFT816补零倍数半功率宽度乘数1.81.41.3带宽安全系数最小带宽下限6*df4*df3*dfdf为原始频率分辨率过渡带宽度20%通带宽或≥5条谱线同左同左抑制Gibbs振铃峰值/噪声底阈值不检查≥10≥8防噪声误判残差能量停止阈值--≤2%~5%噪声水平决定6.5 边界效应短数据要注意首尾失真频域掩码虽然零相位但在数据很短时掩码截断引起的边界效应依然存在。一个实用技巧分解完成后把各分量首尾各1%长度的采样点用淡入淡出窗口处理再重新计算残差防止边界失真传播到下一级。如果数据长度少于几百个点建议直接在边缘丢弃1%到2%的样本再做分析因为边界信息本来就不可靠留着只会误导后续的瞬时频率计算。7. 适用范围与扩展方向7.1 适合的信号类型这套方法最舒服的场景是窄带分量数量少2到5个、频率间隔大于三四倍傅里叶分辨率、各分量幅度差距可以很大比如差20dB以上、噪声是宽带平滑背景。典型例子包括齿轮箱振动信号轴频、啮合频率及其边带的分离、电力谐波分析基波加3次、5次谐波的提取、含呼吸调制的心电信号分离、水声信号中的单频成分提取等。在这些场景下它比固定带通滤波器省去了频率先验比EMD更不容易模态混叠比VMD少了反复调K的麻烦。7.2 不适合的情况如果信号本身是宽频瞬态冲击比如故障轴承初期那种周期性冲击逐级剥离窄带成分的意义就不大。因为冲击在频谱上是连续宽带结构用包络谱或者直接做时域冲击提取更合适。如果信号是强非平稳调频信号频率随时间快速变化单次FFT的主频概念本身就会失效应该用短时傅里叶或时频分析来指导分解。7.3 可扩展的几个方向第一提取出的c1、c2、c3可以继续做Hilbert变换直接得到各分量的瞬时幅值和瞬时频率构造完整的时频特征图这在机械故障诊断和生理信号分析里都很有用。第二带宽收缩系数和过渡带宽度可以通过优化算法自动搜索比如用遗传算法找使残差能量最小化的参数组合。但这类算法计算量大一般只在批量离线分析时值得做在线实时处理还是用固定经验值稳妥。第三配合MUSIC或ESPRIT估计频率可以突破傅里叶分辨率的限制处理频率非常接近的工程信号。我自己的体会是先把这套三级傅里叶剥离流程作为第一遍快速侦察筛查大多数信号如果发现确实存在频谱上挤在一起的峰再针对性地用子空间方法去精细估计这样多数问题已经能解决一大半。这套方法最初只是我处理振动信号时临时搭的方案B后来在好几个项目里反而成了首选因为它流程简单、参数透明、每一步的频谱图都能直接拿来做人工检查。如果让我总结一条最重要的经验就是永远不要只看最终分解结果要把每一级的残差频谱和提取分量的时域波形打印出来看一眼大部分问题都藏在中间级里。希望这篇实现记录能帮你少走几步弯路如果你在自己的数据上试出了不同的系数规律欢迎一起交流。