
做齿轮箱振动分析的时候我第一次被包络谱救回来。当时直接拿原始振动信号做 FFT频谱图上载波分量占了绝对优势调制边带缩在旁边要么被噪声埋住要么根本分不清主次。真正让我把故障频率从一堆边带里拎出来的就是包络谱分析这条链路目标信号→希尔伯特变换→得到解析信号→求解析信号的模→得到包络信号再对包络信号做 FFT。很多教材只写公式真在 Matlab 里落地的时候会碰到不少细节比如hilbert()返回的到底是解析信号还是纯希尔伯特变换、包络求完之后要不要去直流、采样率和数据长度怎么配合才不会被频率分辨率坑到。这篇就按我自己的实操经验把整条链路从原理讲到代码再给两个验证案例最后聊一下那些脚本里不写、但做项目一定会踩到的坑。这篇内容不只是给故障诊断工程师看的。凡是需要从调制信号里提取低频特征的人比如做通信解调、语音共振峰分析、电力系统谐波检测的都能直接复用这套流程。要不是你只想知道“Matlab 里包络谱怎么算”直接看第 2 节和第 3 节的完整代码就行如果你想把算法真正用准、看谱的时候不被假象带偏那从第 1 节开始读更值。1. 常规频谱为什么会被调制信号绕晕1.1 幅值调制信号在工程里到处都是先建立一个前提工程测到的很多信号本质上都是“载波 调制”的形态。旋转机械里最常见的一类故障比如轴承滚道出现剥落、齿轮齿面出现点蚀每当滚动体或轮齿经过缺陷位置就会产生一个周期性的冲击。这个冲击能量会激励起机械结构本身的固有频率振动相当于给一个高频振荡叠加了一个低频的周期性包络。故障越明显冲击越尖锐包络的调制深度越大。通信里的 AM 调制、DSB 信号语音里声门脉冲序列对声道共振峰的激励本质上也都属于同一类结构。所以把“调幅信号的解调技术”学透适用范围远比想象中广。写成数学表达式就是fs 4096; % 采样率 t (0:4096*2-1)/fs; % 2 秒时长 fm 30; % 调制频率也就是要提取的特征频率 fc 1000; % 载波频率 / 结构共振中心频率 sig (1 0.6*cos(2*pi*fm*t)) .* cos(2*pi*fc*t);这里fc1000相当于高频载波fm30就是我们想从包络里提取的“低频特征”0.6是调幅深度。对这个信号展开一下就变成了三条谱线的组合载波分量幅度 1频率fc上边带幅度 0.3频率fc fm下边带幅度 0.3频率fc - fm为什么这对实际诊断是致命伤因为齿轮箱或者轴承座上往往同时存在多个相互独立的调制源。比如一个齿轮箱里有三对齿轮啮合每个啮合频率都会变成一组载波轴承故障特征频率、转频谐波又会调制到这些载波旁边。所有边带叠加起来直接 FFT 的结果就是一大片密集的“草丛”峰值高度接近间隔混乱很难判断哪个峰才是真正的故障特征频率。即便只有一个载波也存在问题边带幅度只有载波的 β/2一旦调制深度小边带很容易被噪声淹没。1.2 包络谱的思路先解调再谱分析换个角度看这个搅在一起的频率成分如果不去管高频振荡本身而是先把信号的“外壳轮廓”恢复出来也就是提取调制信号那么高频载波带来的边带问题就消失了低频特征直接变成了包络信号里的主分量。包络信号再经过一次 FFT原来的边带对应位置会变成独立的低频峰值这就是包络谱分析的核心逻辑。包络谱相比直接频谱分析的优势非常明显调幅信息从“边带”变成了“主峰”抗噪声和抗泄露能力更强特征频率通常远低于载波频率包络谱只关心低频段同样数据长度下频率分辨率可以做得更高多故障叠加时虽然包络谱会混入多个调制分量但配合带通滤波可以把不同频带激励出来的调制分量分离开。你可能会问提取包络是不是把波形波峰连接起来就行理论方向对实际操作不可靠因为离散采样点不一定每个周期都能恰好采到真正的峰值而且噪声会让波峰点上下抖动。正规做法是构造解析信号通过希尔伯特变换得到瞬时包络。这就是标题里“目标信号→希尔伯特变换→得到解析信号→求解析信号的模→得到包络信号”这条链路的来源。简单说包络谱分析就是“解调 FFT”先解调再把解调结果送入频域分析。2. 希尔伯特变换在 Matlab 里的正确打开方式2.1 90 度相移给信号配一个“同步副本”希尔伯特变换的物理意义不难理解它是一个全通滤波器对所有频率分量的幅度不变但相位移动 90 度。也就是说如果输入信号里有一个余弦分量经过希尔伯特变换后会对应变成正弦如果输入是正弦希尔伯特变换就变成负余弦。连续时间域里的定义是一个带 Cauchy 主值的卷积[ \hat{x}(t) H[x(t)] \frac{1}{\pi} \int_{-\infty}^{\infty} \frac{x(\tau)}{t-\tau} d\tau ]不用在这条积分式上死磕。更直观的方式是把它看作给原信号配了一个“完全同步、只是相位错开 90 度”的副本。有了这个副本就可以构造解析信号[ z(t) x(t) j \cdot H[x(t)] ]这个复数信号的实部是原始信号虚部是希尔伯特变换结果。为什么要构造一个复数因为把信号搬进复平面后可以用旋转向量的模长来表示瞬时幅度也就是包络。这里就要说一个极其常见、我自己也栽过的坑Matlab 的hilbert(x)返回的不是纯希尔伯特变换而是解析信号 z(t)。如果你在命令行运行z hilbert(x);那real(z)等于原始信号ximag(z)才是希尔伯特变换结果。很多人第一次用的时候以为函数返回的只有虚部写了imag(hilbert(x))拿去当解析信号用包络当然算不出来。正确获取包络的标准写法是z hilbert(x); % z 就是解析信号实部为原始信号虚部为希尔伯特变换 envelope abs(z); % 取模得到包络信号hilbert在 Matlab 里基于 FFT 实现先把信号变换到频域屏蔽负频率分量再对正频率分量加倍最后做逆 FFT 得到解析信号。这个过程对窄带信号效果非常好但对所有信号来说边缘样本始终会有一定的端点效应这一点后文单独说。2.2 为什么取模之后还要去直流解析信号的模长abs(z)其实就是瞬时幅度包络。对于幅度调制信号例如我在第 1 节里写的sig (1 0.6*cos(2*pi*fm*t)) .* cos(2*pi*fc*t)理论上包络就是1 0.6*cos(2*pi*fm*t)其中静态的“1”会变成包络信号里的直流分量也就是 0 Hz 处一个很大的谱峰。做包络谱之前如果不先把包络信号的直流分量去掉FFT 之后 0 Hz 这个大峰的能量泄漏会污染旁边的低频区段。去除直流非常简单envelope envelope - mean(envelope);这一步为什么重要看看 FFT 的公式就明白了直流分量的能量全部集中在第 1 个频点但加窗和有限长度数据会导致谱泄漏0 Hz 大峰会向相邻几个频点扩散。故障特征频率往往只有几十赫兹很容易被这个泄漏淹掉。所以只要是工程信号我都会在算包络谱之前做一次去均值。补充一点包络谱幅值的理解包络谱里的峰值幅度表示“调制深度”不是原始信号在载波频率处的幅值。比如sig (1 0.6*cos(2*pi*fm*t)) .* cos(2*pi*fc*t)包络谱在fm处的峰值理论值是 0.6单边。所以你用包络谱看实测信号时如果某个特征频率的幅值随载荷增加明显上涨说明这个频点的调制能量在加强故障冲击更剧烈。3. 从目标信号到包络谱的完整 Matlab 过程3.1 可以直接抄的完整脚本我一般把流程写成一个可复用脚本方便后面换数据路径。以下是仿真信号的完整版fs 4096; % 采样率 dur 2; % 数据时长2 秒 t (0:fs*dur-1)/fs; fc 1000; % 载波频率 fm 30; % 调制频率 beta 0.8; % 调幅深度 % 构造调幅信号并加噪声 x (1 beta*cos(2*pi*fm*t)) .* cos(2*pi*fc*t); x x 0.1*randn(size(x)); % 希尔伯特变换求解析信号 z hilbert(x); % 取模得到包络信号并去除直流 env abs(z); env env - mean(env); % 对包络信号做 FFT N length(env); EnvSpec abs(fft(env)) / N; EnvSpec EnvSpec(1:N/21); EnvSpec(2:end-1) EnvSpec(2:end-1) * 2; % 单边谱修正 f_axis (0:N/2) * fs / N; plot(f_axis, EnvSpec); xlim([0 200]); grid on; xlabel(频率 / Hz); ylabel(幅度);运行之后你会看到 30 Hz 处出现一个非常明显的峰值而原始信号直接 FFT 的频谱却是一堆集中在 1000 Hz 附近的边带。这就直观体现了包络谱的解调能力。关于谱线幅值的处理我再多说两句。abs(fft(env))/N得到的是双边谱双边谱里 30 Hz 的能量会被平分到 30 Hz 和 -30 Hz 两个频点所以单边显示时要把除 DC 和 Nyquist 外的所有频点幅值乘 2。这套规则和普通实信号频谱分析完全一样很多人在这个细节上出错导致幅值总觉得小了一半。3.2 加窗的取舍诊断中要先保证峰位准确如果你发现仿真信号的非整周期截断导致谱峰旁边出现了比较明显的拖尾可以给包络信号加窗比如 Hann 窗win hann(N, periodic); envw env .* win; EnvSpec abs(fft(envw)) / N;注意加窗之后峰值幅度会被压低Hann 窗对单频分量幅值的校正因子大约是 2但实际信号往往包含多个频率分量简单乘 2 可能不准确。我的建议是工程诊断先以频率位置为主加窗能把旁瓣压下去、把峰看准这个收益远大于幅度定量上的那点损失。如果真要做严格的幅值定量分析建议用多周期整周期截断或者做窗函数幅值校正不要盲目套一个固定系数。这段加窗逻辑新手最容易忽略但实际项目里也挺常见。不连续截断导致的谱泄漏常常会把一个不大的特征峰“变”成好几个伪峰让你误判为多个故障源并存。3.3 为什么包络谱里偶尔会出现“多出来的峰”给上面的仿真信号加噪声之后你有时会看到 60 Hz、90 Hz 也出现小峰这不一定代表信号里有三个独立的调制源。原因是包络信号1 beta*cos(2*pi*fm*t)虽然是理想正弦形但真实测量信号里存在噪声、转速波动和非整周期截断这些都会让包络不再是完美正弦于是包络频谱里出现了谐波分量。故障冲击越尖锐包络波形越接近连续尖脉冲序列谐波分量会越丰富频谱上会出现 2 倍频、3 倍频甚至更多倍频的峰。这是包络谱分析里一个重要的经验不要一看到谐波就认为是新故障要结合理论特征频率去核对基频在哪。反过来如果包络谱上有明显的基频加一串谐波往往说明冲击非常集中故障程度比较严重。真实故障不再是纯正弦调幅而是周期性冲击调幅包络的尖峰程度直接反映故障冲击的陡峭程度。4. 两个案例帮你建立判读手感4.1 仿真案例原始频谱和包络谱的对比为了让你更直观感受包络谱的优势我把同一段仿真信号分别用两种方式处理然后看结果对比。分析方法数据处理链路频谱主峰位置特征识别难度直接 FFT原始信号 → 加窗 → FFT1000 Hz 载波 970/1030 Hz 边带边带叠在载波附近容易被噪声掩盖包络谱信号 → Hilbert → abs → 去直流 → FFT30 Hz 处一个主峰特征峰独立低频频段干净识别容易直接 FFT 时如果你想通过观察 970 Hz 和 1030 Hz 边带来判断 30 Hz 调制频率需要足够的频率分辨率而且 1000 Hz 附近一旦还有其他载波边带这几个峰就全糊在一起了。包络谱把 30 Hz 的调制成分直接提取出来特征峰周围几乎没有干扰。这就是包络谱在工程里受欢迎的根本原因把问题从“高频边带分辨”转成“低频主峰识别”。4.2 轴承故障案例特征频率理论值怎么和谱峰对上滚动轴承故障诊断是包络谱的高频应用场景。以深沟球轴承 6205-2RS 为例在 1770 r/min 转速下转频是 29.5 Hz常见故障特征频率可以按下表估算故障位置特征频率系数1770 r/min 下的频率外圈 BPFO3.585 × fr105.8 Hz内圈 BPFI5.415 × fr159.7 Hz滚动体 BSF2.322 × fr68.5 Hz保持架 FTF0.398 × fr11.7 Hz拿到实测数据后先用带通滤波截取结构与轴承共振的频带再做包络谱最后把谱峰位置与上面这些理论值做比对。如果外圈出现缺陷包络谱上会在 105.8 Hz 及其倍频处出现明显峰值内圈缺陷谱峰则会在 159.7 Hz 附近并伴随转频 29.5 Hz 的边带。为什么会这样因为外圈固定安装缺陷位置相对传感器稳定冲击间隔非常规则内圈随轴一起转动缺陷进入载荷区的角度不断变化导致冲击幅值被转频调制包络谱上自然浮出转频边带。这个差异在实际谱图里特别有用。两组特征可能频率很接近但边带的构型不同完全可以作为区分内圈和外圈故障的依据。4.3 多故障叠加时包络谱分不开怎么办工程现场不像仿真轴承和齿轮故障可能同时存在所有冲击都会作用到同一个测点上。如果你直接把原始信号做希尔伯特变换再求包络谱得到的结果是所有调制信号的矢量叠加包络谱会同时出现多个家族的特征峰谱图会再次变成“草丛”。这里有一个减少不确定性的关键步骤带通滤波。在做 Hilbert 变换之前先针对某个共振频带做带通滤波只保留与该故障源相关的共振成分再做包络谱。例如fc_band 1500; % 带通中心频率 bw 1000; % 带宽 [b, a] butter(4, [(fc_band-bw/2)/(fs/2), (fc_bandbw/2)/(fs/2)], bandpass); % 零相位滤波避免相位失真破坏包络形状 x_filtered filtfilt(b, a, x); env abs(hilbert(x_filtered));为什么必须用零相位滤波器而不是普通filter普通 IIR 滤波器会产生相位延迟信号通过滤波器后包络会“扭”变形包络谱峰的位置虽然通常还在但幅值和边带结构可能受影响。filtfilt对信号做正反向两次滤波相位延迟相互抵消包络形状保持正确。这是我在实际处理里反复对比过的结果值得特别注意。5. 实操中踩过的坑边界效应、频率分辨率和滤波频带5.1 端点效应处理不当谱图首尾会“翘起来”基于 FFT 的 Hilbert 变换是对整段数据做周期延拓假设的如果数据首尾不连续解析信号两端会剧烈振荡包络谱两端会出现一大坨虚假能量严重时能盖住几十赫兹内的真实低频特征。我常用的对策有三个数据长度尽量多采集至少包含 30 个以上的目标特征周期后再截取在数据前后各延长 10% 样本做过渡段Hilbert 变换算完再丢弃边缘部分对包络信号加窗把边界的虚假振荡压到不显眼的水平。这三种方法可以叠加使用。尤其第一个对任何故障诊断来说都是基础中的基础数据太短意味着频率分辨率不足后面再怎么花样翻新也救不回来。5.2 采样率和数据长度要同时满足不能只盯采样率包络谱关心的频率往往很低比如保持架故障特征可能只有 11.7 Hz。很多人看到传感器采样率是 10 kHz就觉得分辨率肯定够了其实分辨率的公式是[ \Delta f \frac{f_s}{N} ]如果只采集了 0.1 秒数据N1000fs10000频率分辨率就是 10 Hz11.7 Hz 和 30 Hz 根本分不开谱峰全部糊在一起。提高频率分辨率的有效办法是增加采样时长而不是盲目提高采样率。这也是包络谱分析的隐形门槛你不需要很高的采样率去分辨几十赫兹的低频特征但一定需要足够长的数据窗口。我的经验值是目标峰值之间最小频率间隔除以频率分辨率要大于 3也就是 (\Delta f \Delta f_{min}/3)这样才能把两个相邻峰干净地分开。例如要区分 105.8 Hz 外圈和 159.7 Hz 内圈最小间隔约 54 Hz频率分辨率做到 18 Hz 以内就够但如果要分辨 11.7 Hz 保持架频率和 29.5 Hz 转频间隔不到 18 Hz频率分辨率必须做到 6 Hz 以下也就是数据时长至少要 0.17 秒实际我建议留一倍余量到 0.5 秒以上。5.3 带通滤波频段选错包络谱会变成“噪声大会”带通滤波在包络谱分析里不是随便选一段频率就行的。如果你把带通中心频率放在一个没有共振峰的地方滤波后剩下来的主要是噪声包络谱也只是噪声的包络谱找不出任何特征峰。正确做法是先看原始信号的功率谱或者时频谱找出能量集中、明显抬高的共振带再以这个频带做带通滤波器。选择合适的共振频带时我会按以下顺序排查用pwelch或直接fft看功率谱找到能量集中的频带用spectrogram看时频图确认冲击能量随时间的变化是否规则条件允许时用谱峭度spectral kurtosis选带可以快速定位冲击成分集中的频带选完频带后做包络谱如果特征峰不明显回头重新选带。选带的原则可以通俗理解成既然包络谱是“解调”解调的载波频段就必须是真正承载故障冲击的共振频带。你拿错了载波解调出来的噪音和干扰自然堆满整个谱面。5.4 包络信号本身也可以再做平滑理论推导里解析信号的模已经是最干净的时间包络了但数值实现中由于载波频率不一定是 FFT 频点的整数倍、原始信号带宽不窄等原因包络中还是会残留少量载波附近的高频分量。如果这些残留干扰了低频段的观察可以对包络做一次低通滤波把远高于特征频率的噪声去掉fc_lp 500; % 低通截止按实际特征频率调整 [b_lp, a_lp] butter(4, fc_lp/(fs/2), low); env_smooth filtfilt(b_lp, a_lp, env);这个步骤不是必须的但当我需要从波形上肉眼判断冲击间隔时平滑后的包络会清晰很多尤其是现场数据噪声偏大的场景效果明显。6. 从包络延伸到瞬时频率和其他工程应用解析信号的价值不只是求包络。既然已经得到了 (z(t))你还能从相位里提取瞬时频率inst_phase unwrap(angle(z)); inst_freq diff(inst_phase) / (2*pi) * fs;这个序列代表信号瞬时频率的变化轨迹。在旋转机械里瞬时频率可以用于监测转速波动、识别扭转振动在语音处理里瞬时频率和包络组合起来可以重构语谱图。很多问题如果只看包络等于只用了希尔伯特变换一半的信息。此外如果设备转速不是恒定的直接对包络做频谱分析会出现频率模糊这时可以把包络信号按角度域重采样再做阶比分析。这已经是一套完整的方法论但底层仍然离不开“目标信号→希尔伯特变换→解析信号→包络信号”这条链路。先把这个基础在 Matlab 里完全跑通再去扩展变转速场合的阶比谱会顺畅很多。最后分享一个我自己的真实处理习惯每次拿到振动数据我都先把原始信号、包络信号、包络谱三张图画在同一张 figure 里。原始信号看冲击节拍包络信号看幅值调制轮廓包络谱看特征频率三个视角互为印证。这套流程跑顺了之后故障定位的速度会比单看频谱快上一大截也少了很多误判。