MATLAB希尔伯特变换实现包络谱分析：滚动轴承故障诊断实战指南

简介本资源是一份面向信号处理初学者与工程实践者的MATLAB实战代码包聚焦希尔伯特变换在非平稳信号如机械振动、语音包络谱分析中的核心应用。资源提供完整可运行的MATLAB实现流程从原始信号预处理、hilbert函数构造解析信号、abs提取瞬时包络到FFT获取包络谱并可视化结果覆盖理论落地的关键环节。压缩包共2个文件1个PNG示意图用于结果展示1个.m主程序脚本含详细注释总大小仅1KB轻量易读便于快速理解算法逻辑与代码结构。已有1615人学习下载适合高校学生课程设计、故障诊断入门实践及工程师快速复现包络谱分析流程代码简洁规范无需额外依赖开箱即用是掌握希尔伯特变换工程化应用的高效入门材料。 搞旋转机械故障诊断的人应该都经历过这种尴尬频谱图上转频的谐波一排一排立在那儿旁边围着一圈边带看着热闹可你要找的故障特征频率却像躲猫猫一样怎么也对不上号。尤其滚动轴承早期故障能量都集中在高频共振区低频段那几个毫伏级别的冲击频率早被背景噪声吃了直接在原始频谱上读特征频率基本是碰运气。老师傅丢过来一句“做个包络谱看看”我才算开了窍——原来用MATLAB的希尔伯特变换配合一个取模运算就能把高频调制信号里的低频故障频率干干净净地剥出来。这篇文章就把这一套东西讲透包络谱到底解决什么问题希尔伯特变换和解析信号的数学原理一份可以直接复制运行的MATLAB源程序以及我实际踩过的坑——采样率怎么选、数据截多长、要不要先带通滤波、谱线怎么读。无论你是做毕设、搞故障诊断课题还是在现场处理振动数据都能直接用。1. 包络谱到底在解决什么问题——频谱图上看不见的故障特征1.1 调制故障信息藏在共振频率的“振幅”里先想一个场景。把一辆车开过一段金属检修盖板车轮每转一圈就会“咯噔”一下。这个“咯噔”的节律可能是几十赫兹但车轮和盖板碰撞激起的车身嗡鸣可能是几百甚至上千赫兹。你耳朵听到的是高频的嗡鸣但真正能告诉你“车轮有毛病”的是那个低频的“咯噔咯噔”节律。再把这个场景搬到滚动轴承上轴承某一圈出现点蚀滚动体滚过缺陷时会产生一个冲击脉冲这个冲击的重复频率就是故障特征频率。冲击会激起轴承座、端盖等结构的固有共振于是传感器拾取到的信号变成了一种经典的调幅信号——高频共振频率是“载波”低频故障特征频率是“调制波”。信号写成数学形式就是x(t) a(t) · cos(2π·f_c·t) 其他分量其中 f_c 是结构共振频率a(t) 是一个以故障特征频率为周期的慢变信号。也就是说故障的“线索”不在 f_c 附近而在 a(t) 这个包络里。1.2 为什么原始FFT谱很难直接读出故障频率很多人拿到振动信号第一件事就是做FFT结果频谱图上低频段往往只有转频及其谐波故障特征频率处却看不到明显峰值。原因有两个。第一轴承早期故障产生的冲击能量非常小直接体现在低频段的时候幅度可能比背景噪声还低甚至淹没在转频谐波的旁瓣里。第二故障冲击激起的共振分量在频谱图上表现为一个宽频的“山包”这个山包本身不携带明确的故障频率信息它只是告诉我们“这里有共振”但共振频率谁激起来的、以什么节奏激起来的完全看不出来。于是就有了一个很自然的处理思路先把高频载波去掉把“振幅变化”单独提取出来再对这个振幅变化做频谱分析。这就是包络谱。包络谱的横轴不再是原始信号频率而是“调制频率”——也就是冲击的重复频率。在这个域里故障特征频率会以清晰谱线的形式出现。1.3 包络谱的“解调”思路包络谱本质上是一种解调手段。通信工程里解调是为了从调幅广播信号里还原语音故障诊断里的包络谱则是为了从振动信号里还原冲击的重复节律。两者用的数学工具完全一致核心就是希尔伯特变换。实际工程中包络谱通常和带通滤波配合使用先在共振频带附近做带通滤波把有用的调幅成分单独留出来再求包络然后做FFT。这套流程也叫共振解调是滚动轴承故障诊断最经典的做法。接下来的内容会按这条主线展开。2. 希尔伯特变换与解析信号包络是怎么算出来的2.1 希尔伯特变换的物理含义希尔伯特变换的定义很多但做工程的人只需要记住一个最重要的理解对一个实信号 x(t) 做希尔伯特变换等价于让它的所有频率分量相位偏移 -90 度。一个余弦会变成正弦一个正弦会变成负余弦频率成分的幅值保持不变。用公式表示就是x̂(t) H[x(t)] (1/π) · ∫ x(τ)/(t-τ) dτ这个积分形式看着吓人实际用的时候我们根本不直接算它。在频域希尔伯特变换就是一个简单的滤波器幅频响应恒为1相频响应在正频率是 -90 度在负频率是 90 度。理解到这一层就够了。2.2 从解析信号到包络数学推导与直观理解为什么要搞出个 -90 度相移的副本因为有了它就能构造一个复数信号z(t) x(t) j·x̂(t)这个复数信号叫解析信号。它的实部是原始信号虚部是原始信号的希尔伯特变换。解析信号的一个重要性质是原始信号的“包络”可以直接由它的模得到env(t) |z(t)| sqrt(x²(t) x̂²(t))直观理解是这样假设 x(t) a(t)·cos(2πf_c·t)a(t)变化很慢那么它的希尔伯特变换近似等于 a(t)·sin(2πf_c·t)。于是|z(t)| sqrt(a²(t)·cos² a²(t)·sin²) a(t)瞧三角恒等式把载波消掉了剩下的就是幅度调制函数 a(t)也就是我们想要的包络。这个结论对任意窄带信号都近似成立所以希尔伯特变换是包络解调的标准方案。2.3 MATLAB中hilbert()函数的使用边界MATLAB的Signal Processing Toolbox里提供了hilbert函数但很多人第一次用都会迷路明明函数名叫hilbert为什么返回的结果是复数因为MATLAB的hilbert(x)返回的不是希尔伯特变换的结果而是解析信号本身。换句话说它返回的是 z(t)不是 x̂(t)。取实部是原始信号取虚部才是真正的希尔伯特变换结果。要用它求包络必须配合abs()取模env abs(hilbert(x));这个细节我在不少代码里见过写反的。有人直接plot(imag(hilbert(x)))看着像包络又不像有人plot(real(hilbert(x)))发现跟原信号一模一样。记住hilbert函数返回解析信号包络 abs(解析信号)。另外还有个边界条件hilbert函数要求输入是实信号。如果你不小心传了一个复数进去它会把输入当成解析信号的一部分直接不改变虚部结果会非常怪异。所以使用前最好确认x是实数列向量。2.4 自己动手实现一遍希尔伯特包络如果没有工具箱或者想彻底搞明白原理可以用FFT自己实现。思路是这样的对原始信号做FFT把负频率部分置零正频率部分幅度加倍直流和奈奎斯特频率分量保持不变然后做IFFT得到的复数序列就是解析信号取模就是包络。function env my_hilbert_env(x) % 自实现希尔伯特包络提取 % 输入x为实数列向量输出env为包络信号 N length(x); X fft(x); H zeros(N, 1); if mod(N, 2) 0 % N为偶数 H(1) 1; % 直流分量保持不变 H(2:N/2) 2; % 正频率加倍 H(N/21) 1; % 奈奎斯特频率分量保持不变 else % N为奇数 H(1) 1; H(2:(N1)/2) 2; end z ifft(X .* H); env abs(z); end这段代码和MATLAB内置函数的结果几乎一致。跑一遍这个函数再对比abs(hilbert(x))你就能理解解析信号是怎么构造出来的了。3. 直接可用的MATLAB单信号包络谱源程序3.1 最小可用代码原理说完了直接上代码。下面这个版本没有读取文件而是构造了一段仿真信号来演示完整流程好处是你复制粘贴就能跑不需要额外数据。%% 希尔伯特变换求包络谱——仿真信号演示 clear; clc; close all; % 仿真参数 fs 25600; % 采样频率 25600 Hz N 16 * 1024; % 分析点数 16384 t (0:N-1) / fs; fr 29.5; % 轴频 29.5 Hz BPFO 4.78 * fr; % 外圈故障特征频率约 141 Hz fc 3000; % 共振频带中心频率 3000 Hz % 构造仿真信号转频分量 外圈故障调制分量 随机噪声 x 0.8 * sin(2*pi*fr*t) ... 0.6 * sin(2*pi*BPFO*t) .* sin(2*pi*fc*t) ... 0.05 * randn(1, N); x x(:); % 转成列向量 % 希尔伯特变换求包络 x_analytic hilbert(x); % 解析信号 env abs(x_analytic); % 包络 env env - mean(env); % 去直流否则0Hz处会有一个大尖峰 % 包络谱分析 ENV fft(env); % 对包络做FFT f (0:N/2-1) / N * fs; % 单边频率轴 Amp 2 * abs(ENV(1:N/2)) / N; % 单边幅值谱 % 绘制时域包络和包络谱 figure(Color, w, Position, [100 100 1200 500]); subplot(1, 2, 1); plot(t, x, Color, [0.6 0.6 0.6]); hold on; plot(t, env, r-, LineWidth, 1.2); xlim([0 0.1]); legend({原始信号, 包络}, FontSize, 9); xlabel(时间/s); ylabel(幅值); title(原始信号与包络); subplot(1, 2, 2); plot(f, Amp, b-, LineWidth, 1); xlim([0 500]); xlabel(频率/Hz); ylabel(幅值); title(包络谱); grid on;3.2 关键参数说明先看采样频率fs。这里用的25600 Hz是工业现场比较常见的设置能覆盖大多数机械共振频带。如果你的设备共振频率更高采样率也要跟着提上去保证奈奎斯特频率大于共振频率的2倍以上否则带通滤波和包络解调都会出问题。再看分析点数N。N直接决定频率分辨率分辨率Δf fs / N。这里的N16384对应分辨率约1.56 Hz足够分辨141 Hz和29.5 Hz边带。如果数据量有限可以适当减小N但要注意分辨率会变差。hilbert那句是整段代码的核心。x_analytic是复数解析信号abs之后得到包络再减去均值去除直流。很多初学者在这一步会漏掉去直流导致包络谱0 Hz处的谱线巨大低频段被压制得什么都看不出来。3.3 输出结果如何验证跑完这段代码在时域图里应该能看到红色包络呈现明显的周期性起伏这个起伏的周期对应外圈故障特征频率。在包络谱图里141 Hz附近应该有一根清晰的谱线还可能伴随29.5 Hz间隔的边带这正是仿真里调制的体现。验证代码是否正确的办法很简单把BPFO和fr的值换一换看看谱线是否跟着变。如果谱线位置和设定值对得上说明整个流程没有bug。这个习惯非常有用我建议你先跑一遍仿真再换成自己的数据。4. 批量处理多个数据的工程化代码4.1 工程场景说明现场测试或实验研究很少只测一条信号。最常见的情况是采集了一堆数据文件每个文件对应不同工况、不同测点或不同故障状态需要批量计算包络谱并保存结果。手工一条条导入、分析、导图效率极低还容易出错。下面这段代码就是为这种场景准备的。4.2 批量源码与注释%% 批量计算多个MAT文件的包络谱并保存图片 clear; clc; close all; % 路径与参数设置 filePath ./data; % 数据文件夹路径 outPath ./envelope_results; % 结果输出文件夹 if ~exist(outPath, dir) mkdir(outPath); end fs 25600; % 采样频率按实际修改 maxDuration 10; % 单条信号分析时长单位秒 % 获取所有.mat文件 files dir(fullfile(filePath, *.mat)); for k 1:length(files) % 读取数据 S load(fullfile(filePath, files(k).name)); fn fieldnames(S); % 取第一个变量作为信号 x S.(fn{1})(:); % 如果信号太长截取前maxDuration秒 if length(x) maxDuration * fs x x(1:maxDuration * fs); end N length(x); % 希尔伯特变换求包络 analytic hilbert(x); env abs(analytic); env env - mean(env); % 包络FFT ENV fft(env); f (0:N/2-1) / N * fs; Amp 2 * abs(ENV(1:N/2)) / N; % 绘制包络谱并保存 hFig figure(Visible, off, Color, w, Position, [100 100 900 500]); plot(f, Amp, b-, LineWidth, 1); xlim([0 1000]); xlabel(频率/Hz); ylabel(幅值); title(sprintf(%s 包络谱, files(k).name(1:end-4))); grid on; saveas(hFig, fullfile(outPath, [files(k).name(1:end-4) _envelope.png])); close(hFig); fprintf(已处理: %s\n, files(k).name); end disp(批量包络谱计算完成);4.3 批量处理时的几个注意点读取数据时用fieldnames取了.mat文件里的第一个变量这样即使两个文件的变量名不一致也能处理。但要注意如果.mat文件里存了多个变量务必确认第一个变量确实是振动信号。截取前10秒是经验值。包络谱的频率分辨率是1/TT是分析时长10秒对应0.1 Hz分辨率对绝大多数轴承故障诊断都足够。如果你处理的信号本身很短就不要强行截取直接全段分析即可。保存图片时用了Visible,off这样批处理时不会弹出几百个窗口处理速度也快很多。如果你需要在屏幕上实时查看结果去掉这个参数就行。我建议在批量跑之前先用第3章的单信号代码验证一条数据确认特征谱线清晰、频率轴正确再整批处理。不然几十个文件跑完之后发现参数错了返工成本很高。5. 先带通滤波再做包络谱共振解调的标准操作5.1 为什么滤波之后包络谱更干净直接对原始信号求包络谱在有些场合也能看到特征频率但谱线经常比较脏。原因是原始信号里混杂了大量与故障调制无关的成分——转频谐波、齿轮啮合频率、随机噪声、工频干扰。这些成分在求包络时会被“解调”到低频段变成包络谱里杂乱的背景峰。共振解调的思路是先带通滤波把故障冲击激起的共振频带单独切出来滤掉其他所有成分再对这个窄带信号做包络谱。这样一来包络谱里留下的主要就是故障冲击的调制信息信噪比会大幅提升。这也是为什么现场工程师提取轴承包络谱前几乎都会做带通滤波。5.2 带通滤波器的设计与参数选择滤波频带怎么选最直观的办法是看原始信号的FFT谱找一个能量集中的宽频“山包”这就是共振区。以我的经验这个山包通常在1000 Hz到8000 Hz之间具体位置取决于轴承座结构、传感器安装方式和被测设备刚度。下面是一段完整的滤波包络谱代码%% 带通滤波 希尔伯特包络谱 fs 25600; % 假设原始信号存在变量 x 中 % 设计带通滤波器 f_low 2000; % 带通下限 f_high 4000; % 带通上限 [b, a] butter(4, [f_low/(fs/2), f_high/(fs/2)], bandpass); % 零相位滤波 x_f filtfilt(b, a, x); % 对滤波后信号求包络 analytic hilbert(x_f); env abs(analytic); env env - mean(env); % 包络谱 N length(env); ENV fft(env); f (0:N/2-1) / N * fs; Amp 2 * abs(ENV(1:N/2)) / N; figure(Color, w, Position, [100 100 900 500]); plot(f, Amp, b-, LineWidth, 1); xlim([0 500]); xlabel(频率/Hz); ylabel(幅值); title(带通滤波后的包络谱); grid on;这里用了butter设计4阶带通滤波器配合filtfilt做零相位滤波。零相位意味着滤波不会让波形产生时移和相位畸变对包络的峰值位置和形状影响最小。如果你对相位不敏感也可以用filter但我习惯用filtfilt。f_low和f_high的取值范围需要根据实际频谱调整。如果你不确定共振峰位置可以多试几组频带比如2000-4000、2500-5000、3000-6000对比哪组包络谱的特征频率谱线最突出。实际操作中这是个非常有效的土办法。5.3 滤波前后包络谱对比的经验我拿仿真信号试过滤波后的包络谱比滤波前干净很多。原始信号直接做包络谱时转频分量会被解调出一些杂散的谐波成分这些成分在低频段和特征频率混在一起干扰判断。带通滤波之后转频分量被滤掉了包络谱就只留下载波调制的信息141 Hz处的谱线非常干净。如果是实测信号这种差异会更明显。现场振动信号往往包含大量低频振动和工频干扰这些成分在包络解调后都会在低频段留下痕迹。带通滤波相当于在解调之前做了一道“选区”让后面求出来的包络只反映我们关心的共振频带里的幅值变化。我强烈建议实际项目里都走滤波流程。6. 采样率、窗函数、数据长度包络谱参数与避坑清单6.1 频率分辨率的真正含义包络谱的频率分辨率跟原始FFT完全一样都是Δf fs / N也等于1/T。理解这一点很关键。假设数据只有0.5秒分辨率就是2 Hz两条相差不到2 Hz的谱线根本分不开。而滚动轴承特征频率的边带间隔往往就是转频也就是20到30 Hz的量级2 Hz的分辨率通常够用。但如果转频很低比如5 Hz以下分辨率不够就会导致边带糊成一片无法判断。所以判断数据长度够不够不能只看点数多不多要看总时长T对应的分辨率是否小于你要分辨的最小频率间隔。常见的做法是保证T至少包含20到50个故障冲击周期这样包络谱的谱峰才能稳定。6.2 常见工程坑与处理办法下面这张表是我实际使用中总结的也是很多刚接触包络谱的人反复踩的坑。现象根本原因处理方法包络谱0 Hz处有巨大尖峰包络信号含有直流分量对包络先减均值再FFT特征频率谱线弱低频底噪大信号未经带通滤波先共振带带通滤波再解调特征频率旁边出现莫名杂峰原始信号里有其他强调制源检查是否混入齿轮啮合或工频干扰谱线位置与理论值对不上转频估不准或信号截断非整周期用转速计测准转频增加数据长度图两端出现很高尖峰希尔伯特变换端点效应对包络去掉首尾若干点后再FFT同一组数据处理结果漂移MATLAB版本或工具箱差异确认hilbert来自Signal Processing Toolbox其中端点效应值得多说一句。希尔伯特变换本质上是全信号积分信号不完整时开头和结尾会产生明显畸变。这种畸变在做FFT时会泄漏成低频段的宽峰。解决办法很简单包络算出来之后舍弃开头和结尾各一段数据比如1024点只对中间平稳部分做FFT。6.3 核对特征频率的小技巧算完包络谱第一件事不是急着抄谱线位置而是验算。用轴转速除以60得到转频fr再用轴承参数算出理论特征频率然后到包络谱里找对应谱线。如果谱线落在理论值附近1到2个频率分辨率以内基本可以确认是特征频率。我还习惯把一阶、二阶、三阶特征频率同时标注在图上。轴承故障的特征频率往往有谐波基频弱的时候谐波反而清楚。只盯基频容易漏判把谐波一起看判据就更可靠。7. 包络谱读谱实战从频率尖峰到故障定位7.1 滚动轴承特征频率公式包络谱的谱线只有翻译成物理意义才算真正有用。滚动轴承四大部件的故障特征频率公式如下故障位置特征频率外圈BPFO n/2 · fr · (1 - d/D · cosα)内圈BPFI n/2 · fr · (1 d/D · cosα)滚动体BSF D/(2d) · fr · (1 - (d/D · cosα)²)保持架FTF fr/2 · (1 - d/D · cosα)其中n是滚动体数量d是滚动体直径D是轴承节径α是接触角fr是轴转频。这些参数通常可以在轴承型号手册或厂家资料里查到。算出来之后把包络谱里的峰值谱线和这几个理论值对照。比如外圈特征频率大约是3.2倍的转频那就在3.2·fr附近找谱线。如果找到的谱线正好落在BPFO位置且和理论值的偏差在2 Hz以内外圈故障的可能性就很大。7.2 边带的含义与故障程度判断包络谱里除了特征频率本身边带的形态也很有价值。内圈故障的调制过程往往混入转频成分所以BPFI周围会出现fr间隔的边带。外圈故障的调制源相对固定边带通常没有内圈明显。滚动体故障则可能在BSF周围出现保持架转频FTF间隔的边带。如果特征频率周围几乎没有边带说明调制很纯故障冲击的重复性好如果边带丰富、谱峰平坦说明转速波动或冲击能量不稳定。结合边带形态和特征频率幅值随时间的增长趋势可以粗略判断故障是否在扩展。当然精确的故障程度评估还需要结合加速度峰值、峭度等多个指标不止包络谱一个维度。7.3 一个建议的习惯在实际处理数据时我的流程是固定的先看时域波形有没有周期性冲击再做原始频谱确定共振频带然后带通滤波再做希尔伯特包络谱最后将理论特征频率标注在谱图上核对。这个方法的好处是每个环节都能交叉验证。时域冲击能对得上共振频带能对得上包络谱特征频率还能对得上三条线互相印证误判的概率就很小。反之如果只算一个包络谱就下结论万一谱线位置理解错了后面全盘皆输。按照这个流程走下来一套可靠的包络谱分析流程就建立了。我在实际使用中有一个习惯先跑通仿真信号确认算法没问题再上实测数据实测数据如果谱线不够清楚先检查数据长度和共振频带的选择不要急着改算法。包络谱不是越复杂越好而是要让物理意义清晰的谱线说话。你拿自己的数据跑出来的谱线如果有疑问建议回到前面几个环节逐项排查大部分问题都出在参数设置上算法本身反而不容易出错。本文还有配套的精品资源点击获取