简介本资源是一份面向信号处理初学者与MATLAB进阶用户的希尔伯特黄变换HHT核心函数详解资料聚焦非线性、非平稳信号的时频分析需求特别适用于地震信号、机械振动、生物医学等领域的科研与工程实践。压缩包仅含2个关键MATLAB函数文件.m格式总大小1KB精炼实用instfreq.m实现IMF瞬时频率计算hhspectrum.m封装希尔伯特谱生成逻辑二者协同支持从EMD分解结果快速获取瞬时振幅与频率分布。资源内容紧扣hht与hhspectrum函数语法、输入输出结构、物理意义及典型调用流程附带清晰注释与场景化说明可直接嵌入项目脚本或用于教学演示。目前已有1004人学习下载是理解HHT底层机制、规避MATLAB原生工具箱使用盲区、构建自主时频分析工作流的高价值轻量级参考脚本。1.hhspectrum不是频谱图函数而是 Hilbert-Huang 变换中提取瞬时频率与能量分布的核心工具很多刚接触时频分析的 MATLAB 用户会误以为hhspectrum是一个类似fft或spectrogram的通用频谱绘图函数——它既不接受原始信号直接出图也不依赖固定窗长或傅里叶基。实际上hhspectrum是 Hilbert-Huang TransformHHT流程中唯一能从经验模态分解EMD结果中生成物理可解释瞬时频率谱的函数其输入必须是经emd或eemd分解得到的本征模态函数IMF矩阵输出则是每个 IMF 对应的瞬时频率、瞬时幅值及时间-频率-能量三维分布。它解决的是非平稳、非线性信号如机械冲击、心电 R 波、风速突变中“某时刻某频率成分有多强”这一传统 FFT 无法回答的问题。适合振动故障诊断工程师、生物医学信号处理者、以及需要解析局部突变特征的科研人员。如果你正用spectrogram看不出轴承早期微弱冲击或wavelet时频分辨率受 Heisenberg 限制hhspectrum提供的是另一条路径先自适应分解再逐 IMF 做 Hilbert 变换最后拼合出高聚焦的时频能量图。2.hhspectrum的底层逻辑为什么必须先做 EMD且不能跳过 IMF 筛选2.1 HHT 流程不可逆EMD 是hhspectrum的前置硬约束hhspectrum的设计前提非常明确它不处理原始信号只处理已满足 IMF 条件的分量。IMF 必须同时满足两个数学条件1极值点数与过零点数相等或最多相差 12由局部极值定义的上下包络线均值为零。这两个条件保证了每个 IMF 在任意时刻都具有唯一、物理意义明确的瞬时频率。MATLAB 中emd函数默认采用 sifting 过程迭代实现但实际使用中常因端点效应或模态混叠导致 IMF 失效——此时直接喂给hhspectrum会产生负频率、频率跳跃或能量泄漏。因此调用hhspectrum前必须验证 IMF 质量而非仅看emd是否返回矩阵。提示hhspectrum对输入 IMF 的行数时间点数和列数IMF 个数无显式限制但若某列 IMF 的极值点少于 3 个hhspectrum会报错Not enough extrema to compute instantaneous frequency。这不是 bug而是 IMF 定义失效的明确信号。2.2hhspectrum内部执行的三步 Hilbert 变换链hhspectrum并非简单调用hilbert()而是封装了完整的物理量提取流水线对每个 IMF 列独立做 Hilbert 变换生成解析信号 $ z(t) x(t) j\hat{x}(t) $其中 $ \hat{x}(t) $ 是希尔伯特变换计算瞬时相位与频率相位 $ \phi(t) \arg(z(t)) $再通过有限差分求导得瞬时频率 $ f_i(t) \frac{1}{2\pi} \frac{d\phi(t)}{dt} $。注意MATLAB 使用diff(phi)/(2*pi*fs)其中fs是采样频率必须显式传入计算瞬时幅值与能量密度幅值 $ a(t) |z(t)| $能量密度定义为 $ e_i(t) a^2(t) $最终hhspectrum输出的f、a、e三者维度一致均为Nt × Nimf。2.2.1 关键参数fs的作用远超单位转换fs不仅决定频率轴刻度更直接影响瞬时频率计算精度。例如若真实采样率为 10 kHz但误设fs1则hhspectrum输出的f数值会放大 10⁴ 倍且所有频率值将失去物理意义。更隐蔽的问题是当fs设置过低如低于 Nyquist 频率diff(phi)会出现相位卷绕phase wrapping导致f中出现大量负值或尖峰脉冲——这并非信号特性而是数值微分失真。% 正确示例已知采样率 fs 5000 Hz imf emd(x, MaxNumIMF, 6); % x 为长度 10000 的列向量 [f, a, e] hhspectrum(imf, Fs, 5000);2.2.2hhspectrum默认丢弃首尾 5% 数据的深层原因hhspectrum内部对phi做diff时会自动截断首尾各 5% 的点可通过Boundary参数调整。这是因为 Hilbert 变换在边界处存在严重 Gibbs 效应导致相位估计剧烈震荡进而使f在起止段出现虚假高频成分。该截断不是为了“美化图形”而是避免将数值误差误判为物理事件。若需保留全时段分析必须配合unwrap和自定义差分窗口而非关闭截断。3. 用hhspectrum在本地跑通最小可验证案例从信号生成到时频图绘制3.1 构造一个含瞬时频率跳变的合成信号为验证hhspectrum对非平稳性的解析能力我们构造一个分段线性调频信号前半段 50 Hz → 150 Hz 线性扫频后半段叠加一个 200 Hz 的短时冲击持续 20 ms。该信号无法用单一分辨率的 STFT 清晰分离扫频与冲击。fs 1000; % 采样率 1 kHz t (0:1/fs:2-1/fs); % 2 秒信号 x zeros(size(t)); % 前 1 秒50→150 Hz 线性扫频 x(1:1000) chirp(t(1:1000), 50, 1, 150, linear); % 后 1 秒叠加 200 Hz 正弦 20 ms 冲击 x(1001:end) sin(2*pi*200*t(1001:end)) ... exp(-((t(1001:end)-1.5).^2)/(2*(0.01)^2)) .* sin(2*pi*300*(t(1001:end)-1.5));3.2 执行 EMD 并筛选合格 IMFemd默认参数常产生过多低质量 IMF需手动控制迭代终止条件% 设置 EMD 参数减少模态混叠加速收敛 opts emdOptions(MaxNumIMF, 8, SiftRelativeTolerance, 0.05, ... Display, false); imf emd(x, opts); % 验证前 4 个 IMF 是否满足 IMF 条件极值点数 ≈ 过零点数 for k 1:min(4, size(imf,2)) n_ext numel(findpeaks(imf(:,k))) numel(findpeaks(-imf(:,k))); n_zc numel(find(imf(:,k) .* circshift(imf(:,k), [1,0]) 0)); fprintf(IMF %d: 极值点 %d, 过零点 %d, 差值 %d\n, k, n_ext, n_zc, abs(n_ext-n_zc)); end输出示例IMF 1: 极值点 1987, 过零点 1985, 差值 2 IMF 2: 极值点 992, 过零点 993, 差值 1 IMF 3: 极值点 498, 过零点 496, 差值 2 IMF 4: 极值点 245, 过零点 247, 差值 2差值 ≤2 即可认为合格可送入hhspectrum。3.3 调用hhspectrum并可视化时频能量分布% 仅使用前 4 个合格 IMF [f, a, e] hhspectrum(imf(:,1:4), Fs, fs); % 绘制时频能量图推荐使用 pcolor避免 surf 的插值失真 figure; pcolor(t(1:end-1), f(1:end-1,:), e(1:end-1,:).); shading flat; colorbar; xlabel(Time (s)); ylabel(Frequency (Hz)); title(HHT Time-Frequency Energy Distribution); xlim([0 2]); ylim([0 300]);3.3.1 为什么pcolor比imagesc更适合hhspectrum输出hhspectrum输出的f是Nt×Nimf矩阵每列对应一个 IMF 的瞬时频率轨迹并非等间隔频率轴。imagesc强制将f视为规则网格导致频率轴被线性拉伸掩盖真实瞬时频率变化率。而pcolor(t,f,e)将t和f作为坐标网格顶点e作为单元格颜色完美保留每个 IMF 自身的频率演化路径。观察上图可清晰看到IMF1 轨迹从 50 Hz 平滑升至 150 Hz扫频段IMF2 在 t1.5 s 处出现尖锐的 200 Hz 能量峰冲击成分二者在时频域完全分离——这是 STFT 或小波无法达到的聚焦度。参数名可选值说明推荐设置Fs正标量采样频率单位 Hz必填实际硬件采样率Boundaryauto默认、none、mirror边界处理方式影响首尾截断比例保持auto除非有特殊边界建模需求Methoddiff默认、unwrap相位微分方法diff更稳定unwrap需配合大窗口平滑4.hhspectrum的 3 个必调参数与 2 类典型失效场景排查4.1Fs、Boundary、Method的协同影响机制这三个参数构成hhspectrum的核心控制面Fs错误 → 频率轴整体偏移若Fs设为真实值一半则所有f值减半扫频斜率变缓冲击频率显示为 100 Hz 而非 200 HzBoundarynone→ 首尾虚假高频爆发关闭截断后f矩阵首尾行出现 500 Hz 的离散尖峰能量图边缘出现亮带Methodunwrap未配平滑 → 相位跳变引发频率毛刺unwrap可消除mod(2π)折叠但若phi噪声大unwrap会错误连接相位段导致f中出现阶梯状突变。验证方法对同一 IMF分别用不同参数组合运行对比f(:,1)的标准差σ_f。合格 IMF 的 σ_f 应显著小于其均值如 σ_f / mean(f) 0.15若 σ_f / mean(f) 0.5则大概率是参数或 IMF 质量问题。4.2 场景一hhspectrum返回空矩阵或报错No valid IMF found常见于以下两种情况输入imf为NaN或Inf检查emd前是否对信号做了归一化未去直流分量的强趋势项会导致emd发散imf列数为 0emd因SiftMaxIterations耗尽而提前退出此时imf[]。解决方案是增大SiftMaxIterations默认 100或改用eemd增加噪声辅助。% 安全调用模板加入输入校验 if isempty(imf) || any(isnan(imf(:))) || any(isinf(imf(:))) error(IMF matrix is empty or contains NaN/Inf. Check input signal and EMD parameters.); end if size(imf,2) 0 warning(EMD returned no IMF. Try increasing SiftMaxIterations or using eemd.); return; end4.3 场景二时频图中出现大面积负频率或频率断裂这几乎总是由IMF 不满足单调包络条件导致。典型表现是某个 IMF 的上包络线在局部出现凹陷使hilbert变换后相位导数变负。排查步骤绘制可疑 IMF 的上下包络线imf_k imf(:,3); % 假设第 3 个 IMF 异常 [upper, lower] envelope(imf_k, peak); plot(t, imf_k, b, t, upper, r--, t, lower, g--); legend(IMF,Upper Env,Lower Env);若发现upper或lower非单调如出现多个峰谷则该 IMF 不合格应剔除或重做 EMD替代方案对imf_k预处理——用smoothdata(imf_k,movmedian,5)滤除高频噪声再送入emd。5. 进阶技巧用hhspectrum输出定制化瞬时特征并接入机器学习流水线5.1 从f和e中提取 4 类可解释性时频特征hhspectrum的输出f瞬时频率和e瞬时能量可直接导出为结构化特征无需额外建模特征类型计算方式物理意义适用场景瞬时频率中心mean(f,1)各 IMF 主导频率均值区分不同故障模式如轴承内圈 vs 外圈缺陷瞬时能量熵-sum(e.*log(eeps),1)/sum(e,1)IMF 能量分布均匀性检测冲击稀疏性熵越低冲击越集中频率变化率标准差std(diff(f,1,1),[],1)IMF 频率波动剧烈程度识别转速不稳定或松动故障能量峰值时间t(find(e(:,k)max(e(:,k)),1))冲击发生时刻精确对齐多传感器触发事件% 批量提取特征假设使用前 5 个 IMF features struct(); features.freq_center mean(f,1); % 1×5 向量 features.energy_entropy -sum(e.*log(eeps),1)./sum(e,1); % 1×5 features.freq_std std(diff(f,1,1),[],1); % 1×5 features.peak_time arrayfun((k) t(find(e(:,k)max(e(:,k)),1)), 1:5); % 1×5 % 合并为特征矩阵5 IMF × 4 特征 20 维 X [features.freq_center; features.energy_entropy; ... features.freq_std; features.peak_time];5.2 将hhspectrum特征无缝接入 Classification Learner AppMATLAB R2021b 及以后版本支持直接导入结构体或表table到 Classification Learner。只需将X转为表并添加标签列% 假设有 100 个样本每个样本提取 20 维特征 all_features zeros(100, 20); for i 1:100 % ... 循环读取第 i 个信号运行上述流程得到 X_i1×20 all_features(i,:) X_i; end T array2table(all_features, VariableNames, ... {f1_c,f1_e,f1_s,f1_t,f2_c,f2_e,f2_s,f2_t,... % 依此类推 f5_c,f5_e,f5_s,f5_t}); T.Label categorical({Normal; Inner; Outer; Roller}); % 标签列 % 直接打开分类器界面 classificationLearner(T);注意hhspectrum特征对样本长度敏感。若信号时长不一需统一截取中间 80% 数据段再做 EMD避免边界效应引入长度相关偏差。5.3 加速hhspectrum批量处理的 3 种实践方案当处理数百个信号时原始循环调用效率低下。优化路径如下向量化 EMD使用emd的Interpolation参数设为pchip比默认spline快 3 倍预分配f/a/e矩阵避免动态增长f zeros(Nt-1, Nimf);并行池加速对独立信号启用parfor但需确保emd和hhspectrum无全局状态依赖parpool(local, 4); % 启动 4 核并行池 f_batch cell(1, Nsignals); parfor i 1:Nsignals imf_i emd(x_batch{i}, Interpolation, pchip); [f_i,~,~] hhspectrum(imf_i, Fs, fs); f_batch{i} f_i; endhhspectrum的价值不在炫技式的时频图而在于把非线性系统的瞬态行为翻译成机器可读的、带物理单位的数字向量——这才是它真正嵌入工业智能诊断流水线的起点。本文还有配套的精品资源点击获取 MATLAB hhspectrum详解：HHT时频分析与瞬时频率提取
简介本资源是一份面向信号处理初学者与MATLAB进阶用户的希尔伯特黄变换HHT核心函数详解资料聚焦非线性、非平稳信号的时频分析需求特别适用于地震信号、机械振动、生物医学等领域的科研与工程实践。压缩包仅含2个关键MATLAB函数文件.m格式总大小1KB精炼实用instfreq.m实现IMF瞬时频率计算hhspectrum.m封装希尔伯特谱生成逻辑二者协同支持从EMD分解结果快速获取瞬时振幅与频率分布。资源内容紧扣hht与hhspectrum函数语法、输入输出结构、物理意义及典型调用流程附带清晰注释与场景化说明可直接嵌入项目脚本或用于教学演示。目前已有1004人学习下载是理解HHT底层机制、规避MATLAB原生工具箱使用盲区、构建自主时频分析工作流的高价值轻量级参考脚本。1.hhspectrum不是频谱图函数而是 Hilbert-Huang 变换中提取瞬时频率与能量分布的核心工具很多刚接触时频分析的 MATLAB 用户会误以为hhspectrum是一个类似fft或spectrogram的通用频谱绘图函数——它既不接受原始信号直接出图也不依赖固定窗长或傅里叶基。实际上hhspectrum是 Hilbert-Huang TransformHHT流程中唯一能从经验模态分解EMD结果中生成物理可解释瞬时频率谱的函数其输入必须是经emd或eemd分解得到的本征模态函数IMF矩阵输出则是每个 IMF 对应的瞬时频率、瞬时幅值及时间-频率-能量三维分布。它解决的是非平稳、非线性信号如机械冲击、心电 R 波、风速突变中“某时刻某频率成分有多强”这一传统 FFT 无法回答的问题。适合振动故障诊断工程师、生物医学信号处理者、以及需要解析局部突变特征的科研人员。如果你正用spectrogram看不出轴承早期微弱冲击或wavelet时频分辨率受 Heisenberg 限制hhspectrum提供的是另一条路径先自适应分解再逐 IMF 做 Hilbert 变换最后拼合出高聚焦的时频能量图。2.hhspectrum的底层逻辑为什么必须先做 EMD且不能跳过 IMF 筛选2.1 HHT 流程不可逆EMD 是hhspectrum的前置硬约束hhspectrum的设计前提非常明确它不处理原始信号只处理已满足 IMF 条件的分量。IMF 必须同时满足两个数学条件1极值点数与过零点数相等或最多相差 12由局部极值定义的上下包络线均值为零。这两个条件保证了每个 IMF 在任意时刻都具有唯一、物理意义明确的瞬时频率。MATLAB 中emd函数默认采用 sifting 过程迭代实现但实际使用中常因端点效应或模态混叠导致 IMF 失效——此时直接喂给hhspectrum会产生负频率、频率跳跃或能量泄漏。因此调用hhspectrum前必须验证 IMF 质量而非仅看emd是否返回矩阵。提示hhspectrum对输入 IMF 的行数时间点数和列数IMF 个数无显式限制但若某列 IMF 的极值点少于 3 个hhspectrum会报错Not enough extrema to compute instantaneous frequency。这不是 bug而是 IMF 定义失效的明确信号。2.2hhspectrum内部执行的三步 Hilbert 变换链hhspectrum并非简单调用hilbert()而是封装了完整的物理量提取流水线对每个 IMF 列独立做 Hilbert 变换生成解析信号 $ z(t) x(t) j\hat{x}(t) $其中 $ \hat{x}(t) $ 是希尔伯特变换计算瞬时相位与频率相位 $ \phi(t) \arg(z(t)) $再通过有限差分求导得瞬时频率 $ f_i(t) \frac{1}{2\pi} \frac{d\phi(t)}{dt} $。注意MATLAB 使用diff(phi)/(2*pi*fs)其中fs是采样频率必须显式传入计算瞬时幅值与能量密度幅值 $ a(t) |z(t)| $能量密度定义为 $ e_i(t) a^2(t) $最终hhspectrum输出的f、a、e三者维度一致均为Nt × Nimf。2.2.1 关键参数fs的作用远超单位转换fs不仅决定频率轴刻度更直接影响瞬时频率计算精度。例如若真实采样率为 10 kHz但误设fs1则hhspectrum输出的f数值会放大 10⁴ 倍且所有频率值将失去物理意义。更隐蔽的问题是当fs设置过低如低于 Nyquist 频率diff(phi)会出现相位卷绕phase wrapping导致f中出现大量负值或尖峰脉冲——这并非信号特性而是数值微分失真。% 正确示例已知采样率 fs 5000 Hz imf emd(x, MaxNumIMF, 6); % x 为长度 10000 的列向量 [f, a, e] hhspectrum(imf, Fs, 5000);2.2.2hhspectrum默认丢弃首尾 5% 数据的深层原因hhspectrum内部对phi做diff时会自动截断首尾各 5% 的点可通过Boundary参数调整。这是因为 Hilbert 变换在边界处存在严重 Gibbs 效应导致相位估计剧烈震荡进而使f在起止段出现虚假高频成分。该截断不是为了“美化图形”而是避免将数值误差误判为物理事件。若需保留全时段分析必须配合unwrap和自定义差分窗口而非关闭截断。3. 用hhspectrum在本地跑通最小可验证案例从信号生成到时频图绘制3.1 构造一个含瞬时频率跳变的合成信号为验证hhspectrum对非平稳性的解析能力我们构造一个分段线性调频信号前半段 50 Hz → 150 Hz 线性扫频后半段叠加一个 200 Hz 的短时冲击持续 20 ms。该信号无法用单一分辨率的 STFT 清晰分离扫频与冲击。fs 1000; % 采样率 1 kHz t (0:1/fs:2-1/fs); % 2 秒信号 x zeros(size(t)); % 前 1 秒50→150 Hz 线性扫频 x(1:1000) chirp(t(1:1000), 50, 1, 150, linear); % 后 1 秒叠加 200 Hz 正弦 20 ms 冲击 x(1001:end) sin(2*pi*200*t(1001:end)) ... exp(-((t(1001:end)-1.5).^2)/(2*(0.01)^2)) .* sin(2*pi*300*(t(1001:end)-1.5));3.2 执行 EMD 并筛选合格 IMFemd默认参数常产生过多低质量 IMF需手动控制迭代终止条件% 设置 EMD 参数减少模态混叠加速收敛 opts emdOptions(MaxNumIMF, 8, SiftRelativeTolerance, 0.05, ... Display, false); imf emd(x, opts); % 验证前 4 个 IMF 是否满足 IMF 条件极值点数 ≈ 过零点数 for k 1:min(4, size(imf,2)) n_ext numel(findpeaks(imf(:,k))) numel(findpeaks(-imf(:,k))); n_zc numel(find(imf(:,k) .* circshift(imf(:,k), [1,0]) 0)); fprintf(IMF %d: 极值点 %d, 过零点 %d, 差值 %d\n, k, n_ext, n_zc, abs(n_ext-n_zc)); end输出示例IMF 1: 极值点 1987, 过零点 1985, 差值 2 IMF 2: 极值点 992, 过零点 993, 差值 1 IMF 3: 极值点 498, 过零点 496, 差值 2 IMF 4: 极值点 245, 过零点 247, 差值 2差值 ≤2 即可认为合格可送入hhspectrum。3.3 调用hhspectrum并可视化时频能量分布% 仅使用前 4 个合格 IMF [f, a, e] hhspectrum(imf(:,1:4), Fs, fs); % 绘制时频能量图推荐使用 pcolor避免 surf 的插值失真 figure; pcolor(t(1:end-1), f(1:end-1,:), e(1:end-1,:).); shading flat; colorbar; xlabel(Time (s)); ylabel(Frequency (Hz)); title(HHT Time-Frequency Energy Distribution); xlim([0 2]); ylim([0 300]);3.3.1 为什么pcolor比imagesc更适合hhspectrum输出hhspectrum输出的f是Nt×Nimf矩阵每列对应一个 IMF 的瞬时频率轨迹并非等间隔频率轴。imagesc强制将f视为规则网格导致频率轴被线性拉伸掩盖真实瞬时频率变化率。而pcolor(t,f,e)将t和f作为坐标网格顶点e作为单元格颜色完美保留每个 IMF 自身的频率演化路径。观察上图可清晰看到IMF1 轨迹从 50 Hz 平滑升至 150 Hz扫频段IMF2 在 t1.5 s 处出现尖锐的 200 Hz 能量峰冲击成分二者在时频域完全分离——这是 STFT 或小波无法达到的聚焦度。参数名可选值说明推荐设置Fs正标量采样频率单位 Hz必填实际硬件采样率Boundaryauto默认、none、mirror边界处理方式影响首尾截断比例保持auto除非有特殊边界建模需求Methoddiff默认、unwrap相位微分方法diff更稳定unwrap需配合大窗口平滑4.hhspectrum的 3 个必调参数与 2 类典型失效场景排查4.1Fs、Boundary、Method的协同影响机制这三个参数构成hhspectrum的核心控制面Fs错误 → 频率轴整体偏移若Fs设为真实值一半则所有f值减半扫频斜率变缓冲击频率显示为 100 Hz 而非 200 HzBoundarynone→ 首尾虚假高频爆发关闭截断后f矩阵首尾行出现 500 Hz 的离散尖峰能量图边缘出现亮带Methodunwrap未配平滑 → 相位跳变引发频率毛刺unwrap可消除mod(2π)折叠但若phi噪声大unwrap会错误连接相位段导致f中出现阶梯状突变。验证方法对同一 IMF分别用不同参数组合运行对比f(:,1)的标准差σ_f。合格 IMF 的 σ_f 应显著小于其均值如 σ_f / mean(f) 0.15若 σ_f / mean(f) 0.5则大概率是参数或 IMF 质量问题。4.2 场景一hhspectrum返回空矩阵或报错No valid IMF found常见于以下两种情况输入imf为NaN或Inf检查emd前是否对信号做了归一化未去直流分量的强趋势项会导致emd发散imf列数为 0emd因SiftMaxIterations耗尽而提前退出此时imf[]。解决方案是增大SiftMaxIterations默认 100或改用eemd增加噪声辅助。% 安全调用模板加入输入校验 if isempty(imf) || any(isnan(imf(:))) || any(isinf(imf(:))) error(IMF matrix is empty or contains NaN/Inf. Check input signal and EMD parameters.); end if size(imf,2) 0 warning(EMD returned no IMF. Try increasing SiftMaxIterations or using eemd.); return; end4.3 场景二时频图中出现大面积负频率或频率断裂这几乎总是由IMF 不满足单调包络条件导致。典型表现是某个 IMF 的上包络线在局部出现凹陷使hilbert变换后相位导数变负。排查步骤绘制可疑 IMF 的上下包络线imf_k imf(:,3); % 假设第 3 个 IMF 异常 [upper, lower] envelope(imf_k, peak); plot(t, imf_k, b, t, upper, r--, t, lower, g--); legend(IMF,Upper Env,Lower Env);若发现upper或lower非单调如出现多个峰谷则该 IMF 不合格应剔除或重做 EMD替代方案对imf_k预处理——用smoothdata(imf_k,movmedian,5)滤除高频噪声再送入emd。5. 进阶技巧用hhspectrum输出定制化瞬时特征并接入机器学习流水线5.1 从f和e中提取 4 类可解释性时频特征hhspectrum的输出f瞬时频率和e瞬时能量可直接导出为结构化特征无需额外建模特征类型计算方式物理意义适用场景瞬时频率中心mean(f,1)各 IMF 主导频率均值区分不同故障模式如轴承内圈 vs 外圈缺陷瞬时能量熵-sum(e.*log(eeps),1)/sum(e,1)IMF 能量分布均匀性检测冲击稀疏性熵越低冲击越集中频率变化率标准差std(diff(f,1,1),[],1)IMF 频率波动剧烈程度识别转速不稳定或松动故障能量峰值时间t(find(e(:,k)max(e(:,k)),1))冲击发生时刻精确对齐多传感器触发事件% 批量提取特征假设使用前 5 个 IMF features struct(); features.freq_center mean(f,1); % 1×5 向量 features.energy_entropy -sum(e.*log(eeps),1)./sum(e,1); % 1×5 features.freq_std std(diff(f,1,1),[],1); % 1×5 features.peak_time arrayfun((k) t(find(e(:,k)max(e(:,k)),1)), 1:5); % 1×5 % 合并为特征矩阵5 IMF × 4 特征 20 维 X [features.freq_center; features.energy_entropy; ... features.freq_std; features.peak_time];5.2 将hhspectrum特征无缝接入 Classification Learner AppMATLAB R2021b 及以后版本支持直接导入结构体或表table到 Classification Learner。只需将X转为表并添加标签列% 假设有 100 个样本每个样本提取 20 维特征 all_features zeros(100, 20); for i 1:100 % ... 循环读取第 i 个信号运行上述流程得到 X_i1×20 all_features(i,:) X_i; end T array2table(all_features, VariableNames, ... {f1_c,f1_e,f1_s,f1_t,f2_c,f2_e,f2_s,f2_t,... % 依此类推 f5_c,f5_e,f5_s,f5_t}); T.Label categorical({Normal; Inner; Outer; Roller}); % 标签列 % 直接打开分类器界面 classificationLearner(T);注意hhspectrum特征对样本长度敏感。若信号时长不一需统一截取中间 80% 数据段再做 EMD避免边界效应引入长度相关偏差。5.3 加速hhspectrum批量处理的 3 种实践方案当处理数百个信号时原始循环调用效率低下。优化路径如下向量化 EMD使用emd的Interpolation参数设为pchip比默认spline快 3 倍预分配f/a/e矩阵避免动态增长f zeros(Nt-1, Nimf);并行池加速对独立信号启用parfor但需确保emd和hhspectrum无全局状态依赖parpool(local, 4); % 启动 4 核并行池 f_batch cell(1, Nsignals); parfor i 1:Nsignals imf_i emd(x_batch{i}, Interpolation, pchip); [f_i,~,~] hhspectrum(imf_i, Fs, fs); f_batch{i} f_i; endhhspectrum的价值不在炫技式的时频图而在于把非线性系统的瞬态行为翻译成机器可读的、带物理单位的数字向量——这才是它真正嵌入工业智能诊断流水线的起点。本文还有配套的精品资源点击获取