希尔伯特谱（Hilbert Spectrum）详解

DuHz

已于 2025-02-10 14:51:41 修改

阅读量1.3k

点赞数 41

文章标签：算法信息与通信信号处理数学建模傅里叶分析

于 2025-02-10 14:49:03 首次发布

本文链接：https://blog.csdn.net/qq_44648285/article/details/145550227

版权

希尔伯特谱（Hilbert Spectrum）详解

背景简介
希尔伯特变换与解析信号
瞬时幅值、瞬时相位和瞬时频率
本征模态分解（EMD）步骤与数学推导
希尔伯特变换在EMD中的应用
希尔伯特谱的定义与构建
希尔伯特边际谱
希尔伯特谱的特征与优势
实际应用举例
数学公式补充与整理
代码示例
代码简要说明

1. 背景简介

对于很多在实际生活、工业生产、医学检测中遇到的非平稳、非线性信号，传统的时频分析方法（如傅里叶变换、短时傅里叶变换STFT、小波变换等）有一定局限：往往需要事先假设或选择某种基函数，或者频率分辨率与时间分辨率无法同时兼顾。

希尔伯特-黄变换（Hilbert-Huang Transform, HHT） 是由NASA的黄锷（N. E. Huang）等人提出，主要包括两大步骤：

本征模态分解（EMD, Empirical Mode Decomposition）：将原信号分解为若干个自适应的、本征模态函数（IMF）。
希尔伯特谱（Hilbert Spectrum）：对每个IMF做希尔伯特变换，提取瞬时幅值和瞬时频率，并将它们在时频平面中以幅值或能量的形式可视化。

核心目标：更好地分析和可视化信号在时间和频率上的局部细节，且不需要先验的基函数。

2. 希尔伯特变换与解析信号

2.1 希尔伯特变换的定义

给定一个实值信号 $x (t)$ ，它的希尔伯特变换（Hilbert Transform）记作：
$\mathcal{H}\{x(t)\} = \frac{1}{\pi} \text{P.V.} \int_{-\infty}^{+\infty} \frac{x(\tau)}{t - \tau} \, d\tau,$
其中，“ $\text{P.V.}$ ”表示柯西主值积分（Principal Value）。

例如，当我们在离散或数值计算时，往往使用基于卷积或快速傅里叶变换(FFT)的方式来实现希尔伯特变换。而在理论层面，这个积分揭示了一个与原信号相位差为 $\pm \frac{\pi}{2}$ 的正交分量。

2.2 解析信号

将原信号 $x (t)$ 与其希尔伯特变换 $y (t)$ 组合起来，可以得到解析信号（analytic signal） $z (t)$ :
$\, y(t).$
此时 $x (t)$ 为实部， $y (t)$ 为虚部，使得 $z (t)$ 成为一个在正频域具有能量分布的复信号。

3. 瞬时幅值、瞬时相位和瞬时频率

解析信号 $z (t)$ 的极坐标形式可以写作：
$e^{j\theta(t)},$
其中：

瞬时幅值（Instantaneous Amplitude）:
$\sqrt{x^2(t) + y^2(t)},$
瞬时相位（Instantaneous Phase）:
$\theta(t) = \arctan\!\Big(\frac{y(t)}{x(t)}\Big) \quad (\text{通常实际计算中做相位展开}).$
瞬时频率（Instantaneous Frequency）:
$\omega(t) = \frac{d \theta(t)}{d t},$
或者以赫兹表示则为 $\frac{\omega(t)}{2\pi}$ 。

这些量能够描述信号在局部时间（“瞬时”）上的幅值和频率变化，这对于分析非平稳信号具有极大意义。

4. 本征模态分解（EMD）步骤与数学推导

本征模态分解（EMD）是希尔伯特-黄变换的第一步，用于将原信号拆分成若干本征模态函数（IMF, Intrinsic Mode Function）。每个IMF满足如下两条性质：

在全局范围内，极值点数与零交叉点数（或相差不超过1）相近或相等。
在任何时刻，IMF的上包络线与下包络线的平均值为零。

4.1 EMD算法的核心“筛选（Sifting）”过程

设原信号为 $x (t)$ 。EMD的“筛选”循环通常这样进行：

找极值：在 $x (t)$ 中识别所有极大值点 ${t_i^{\max}, x_i^{\max}\}$ 和极小值点 ${t_j^{\min}, x_j^{\min}\}$ 。
插值构造上、下包络：
- 通过这些极大值点，用三次样条或其他平滑插值法得到一个上包络函数 $e_{\max}(t)$ 。
- 通过这些极小值点，同理得到一个下包络函数 $e_{\min}(t)$ 。
计算局部平均：
$m_1(t) = \frac{e_{\max}(t) + e_{\min}(t)}{2}.$
提取候选分量：
$h_1(t) = x(t) - m_1(t).$
这步旨在“去除”信号的低频或高频趋势，保留相对“振荡更快”的部分。
重复筛选：
- 检查 $h_1(t)$ 是否满足IMF条件（极值点数与过零点数接近，且上/下包络平均值近乎零）。
- 若不满足，则令 $\leftarrow h_1(t)$ ，再次执行 步骤1-4，得 $h_2(t)$ 。
- 继续重复，直到获得一个满足IMF定义的函数 $\tilde{c}_1(t)$ 。把它称作第一IMF。
分离并迭代：
$r_1(t) = x(t) - \tilde{c}_1(t).$
将第一IMF从信号中移除，得到剩余信号 $r_1(t)$ 。然后对 $r_1(t)$ 做同样的过程，依次提取 第二IMF、第三IMF……直到剩余信号 $r_n(t)$ 不再含有明显振荡或达到某个停止准则。

最后，原信号被拆分为：
$\sum_{k=1}^{n} c_k(t) + r_n(t),$
其中 $c_k(t)$ 为第 $k$ 个IMF， $r_n(t)$ 是剩余项（可以看作全局趋势或误差）。

4.2 数学层面的简要说明

(极值-包络插值) 让我们在每一层级分离时，充分利用信号自身的极值结构去定义局部振荡成分；
(上/下包络平均) 使得提取出的分量保证“振荡对称性”；
(多次筛选迭代) 则是自适应地“逼近”一个符合IMF定义的分量。

这个过程十分自适应，不需要任何先验基函数（比如傅里叶基或小波基），它直接利用数据本身的局部特征完成分解。

5. 希尔伯特变换在EMD中的应用

当我们利用EMD得到了一系列IMF之后，每个IMF往往被视为一个“近似单分量”或“本征振荡”，具有更明确的瞬时频率意义。

对每个IMF（记作 $c_k(t)$ ）进行希尔伯特变换，可得到其对应的瞬时幅值、瞬时相位和瞬时频率：

解析信号：
$z_k(t) = c_k(t) + j\, \mathcal{H}\{c_k(t)\}.$
瞬时幅值：
$A_k(t) = \sqrt{ c_k^2(t) + \mathcal{H}\{c_k(t)\}^2 }.$
瞬时相位：
$\theta_k(t) = \arctan\!\Big(\frac{\mathcal{H}\{c_k(t)\}}{c_k(t)}\Big).$
瞬时角频率：
$\omega_k(t) = \frac{d\theta_k(t)}{dt}.$

这样，我们就能得到不同时间点上每个本征模态分量的瞬时频率与幅值。

6. 希尔伯特谱的定义与构建

希尔伯特谱（Hilbert Spectrum） 是将所有IMF的瞬时幅值及瞬时频率信息集中可视化的三维表示（时间-频率-幅值）。在理论上，我们可以定义一个时频分布函数 $H(\omega, t)$ ：
$H(\omega, t) = \sum_{k=1}^{n} A_k(t)\,\delta\bigl(\omega - \omega_k(t)\bigr),$
其中 $\omega_k(t)$ 是第 $k$ 个IMF在时刻 $t$ 的瞬时角频率， $A_k(t)$ 为对应的瞬时幅值， $\delta(\cdot)$ 为狄拉克 $\delta$ -函数。

在离散或数值计算中，我们通常会将每个IMF在每个采样时间点对应的瞬时频率 $\omega_k(t_i)$ 映射到频率轴，并用幅值 $A_k(t_i)$ 做加权。然后通过插值或散点绘图的方式，构建时间-频率二维坐标，再用颜色深浅/亮度等方式表达幅值或能量：

横轴: 时间 $t$
纵轴: 频率 $\omega / (2\pi)$
颜色/亮度: 幅值 $A_k(t)$ 或能量 $A_k^2(t)$

这种时频图就是希尔伯特谱，可以清晰显示信号在不同时刻的频率成分及其强度。

7. 希尔伯特边际谱

希尔伯特边际谱（Hilbert Marginal Spectrum）是对希尔伯特谱在时间上的积分或某种统计汇总。可写作：
$h(\omega) = \int_{0}^{T} H(\omega, t) \, dt,$
也就是把每一时刻上在频率 $\omega$ 处的幅值（或能量）做积分，得到一个随频率变化的整体分布。形式上可以展开为：
$h(\omega) = \int_{0}^{T} \sum_{k=1}^{n} A_k(t)\,\delta(\omega - \omega_k(t)) \, dt = \sum_{k=1}^{n}\int_{0}^{T} A_k(t)\,\delta(\omega - \omega_k(t)) \, dt.$
它类似于傅里叶幅值谱，但又是基于瞬时频率在时间上的变化进行统计得到的。它更能体现非平稳信号在整个观察区间里各频率所占的“权重”。

8. 希尔伯特谱的特征与优势

自适应性：无需先验基函数，EMD过程完全利用信号本身特征，适用于非线性、非平稳信号。
高时频分辨率：与STFT、小波相比，希尔伯特谱在某些情况下能够更清晰地展现瞬时频率。
物理意义直观：瞬时频率和瞬时幅值常能直接映射到系统或设备的真实运动、振动或活动模式中，更易解释。
可视化强：在时-频平面，用颜色/灰度显示幅值或能量，能直观看到频率成分的随时间变化轨迹。

9. 实际应用举例

机械故障检测：旋转机械、齿轮箱等的振动信号通常是非平稳的，如转速不恒定时的故障特征频率会随时间变化；希尔伯特谱可帮助识别故障发生的时间和频率区间。
生物医学信号分析：对脑电(EEG)、心电(ECG)等复杂信号做EMD + 希尔伯特谱，可获得不同波段的神经活动或心律变化随时间的细节。
语音信号：语音往往包含多种共振峰且时变性强，利用瞬时频率可更好地捕捉说话人发音动态。
金融数据：研究非平稳金融时间序列的局部波动模式，也能用EMD/HHT进行探索性分析。

10. 数学公式补充与整理

下表给出本文中涉及的主要数学公式总结，便于快速查阅：

内容	公式
希尔伯特变换定义	$\displaystyle \mathcal{H}\{x(t)\} \;=\; \frac{1}{\pi} \,\text{P.V.}\int_{-\infty}^{+\infty}\! \frac{x(\tau)}{\,t-\tau\,}\, d\tau$ $
解析信号	$\displaystyle z(t) = x(t) + j\,\mathcal{H}\{x(t)\}$
瞬时幅值	$\displaystyle A(t) =
瞬时相位	$\displaystyle \theta(t) = \arctan\!\Bigl(\frac{\mathcal{H}\{x(t)\}}{x(t)}\Bigr)$
瞬时角频率	$\displaystyle \omega(t) = \frac{d \theta(t)}{d t}$
EMD分解形式	$\displaystyle x(t) = \sum_{k=1}^{n} c_k(t) + r_n(t)$
EMD求第 $k$ 个IMF	1) 识别极值并插值构造包络 2) 求平均 $m_1(t)$ 3) $h_1(t) = x(t) - m_1(t)$ 4) 反复筛选以得到 $c_k(t)$
对IMF的希尔伯特变换	$\displaystyle z_k(t) = c_k(t) + j\,\mathcal{H}\{c_k(t)\},\;\; A_k(t),\,\omega_k(t)\,\text{同理可得}$
希尔伯特谱	$\displaystyle H(\omega, t) = \sum_{k=1}^{n} A_k(t)\,\delta\bigl(\omega - \omega_k(t)\bigr)$
希尔伯特边际谱	$\displaystyle h(\omega) = \int_{0}^{T}\! H(\omega, t)\,dt \;=\; \sum_{k=1}^{n}\!\int_{0}^{T}\! A_k(t)\,\delta(\omega - \omega_k(t))\, dt$

11. 代码示例

下面示例使用Python演示如何对一个简单的三频正弦叠加信号进行EMD分解并计算希尔伯特谱。这里使用第三方库 PyEMD 完成 EMD，使用 scipy.signal.hilbert 完成希尔伯特变换，并用散点方式画出希尔伯特谱的雏形。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
from PyEMD import EMD

# 1. 生成测试信号（混合三个不同频率的正弦波）
fs = 1000      # 采样率
t = np.linspace(0, 1, fs, endpoint=False)

s1 = np.sin(2 * np.pi * 50 * t)        # 50Hz, 幅值1
s2 = 0.5 * np.sin(2 * np.pi * 120 * t) # 120Hz, 幅值0.5
s3 = 0.2 * np.sin(2 * np.pi * 180 * t) # 180Hz, 幅值0.2
s = s1 + s2 + s3

# 2. EMD分解
emd = EMD()
imfs = emd(s)  # imfs.shape -> (IMF个数, 信号长度)

# 3. 对每个IMF做希尔伯特变换，得到瞬时幅值与瞬时频率
num_imfs = imfs.shape[0]
instantaneous_amps = []
instantaneous_freqs = []

for i in range(num_imfs):
    # (a) 解析信号
    analytic_signal = hilbert(imfs[i])
    # (b) 瞬时幅值
    amp = np.abs(analytic_signal)
    # (c) 相位
    phase = np.unwrap(np.angle(analytic_signal))
    # (d) 瞬时频率 (rad/s -> Hz)
    inst_freq = np.diff(phase) * fs / (2.0 * np.pi)
    
    # 与 inst_freq 对齐，幅值少最后一点
    instantaneous_amps.append(amp[:-1])
    instantaneous_freqs.append(inst_freq)

# 4. 简易散点画希尔伯特谱 (time, freq, amp)
time = t[:-1]  # 配合freq长度
plt.figure(figsize=(10,6))

for i in range(num_imfs):
    # 用散点图表达：(time, freq)，颜色代表瞬时幅值
    plt.scatter(time, instantaneous_freqs[i], 
                s=5, c=instantaneous_amps[i], cmap='jet', alpha=0.7)

plt.xlabel('Time (s)')
plt.ylabel('Frequency (Hz)')
plt.title('Hilbert Spectrum (scatter demo)')
plt.colorbar(label='Amplitude')
plt.ylim(0, 250)
plt.tight_layout()
plt.show()