别再只盯着傅里叶了！用Python实战HHT（希尔伯特-黄变换）分析你的非平稳信号
别再只盯着傅里叶了！用Python实战HHT分析你的非平稳信号
当股票价格突然剧烈波动、工业设备发出异常振动、或是语音信号中出现急促的音调变化时，传统傅里叶变换给出的"平均化"频谱图往往让人束手无策。这就是为什么NASA科学家黄锷华教授在1998年提出的希尔伯特-黄变换（HHT）会成为分析非平稳信号的利器——它能捕捉到信号中瞬息万变的频率特征，就像给信号装上了"显微镜"。
本文将用Python带你完整实现HHT全流程，从一段模拟的变频信号开始，通过PyEMD库进行经验模态分解（EMD），再用希尔伯特变换提取各分量的瞬时频率，最后与短时傅里叶变换（STFT）结果对比。你会发现，当信号频率随时间快速变化时，HHT展现的时频分辨率优势令人惊艳。
1. 为什么傅里叶变换会"失明"？
傅里叶变换假设信号是平稳的——即频率成分在整个时间段内保持不变。这就像用一张长时间曝光的照片拍摄飞鸟：最终只能得到模糊的轨迹，而无法看清翅膀扇动的每个瞬间。具体到代码中，我们生成一个频率从5Hz线性增加到25Hz的测试信号：
import numpy as np import matplotlib.pyplot as plt fs = 1000 # 采样率 t = np.linspace(0, 2, 2*fs) # 2秒时长 freq = 5 + 10*t # 线性变化的频率 signal = np.sin(2*np.pi*freq*t) # 变频信号 plt.figure(figsize=(10,4)) plt.plot(t, signal) plt.title('频率从5Hz增加到25Hz的测试信号') plt.xlabel('时间(s)'); plt.ylabel('振幅') plt.show()用FFT分析这个信号会得到什么？下图显示频谱在5-25Hz区间形成"一片"能量分布，完全无法反映频率随时间变化的细节：
from scipy.fft import fft fft_result = fft(signal) freqs = np.linspace(0, fs, len(fft_result)) plt.figure(figsize=(10,4)) plt.plot(freqs[:100], np.abs(fft_result)[:100]) # 显示前100Hz plt.title('FFT频谱分析结果') plt.xlabel('频率(Hz)'); plt.ylabel('幅值') plt.show()关键局限：傅里叶变换的"全局性"使其无法定位频率成分出现的时间点，这正是非平稳信号分析需要突破的核心问题。
2. HHT双剑客：EMD分解与希尔伯特变换
HHT的独特之处在于它的两阶段处理流程，先通过EMD将信号拆解成本征模态函数（IMF），再对每个IMF进行希尔伯特变换获得瞬时频率。
2.1 经验模态分解（EMD）
EMD就像一台智能滤波器，能自适应地将信号分解为不同时间尺度的振动分量。安装PyEMD库后，分解过程只需几行代码：
pip install EMD-signalfrom PyEMD import EMD emd = EMD() IMFs = emd(signal) plt.figure(figsize=(10,8)) for i, imf in enumerate(IMFs): plt.subplot(len(IMFs), 1, i+1) plt.plot(t, imf) plt.ylabel(f'IMF {i+1}') plt.xlabel('时间(s)') plt.tight_layout() plt.show()EMD分解出的IMF具有以下特征：
- 每个IMF的过零点数与极值点数相差不超过1
- 局部均值为零（对称性）
- 按频率从高到低顺序排列
2.2 希尔伯特谱分析
对每个IMF进行希尔伯特变换，可以构造解析信号并计算瞬时频率。PyEMD提供了直接生成希尔伯特谱的方法：
from PyEMD import Visualisation # 创建可视化对象 vis = Visualisation() # 计算希尔伯特谱 hilbert_spectrum = vis._Hilbert_spectrum(IMFs, t, freq_res=0.5) # 绘制时频图 plt.figure(figsize=(10,6)) vis.plot_instant_freq(t, hilbert_spectrum) plt.title('HHT时频分析结果') plt.ylabel('频率(Hz)') plt.show()与STFT的对比实验最能体现HHT的优势。我们使用相同的信号进行短时傅里叶变换：
from scipy.signal import stft f, t_stft, Zxx = stft(signal, fs, nperseg=256) plt.figure(figsize=(10,6)) plt.pcolormesh(t_stft, f, np.abs(Zxx), shading='gouraud') plt.title('STFT时频分析') plt.ylabel('频率(Hz)'); plt.xlabel('时间(s)') plt.colorbar(label='幅值') plt.show()| 对比指标 | HHT | STFT |
|---|---|---|
| 时间分辨率 | 自适应（高） | 固定（依赖窗长） |
| 频率分辨率 | 瞬时频率计算 | 受限于不确定性原理 |
| 边界效应 | 无 | 明显 |
| 计算复杂度 | 较高 | 较低 |
3. 实战：轴承故障振动信号分析
让我们用HHT处理真实的工业数据——来自凯斯西储大学轴承数据中心的故障振动信号。这个数据集包含内圈故障、外圈故障等典型机械故障信号。
import pandas as pd # 加载内圈故障数据 df = pd.read_csv('bearing_fault.csv') vibration = df['vibration'].values t = np.arange(len(vibration))/12000 # 12kHz采样率 # EMD分解 IMFs = emd(vibration) # 计算希尔伯特谱 hilbert_spectrum = vis._Hilbert_spectrum(IMFs, t[:len(IMFs[0])], freq_res=1) # 故障特征频率标记 fault_freq = 157 # 内圈故障特征频率(Hz) plt.figure(figsize=(12,6)) vis.plot_instant_freq(t[:len(IMFs[0])], hilbert_spectrum) plt.axhline(y=fault_freq, color='r', linestyle='--', label='故障特征频率') plt.legend() plt.show()在这个案例中，HHT清晰地捕捉到了157Hz附近的持续高频成分（红色虚线），这正是轴承内圈故障的典型特征。而传统的FFT分析由于频率混叠效应，很难如此直观地识别这一特征频率。
4. 避开HHT的"雷区"：实用技巧与优化
虽然HHT功能强大，但实际应用中需要注意以下几个关键点：
端点效应抑制
EMD在信号两端会产生失真，可通过镜像延拓缓解：emd = EMD(spline_kind='akima', mirror_extremes=True)IMF筛选策略
并非所有IMF都有物理意义，通常保留前几个高频分量：meaningful_IMFs = IMFs[:3] # 取前三个IMF采样率选择
根据Nyquist定理，采样率应至少是最高分析频率的2倍：min_sample_rate = 2 * max_expected_freq计算加速技巧
对于长信号，可先降采样处理：from scipy import signal vibration_ds = signal.decimate(vibration, q=2) # 降采样2倍
专业建议：在金融时间序列分析中，HHT对突发性波动极为敏感。但要注意，IMF分量可能对应着不同时间尺度的市场行为（高频噪声、主力操作、长期趋势等），需要结合领域知识解读。