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/2+1) = 1; % 奈奎斯特频率分量保持不变 else % N为奇数 H(1) = 1; H(2:(N+1)/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。这里的N=16384，对应分辨率约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/T，T是分析时长，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 常见工程坑与处理办法

下面这张表是我实际使用中总结的，也是很多刚接触包络谱的人反复踩的坑。

其中端点效应值得多说一句。希尔伯特变换本质上是全信号积分，信号不完整时开头和结尾会产生明显畸变。这种畸变在做FFT时会泄漏成低频段的宽峰。解决办法很简单：包络算出来之后，舍弃开头和结尾各一段数据，比如1024点，只对中间平稳部分做FFT。

6.3 核对特征频率的小技巧

算完包络谱，第一件事不是急着抄谱线位置，而是验算。用轴转速除以60得到转频fr，再用轴承参数算出理论特征频率，然后到包络谱里找对应谱线。如果谱线落在理论值附近1到2个频率分辨率以内，基本可以确认是特征频率。

我还习惯把一阶、二阶、三阶特征频率同时标注在图上。轴承故障的特征频率往往有谐波，基频弱的时候谐波反而清楚。只盯基频容易漏判，把谐波一起看，判据就更可靠。

7. 包络谱读谱实战：从频率尖峰到故障定位

7.1 滚动轴承特征频率公式

包络谱的谱线只有翻译成物理意义才算真正有用。滚动轴承四大部件的故障特征频率公式如下：

其中n是滚动体数量，d是滚动体直径，D是轴承节径，α是接触角，fr是轴转频。这些参数通常可以在轴承型号手册或厂家资料里查到。

算出来之后，把包络谱里的峰值谱线和这几个理论值对照。比如外圈特征频率大约是3.2倍的转频，那就在3.2·fr附近找谱线。如果找到的谱线正好落在BPFO位置，且和理论值的偏差在2 Hz以内，外圈故障的可能性就很大。

7.2 边带的含义与故障程度判断

包络谱里除了特征频率本身，边带的形态也很有价值。内圈故障的调制过程往往混入转频成分，所以BPFI周围会出现fr间隔的边带。外圈故障的调制源相对固定，边带通常没有内圈明显。滚动体故障则可能在BSF周围出现保持架转频FTF间隔的边带。

如果特征频率周围几乎没有边带，说明调制很纯，故障冲击的重复性好；如果边带丰富、谱峰平坦，说明转速波动或冲击能量不稳定。结合边带形态和特征频率幅值随时间的增长趋势，可以粗略判断故障是否在扩展。当然，精确的故障程度评估还需要结合加速度峰值、峭度等多个指标，不止包络谱一个维度。

7.3 一个建议的习惯

在实际处理数据时，我的流程是固定的：先看时域波形有没有周期性冲击，再做原始频谱确定共振频带，然后带通滤波，再做希尔伯特包络谱，最后将理论特征频率标注在谱图上核对。

这个方法的好处是每个环节都能交叉验证。时域冲击能对得上，共振频带能对得上，包络谱特征频率还能对得上，三条线互相印证，误判的概率就很小。反之，如果只算一个包络谱就下结论，万一谱线位置理解错了，后面全盘皆输。

按照这个流程走下来，一套可靠的包络谱分析流程就建立了。我在实际使用中有一个习惯：先跑通仿真信号确认算法没问题，再上实测数据；实测数据如果谱线不够清楚，先检查数据长度和共振频带的选择，不要急着改算法。包络谱不是越复杂越好，而是要让物理意义清晰的谱线说话。你拿自己的数据跑出来的谱线如果有疑问，建议回到前面几个环节逐项排查，大部分问题都出在参数设置上，算法本身反而不容易出错。

本文还有配套的精品资源，点击获取