企业官网建设流程全解析

简介：中央大学数据研究中心提供的信号处理算法源码包，围绕经验模态分解与希尔伯特-黄变换两大核心方法展开，并同时提供标准分解与集合分解两种算法版本，主要面向从事振动分析、故障诊断、生物医学信号处理等领域研究的高校师生和科研工程师，也适合作为相关课程与论文实验的对照实现。压缩包总计45个文件，容量仅188KB，以41个函数脚本为主体，分别承担极值点检测、包络拟合、希尔伯特谱计算等环节；另附2个数据文件用于测试运行，1个加密函数及1个说明文档便于辅助使用。迄今已有2199人学习下载。代码模块划分清晰，入口脚本可直接运行，配合示例数据能够快速复现完整分解流程，也便于对照其他开源实现版本展开算法对比与二次开发。

1. 这套EMD-HHT的MATLAB源码，为什么值得花半小时跑一遍

拿到这套流传很广的“EMD+HHT Matlab源代码”，第一步不是改代码，而是找一份非平稳信号把它完整跑一遍，比如低速轴承故障数据或某段结构振动实测序列。EMD（经验模态分解）不预设备基函数，而是从数据自身一层层抽出固有模态函数（IMF），再对每个IMF做希尔伯特变换，得到瞬时频率与幅值，这就是HHT的全过程。它的直接价值是把非平稳、非线性信号里“什么时候、以多大能量、在哪个频率附近振动”完整拆出来，这在故障诊断、生物医学数据处理和结构健康监测里比FFT更贴近物理真实。适合手里有实测信号、想做时频特征提取的工程师，也适合正在复现论文里EMD/HHT对比实验的研究者。

2. EMD和HHT的底层逻辑：从筛分迭代到瞬时频率

2.1 EMD的筛分过程：包络、均值、IMF三个概念

FFT把信号看成若干固定频率正弦分量的叠加，这个假设对平稳信号成立，但现场采集的振动和生理信号往往是“频率随时间漂移”的。EMD的思路完全反过来：它不要求预先给定基函数，而是用信号自身的极值点构造上下包络，从数据中一层层提取波动模式，这个模式就是IMF（固有模态函数）。IMF必须满足两个条件：整个时间范围内过零点次数与极值点次数之差最多不超过1；局部上下包络的均值趋近于零。

实际筛分过程是一个反复迭代的剥层操作，常见实现是这样的：

1. 找出原始信号 x(t) 的全部局部极大值与局部极小值；
2. 用三次样条函数分别插值出上包络线与下包络线；
3. 两条包络求平均，得到均值线 m1(t)；
4. 原始信号减去均值线，得到候选分量 h1(t) = x(t) - m1(t)；
5. 检查 h1(t) 是否满足IMF的两个条件，不满足就以 h1(t) 为新输入，重复第1到第4步，直到满足条件，得到 imf1；
6. 残差 r(t) = x(t) - imf1，再对残差重复第1到第5步，得到 imf2、imf3，直到残差变成单调序列或幅度小于设定阈值。

这个内层循环称为“筛分”（sifting）。停止筛选时用到两个控制参数：标准差阈值SD和最大迭代次数。经验上SD取0.05到0.2之间比较稳妥，设太小会让筛分次数快速增长，设太大会提前结束导致IMF不纯；最大迭代次数常见默认在几百到几千之间，防止某些端点畸变信号把循环拖死。

工程选型方面，既然有小波变换为什么还要用EMD？很多做故障诊断的人对比过：小波需要人为选择小波基函数和分解层数，换一种信号类型往往就得重新调参。EMD只有停止条件、端点处理两处可调，对不熟悉的信号首次跑出来的结果通常也能解释通。它属于“数据自适应型分解”，好处是比较通用，代价是迭代多、速度比FFT慢，而且数学上的严格收敛性到现在还在不断补证明，工程上重点是看实际分解效果。

2.2 HHT谱的构成：瞬时幅值与瞬时频率的由来

分解只是前半段。HHT的第二步是对每个IMF做希尔伯特变换，把实信号变成复解析信号。希尔伯特变换可以理解为对信号做90度相移，公式是卷积形式，MATLAB里直接用hilbert函数一步完成。

每个IMF对应的解析信号写成：

z_i(t) = imf_i(t) + j * H[imf_i(t)]

从这个复信号里可以提取三个物理量：

• 瞬时幅值 a_i(t) = sqrt(imf_i^2 + H[imf_i]^2)，反映该模态在某一刻的能量大小；
• 瞬时相位 φ_i(t) = arctan( H[imf_i] / imf_i )；
• 瞬时频率 f_i(t) = (1/2π) * dφ_i(t)/dt。

把时间、瞬时频率、瞬时幅值三个维度画在一张图上，就得到HHT时频谱，也就是常说的Hilbert谱。和短时傅里叶变换的最大区别是，HHT不受海森堡测不准原理那种时间分辨率与频率分辨率互相制约的限制，瞬时频率的定位是逐点计算的，所以对频率突变和瞬态成分更敏感。

举个直观例子：如果第三个IMF对应旋转机械的工频分量，那么它的瞬时频率会在60Hz附近微小波动，瞬时幅值则代表该时刻的振动能量强弱。分析边频带和启停车过程时，FFT看的是一个时间段内的平均频谱，HHT则能把转子升速过程里频率慢慢爬升的轨迹画出来。这套源码的价值就是把这些步骤封装成可调用函数，不必每次从零写筛分循环。

3. 把源码装进环境：路径配置、首次调用与参数初读

3.1 下载后第一件事：看清文件结构，把目录加进MATLAB路径

压缩包解压后通常能看到一批.m文件，核心的几个是emd.m、hht.m、toimage.m、extrema.m以及演示脚本。我一般先打开主目录，确认没有嵌套过深再动手配路径。直接双击某个demo文件运行不一定能成功，因为MATLAB只认当前工作目录，找不到其他函数就会报错，所以第一步是把整个工具箱目录永久加进搜索路径。

cd('你的EMD工具箱绝对路径'); addpath(genpath(pwd)); savepath;

逻辑说明：cd切换当前目录；addpath(genpath(pwd))把当前路径下的所有子目录递归加入搜索路径，保证各个.m文件之间相互调用时不会被“未定义函数”卡住；savepath把路径保存到MATLAB的pathdef.m里，重启后依然有效。

参数说明：pwd是当前工作目录的路径字符串，genpath负责把子目录拼成完整路径列表。如果工具箱目录名带中文或者空格，老版本MATLAB偶尔会路径解析异常，此时可以先改目录名或把工具箱挪到纯英文路径再执行。

3.2 用一段合成信号验证安装是否成功

配置好路径后，最稳妥的验证方式是生成一段双频正弦信号，调用一次emd看看能不能正常返回结果。

fs = 200; t = (0:1/fs:1)'; x = sin(2*pi*20*t) + 0.5*sin(2*pi*60*t); imf = emd(x); disp(size(imf));

逻辑说明：这段代码构造了20Hz和60Hz两个正弦叠加的测试信号，采样率200Hz，一共201个采样点。调用emd后，输出矩阵imf的行数代表分解得到的固有模态函数数量，列数等于输入信号长度。正常情况下会分解出两个主要IMF加一条残差。

参数说明：采样率200Hz对20Hz和60Hz信号是足够的。如果代码报错提示emd未定义，说明路径没有加成功；如果输出IMF数量异常多或第一个IMF混着两个频率，就需要往下看筛选停止参数设置。

3.3 返回值结构和常用参数初读

这套工具箱里最常见的函数是emd，语法通常写成这样：

[imf, residual, info] = emd(x, 'MaxNumIMF', 4, 'SiftLog', 0);

逻辑说明：imf是主输出矩阵，每一行是一个固有模态函数，第一行对应最高频分量，越往下频率越低；residual是最后剩下的趋势项；info是结构体，包含迭代次数、停止条件等信息，调试时可从这里确认筛分收敛情况。

参数说明：MaxNumIMF限制最多分解出多少个模态，实际工程中有时不想让残差分解得过碎，直接写4或5能加快运行；SiftLog控制是否在命令行打印筛分日志，0表示关闭，1表示打开。另一个常见参数是SiftMaxIterations，限制单个IMF的最大迭代次数，默认值在1000到2000之间，信号很长时可以适当降低。

4. 用振动信号做一次实战：完整调用链与结果解读

4.1 先构造一个贴近现场的非平稳信号

纯正弦测试只能证明程序在跑，真正验证算法能力需要非平稳成分。现场采集的振动信号通常由工频正弦、低频缓变、瞬态冲击和随机噪声叠加而成，我习惯先用仿真信号复现这种结构：

fs = 1000; t = (0:1999)/fs; x1 = sin(2*pi*50*t); x2 = 0.6*sin(2*pi*8*t); x3 = sin(2*pi*120*t) .* exp(-((t-0.6)/0.05).^2); x = x1 + x2 + x3 + 0.1*randn(size(t));

逻辑说明：x1模拟工频干扰，x2模拟低频结构振动，x3是一个高斯包络调制的冲击成分，冲击中心在0.6秒附近，频率120Hz，幅值从零涨到1再衰减到零，最后叠加白噪声模拟采集环境。

参数说明：fs=1000Hz决定了可分析的最高频率到500Hz，120Hz的冲击完全落在有效范围内；0.05是高斯包络的宽度系数，决定冲击持续时间；0.1是噪声标准差，用来模拟信号中的随机干扰。

4.2 完整调用链：emd分解加希尔伯特求瞬时频率

仿真信号构造好后，直接进入完整分析流程：

[imf, residual, info] = emd(x, 'SiftLog', 0); figure; for k = 1:size(imf, 1) subplot(size(imf, 1)+1, 1, k); plot(t, imf(k, :)); title(['IMF ', num2str(k)]); end subplot(size(imf,1)+1, 1, size(imf,1)+1); plot(t, residual); title('Residual');

逻辑说明：分解得到的每个IMF分别画在子图里，最下方画残差趋势项。这样做的目的是先肉眼确认分解质量：第一个IMF应该接近冲击和噪声，第二个IMF对应50Hz工频，再往下对应8Hz低频缓变成分，如果有模态混叠，图里会出现“这个IMF里明显夹着另一个频率”的情况。

接下来对每个IMF计算瞬时频率和瞬时幅值：

hilbert_amp = zeros(size(imf)); inst_freq = zeros(size(imf)); for k = 1:size(imf, 1) z = hilbert(imf(k, :)); hilbert_amp(k, :) = abs(z); inst_phase = angle(z); inst_freq(k, 1:end-1) = diff(unwrap(inst_phase)) * fs / (2*pi); end

逻辑说明：hilbert把IMF变成解析信号，abs得到瞬时幅值包络，angle得到瞬时相位，unwrap把相位展开避免跳变，差分后乘以fs/(2π)就得到物理频率值。瞬时频率数组比原信号少一个点，因为diff会损失一位。

参数说明：这里的fs=1000是采样率，直接决定diff结果换算成Hz时的比例系数。如果信号本身有趋势项导致hilbert变换输出异常，可以只对前几个IMF计算瞬时频率，残差通常不做这一步。

4.3 结果怎么读：IMF、瞬时频率和原始频率成分的对应关系

跑完上面的代码，通常第一个IMF的瞬时频率集中在120Hz附近，而且冲击发生的0.6秒前后幅值包络明显增大；第二个IMF的瞬时频率在50Hz上下波动；第三个IMF在8Hz附近。如果看到某条IMF的瞬时频率曲线剧烈跳变到负值，说明该IMF分解不纯，混入了另一个频带成分。

实际工程里，我一般先看HHT谱，再看某条窄带分量的瞬时幅值包络特征，比如轴承局部损伤时，冲击引起的IMF包络会有周期性尖峰，包络谱上能看到故障特征频率的倍频。这套源码里的HHT绘图函数会直接生成三维或二维时频图，但更实用的做法是自己提取瞬时幅值交给后续分类器，而不是停留在看图上。

5. 常见问题与避坑：端点效应、模态混叠和版本差异

5.1 端点效应：两端出现大幅发散毛刺

现象：分解出的第一个IMF在信号两端出现幅度远超中间的振荡，看起来像伪造的假信号，导致瞬时频率也乱跳。

原因：EMD做三次样条包络时，采样信号首尾各只有一个极值点，样条在端点处的斜率不受约束，包络线在两端经常被抬高或压低，筛分过程就会被这种错误包络带偏。

解决：常用做法是端点延拓，即分解前给信号两端各延长若干个极值点，延拓长度取信号总长的2%到5%比较常见。实操中我会先调用extrema函数看首尾极值分布，再用镜像延拓方法复制出去一段，分解完成后再把延拓部分裁掉。工具箱自带端点处理选项时，优先选镜像延拓而不是零延拓。

5.2 模态混叠：一个IMF里面混着两个频率段的成分

现象：某一条IMF的波形和频谱同时包含高频冲击和低频工频，或者两个相邻IMF频谱重叠严重，瞬时频率曲线在多个值之间来回跳变。

原因：信号里存在间歇性小幅度高频成分时，筛选过程会在局部被高频干扰带偏，导致本应属于不同时间尺度的分量被并进同一个IMF。这种现象在含冲击和启停车信号的实测数据里最常见。

解决：先对信号做带通预处理，把明显的高频噪声在分解前滤掉一部分；或者改用集合经验模态分解思路，给原信号叠加多次白噪声再取平均，代价是耗时成倍增加。单独跑EMD时，我更倾向于先观察原始信号频谱，确定感兴趣频带边界，再用滤波把强干扰频段压下去。

5.3 工具箱hht函数与MATLAB自带hht函数混淆

现象：调用hht时报参数错误，或者hht返回的尺寸和手动算的瞬时频率对不上。

原因：MATLAB在较新版本里也内置了hht函数，用于对IMF做希尔伯特-黄变换生成谱图，而这份源码包里同样有一个hht.m。路径搜索顺序的问题在于，如果工具箱目录排在MATLAB内置函数之后，实际调用的就会是内置版本，语法和数据封装方式完全不同。

解决：用which hht查看当前解析到哪个文件，确认返回路径是否指向工具箱目录。如果指向内置函数，就把工具箱目录在路径管理器里上移，或者直接调用define路径下的完整相对路径写法，例如EMD工具箱目录下的hht（此处省略具体目录）。更稳妥的做法是不用高频重名的hht，自己按第4章的hilbert流程算瞬时频率，控制权更大。

5.4 长信号分解性能差：几分钟跑不完或内存耗尽

现象：几百万个采样点的信号，调用emd后运行速度极慢，甚至长时间不响应，或报内存不足。

原因：EMD每次筛分都要对整段信号做极值查找、样条插值和迭代，算法复杂度随数据点数呈超线性增长，点数过大时每次循环都会产生大量中间矩阵。

解决：先把数据降采样到满足奈奎斯特要求的采样率，再做一次去趋势预处理。对超长信号，常见的做法是把信号按时间滑窗切块，每块独立做EMD再拼接，每个分块长度在数千到数万点之间比较合理。分块之间要保留10%左右的重叠，可以有效减轻端点效应跨块影响。另外，循环里做好内存清理，大数组用完后及时赋空。

6. 和FFT对照验证，把EMD用在批处理里

6.1 用FFT校准EMD结果：判断分解是否可靠

HHT谱只是看起来漂亮不行，频率对不对要用FFT交叉验证。做法很简单：对每个IMF做FFT，找峰值频率，再和原始信号的FFT峰值比对。如果某个IMF的FFT峰值明显偏离原始信号中存在的频率成分，这个IMF多半分解无效。

for k = 1:size(imf, 1) f = (0:length(imf)-1) * fs / length(imf); spec = abs(fft(imf(k, :))); [~, idx] = max(spec(2:floor(end/2))); fprintf('IMF%d peak: %.2f Hz\n', k, f(idx+1)); end

逻辑说明：对每个IMF做FFT并打印主峰频率。实际使用时我会再叠加一个阈值判断，如果主峰频率落在已知成分的±5%范围内，就认为这个IMF是可信的。

参数说明：f是频率轴，idx在频谱前半段找最大峰值，因为第二半是共轭对称部分，不存在物理信息。fft长度默认等于信号长度，点数过少时频率分辨率不够，此时可以用fft(imf, 4096)补零加密。

6.2 批处理几十个文件时的参数习惯

EMD不像FFT那样一次算完，参数设置对批处理效率影响很大。我处理一批同类信号时会固定一个MaxNumIMF，比如4或者5，再固定SiftMaxIterations，保证每个文件的分解深度一致。这样得到的IMF矩阵尺寸一致，后续特征提取代码才能无差别处理。

批量脚本里我还会把采样率、分解参数、每个IMF的瞬时频率均值等写成一行一行的表格保存，而不是每跑一个文件就弹一张图。数据的优势就是在大量样本上都稳定时才显现的，单个样本结果不错并不代表批量可用。

在那之后，我每次拿到陌生信号都会强制走一遍固定流程：先看原始波形和频谱，再降采样、去趋势，然后固定参数跑批量，抽出几个样本做FFT对照，确认无误才把结果交给分类模型，这套习惯帮我挡掉了不少“看起来像缺陷实际是端点效应作祟”的翻车现场。希望帮到你。

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