# emd
经验模式分解
函数库: TySignalProcessing
# 语法
imf,residual = emd(x)
imf,residual,info = emd(x)
___ = emd(___;Name=Value)
emd(___;plotfig=true)
# 说明
imf,residual = emd(x) 使用 emd 可以将复杂信号分解和简化为有限个本征模态函数，以便进行希尔伯特频谱分析。
imf,residual,info = emd(x) 返回 IMFs 和残差信号的额外信息，用于诊断目的。 示例
___ = emd(___;Name=Value) 使用一个或多个 Name=Value 对参数指定的附加选项执行经验模式分解。 示例
emd(___;plotfig=true) 在同一图中以子图的形式绘制原始信号、IMFs 和残差信号。
# 示例
正弦信号的本征模函数的过零点和极值点
这个三角恒等式展现了同一个物理信号的两种不同视角：
生成正弦波信号 s，使用 emd 来计算信号的本征模态函数（IMFs）和额外的诊断信息。该函数默认输出一个表格，指示每个 IMF 的筛选迭代次数、相对容差和筛选停止标准。
using TySignalProcessing
using TyMath
using TyPlot
using TyBase
using Printf
t = collect(0:1e-3:10)'
omega1 = 2 * pi * 100
omega2 = 2 * pi * 20
s =
0.25 * cos.((omega1 - omega2) .* t) +
2.5 * cos.(omega1 .* t) +
0.25 * cos.((omega1 + omega2) .* t)
imf, ~, info = emd(s)
绘制 IMF，并选择从 2 s 开始的 0.5 s 间隔。IMF 是一种 AM 信号，因为 emd 将信号视为幅度调制。
plot(t, imf)
xlim([2 2.5])
xlabel("Time (s)")
ylabel("IMF")
计算振动信号的本征模函数
模拟损坏轴承的振动信号。计算该信号的希尔伯特频谱并查找缺陷。
一个节圆直径为 12 cm 的轴承有八个滚动体。每个滚动体的直径为 2 cm。当内滚道以每秒 25 个周期的速度被驱动时，外滚道保持静止。加速度计以 10 kHz 的频率对轴承振动进行采样。
using TySignalProcessing
using TyPlot
using TyMath
using TyBase
fs = 10000
f0 = 25
n = 8
d = 0.02
p = 0.12
健康轴承的振动信号包括几个数量级的驱动频率。
t = collect(0:(1 / fs):(10 - 1 / fs))'
yHealthy = [1 0.5 0.2 0.1 0.05] * sin.(2 * pi * f0 .* [1 2 3 4 5]' .* t) / 5
在测量过程的中途，轴承振动会产生共振。
yHealthy =
(1 .+ 1 ./ (1 .+ collect(-10:(20 / (length(yHealthy) - 1)):10)' .^ 4)) .* yHealthy
共振在轴承外圈中产生缺陷，导致轴承逐渐磨损。该缺陷会导致轴承的滚珠通过频率外圈（BPFO）重复发生一系列撞击：
其中
ca = 15
bpfo = n * f0 / 2 * (1 - d / p * cosd(ca))
使用 pulstran 函数将撞击模拟为一列周期性的 5 ms 正弦波。每个 3 kHz 的正弦波都有一个平顶窗口。使用幂律在轴承振动信号中引入渐进磨损。
fImpact = 3000
tImpact = collect(0:(1 / fs):(5e-3 - 1 / fs))'
wImpact = flattopwin(length(tImpact))' / 10
xImpact = sin.(2 * pi * fImpact .* tImpact) .* wImpact
tx = collect(0:(1 / bpfo):t[end])'
tx = [tx; 1.3 .^ tx .- 2]
nWear = 49000
nSamples = 100000
yImpact = pulstran(t, tx', xImpact, fs) ./ 5
yImpact = [zeros(1, nWear) yImpact[1, (nWear + 1):nSamples]']
将冲击信号与健康轴承信号相加，生成 BPFO 振动信号。绘制信号图，选择从 5.0 s 开始的 0.3 s 间隔。
yBPFO = yImpact + yHealthy
xLimLeft = 5.0
xLimRight = 5.3
yMin = -0.6
yMax = 0.6
plot(t, yBPFO)
hold("on")
limLeft, limRight = meshgrid2([xLimLeft; xLimRight], [yMin; yMax])
plot(limLeft, limRight, "--")
hold("off")
放大所选区间，直观显示影响效果。
xlim([xLimLeft xLimRight])
为信号添加高斯白噪声。指定噪声方差为 1/150^2。
rn = 150
rng = mt19937ar(1234)
yGood = yHealthy + randn(rng, size(yHealthy)) / rn
rng1 = mt19937ar(1234)
yBad = yBPFO + randn(rng1, size(yHealthy)) / rn
plot(t, yGood, t, yBad)
xlim([xLimLeft xLimRight])
legend("Healthy", "Damaged")
使用 emd 对健康轴承信号进行经验模态分解。计算前五个本征模式函数（IMF）。使用 Display 名称-值参数输出一个表格，显示每个 IMF 的筛选迭代次数、相对容差和筛选停止标准。
imfGood, = emd(yGood; MaxNumIMF=5, Display=true)
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
1 | 3 | 0.016506 | SiftMaxRelativeTolerance
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
2 | 3 | 0.1098 | SiftMaxRelativeTolerance
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
3 | 7 | 0.12073 | SiftMaxRelativeTolerance
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
4 | 1 | 0.0075983 | SiftMaxRelativeTolerance
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
5 | 2 | 0.020977 | SiftMaxRelativeTolerance
shared_signalwavelet:emd:general:MaxNumIMFHit
使用 emd（不带输出参数）可视化前三个 IMF 和残差。
emd(yGood, MaxNumIMF=5; plotfig=true)
计算并显示故障轴承信号的 IMF。第一个经验模式显示了高频冲击。这种高频模式的能量随着磨损的加剧而增加。
imfBad, = emd(yBad; MaxNumIMF=5, Display=true)
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
1 | 2 | 0.041513 | SiftMaxRelativeTolerance
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
2 | 3 | 0.13395 | SiftMaxRelativeTolerance
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
3 | 6 | 0.18922 | SiftMaxRelativeTolerance
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
4 | 1 | 0.015403 | SiftMaxRelativeTolerance
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
5 | 2 | 0.027262 | SiftMaxRelativeTolerance
shared_signalwavelet:emd:general:MaxNumIMFHit
emd(yBad, MaxNumIMF=5; plotfig=true)
# 输入参数
x - 输入信号向量
时间域信号，指定为实值向量。
数据类型： Float | Int
# 名称-值对参数
指定可选的参数对为 Name1=Value1,...,NameN=ValueN，其中 Name 是参数名称，Value 是相应的值。Name = Value 参数必须出现在其他参数之后，但参数对的顺序并不重要。
示例： "MaxNumIMF"=5
SiftRelativeTolerance - 柯西收敛准则0.2（默认值）| 正标量
收敛准则，指定为由 "SiftRelativeTolerance" 和正标量组成的名称参数对。SiftRelativeTolerance 是筛选停止准则之一，也就是说，当当前相对容差小于 SiftRelativeTolerance 时，筛选停止。有关更多信息，请参见筛分相对容差。
数据类型： Float | Int
SiftMaxIterations - 筛选迭代的最大次数100（默认值）| 正标量整数
筛选迭代的最大次数，指定为由 "SiftMaxIterations" 和正标量整数组成的名称参数对。SiftMaxIterations 是筛选停止标准之一，即当当前迭代次数大于 SiftMaxIterations 时，筛选停止。
SiftMaxIterations 只能指定为正整数。
数据类型： Int
MaxNumIMF - 提取的 IMF 的最大数量10（默认值）| 正值标量整数
提取的 IMF 的最大数量，指定为由 "MaxNumIMF" 和正值标量整数组成的名称参数对。MaxNumIMF 是分解停止准则之一，即当生成的 IMF 数量等于 MaxNumIMF 时，分解停止。
MaxNumIMF 只能指定为正整数。
数据类型： Int
MaxNumExtrema - 残差信号中的极值的最大数量1（默认）| 正标量整数
残差信号中的极值的最大数量，指定为由 "MaxNumExtrema" 和正标量整数组成的名称参数对。MaxNumExtrema 是分解停止标准之一，即当极值数量小于 MaxNumExtrema 时停止分解。
MaxNumExtrema 只能指定为正整数。
数据类型： Int
MaxEnergyRatio - 信号与残余能量比20（默认值）| 标量
信号到残余能量比。MaxEnergyRatio 是信号在筛选开始时与平均包络能量之比。MaxEnergyRatio 是分解停止标准之一，即当当前能量比大于MaxEnergyRatio 时，分解停止。有关更多信息，请参见能量比。
数据类型： Float | Int
Interpolation - 用于包络构建的插值方法"spline"（默认值）| "pchip"
用于包络构建的插值方法，指定为名称参数对，例如 "Interpolation"="spline" 或 "Interpolation"="pchip"。
将 Interpolation 的指定值有两种情况:
如果 x 是光滑信号，则为 "spline"；
如果 x 是非光滑信号，则为 "pchip"。
"spline" 插值方法使用三次样条，而 "pchip" 使用分段三次 Hermite 插值多项式。
数据类型： String
Display - 在命令窗口中切换信息显示false（默认）| true
在命令窗口中切换信息显示，指定为以 "Display" 和 false 或 true 组成的名称参数对。在命令窗口生成的表格指示每个生成的 IMF 的 sift 迭代次数、相对容限和 sift 停止准则。将 Display 指定为 true 以显示表格，或者指定为 false 以隐藏表格。
数据类型： Bool
# 输出参数
IMF - 内在模态函数矩阵
内在模态函数（IMF），返回为一个矩阵。每个 IMF 都是一个振幅和频率调制的信号，具有正向和缓慢变化的包络。为了对信号进行谱分析，您可以将其 IMF 应用希尔伯特-黄变换。
数据类型： Float
residual - 信号的残差列向量
信号的残差，以列向量的形式返回。残差表示原始信号 x 中未被 emd 分解的部分。
数据类型： Float | Int
info - 诊断的附加信息结构体
诊断的附加信息，以包含以下字段的结构形式返回：
NumIMF - 提取的 IMF 数量
NumIMF 是一个从 1 到 N 的向量，其中 N 是 IMF 的数量。如果没有提取 IMF，则 NumIMF 为空；
NumExtrema - 每个 IMF 中的极值个数
NumExtrema 是一个长度与 IMF 数量相等的向量。NumExtrema 的第 k 个元素是在第 k 个 IMF 中发现的极值个数。如果没有提取 IMF，则 NumExtrema 为空；
NumZerocrossing - 每个 IMF 中的过零点个数
每个 IMF 中的零交叉次数。NumZerocrossing 是一个长度与 IMF 数量相等的向量。NumZerocrossing 的第 k 个元素是第 k 个 IMF 中的过零点个数。如果没有提取 IMF，则 NumZerocrossing 为空；
NumSifting - 用于提取每个 IMF 的筛选迭代次数
NumSifting 是一个长度与 IMF 数量相等的向量。NumSifting 的第 k 个元素是用于提取第 k 个 IMF 的筛选迭代次数。如果没有提取 IMF，则 NumSifting 为空；
MeanEnvelopeEnergy - 为每个 IMF 提取的上下包络平均值的能量
如果 UE 是上包络线，LE 是下包络线，则 MeanEnvelopeEnergy 为 mean((LE+UL)/2).^2)。MeanEnvelopeEnergy 是一个长度与 IMF 数量相等的向量。MeanEnvelopeEnergy 的第 k 个元素是第 k 个 IMF 的平均包络能。如果没有提取 IMF，则 MeanEnvelopeEnergy 为空；
RelativeTolerance - 每个 IMF 的最终残差相对容差
相对容差的定义是前一筛分步骤的残差与当前筛分步骤的残差之差的平方 2-norm 与第 i 个筛分步骤的残差的平方 2-norm 之比。当相对容差小于 SiftRelativeTolerance 时，筛分过程停止。更多信息，请参见筛选相对容差。相对容差是一个长度与 IMF 数量相等的向量。RelativeTolerance 的第 k 个元素是第 k 个 IMF 的最终相对容差。如果没有提取 IMF，则 RelativeTolerance 为空。
# 更多
经验模式分解
经验模态分解（EMD）算法通过迭代过程将信号
首先找到
的局部最小值和最大值； 然后利用局部极值分别构建
的下包络 和上包络 。得出包络的平均值 ； 从
减去均值，得到残差： 。
分解概述如下：
首先，设
，其中 为初始信号，设 ； 在筛选之前，先检查
： (1) 求
的局部极值总数 (TN)； (2) 找出
的能量比 (ER)。 如果 (ER > MaxEnergyRatio) 或 (TN < MaxNumExtrema) 或 (IMF 数量 > MaxNumIMF) 则停止分解；
让
； 筛选
，得到 ； 检查
： (1) 查找
的相对容差 (RT)； (2) 获取当前的筛分迭代次数 (IN)。
如果 (RT < SiftRelativeTolerance) 或 (IN > SiftMaxIterations) 则停止筛分。找到 IMF：
。否则，让 并转到步骤 5； 让
； 让
。返回步骤 2。
更多信息，请参见 [1] 和 [3]。
本征模式函数
EMD 算法通过迭代筛选过程，将信号 x(t) 分解为 IMFs
在 Huang 等人首次提出 IMF 时[1]，IMF 被定义为具有两个特征的函数：
局部极值的数量（局部最小值和局部最大值的总数）和过零点的数量最多相差一个；
由局部极值构建的上下包络的平均值为零。
然而，正如文献[4]所指出的，筛选直到获得严格的 IMF 可能会导致没有物理意义的 IMF。具体来说，筛选到过零点和局部极值的数量最多相差一个时，就会得到类似纯音的 IMF，换句话说，与在傅立叶基础上投影得到的函数非常相似。这种情况正是 EMD 所要避免的，它更倾向于使用具有物理意义的 AM-FM 调制成分。
参考文献 [4] 提出了获得有物理意义的结果的方案。emd 函数通过使用 Sift Relative Tolerance（一种 Cauchy 类型的停止准则）来放宽原始 IMF 定义。emd 函数迭代提取自然的 AM-FM 模式。生成的 IMF 可能无法满足局部极值-零交叉标准。请参见正弦波内在模式函数中的零交叉和极值。
筛分相对容差
Sift 相对容差是 [4] 中提出的一种 Cauchy 型停止准则。当当前相对容差小于 SiftRelativeTolerance 时，筛分停止。当前相对容差定义为：
Relative Tolerance
由于柯西准则并不直接计算过零点和局部极值的数量，因此分解返回的 IMF 有可能不符合本征模态函数的严格定义。在这种情况下，可以尝试降低 SiftRelativeTolerance 的默认值。有关停止标准的详细讨论，请参见 [4]。该参考文献还讨论了在经验模式分解中坚持使用严格定义的 IMF 的优缺点。
能量比
能量比是筛选开始时的信号能量与平均包络能量之比 [2]。当当前能量比大于最大能量比时，分解停止。对于第 i 个 IMF，能量比定义为：
Energy Ratio
# 参考文献
[1] Huang, Norden E., Zheng Shen, Steven R. Long, Manli C. Wu, Hsing H. Shih, Quanan Zheng, Nai-Chyuan Yen, Chi Chao Tung, and Henry H. Liu. “The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis.” Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, no. 1971 (March 8, 1998): 903–95. https://doi.org/10.1098/rspa.1998.0193.
[2] Rato, R.T., M.D. Ortigueira, and A.G. Batista. “On the HHT, Its Problems, and Some Solutions.” Mechanical Systems and Signal Processing 22, no. 6 (August 2008): 1374–94. https://doi.org/10.1016/j.ymssp.2007.11.028.
[3] Rilling, Gabriel, Patrick Flandrin, and Paulo Gonçalves. "On Empirical Mode Decomposition and Its Algorithms." IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing 2003. NSIP-03. Grado, Italy. 8–11.
[4] Wang, Gang, Xian-Yao Chen, Fang-Li Qiao, Zhaohua Wu, and Norden E. Huang. “On Intrinsic Mode Function.” Advances in Adaptive Data Analysis 02, no. 03 (July 2010): 277–93. https://doi.org/10.1142/S1793536910000549.