基于MATLAB的HHT变换完整实现（含EMD分解与三维时频谱生成）

一、核心流程设计

% 主程序流程
[原始信号] → [EMD分解] → [IMF分量] → [希尔伯特变换] → [幅值/相位/瞬时频率] → [三维时频谱]

二、完整MATLAB代码实现

%% 参数设置
fs = 1000;          % 采样频率
t = 0:1/fs:1-1/fs;  % 时间轴
f1 = 50; f2 = 120;  % 信号频率成分
x = sin(2*pi*f1*t) + 0.5*sin(2*pi*f2*t) + 0.3*randn(size(t)); % 含噪信号

%% EMD分解（关键步骤）
imf = emd(x, 'MaxNumIMF', 5, 'Display', 1); % 分解为5个IMF分量
residual = emd(x, 'Interpolation', 'pchip'); % 残差分量

%% 希尔伯特变换（特征提取）
analytic = cell(size(imf));
for i = 1:size(imf,2)
    analytic{
   i} = hilbert(imf(:,i)); % 解析信号计算
end

% 瞬时特征提取
[inst_amp, inst_phase] = deal(zeros(size(imf)));
inst_freq = cell(size(imf));
for i = 1:size(imf,2)
    phase = angle(analytic{
   i}); 
    inst_amp(:,i) = abs(analytic{
   i}); % 幅值谱
    inst_phase(:,i) = unwrap(phase); % 相位谱

    % 瞬时频率计算（改进方法）
    delta_t = diff(t);
    delta_phase = diff(inst_phase(:,i));
    inst_freq{
   i} = [inst_phase(1,i); delta_phase./(2*pi*delta_t)]; % 弧度/秒
end

%% 三维时频谱构建
[~,f_idx] = min(abs(fft_freq(fs,1/fs)-50:1:120)); % 频率索引
t_res = 0.01; % 时间分辨率
f_res = 1; % 频率分辨率

% 构建时频矩阵（Ridge提取算法）
Z = zeros(length(t_res), length(f_idx));
for i = 1:size(imf,2)
    [t_f, f_f] = tfridge(inst_amp(:,i), inst_freq{
   i}, fs);
    Z(:,f_idx) = Z(:,f_idx) + interp1(t_f, f_f, t, 'linear', 0);
end

%% 可视化分析
figure;

% 子图1：IMF分量与原始信号
subplot(3,2,1);
plot(t, x, 'k', t, imf(:,1), 'r', t, imf(:,2), 'b');
title('原始信号与IMF分量'); 
xlabel('时间(s)'); ylabel('幅值');

% 子图2：瞬时频率分布
subplot(3,2,2);
imagesc(t, 1:size(imf,2), inst_freq);
set(gca,'YDir','normal');
title('瞬时频率分布');
xlabel('时间(s)'); ylabel('IMF序号');
colorbar;

% 子图3：三维时频谱
subplot(3,2,3);
surf(linspace(0,1,length(t)), linspace(1,size(imf,2),size(imf,2)), Z);
shading interp;
title('三维时频谱');
xlabel('时间(s)'); ylabel('IMF序号'); zlabel('幅值');

% 子图4：Hilbert时频谱
subplot(3,2,4);
hs = hht(imf,fs); 
imagesc(hs.time, hs.f, hs.power);
title('Hilbert时频谱');
xlabel('时间(s)'); ylabel('频率(Hz)');
colorbar;

三、参数优化

四、结果分析要点

1. IMF分量验证

   • 检查每个IMF是否满足窄带条件（带宽<0.3倍采样率/长度）

   • 计算重构误差：

     reconstruction_error = norm(x - sum(imf,2))/norm(x);
     

2. 时频特征提取 幅值谱：反映信号能量随时间和频率的分布 相位谱：显示信号相位随时间的变化规律 瞬时频率：通过inst_freq = diff(phase)/(2*pi*dt)计算

3. 三维时频谱优化 使用Ridge提取算法增强主频线 采用对数频率轴提升低频分辨率

参考代码 HHT变换，通过对实测波进行EMD分解为IMF，生成幅值，相位谱，瞬时频率，三维时频谱 www.youwenfan.com/contentalg/84545.html

五、工程应用案例

1. 机械故障诊断 轴承故障信号分解后，IMF2出现120Hz特征频率（对应外圈损伤）
2. 生物医学信号分析 ECG信号分解后，IMF1对应QRS波群，IMF2对应P波
3. 地球物理勘探 地震信号分解后，高频IMF对应微破裂事件

六、常见问题解决方案

1. 模态混叠 改用EEMD（集合经验模态分解） 增加停止阈值至0.3-0.5

2. 端点效应

   • 采用镜像延拓法：

     function x_pad = boundary_extension(x,n_pad)
         x_pad = [flipud(x(1:n_pad)); x; flipud(x(end-n_pad+1:end))];
     end
     

3. 频率漂移

   • 对长时间信号分段处理（每段1024点）

   • 使用自适应窗长：

     window = 2^nextpow2(len(segment));
     

七、扩展功能实现

1. 边际谱计算

   marginal_spectrum = sum(inst_freq, axis=0) * fs/length(inst_freq);
   

2. 时频分辨率优化

   opts.t = 0.1;    % 时间分辨率
   opts.f = 1;      % 频率分辨率
   [hs, f, t] = hht(IMFs, fs, opts);
   

3. 交互式可视化

   import plotly.graph_objects as go
   fig = go.Figure(data=[go.Surface(z=Z, x=time_grid, y=freq_grid)])
   fig.update_layout(title='Interactive 3D HHT Spectrum')
   fig.show()