基于Hilbert-Huang变换（HHT）的模态分解与瞬时频率计算

一、HHT算法原理框架

graph TD A[原始信号] --> B[经验模态分解] B --> C{IMF条件验证} C -->|满足| D[保存IMF] C -->|不满足| B D --> E[希尔伯特变换] E --> F[瞬时频率计算] F --> G[时频谱重构]

---

二、核心实现步骤

2.1 经验模态分解（EMD）

算法流程：

1. 极值点检测：找出信号所有局部极大/极小值点
2. 包络线拟合：三次样条插值得到上下包络线
3. 均值计算：包络线均值作为本征模态函数(IMF)的基准
4. 迭代筛选：原始信号减去均值后重复检测，直至满足IMF条件

MATLAB实现：

function imfs = emd(signal)
    imfs = {};
    residual = signal;
    while true
        % 极值点检测
        [max_peaks, min_peaks] = find_extrema(residual);
        % 包络线拟合
        upper_env = spline(max_peaks(:,1), max_peaks(:,2), 1:length(residual));
        lower_env = spline(min_peaks(:,1), min_peaks(:,2), 1:length(residual));
        % 计算均值
        mean_env = (upper_env + lower_env)/2;
        % 筛选IMF
        imf_candidate = residual - mean_env;
        % IMF条件验证
        if is_imf(imf_candidate)
            imfs{end+1} = imf_candidate;
            residual = residual - imf_candidate;
            if norm(residual) < 0.1*norm(signal)
                break;
            end
        else
            residual = imf_candidate;
        end
    end
end

2.2 希尔伯特变换

数学表达：
\(z(t)=x(t)+jx^(t)\)
其中\(x^(t)\)为希尔伯特变换结果，解析信号z(t)的相位导数即为瞬时频率。

MATLAB实现：

function [amp, freq] = hilbert_analysis(imf)
    analytic = hilbert(imf);
    amp = abs(analytic);
    phase = angle(analytic);
    freq = diff(phase)/(2*pi*0.001); % 假设采样间隔0.001s
end

---

三、关键参数优化

参数	推荐值	作用说明
停止阈值	0.1-0.3	控制IMF分解精度
包络插值方法	三次样条	减少端点效应
端点处理	镜像延拓	延长信号长度至2倍原始长度
筛选迭代次数	10-15次	平衡分解效率与精度

---

四、典型应用案例

4.1 机械故障诊断

案例描述：轴承故障信号分析

% 加载振动信号
[fs, data] = load_vibration_data('bearing_fault.mat');

% HHT分解
imfs = emd(data);
num_imfs = length(imfs);

% 瞬时频率计算
[amp, freq] = deal(cell(1,num_imfs));
for i = 1:num_imfs
    [amp{i}, freq{i}] = hilbert_analysis(imfs{i});
end

% 特征频率提取
fault_freq = 1200; % 理论故障频率Hz
[~, idx] = min(abs(mean(freq{2}) - fault_freq));
diagnosis_result = ['外圈损伤 (IMF' num2str(idx) ')'];

4.2 生物医学信号处理

ECG信号分析：

% 加载ECG信号
[fs, ecg] = load_ecg_signal('mitdb_100.mat');

% HHT分解
imfs = emd(ecg);
num_imfs = length(imfs);

% 时频谱重构
t = (1:length(ecg))/fs;
figure;
imagesc(t, 1:num_imfs, amp);
xlabel('时间(s)');
ylabel('IMF分量');
title('Hilbert时频谱');

---

五、结果验证方法

1. 重构误差计算：

   reconstruction_error = norm(signal - sum(imfs, 2))/norm(signal);
   

2. 模态混叠检测：

   function mixing_degree = detect_mode_mixing(imfs)
       cross_corr = xcorr(imfs{1}, imfs{2});
       [~, peak] = max(cross_corr);
       mixing_degree = 1 - abs(peak)/length(imfs{1});
   end
   

3. 物理意义验证：

   % 检查IMF是否符合窄带条件
   bandwidth = bandwidth_estimator(imfs{1});
   is_physical = bandwidth < 0.3*fs/length(imfs{1});
   

参考代码 HHT变换实现模态分解和瞬时频率计算 www.youwenfan.com/contentcnh/64429.html

六、扩展功能

1. 多尺度分析：

   function multi_scale = multi_scale_analysis(signal, scales)
       multi_scale = cell(1,length(scales));
       for i = 1:length(scales)
           scaled_signal = resample(signal, scales(i), 1);
           multi_scale{i} = emd(scaled_signal);
       end
   end
   

2. 特征融合：

   function fused_feat = feature_fusion(imfs)
       num_imfs = length(imfs);
       for i = 1:num_imfs
           [amp, freq] = hilbert_analysis(imfs{i});
           features(i,:) = [mean(amp), std(freq), skewness(amp)];
       end
       fused_feat = mean(features);
   end
   

3. 实时处理框架：

   function real_time_hht(data_stream)
       persistent emd_obj;
       if isempty(emd_obj)
           emd_obj = emd('Init');
       end
       imf = process_block(emd_obj, data_stream);
       [amp, freq] = hilbert_analysis(imf);
       update_display(amp, freq);
   end
   

---

该方法通过结合经验模态分解与希尔伯特变换，实现了非线性非平稳信号的时频特征提取。实际应用中建议根据信号特性调整分解层数，并配合显著性检验筛选有效IMF分量。