基于希尔伯特-黄变换的EMD算法实现

资 源 简 介

本项目旨在MATLAB环境下实现经验模态分解（Empirical Mode Decomposition, EMD）算法。EMD是一种由Norden E. Huang提出的针对非线性、非平稳信号的基础处理方法，是希尔伯特-黄变换（HHT）的核心组成部分。该算法的核心逻辑是将复杂的原始信号根据其自身的时间尺度特性，逐步分解为一系列具有不同特征尺度的固有模态函数（Intrinsic Mode Functions, IMF）。每个IMF分量必须满足两个基本条件：其极值点的数量与零点数量相等或最多差一个，且由局部极大值定义的上包络线和局部极小值定义的下包络线的局部平均值为零。在具体实现过程中，程序通过识别信号的所有局部极值点并利用三次样条插值（Cubic Spline Interpolation）拟合出信号的上、下包络线，计算包络线的均值。通过不断的迭代筛选（Sifting Process），将均值从信号中扣除，提取出符合标准的IMF，直到剩余部分成为一个单调函数或不再具备分解价值的残余分量。该程序具备高度自适应性，能够有效捕捉非平稳信号的瞬时跳变和局部细节，广泛应用于地震勘探、生物电信号（如ECG/EEG）处理、机械故障诊断以及水文气象数据分析等多个科研与工程领域。

详 情 说 明

基于希尔伯特-黄变换的经验模态分解(EMD)算法实现

项目介绍

本项目在MATLAB环境下实现了一种非线性、非平稳信号处理的核心算法——经验模态分解（Empirical Mode Decomposition, EMD）。该方法属于希尔伯特-黄变换（HHT）的基础框架，其核心价值在于不依赖于预定义的基函数（如傅里叶变换中的正弦基或小波变换中的小波基），而是根据信号自身的时间尺度特征自适应地提取出各阶固有模态函数（IMF）。

功能特性

1. 自适应信号分解：能够自动将复合信号分解为从高频到低频排列的一系列IMF分量。
2. 完整的筛选过程：内置Sifting算法，通过多次迭代确保每个IMF分量满足局部均值为零及极值点与零点数基本相等的物理约束。
3. 包络拟合技术：采用三次样条插值（Cubic Spline Interpolation）拟合极值点，构建精确的上、下包络线。
4. 边界处理机制：实现了基础的端点延拓策略，以缓解样条插值在信号边缘产生的发散效应。
5. 多维度结果评价：程序不仅提供时域波形图，还同步计算并展示了各IMF分量的频谱分析结果（FFT）。

系统要求

1. 软件环境：MATLAB R2016a 或更高版本。
2. 基础功能需求：需要MATLAB核心函数库支持，特别是样条插值（spline）与快速傅里叶变换（fft）相关函数。

核心实现逻辑

程序遵循标准EMD算法流程，具体步骤如下：

1. 信号模拟与生成

程序首先构建了一个复杂的合成信号，包含15Hz和3Hz的正弦信号、一个二次函数趋势项（代表非平稳趋势）以及随机的高斯白噪声，用于验证算法的分解性能。

1. 算法初始化

设定最大IMF提取上限为10。 定义筛选停止的SD阈值（标准差准则）为0.3。 限制单个IMF提取的最大迭代次数为50次，防止程序无限循环。

1. 提取固有模态函数 (Intrinsic Mode Functions)

通过外层循环控制分解层数，内层循环执行筛选过程： 寻找局部极值：通过一阶差分的正负号切换识别信号的所有局部极大值和极小值。 端点延拓：将信号的起始点和结束点加入极值序列，作为边界控制点。 包络构建：利用三次样条插值，根据极大值序列生成上包络，根据极小值序列生成下包络。 均值扣除：计算上下包络的平均值，并从当前分量中减去此均值，得到更接近IMF定义的中间信号。 收敛判定：计算两次迭代之间的标准差（SD），当SD小于设定阈值时，认为当前信号已达到IMF标准，存入结果矩阵。

1. 残余分量处理

在提取完一个IMF后，从原始信号中扣除该分量。如果剩余部分的极值点数量少于3个（无法再形成包络线），则判定分解结束。

1. 结果可视化

时域分析：垂直排列展示原始信号、各个提取出的IMF分量以及最后的残余分量（Residue），直观反映信号从高频到低频的演变。 频谱分析：对每个IMF进行单边振幅谱计算，帮助用户识别不同模态对应的中心频率，通过频谱的独立性验证分解效果。

关键函数与算法细节说明

1. 极值检测 (get_extrema)

该功能通过计算序列的一阶差分来定位斜率变号的位置。差分值由正变负代表极大值，由负变正代表极小值。

1. 标准差停止准则 (SD Criterion)

程序通过计算相邻两次迭代结果的平方差之和与序列能量的比值（SD）来控制筛选质量。公式为：SD = sum((h_prev - h)^2) / (sum(h_prev^2) + eps)。当SD足够小时，表明信号已趋于对称稳定。

1. 频谱构建 (compute_fft)

利用快速傅里叶变换将提取出的IMF转换至频域。为提高计算效率，程序使用了2的幂次补零（nextpow2）策略，并对幅值进行了归一化和双边对称修正，展示0至50Hz的低频段特征。

1. 辅助工具 (num_2_str_custom)

内置了简单的数值转字符串工具，用于在图形绘制时自动生成自适应的坐标轴标签，增强了结果展示的自动化程度。