简介本资源是一份面向信号处理初学者与工程实践者的MATLAB轻量级教学代码包聚焦希尔伯特变换在非平稳信号如机械振动、语音包络谱分析中的核心应用。资源共2个文件1个MATLAB源程序.m完整实现数据预处理、hilbert()变换、abs()求包络、FFT计算包络谱及可视化全流程1张PNG示意图直观展示关键结果对比便于理解瞬时幅度提取原理。压缩包仅1KB结构精炼无冗余适合嵌入课程实验、故障诊断项目或自学验证。已有1615人学习下载代码注释清晰、步骤对应理论严谨既可作为入门者掌握解析信号构建与包络解调的实操范例也便于工程师快速复用核心逻辑至轴承故障特征提取等实际场景。 做信号处理的朋友应该都有过这种经历手里拿到一段振动信号频谱图上密密麻麻全是峰值但你想看的低频调制成分却根本看不出来。比如滚动轴承出现局部缺陷时故障特征频率往往以调制边频带的形式出现在高频固有共振附近如果直接对原始信号做FFT低频段的故障频率会被“藏”在载波旁边根本没法直接读出来。这时候就该上希尔伯特变换求包络谱了。这篇文章围绕“MATLAB希尔伯特变换求包络谱源程序代码”这个主题把从原理到代码、从参数设置到避坑经验的完整链路都过一遍。适合正在做机械故障诊断、振动信号分析、通信解调或生物医学信号处理的朋友尤其是那些已经会基本FFT分析、但面对调制信号不知道如何下手的入门到进阶用户。我会直接给出可运行的MATLAB源码逐段拆解每一行的作用并解释采样率、滤波带宽、FFT点数这些关键参数到底怎么定最后再补充几个我实际调试中踩过的坑。1. 包络谱到底在做什么一个必须理解的前提1.1 为什么直接做FFT看不到调制信息先想一个问题一个1000Hz的正弦载波被60Hz的低频信号调制得到的信号在频谱上会长什么样学过通信原理的都知道频谱上出现的不是1000Hz一根谱线而是1000Hz两侧各出现一根边频分别是940Hz和1060Hz。如果调制信号不是一个单频而是一段频率范围那频谱上就会出现以1000Hz为中心的一簇边频带。问题就在这里。机械故障信号里滚动轴承外圈故障特征频率BPFO一般在几十赫兹到几百赫兹能量很小而它调制的高频固有振动能量很大。如果直接做FFT故障信息在频谱图上只是载频附近的两根小边带很容易被噪声淹没肉眼根本分辨不出来。这时候包络谱就派上用场了——它把信号“解调”回基带把那些藏在边带里的低频成分重新拉回到频谱的低频段让故障特征频率直接暴露在谱图上。1.2 希尔伯特变换为什么是求包络的标准方法求包络不是只有希尔伯特变换一种办法比如二极管检波、绝对值整流加低通滤波也能得到类似效果。但这些方法要么对信噪比要求高要么会引入非线性失真。希尔伯特变换的好处在于它是在频域上精确构造解析信号从数学上给出了“包络”这个概念的严格定义。具体来说对一个实信号x(t)做希尔伯特变换得到x̂(t)再构造复信号z(t)x(t)j·x̂(t)这个z(t)叫做解析信号。解析信号在复平面上的幅值|z(t)|就是信号的包络相位角就是瞬时相位对瞬时相位求导还能得到瞬时频率。所以希尔伯特变换一次给了三个东西包络、瞬时相位、瞬时频率。我们在故障诊断里最常用的是包络也就是|z(t)|在通信解调里瞬时相位和瞬时频率又特别重要。MATLAB里的hilbert函数返回的就是完整解析信号注意它返回的是一整个复数数组不是直接返回包络。很多人第一次用的时候直接plot(hilbert(x))画出来发现是一堆乱七八糟的复数然后懵了。实际上你要做的是先调用hilbert得到z再用abs(z)才是包络再用angle(z)才是瞬时相位。这个流程说起来简单但我在不少咨询里看到有人卡在这一步上。1.3 适用场景与不适用场景包络谱最经典的场景是滚动轴承故障诊断。滚动体滚过局部缺陷时会产生周期冲击激起轴承结构的高频固有振动形成“高频载波低频重复频率”的调制信号。包络谱能把低频重复频率提取出来对应BPFO、BPFI、BSF、FTF这些特征频率。齿轮箱故障也常用类似思路齿面剥落、断齿都会产生调制成分。通信领域更不用说AM解调最直接的方法就是希尔伯特变换求包络。但要注意通信里的AM信号包络和机械故障信号包络有区别AM信号调制深度不大、波形规则而机械冲击信号是非平稳的、带噪的实际处理起来对参数设置更敏感。包络谱也有不适用的时候。如果信号本身没有调制结构比如纯正弦叠加或纯随机噪声包络谱反而会把高频能量“搬”到低频段造成假峰。另外如果目标故障频率低于采样频率分辨率所对应的最小频率间隔分不清也正常。所以包络谱不是万能的它解决的是“调制信号解调”这一类问题先用频谱图判断有没有调制现象再决定要不要做包络谱这个逻辑顺序不能乱。2. 从原理到代码希尔伯特包络谱的完整实现2.1 先构造一个仿真信号验证思路写代码之前我习惯先做一个仿真信号把已知的调制频率注入进去然后把完整的解调流程跑通确认每个环节输出正确了再拿去处理真实数据。这个习惯帮我省了很多调试时间。这里我们构造一个采样率fs5000Hz、时长1秒的信号载波频率1000Hz调制信号60Hz调制深度0.5再加入一点随机噪声。MATLAB代码很简单fs 5000; % 采样率 5000 Hz t (0:fs-1)/fs; % 1 秒时间轴 fc 1000; % 载波频率 1000 Hz fm 60; % 调制频率 60 Hz carrier sin(2*pi*fc*t); % 载波 modulating 1 0.5*sin(2*pi*fm*t); % 调制信号直流偏置1保证包络恒正 y modulating .* carrier 0.1*randn(size(t)); % 调制加噪声注意modulating里加的“1”很重要。如果调制信号是纯交流、没有直流偏置包络会在零点上下波动取绝对值后会多出两倍频分量。加1是为了让包络始终为正这样后续做FFT时才能干净地解出60Hz这个频率。调制深度0.5表示包络在0.5到1.5之间摆动不会过调制。2.2 带通滤波先把干扰剔掉再解调真实信号直接做希尔伯特变换不是不行但效果往往很差。原因在于希尔伯特变换是全局操作信号中所有频率分量都会参与构造解析信号如果存在明显的带外噪声或无关频率强峰包络里会混入大量干扰成分。所以标准流程是先带通滤波把载波及其边带所在的频率范围提出来再做解调。载波在1000Hz调制频率60Hz所以边带在940Hz和1060Hz。滤波器通带设在800Hz到1200Hz足够覆盖同时滤掉低频和高频噪声。这里我推荐用filtfilt做零相位滤波虽然它是MATLAB的Signal Processing Toolbox函数但用起来真的很省心——它前后各做一次滤波相位偏移相互抵消不产生信号畸变。% 设计带通滤波器800-1200 Hz巴特沃斯4阶 [b, a] butter(4, [800 1200]/(fs/2), bandpass); y_filtered filtfilt(b, a, y);有人会问为什么不用filter而用filtfilt。filter是普通IIR滤波会引入非线性相位导致包络波形前后错位filtfilt虽然计算量翻倍但零相位特性对包络分析特别关键。尤其后面还要对包络做FFT任何相位畸变都会直接影响谱峰的位置和幅值。2.3 核心三行代码希尔伯特、取包络、FFT接下来是最核心的三行analytic_signal hilbert(y_filtered); % 解析信号 envelope abs(analytic_signal); % 包络 envelope_spectrum abs(fft(envelope - mean(envelope))); % 包络减直流再FFT这里有个细节很多新手会忽略先对包络去均值再做FFT。原因很简单包络信号里有一个很大的直流分量包络的正向偏置如果不去掉FFT结果在0Hz处会出现一个巨大的谱峰旁边频段的细节全被压低。减去均值之后直流分量被归零60Hz的谱峰才能以正常的幅度显现出来。频率轴计算也不能出错。包络信号的采样率还是fs5000HzFFT点数N5000频率分辨率fs/N1Hz。如果只看单边谱只取前N/2个点N length(envelope); f (0:N/2-1) * fs / N; single_side envelope_spectrum(1:N/2) * 2 / N;画出来之后在f60Hz处应该能看到一个明显的谱峰这就是被解调出来的调制频率。实测频谱图上940Hz和1060Hz两根边带对应到包络谱就是60Hz一根线这说明解调成功。2.4 完整可运行的源码把上面所有代码整合到一起就是一个完整的包络谱分析脚本%% 希尔伯特变换求包络谱示例 clear; clc; close all; %% 1. 生成仿真信号 fs 5000; t (0:fs-1)/fs; fc 1000; fm 60; carrier sin(2*pi*fc*t); modulating 1 0.5*sin(2*pi*fm*t); y modulating .* carrier 0.1*randn(size(t)); %% 2. 带通滤波 [b, a] butter(4, [800 1200]/(fs/2), bandpass); y_filtered filtfilt(b, a, y); %% 3. 希尔伯特变换求包络谱 analytic_signal hilbert(y_filtered); envelope abs(analytic_signal); envelope_spectrum abs(fft(envelope - mean(envelope))); N length(envelope); f (0:N/2-1) * fs / N; single_side envelope_spectrum(1:N/2) * 2 / N; %% 4. 画图对比 figure; subplot(3,1,1); plot(t, y); xlim([0 0.05]); title(原始信号); subplot(3,1,2); plot(t, y_filtered, b); hold on; plot(t, envelope, r, LineWidth, 1.5); xlim([0 0.05]); title(滤波后信号与包络); subplot(3,1,3); plot(f, single_side); xlim([0 200]); title(包络谱); xlabel(频率 (Hz)); ylabel(幅值);运行这个脚本第一张图能看到高频载波被低频包络调制的效果第二张图红线勾勒出包络形状第三张图在60Hz处出现清晰谱峰。整个过程跑通之后你可以把fm换成别的频率比如25Hz、120Hz看看包络谱峰值是否跟随变化。3. 高频参数的决定细节采样率、滤波带宽与FFT点数3.1 采样率不能只看最高频率还要看你要解调的包络带宽选择采样率的时候有个常见的认知误区有人觉得载波是1000Hz那采样率只要大于2000Hz就够了因为奈奎斯特定理就是这么说的。这个理解本身没错但忽略了包络分析的特殊性。包络信号本身频率很低比如60Hz但它必须在时域上由原始高频信号经过计算得到采样率太低会导致每个载波周期只有少数几个采样点包络的幅值估计精度会严重下降。想象一下一个正弦波一个周期只采2个点你怎么准确判断它的幅值变化所以工程实践中一般要求载波频率至少被10到20个点覆盖。对1000Hz载波采样率至少要5000到10000Hz才比较稳。另外带通滤波器的过渡带也需要采样率支撑。如果采样率刚好卡在奈奎斯特频率附近滤波器没有空间设计过渡带会和镜像频率发生混叠。我的经验是采样率设置成目标载波频率的5到10倍既保证包络精度又不至于数据量太大。3.2 滤波带宽宽了不行窄了更不行带通滤波器的带宽选取是包络谱分析里最需要“感觉”的参数。带宽太宽把无关频率也放进来了包络里会混入背景噪声带宽太窄可能把边带滤掉一部分导致解调出来的包络幅值偏小、频率成分失真。实际工程中如果已经知道载波频率fc和可能的最大调制频率fm_max那滤波器带宽至少要覆盖fc ± fm_max这个范围。拿轴承诊断举例外圈故障特征频率如果是107Hz再考虑到可能存在2倍频、3倍频调制频率范围可能到300Hz甚至更高那滤波带宽就应该是fc ± 300Hz以上才不至于丢信息。滤波器阶数也要注意。阶数越高过渡带越窄、截止特性越陡峭但相位畸变越大。我一般用butter(4, ...)或butter(6, ...)作为折中四阶巴特沃斯在实际信号处理里已经很常用。如果用的是filtfilt4阶的等效滤波效果实际上相当于8阶但零相位足够应付大多数场景。3.3 FFT点数与频率分辨率想分辨多细取决于你采了多长包络谱的频率分辨率由fs/N决定N是参与FFT的数据长度。如果采集了T秒信号Nfs·T分辨率就等于1/T。想看60Hz旁边的59Hz和61Hz两条谱线分辨率至少要1Hz那数据长度至少要1秒。想看0.5Hz的频率间隔就得采2秒数据。这个逻辑在仿真信号里不明显因为数据随便生成多长都行但真实采集场景里数据长度受硬件、存储和实时性限制。如果采集时间不够可以通过补零来让FFT曲线更平滑但注意补零不能提高真实频率分辨率——它只是插值不会让你分辨出原本分不开的两个频率峰。我调试时遇到过一种情况包络谱上某个峰看着像70Hz但旁边还有一个比较宽的凸起一开始以为是两个频率后来把数据长度从1秒增加到4秒才发现原来是一个单峰被矩形窗的旁瓣“拉宽”了。数据长度足够之后谱线变得尖锐问题自然消失。所以遇到谱线过宽、无法辨认的现象第一反应不是改算法而是加长采样时间。4. 实操避坑包络谱最容易翻车的几个细节4.1 0Hz处巨大谱峰去直流是必须做的前面代码里有一行envelope - mean(envelope)这行看似不起眼实际上决定了包络谱的图面质量。如果不做这一步0Hz峰值的幅值可能是60Hz峰值的几十倍图面自动缩放后整个低频段都被“压”成一条直线。你可能还会想“是不是调制频率太小了”其实只是直流没去掉。真实信号里包络的直流分量不仅仅来自调制信号的直流偏置还包括传感器自身偏置、信号调理电路引入的基线漂移。所以就算原始信号可能没有明显偏置包络做出之后也建议主动减一次均值。这是包络谱分析的基本卫生习惯。4.2 包络波形首尾翘起边界效应要留个心眼filtfilt的零相位滤波在大多数情况下是好事但它会带来一个副作用——信号首尾的瞬态响应比较明显滤波后的波形开头和结尾会有一段不真实的“翘起”或“塌陷”。如果这段瞬态进入后续的FFT会在包络谱中产生虚假成分。处理办法有两种。一种简单粗暴但有效截掉首尾各0.1到0.2秒的数据再做包络谱。另一种是先用适当长度的信号做滤波然后在边缘处加窗平滑但操作起来要细心窗函数选不好反而引入新问题。我个人更偏好截断因为简单可控。代码示例start_idx round(0.1*fs) 1; end_idx round(0.9*fs); y_truncated y_filtered(start_idx:end_idx); envelope abs(hilbert(y_truncated));4.3 矩形窗泄漏包络谱上加窗的意义对包络做FFT时不显式加窗其实等于用了矩形窗。矩形窗的旁瓣泄漏是最严重的如果你的目标谱峰旁边有强干扰旁瓣会拖尾把旁边的小谱峰淹没。特别是在轴承故障诊断中故障特征频率的倍频一次比一次弱如果没有加窗抑制旁瓣第二倍频、第三倍频很容易被第一个峰的旁瓣掩盖。我推荐对包络信号做FFT前加一叶汉宁窗Hann窗。汉宁窗旁瓣衰减比矩形窗好很多主瓣稍微宽一点但损失可以接受。唯一要注意的是加窗后信号总能量变了如果要比较绝对幅值需要做幅度修正。幅值修正系数大约为2因为汉宁窗的相干增益在0.5附近。win hann(N, periodic); envelope_windowed (envelope - mean(envelope)) .* win; envelope_spectrum abs(fft(envelope_windowed)); single_side envelope_spectrum(1:N/2) * 4 / N; % 2/0.54这个4的来源是FFT单边谱本来要乘2加上汉宁窗能量增益约0.5所以要乘4才能恢复真实幅值。注意应用场景如果你只关心频率位置和相对幅值不关心绝对幅值可以不做修正但做绝对标定的时候一定不能漏掉这个系数。4.4 包络谱怎么看只在低频段找峰包络谱分析完成后重点关注的频段一般是0到几百赫兹。载波频率本身在包络谱中应该几乎消失如果能完全解调的话大部分能量都集中在调制频率及其倍频上。所以画图时不要一上来就把整个频段都展示出来直接把xlim限制在0到500Hz或0到1000Hz信息密度高得多。如果你发现包络谱在载波频率附近还有一个大峰说明滤波没有把载波边带滤干净或者解调不彻底。这时候回头检查带通滤波器参数看看是不是带宽设置太宽了。5. 一个轴承故障仿真案例从信号到诊断的完整闭环5.1 故障特征频率的计算滚动轴承外圈故障特征频率BPFO的近似计算公式为BPFO ≈ n/2 × fr × (1 - d/D × cosα)其中n是滚动体数量fr是转频d是滚动体直径D是节圆直径α是接触角。实际应用里很多人直接取近似值就可以了因为转速波动和轴承尺寸误差本来就会带来几个赫兹的偏差。拿一个典型参数举例滚动体个数n9转频fr30Hz相当于1800rpmd/D0.2接触角很小可以近似cosα≈1。那么BPFO≈9/2×30×(1-0.2)108Hz。5.2 故障信号建模与包络谱分析仿真轴承外圈故障信号思路和前面一致高频固有振动作为载波频率取3000Hz外圈故障脉冲的重复频率作为调制频率108Hz。同时可以加入200Hz、300Hz的微弱倍频模拟实际故障信号中的谐波成分。采样率取10000Hz时长取1秒保证频率分辨率到1Hz。完整代码是clear; clc; close all; fs 10000; t (0:fs-1)/fs; fr 30; % 转频 30 Hz bpfo 9/2*fr*0.8; % 约 108 Hz f_carrier 3000; % 轴承固有频率 impact 1 0.8*sin(2*pi*bpfo*t) 0.3*sin(2*pi*2*bpfo*t) 0.15*sin(2*pi*3*bpfo*t); y impact .* sin(2*pi*f_carrier*t) 0.2*randn(size(t)); % 带通滤波中心频率 3000 Hz带宽 ±400 Hz [b, a] butter(4, [2600 3400]/(fs/2), bandpass); y_filtered filtfilt(b, a, y); % 包络谱 start_idx round(0.1*fs) 1; end_idx round(0.9*fs); y_truncated y_filtered(start_idx:end_idx); envelope abs(hilbert(y_truncated)); envelope envelope - mean(envelope); N length(envelope); win hann(N, periodic); envelope_spectrum abs(fft(envelope .* win)); f (0:N/2-1) * fs / N; single_side envelope_spectrum(1:N/2) * 4 / N; figure; subplot(2,1,1); plot((start_idx:end_idx)/fs, envelope); title(包络波形); xlabel(时间 (s)); subplot(2,1,2); plot(f, single_side); xlim([0 500]); xlabel(频率 (Hz)); ylabel(幅值); title(包络谱);运行后的包络谱应该能看到108Hz附近的主峰以及216Hz、324Hz附近逐渐衰减的倍频峰。这个谱图特征和真实轴承外圈故障的包络谱形态非常相似所以用来理解理论到实践的映射很合适。5.3 频谱图与包络谱图为什么包络谱是“放大镜”如果你把原始信号的FFT结果画出来在3000Hz附近会看到一簇密集的边频带但很难分清边带的间隔到底是108Hz还是别的值。而包络谱把高频边带“折叠”回低频段后108Hz主峰直接耸立在低频段上旁边倍频排列清晰诊断结论一目了然。这就是包络谱的意义它不是替代频谱分析而是在频谱分析基础上做二次提取。遇到疑似调制信号的问题先频谱图粗看有没有边带再包络谱精读边带间隔。6. 常见问题排查实录与参数速查6.1 常见问题排查表下面这张表是我在实际调试中频繁遇到的情况整理成速查形式方便你遇到问题的时候对照。现象可能原因解决办法包络谱0Hz处巨大谱峰包络未去直流对包络信号减均值后再FFT包络谱在载波频率附近仍有大峰带通滤波器带宽过宽或滤波无效收窄滤波带宽重新滤波包络波形首尾明显异常filtfilt边界效应截掉首尾0.1~0.2秒数据目标谱峰不尖锐像宽包数据长度太短分辨率不够加长采样时间或检查是否加了不合适的窗包络谱出现两个相邻相等峰可能是载波泄漏或真实频率成分对比不同数据段确认希尔伯特变换结果无法plot忘记取abs或angle记住hilbert返回复数解析信号包络谱幅值异常偏大或偏小加窗后没做幅度修正加汉宁窗后修正系数4matlab单边谱窗增益最后一个问题多说一句用hilbert函数时很多人第一次运行会对结果直接plot画出来是复数不能直接显示然后卡住。其实只需要记住hilbert是复杂数构造工具后续取abs、取angle各取所需就不会在这个基础问题上浪费时间。6.2 参数选择的经验速查表参数建议范围/原则备注采样率fs载波频率的5~10倍保证包络幅值精度带通滤波器通带载波频率±预期最大调制频率边带覆盖要完整滤波器阶数4~6阶巴特沃斯阶数过高会引入相位失真数据长度T至少1/期望频率分辨率想分0.5Hz就采2秒FFT点数N通常等于数据长度补零不能提高真实分辨率窗函数汉宁窗旁瓣抑制好需幅度修正6.3 关于MATLAB版本与工具箱的说明上面代码用到了butter、filtfilt、hilbert和hann这些都来自Signal Processing Toolbox。如果你的MATLAB没有装这个工具箱运行会直接报错。没有工具箱的话可以用最原始的fft和ifft自己实现希尔伯特变换function analytic my_hilbert_analytic(x) n length(x); X fft(x); h zeros(n, 1); if mod(n, 2) 0 h([1, n/21]) 1; h(2:n/2) 2; else h(1) 1; h(2:(n1)/2) 2; end analytic ifft(X .* h); end这段代码的原理是构造一个频域掩蔽数组h把原始频谱的负频率部分置零正频率部分加倍再反变换得到解析信号。和MATLAB内置hilbert函数的功能完全一致也能用于教学理解。当然能用内置函数还是用内置函数性能和边界处理都更可靠。最后再分享一个小技巧我每次做包络谱分析都会同时把原始信号频谱和包络谱画在同一张图上对比即使最终只需要包络谱。这个习惯帮我在很多项目里快速定位问题——如果原始频谱上载波频率周边有明显的边带结构但包络谱里对应频率没出现谱峰那八成是滤波带宽或滤波器参数有问题反过来如果原始频谱上看着很干净包络谱却冒出一个莫名其妙的峰那就要怀疑是不是把噪声当信号解调了。包络谱分析这件事说难不难说简单也不简单。难在它对参数和细节敏感一个滤波带宽没调好、一个直流没去掉结果就可能完全跑偏。我的经验是先把仿真信号跑通理解每条代码在做什么再带着清楚的目的去处理真实信号遇到问题从“采样率-滤波带宽-数据长度-加窗”这套参数链路里逐项排查基本都能定位到原因。希望这篇文章能帮你少走一些弯路。本文还有配套的精品资源点击获取