# hht
Hilbert-Huang 变换
函数库: TySignalProcessing
# 语法
hs = hht(imf)
hs = hht(imf, fs)
hs, f, t = hht(___)
hs, f, t, imfinsf, imfinse = hht(___)
___ = hht(___; Name=Value)
# 说明
hs = hht(imf) 返回由本征模式函数 imf 指定的信号的希尔伯特频谱 hs。hs 可用于分析由光谱内容随时间变化的混合信号组成的信号。使用 hht 对信号进行希尔伯特频谱分析，以识别局部特征。
hs = hht(imf, fs) 返回以 fs 速率采样的信号的希尔伯特频谱 hs。
hs, f, t = hht(___) 除 hs 外，还会返回频率向量 f 和时间向量 t。这些输出参数可与前面的任一输入语法配合使用。
hs, f, t, imfinsf, imfinse = hht(___) 还返回本征模式函数的瞬时频率 imfinsf 和瞬时能量 imfinse，用于信号诊断。
___ = hht(___; Name=Value) 估算希尔伯特频谱参数，附加选项由一个或多个名称-值参数指定。
# 示例
二次 chirp 信号的希尔伯特频谱
生成高斯调制的二次 chirp 信号。指定采样率为 2 kHz，信号持续时间为 2 秒。
using TySignalProcessing
using TyMath
using TyPlot
fs = 2000
t = collect(0:(1 / fs):(2 - 1 / fs))'
q = chirp(t .- 2, 4, 1 / 2, 6, "quadratic", 100) .* exp.(-4 * (t .- 1) .^ 2)
plot(t, q)
使用 emd 直观显示本征模式函数（IMF）和残差。
emd(q; plotfig=true)
计算信号的 IMF。使用 Display 名称-值对输出一个表格，显示每个 IMF 的筛选迭代次数、相对容差和筛选停止标准。
imf, = emd(q; Display=true)
Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit
1 | 1 | 9.1977e-06 | SiftMaxRelativeTolerance
shared_signalwavelet:emd:general:MaxEnergyRatioHitMaxEnergyRatio
使用计算出的 IMF 计算二次 chirp 信号的希尔伯特频谱。将频率范围限制在 0 Hz 至 20 Hz 之间。
hs, f, t, imfinsf, imfinse = hht(imf, fs; FrequencyLimits=[0 20])
多分量信号的 VMD
生成由三个频率分别为 2 Hz、10 Hz 和 30 Hz 的正弦波组成的多分量信号。正弦波的采样频率为 1 kHz，持续 2 秒钟。将信号嵌入方差为 0.01² 的高斯白噪声中。
using TySignalProcessing
using TyPlot
using TyMath
using TyBase
rng = MT19937ar(1234)
fs = 1e3
t = collect(1:(1 / fs):(2 - 1 / fs))'
x =
cos.(2 * pi * 2 .* t) +
2 .* cos.(2 * pi * 10 .* t) +
4 .* cos.(2 * pi * 30 .* t) +
0.01 * randn(rng, 1, length(t))
计算噪声信号的 IMFs，并将其可视化为三维图。
imf = vmd(x)
p, q = ndgrid(t, 1:size(imf[1], 2))
plot3(p, q, imf[1])
grid("on")
xlabel("Time Values")
ylabel("Mode Number")
zlabel("Mode Amplitude")
使用计算出的 IMF 计算多分量信号的希尔伯特频谱。将频率范围限制在 [0, 40] Hz。
hs, f, t, imfinsf, imfinse = hht(imf[1], fs; FrequencyLimits=[0, 40])
计算振动信号的希尔伯特频谱
模拟损坏轴承的振动信号。计算该信号的希尔伯特频谱并查找缺陷。
一个节圆直径为 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)
计算故障轴承信号第一个经验模式的希尔伯特频谱。第一个模式捕捉到了高频冲击的影响。随着轴承磨损的加剧，冲击的能量也在增加。
hs1, f1, t1, imfinsf1, imfinse1 = hht(imfBad[:, 1], fs)
第三模的希尔伯特频谱显示了振动信号中的共振。将频率范围限制在 0 赫兹到 100 赫兹之间。
hs2, f2, t2, imfinsf2, imfinse2 = hht(imfBad[:, 3], fs; FrequencyLimits=[0 100])
计算健康轴承信号的第一和第三模态的希尔伯特频谱。
hs3, f3, t3, imfinsf3, imfinse3 = hht(imfGood[:, 1], fs)
hs4, f4, t4, imfinsf4, imfinse4 = hht(imfGood[:, 3], fs; FrequencyLimits=[0 100])
# 输入参数
imf - 本征模式函数矩阵
本征模式函数，以矩阵的形式指定。imf 是包络相对于零对称的任何信号，其极值和过零点的数量最多相差一个。emd 用于将复杂信号分解和简化为执行希尔伯特频谱分析所需的有限个本征模式函数。
hht 将 imf 中的每一列都视为一个本征模式函数。有关计算 imf 的更多信息，请参阅 emd。
fs - 采样率2π（默认） | 正标量
采样率，指定为正标量。如果没有提供 fs，则使用 2π 的归一化频率计算希尔伯特频谱。如果 imf 被指定为时间表，则会从中推断出采样率。
freqlocation - 频率轴在曲线图上的位置"yaxis"（默认） | "xaxis"
频率轴在曲线图上的位置，指定为 "yaxis" 或 "xaxis"。要在曲线图的 y 轴或 x 轴上显示频率数据，请将 freqlocation 分别指定为 "yaxis" 或 "xaxis"。
# 名称-值参数
以 Name1=Value1,...,NameN=ValueN 的形式指定可选的参数对，其中 Name 是参数名，Value 是相应的值。名称-值参数必须出现在其他参数之后，但参数对的顺序并不重要。
示例： FrequencyResolution=1
FrequencyLimits - 计算希尔伯特频谱的频率限制[0,fs/2]（默认） | 1×2 整数值向量
计算希尔伯特频谱的频率限制，指定为 1×2 整数值向量。FrequencyLimits 的单位为 Hz。
FrequencyResolution - 用于离散频率限制的频率分辨率(f_high-f_low)/100 （默认值） | 正标量
用于离散频率限制的频率分辨率，指定为正标量。频率分辨率的单位为 Hz。如果未指定 FrequencyResolution，则从 FrequencyLimits 推断出 (fhigh-flow)/100 的值。这里，fhigh 是频率上限，flow 是频率下限。
MinThreshold - 希尔伯特频谱的最小阈值-inf（默认） | 标量
希尔伯特谱的最小阈值，以标量形式指定。MinThreshold 将 hs 的元素设置为 0，此时
# 输出参数
hs - 信号的希尔伯特频谱稀疏矩阵
信号的希尔伯特频谱，以稀疏矩阵形式返回。使用 hs 进行时频分析，并识别信号中的局部特征。
f - 频率值向量
hht 使用频率向量 f 和时间向量 t 绘制希尔伯特频谱图。
数学上，f 表示为：
t - 时间值向量
hht 使用时间向量 t 和频率向量 f 绘制希尔伯特频谱图。
imfinsf - 每个 IMF 的瞬时频率向量 | 矩阵
每个 IMF 的瞬时频率，以向量、矩阵的形式返回。
imfinsf 的列数与 imf 相同，并以以下形式返回：
向量（如果 imf 被指定为向量）；
矩阵（如果 imf 被指定为矩阵）。
imfinse - 每个 IMF 的瞬时能量向量 | 矩阵
每个 IMF 的瞬时能量，以向量、矩阵的形式返回。
imfinse 的列数与 imf 相同，并以以下形式返回：
向量（如果 imf 被指定为向量）；
矩阵（如果 imf 被指定为矩阵）。
# 参考文献
[1] Huang, Norden E, and Samuel S P Shen. Hilbert–Huang Transform and Its Applications. 2nd ed. Vol. 16. Interdisciplinary Mathematical Sciences. WORLD SCIENTIFIC, 2014.
[2] Huang, Norden E., Zhaohua Wu, Steven R. Long, Kenneth C. Arnold, Xianyao Chen, and Karin Blank. “ON INSTANTANEOUS FREQUENCY.” Advances in Adaptive Data Analysis 01, no. 02 (April 2009): 177–229.