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

简介: 信号分析中的经验模态分解(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) 等,它们通过引入噪声辅助分析,有效抑制了模态混叠问题。

相关文章
|
并行计算 算法 计算机视觉
【MATLAB 】 CEEMDAN 信号分解+模糊熵(近似熵)算法
【MATLAB 】 CEEMDAN 信号分解+模糊熵(近似熵)算法
871 0
|
算法 计算机视觉
【MATLAB 】 EMD信号分解+希尔伯特黄变换+边际谱算法
【MATLAB 】 EMD信号分解+希尔伯特黄变换+边际谱算法
685 0
|
并行计算 算法 计算机视觉
【MATLAB 】 EMD信号分解+模糊熵(近似熵)算法
【MATLAB 】 EMD信号分解+模糊熵(近似熵)算法
529 0
|
并行计算 TensorFlow 调度
推荐场景GPU优化的探索与实践:CUDA Graph与多流并行的比较与分析
RTP 系统(即 Rank Service),是一个面向搜索和推荐的 ranking 需求,支持多种模型的在线 inference 服务,是阿里智能引擎团队沉淀多年的技术产品。今年,团队在推荐场景的GPU性能优化上又做了新尝试——在RTP上集成了Multi Stream,改变了TensorFlow的单流机制,让多流的执行并行,作为增加GPU并行度的另一种选择。本文详细介绍与比较了CUDA Graph与多流并行这两个方案,以及团队的实践成果与心得。
|
网络协议 Windows
网络连接正常但百度网页打不开显示无法访问此网站解决方案
网络连接正常但百度网页打不开显示无法访问此网站解决方案
4734 0
网络连接正常但百度网页打不开显示无法访问此网站解决方案
|
人工智能 算法 安全
升级版I 算法备案(实操指南)
本文详解算法备案全流程,涵盖法规依据、适用对象、材料准备及操作步骤,助企业合规完成备案,避免罚款与上架风险。
975 0
升级版I 算法备案(实操指南)
|
机器学习/深度学习
利用matlab提取出频域和时域信号的29个特征
利用matlab提取出频域和时域信号的29个特征
|
机器学习/深度学习 算法 数据挖掘
稀疏促进动态模态分解(SPDMD)详细介绍以及应用
稀疏促进动态模态分解(SPDMD)结合了动态模态分解(DMD)的数学优雅性和稀疏优化技术,有效提取高维数据中的关键特征。SPDMD通过稀疏约束自动筛选出最重要模态,去除冗余信息,提升模型的可解释性和计算效率。该方法在流体动力学、图像处理、时间序列分析及金融数据等领域广泛应用,能够识别主要趋势、周期性模式及异常现象。SPDMD不仅提高了数据分析效率,还为各领域研究提供了强有力的工具。通过自动选择最相关的模态,SPDMD尤其适用于大规模数据集和实时应用。
607 4
|
机器学习/深度学习 Python
训练集、测试集与验证集:机器学习模型评估的基石
在机器学习中,数据集通常被划分为训练集、验证集和测试集,以评估模型性能并调整参数。训练集用于拟合模型,验证集用于调整超参数和防止过拟合,测试集则用于评估最终模型性能。本文详细介绍了这三个集合的作用,并通过代码示例展示了如何进行数据集的划分。合理的划分有助于提升模型的泛化能力。