信号分析中的经验模态分解和希尔伯特变换

简介: 信号分析中的经验模态分解(EMD) 和希尔伯特变换(HT) 是处理非线性、非平稳信号的强大工具。它们共同构成了希尔伯特-黄变换(Hilbert-Huang Transform, HHT) 的核心。

信号分析中的经验模态分解(EMD) 和希尔伯特变换(HT) 是处理非线性、非平稳信号的强大工具。它们共同构成了希尔伯特-黄变换(Hilbert-Huang Transform, HHT) 的核心。

1. 经验模态分解 (EMD)

EMD的终极目标是将一个复杂的信号自动分解为一系列本质模态函数(Intrinsic Mode Functions, IMFs)。IMFs是满足以下两个条件的函数:

  1. 极值点与过零点数量相等或最多相差一个
  2. 在任意时间点,由局部极大值定义的包络线和由局部极小值定义的包络线的均值为零

EMD算法步骤(筛分过程 - Sifting Process):

  1. 识别极值点: 找到原始信号 $x(t)$ 的所有局部极大值和局部极小值。
  2. 构造包络线: 分别用三次样条插值法连接所有极大值点形成上包络线 $$e_{\text{max}}(t)$$,连接所有极小值点形成下包络线 $e_{\text{min}}(t)$。
  3. 计算均值包络线: $$m_1(t) = \frac{e_{\text{max}}(t) + e_{\text{min}}(t)}{2}$$
  4. 提取中间信号: $h_1(t) = x(t) - m_1(t)$
  5. 判断是否为IMF: 检查 $h_1(t)$ 是否满足IMF的两个条件。
    • 如果满足,则 $c_1(t) = h_1(t)$ 就是第一个IMF。
    • 如果不满足,则将 $h_1(t)$ 作为新的原始信号,重复步骤1-4,直到满足条件为止。
  6. 分离剩余分量: 从原始信号中分离出第一个IMF,得到剩余信号 $r_1(t) = x(t) - c_1(t)$。
  7. 迭代分解: 将 $r_1(t)$ 作为新的原始信号,重复上述所有步骤,依次提取出 $c_2(t), c_3(t), ..., c_n(t)$。当剩余信号 $r_n(t)$ 成为一个单调函数常量,无法再提取出IMF时,分解过程停止。

最终,原始信号被分解为:
$$x(t) = \sum_{i=1}^{n} c_i(t) + r_n(t)$$


2. 希尔伯特变换 (HT)

希尔伯特变换用于计算一个实信号的解析信号(Analytic Signal),从而得到信号的瞬时幅度瞬时相位(进而得到瞬时频率)。

对于一个实信号 $x(t)$,其希尔伯特变换 $\hat{x}(t)$ 定义为:
$$\hat{x}(t) = \mathcal{H}[x(t)] = \frac{1}{\pi} \text{P.V.} \int_{-\infty}^{\infty} \frac{x(\tau)}{t - \tau} d\tau$$
其中 P.V. 表示柯西主值。

  • 解析信号 (Analytic Signal): $z(t) = x(t) + i\hat{x}(t) = a(t)e^{i\theta(t)}$
  • 瞬时幅度 (Instantaneous Amplitude): $a(t) = \sqrt{x^2(t) + \hat{x}^2(t)}$
  • 瞬时相位 (Instantaneous Phase): $\theta(t) = \arctan\left(\frac{\hat{x}(t)}{x(t)}\right)$
  • 瞬时频率 (Instantaneous Frequency): $f(t) = \frac{1}{2\pi} \frac{d\theta(t)}{dt}$

注意:直接对任意信号求瞬时频率在物理上可能是无意义的。瞬时频率的概念只有在单分量窄带的信号中才有明确的物理意义。这正是EMD存在的价值——它将复杂信号分解为IMF,而每个IMF可以看作是单分量的,从而使得求取的瞬时频率具有物理意义。


3. EMD与希尔伯特变换的结合 (Hilbert-Huang Transform)

HHT就是先通过EMD将信号分解为IMF,再对每个IMF分量进行希尔伯特变换的过程。

对每个IMF分量 $c_i(t)$ 进行希尔伯特变换后,可以得到:
$$c_i(t) = a_i(t) \cos\left( \int \omega_i(t) dt \right)$$

原始信号可以表示为:
$$x(t) = \Re \sum_{i=1}^{n} a_i(t) e^{j \int \omega_i(t) dt}$$

由此可以得到希尔伯特谱 (Hilbert Spectrum)
$$H(\omega, t) = \sum_{i=1}^{n} a_i(t) e^{j \int \omega_i(t) dt}$$

希尔伯特谱精确地描述了幅度(或能量)在时间和频率平面上的分布,非常适合分析非平稳信号。


4. MATLAB代码

MATLAB的信号处理工具箱提供了emdhilbert函数,可以方便地进行计算。

示例1:EMD分解

% 生成一个混合非平稳信号
fs = 1000; % 采样频率 1000Hz
t = 0:1/fs:2; % 2秒时间向量

% 信号成分:一个频率变化的成分 + 一个突变成分 + 噪声
comp1 = chirp(t, 10, 1, 50, 'quadratic'); % 10Hz到50Hz的线性调频信号
comp2 = zeros(size(t));
comp2(t >= 0.5 & t < 0.8) = 1.5 * sin(2*pi*100*t(t >= 0.5 & t < 0.8)); % 100Hz的短时脉冲
comp3 = 0.1 * randn(size(t)); % 高斯白噪声

x = comp1 + comp2 + comp3; % 混合信号

% 执行EMD分解
[imf, residual] = emd(x, 'Interpolation', 'pchip', 'Display', 1);

% 绘制结果
figure;
subplot(length(imf(1,:))+2, 1, 1);
plot(t, x);
title('原始信号');
ylabel('幅度');

for i = 1:size(imf, 2)
    subplot(size(imf,2)+2, 1, i+1);
    plot(t, imf(:, i));
    title(['IMF ', num2str(i)]);
    ylabel('幅度');
end

subplot(size(imf,2)+2, 1, size(imf,2)+2);
plot(t, residual);
title('残差');
xlabel('时间 (s)');
ylabel('幅度');

示例2:希尔伯特变换与希尔伯特谱

% 继续使用上面的EMD结果
% 选择第一个IMF进行希尔伯特变换分析
imf1 = imf(:, 1);

% 计算希尔伯特变换和解析信号
analytic_signal = hilbert(imf1);
inst_amplitude = abs(analytic_signal); % 瞬时幅度
inst_phase = unwrap(angle(analytic_signal)); % 瞬时相位(解卷绕)
inst_freq = diff(inst_phase) / (2*pi) * fs; % 瞬时频率 (需要差分)

% 绘制单个IMF的希尔伯特分析结果
figure;
subplot(4, 1, 1);
plot(t, imf1);
title('IMF 1');
ylabel('幅度');

subplot(4, 1, 2);
plot(t, inst_amplitude);
title('瞬时幅度');
ylabel('幅度');

subplot(4, 1, 3);
plot(t, inst_phase);
title('瞬时相位');
ylabel('弧度');

subplot(4, 1, 4);
plot(t(2:end), inst_freq); % 频率少一个点
title('瞬时频率');
xlabel('时间 (s)');
ylabel('频率 (Hz)');
ylim([0, 100]); % 限制频率范围以便观察

% 使用内置函数绘制所有IMF的希尔伯特谱(需要安装Time-Frequency Toolbox或其他自定义函数)
% 这里提供一个简易版本的希尔伯特谱绘制
figure;
[hsp, f, t] = hilbertspectrum(imf, fs); % 需要自定义hilbertspectrum函数或使用其他工具箱
imagesc(t, f, 10*log10(hsp)); % 用分贝显示
set(gca, 'YDir', 'normal');
xlabel('时间 (s)');
ylabel('频率 (Hz)');
title('希尔伯特谱');
colorbar;

注意: MATLAB没有直接绘制Hilbert Spectrum的内置函数,你可能需要自己编写一个循环来计算每个IMF的瞬时频率和幅度,然后在时间-频率网格上进行累加。

参考代码 信号分析中用到的EMD分解和希尔伯特变换 www.youwenfan.com/contentale/55198.html


5. 优缺点与应用

优点:

  • 自适应性: EMD完全由数据驱动,基函数自适应产生,不需要预先设定基函数。
  • 有效性: 非常适合分析非线性、非平稳信号。
  • 高分辨率: 希尔伯特谱在时间和频率上都具有很高的分辨率。

缺点:

  • 端点效应: EMD在信号端点处的样条插值可能不稳定,导致分解失真。
  • 模态混叠: 不同物理过程的尺度可能出现在同一个IMF中,或者同一尺度出现在不同IMF中。
  • 计算量大: 筛分过程是迭代的,计算成本较高。
  • 理论不完善: 缺乏严格的数学理论基础。

应用:

  • 机械故障诊断: 轴承、齿轮的振动信号分析。
  • 生物医学信号处理: 脑电图(EEG)、心电图(ECG)、肌电图(EMG)分析。
  • 地球物理信号分析: 地震波、海浪、气候变化分析。
  • 金融时间序列分析: 股票价格波动分析。

为了克服EMD的缺点,后续发展出了很多变体算法,如集合经验模态分解(EEMD)互补集合经验模态分解(CEEMDAN) 等,它们通过引入噪声辅助分析,有效抑制了模态混叠问题。

相关文章
|
机器学习/深度学习 移动开发 算法
经验模态分解(Empirical Mode Decomposition ,EMD)特征提取及其原理
经验模态分解(Empirical Mode Decomposition ,EMD)特征提取及其原理
经验模态分解(Empirical Mode Decomposition ,EMD)特征提取及其原理
|
并行计算 算法 计算机视觉
【MATLAB 】 CEEMDAN 信号分解+模糊熵(近似熵)算法
【MATLAB 】 CEEMDAN 信号分解+模糊熵(近似熵)算法
889 0
|
算法 计算机视觉
【MATLAB 】 EMD信号分解+希尔伯特黄变换+边际谱算法
【MATLAB 】 EMD信号分解+希尔伯特黄变换+边际谱算法
703 0
|
并行计算 算法 计算机视觉
【MATLAB 】 EMD信号分解+模糊熵(近似熵)算法
【MATLAB 】 EMD信号分解+模糊熵(近似熵)算法
541 0
|
数据采集 机器学习/深度学习 算法
解密EEMD分析:Rlibeemd包带你玩转信号分解和时间序列预测
EEMD(Ensemble Empirical Mode Decomposition)是一种信号分解方法,它旨在分解非线性、非平稳或非白噪声的信号,以揭示复杂信号的局部特征和周期性成分。EEMD不同于传统的余弦变换、小波变换等线性变换方法,而是基于经验模态分解(EMD)的思想,通过添加噪声和多次重复分解的方式,避免了EMD中出现的模态混叠等问题。
1442 0
|
机器学习/深度学习
【从零开始学习深度学习】21. 卷积神经网络(CNN)之二维卷积层原理介绍、如何用卷积层检测物体边缘
【从零开始学习深度学习】21. 卷积神经网络(CNN)之二维卷积层原理介绍、如何用卷积层检测物体边缘
|
机器学习/深度学习 搜索推荐 PyTorch
【机器学习】图神经网络:深度解析图神经网络的基本构成和原理以及关键技术
【机器学习】图神经网络:深度解析图神经网络的基本构成和原理以及关键技术
4281 3
|
9月前
|
人工智能 架构师 程序员
学历对程序员的深远影响:2025年的现实与思考-优雅草卓伊凡
学历对程序员的深远影响:2025年的现实与思考-优雅草卓伊凡
304 12
学历对程序员的深远影响:2025年的现实与思考-优雅草卓伊凡
|
XML 缓存 Linux
17.2 Linux libxml2安装
libxml2 是一个用来解析XML文档的函数库。它用 C 语言写成,并且能被多种语言所调用,如 C、c++、XSH、C#、Python、Kylix、Delphi、Ruby、PHP 等。它最初是为 GNOME 开发的项目,但是现在可以用在各种各样的项目中。
885 0
17.2 Linux libxml2安装