MATLAB实现经验模态分解（EMD）与希尔伯特变换获取能量谱

一、代码

%% 信号生成与预处理
fs = 1000; % 采样率
t = 0:1/fs:1; % 时间向量
% 生成测试信号（含多分量）
x = sin(2*pi*50*t) + 0.5*sin(2*pi*150*t) + 0.2*randn(size(t)); 

%% EMD分解
imf = emd(x, 'Interpolation', 'pchip', 'Display', 1); % 分解IMF分量
residual = emd(x, 'Interpolation', 'pchip', 'Display', 0); % 残余分量

%% 希尔伯特变换与能量计算
num_IMF = size(imf, 2);
energy_spectrum = zeros(num_IMF, length(t));

for i = 1:num_IMF
    % 希尔伯特变换
    analytic_signal = hilbert(imf(:,i));
    inst_amp = abs(analytic_signal); % 瞬时振幅
    inst_phase = unwrap(angle(analytic_signal)); % 瞬时相位
    inst_freq = (diff(inst_phase)/(2*pi)) * fs; % 瞬时频率
    
    % 能量密度计算
    energy = inst_amp.^2; % 能量密度
    
    % 时频分布
    energy_spectrum(i,:) = energy;
end

%% 能量谱可视化
figure;
for i = 1:num_IMF
    subplot(num_IMF+1,1,i);
    imagesc(t, [1,num_IMF], energy_spectrum');
    set(gca, 'YDir', 'normal');
    title(sprintf('IMF%d 能量分布', i));
    xlabel('时间 (s)');
    ylabel('IMF分量');
end
subplot(num_IMF+1,1,num_IMF+1);
plot(t, sum(energy_spectrum,1));
title('总能量谱');
xlabel('时间 (s)');
ylabel('总能量');

---

二、关键参数

参数	推荐值	作用说明
Interpolation	'pchip'	三次样条插值减少端点效应
Display	1	显示分解进度与IMF数量
窗函数	hann(256)	STFT时频分析窗长选择
重叠率	50%	平衡时频分辨率与计算量

---

三、结果分析要点

1. IMF分量特性 前3个IMF通常包含主要频率成分（示例中50Hz和150Hz） 残余分量反映信号趋势项
2. 能量分布特征 各IMF能量集中在对应频率段（如IMF1能量集中在50Hz附近） 总能量谱呈现多峰结构，反映信号多分量特性
3. 时频分辨率 短时傅里叶变换(STFT)的频率分辨率由窗长决定 建议使用nfft=1024保证频谱细节

---

四、高级功能

1. 边际谱计算

% 计算边际谱（频率能量累积）
marginal_spectrum = sum(energy_spectrum, 1);
figure;
plot(linspace(0,fs/2,length(marginal_spectrum)), marginal_spectrum);
title('边际谱');
xlabel('频率 (Hz)');
ylabel('能量');

2. 三维能量谱可视化

% 创建三维能量谱矩阵
[X,Y] = meshgrid(t, 1:num_IMF);
Z = energy_spectrum';

% 绘制三维曲面
figure;
surf(X, Y, Z, 'EdgeColor', 'none');
xlabel('时间 (s)');
ylabel('IMF分量');
zlabel('能量密度');
shading interp;

---

五、工程应用案例

1. 机械故障诊断（轴承振动信号）

% 加载轴承振动数据
load('bearing_signal.mat');
x = bearing_signal;

% EMD分解与能量分析
imf = emd(x);
energy_spectrum = compute_energy_spectrum(imf);

% 特征频率检测
fault_freq = detect_fault_frequency(energy_spectrum, 1000);

2. 生物医学信号分析（心电ECG）

% 加载ECG信号
load('ecg_signal.mat');
x = ecg_signal;

% 去噪处理
imf = emd(x);
clean_signal = sum(imf(:,1:2), 2);

% 能量特征提取
energy_features = sum(energy_spectrum, 2);

---

六、参考资料

1. 《Hilbert-Huang Transform and Its Applications》（Norden E. Huang）
2. 代码 emd分解之后再进行希尔伯特变换，获得能量谱 www.youwenfan.com/contentcnj/77762.html
3. MATLAB官方文档：emd函数说明
4. 专利CN113267837A：改进的EMD端点处理方法