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

信号分析中的经验模态分解（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的信号处理工具箱提供了emd和hilbert函数，可以方便地进行计算。

示例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） 等，它们通过引入噪声辅助分析，有效抑制了模态混叠问题。