希尔伯特黄变换 (HHT) MATLAB 算法程序

一、算法原理

希尔伯特黄变换 (Hilbert-Huang Transform, HHT) 是一种用于非线性、非平稳信号处理的自适应时频分析方法，由经验模态分解 (EMD) 和 希尔伯特谱分析 (HSA) 两部分组成。

算法流程

原始信号 → EMD分解 → 固有模态函数(IMF) → 希尔伯特变换 → 时频谱

---

二、步骤

2.1 主程序 (HHT_Main.m)

%% =====================================================
%% 希尔伯特黄变换 (HHT) 主程序
%% =====================================================

clc; clear; close all;

%% 1. 生成测试信号
fs = 1000;              % 采样频率 (Hz)
t = 0:1/fs:2;          % 时间向量 (2秒)
N = length(t);         % 信号长度

% 构造非线性非平稳信号
f1 = 5; f2 = 50;       % 频率成分
A1 = 1; A2 = 0.5;      % 幅值

% 信号1: 线性调频信号
signal1 = A1 * chirp(t, f1, max(t), f2);

% 信号2: 幅值调制信号
signal2 = A2 * sin(2*pi*f1*t) .* exp(-2*t);

% 信号3: 低频趋势项
signal3 = 0.3 * sin(2*pi*0.5*t);

% 合成信号
x = signal1 + signal2 + signal3;

%% 2. 经验模态分解 (EMD)
fprintf('正在进行经验模态分解 (EMD)...\n');
[imf, residual] = emd_decomposition(x, fs);

num_imf = size(imf, 1);  % IMF数量
fprintf('EMD分解完成，共得到 %d 个IMF分量\n', num_imf);

%% 3. 希尔伯特谱分析
fprintf('正在进行希尔伯特谱分析...\n');
[H, f, t_hht] = hilbert_spectrum_analysis(imf, fs);

%% 4. 计算边际谱
marginal_spectrum = mean(abs(H), 2);

%% 5. 可视化结果
figure('Position', [100, 100, 1400, 900]);

% 原始信号
subplot(4, 2, 1);
plot(t, x, 'b-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('幅值');
title('原始信号');
grid on;

% EMD分解结果
for i = 1:min(num_imf, 4)
    subplot(4, 2, i+1);
    plot(t, imf(i, :), 'k-', 'LineWidth', 1.2);
    xlabel('时间 (s)'); ylabel('幅值');
    title(sprintf('IMF%d', i));
    grid on;
end

% 残差项
subplot(4, 2, 6);
plot(t, residual, 'r-', 'LineWidth', 1.2);
xlabel('时间 (s)'); ylabel('幅值');
title('残差项');
grid on;

% 希尔伯特时频谱
subplot(4, 2, 7);
imagesc(t_hht, f, abs(H));
set(gca, 'YDir', 'normal');
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('希尔伯特时频谱');
colorbar;

% 边际谱
subplot(4, 2, 8);
plot(f, marginal_spectrum, 'b-', 'LineWidth', 1.5);
xlabel('频率 (Hz)'); ylabel('幅值');
title('边际谱');
grid on;

fprintf('HHT分析完成！\n');

%% 6. 瞬时频率和瞬时幅值分析
figure('Position', [200, 200, 1000, 600]);
for i = 1:min(num_imf, 3)
    subplot(3, 2, 2*i-1);
    inst_amp = abs(hilbert(imf(i, :)));
    plot(t, inst_amp, 'r-', 'LineWidth', 1.2);
    xlabel('时间 (s)'); ylabel('瞬时幅值');
    title(sprintf('IMF%d 瞬时幅值', i));
    grid on;
    
    subplot(3, 2, 2*i);
    inst_freq = instantaneous_frequency(imf(i, :), fs);
    plot(t, inst_freq, 'b-', 'LineWidth', 1.2);
    xlabel('时间 (s)'); ylabel('瞬时频率 (Hz)');
    title(sprintf('IMF%d 瞬时频率', i));
    ylim([0, fs/2]);
    grid on;
end

2.2 经验模态分解函数 (emd_decomposition.m)

function [imf, residual] = emd_decomposition(x, fs)
%% =====================================================
%% 经验模态分解 (Empirical Mode Decomposition)
%% 输入:
%%   x  - 原始信号
%%   fs - 采样频率
%% 输出:
%%   imf      - 固有模态函数矩阵 (每一行为一个IMF)
%%   residual - 残差项
%% =====================================================

% 参数设置
max_iter = 1000;        % 最大迭代次数
stop_thresh = 0.05;    % 停止准则阈值
num_sift = 10;         % 筛选次数

N = length(x);
imf = [];
residual = x;

iter_count = 0;

while iter_count < max_iter
    iter_count = iter_count + 1;
    
    % 当前残差
    h = residual;
    
    % 筛选过程
    for sift_count = 1:num_sift
        % 寻找上下包络
        [upper_env, lower_env] = envelope_extraction(h);
        
        % 计算均值包络
        mean_env = (upper_env + lower_env) / 2;
        
        % 减去均值包络
        h_prev = h;
        h = h - mean_env;
        
        % 检查停止准则
        if check_stop_criterion(h, h_prev, stop_thresh)
            break;
        end
    end
    
    % 检查是否为有效的IMF
    if is_valid_imf(h)
        imf = [imf; h];
        residual = residual - h;
    else
        break;
    end
    
    % 检查残差是否足够小
    if norm(residual) < 0.01 * norm(x)
        break;
    end
end

% 添加最终残差
imf = [imf; residual];

end

%% 包络提取函数
function [upper_env, lower_env] = envelope_extraction(x)
N = length(x);

% 寻找极值点
[peaks, peak_locs] = find_peaks(x);
[troughs, trough_locs] = find_peaks(-x);
troughs = -troughs;

% 如果没有足够的极值点，使用插值
if length(peak_locs) < 3 || length(trough_locs) < 3
    upper_env = spline(1:N, x, 1:N);
    lower_env = spline(1:N, x, 1:N);
else
    % 三次样条插值
    upper_env = spline(peak_locs, peaks, 1:N);
    lower_env = spline(trough_locs, troughs, 1:N);
end
end

%% 寻找峰值函数
function [peaks, locs] = find_peaks(x)
% 寻找局部最大值
diff_x = diff(x);
idx = find(diff_x(1:end-1) > 0 & diff_x(2:end) < 0) + 1;
peaks = x(idx);
locs = idx;
end

%% 停止准则检查
function stop = check_stop_criterion(h, h_prev, thresh)
% 计算标准差准则
sd = sum((h - h_prev).^2) / sum(h_prev.^2);
stop = sd < thresh^2;
end

%% 检查是否为有效IMF
function valid = is_valid_imf(x)
% IMF应满足：极值点数与过零点数相差不超过1
[~, peak_locs] = find_peaks(x);
[~, trough_locs] = find_peaks(-x);

zero_crossings = sum(x(1:end-1) .* x(2:end) < 0);

valid = abs(length(peak_locs) - zero_crossings) <= 1 && ...
        abs(length(trough_locs) - zero_crossings) <= 1;
end

2.3 希尔伯特谱分析函数 (hilbert_spectrum_analysis.m)

function [H, f, t] = hilbert_spectrum_analysis(imf, fs)
%% =====================================================
%% 希尔伯特谱分析
%% 输入:
%%   imf - IMF分量矩阵
%%   fs  - 采样频率
%% 输出:
%%   H - 希尔伯特谱 (频率×时间)
%%   f - 频率向量
%%   t - 时间向量
%% =====================================================

[num_imf, N] = size(imf);

% 频率分辨率
freq_res = 0.1;
f_max = fs/2;
f = 0:freq_res:f_max;
num_freq = length(f);

t = (0:N-1)/fs;

% 初始化希尔伯特谱矩阵
H = zeros(num_freq, N);

% 对每个IMF进行希尔伯特变换
for i = 1:num_imf
    % 希尔伯特变换
    analytic_signal = hilbert(imf(i, :));
    
    % 瞬时幅值和相位
    inst_amplitude = abs(analytic_signal);
    inst_phase = unwrap(angle(analytic_signal));
    
    % 瞬时频率
    inst_freq = diff(inst_phase) * fs / (2*pi);
    inst_freq = [inst_freq(1), inst_freq];  % 补齐长度
    
    % 构建时频谱
    for n = 1:N
        freq_idx = round(inst_freq(n) / freq_res) + 1;
        if freq_idx >= 1 && freq_idx <= num_freq
            H(freq_idx, n) = H(freq_idx, n) + inst_amplitude(n);
        end
    end
end

% 归一化
H = H / max(H(:));
end

2.4 瞬时频率计算函数 (instantaneous_frequency.m)

function inst_freq = instantaneous_frequency(x, fs)
%% =====================================================
%% 计算瞬时频率
%% 输入:
%%   x  - 信号
%%   fs - 采样频率
%% 输出:
%%   inst_freq - 瞬时频率 (Hz)
%% =====================================================

N = length(x);

% 希尔伯特变换
analytic_signal = hilbert(x);

% 瞬时相位
inst_phase = unwrap(angle(analytic_signal));

% 瞬时频率 (Hz)
inst_freq = diff(inst_phase) * fs / (2*pi);

% 补齐长度
inst_freq = [inst_freq(1), inst_freq];

% 去除异常值
inst_freq(inst_freq < 0) = 0;
inst_freq(inst_freq > fs/2) = fs/2;
end

2.5 快速HHT函数（简化版）(fast_hht.m)

function [imf, H, f, t] = fast_hht(x, fs, varargin)
%% =====================================================
%% 快速希尔伯特黄变换（简化版）
%% 适用于快速原型开发和教学演示
%% =====================================================

% 参数解析
p = inputParser;
addParameter(p, 'MaxIMF', 10, @isnumeric);
addParameter(p, 'StopThresh', 0.05, @isnumeric);
parse(p, varargin{:});

max_imf = p.Results.MaxIMF;
stop_thresh = p.Results.StopThresh;

N = length(x);
t = (0:N-1)/fs;

% 简化的EMD分解
imf = zeros(max_imf, N);
residual = x;

for i = 1:max_imf
    % 简单的筛选过程
    h = residual;
    
    for iter = 1:10
        % 上下包络
        [upper_env, lower_env] = simple_envelope(h);
        mean_env = (upper_env + lower_env) / 2;
        h_new = h - mean_env;
        
        % 停止准则
        if norm(h_new - h) < stop_thresh * norm(h)
            break;
        end
        h = h_new;
    end
    
    imf(i, :) = h;
    residual = residual - h;
    
    if norm(residual) < 0.01 * norm(x)
        break;
    end
end

% 希尔伯特谱
H = zeros(fs/2+1, N);
f = 0:1:fs/2;

for i = 1:size(imf, 1)
    analytic_signal = hilbert(imf(i, :));
    inst_amp = abs(analytic_signal);
    inst_phase = unwrap(angle(analytic_signal));
    inst_freq = diff(inst_phase) * fs / (2*pi);
    inst_freq = [inst_freq(1), inst_freq];
    
    for n = 1:N
        freq_idx = round(inst_freq(n)) + 1;
        if freq_idx >= 1 && freq_idx <= length(f)
            H(freq_idx, n) = H(freq_idx, n) + inst_amp(n);
        end
    end
end

H = H / max(H(:));
end

function [upper_env, lower_env] = simple_envelope(x)
% 简单的包络提取
N = length(x);
upper_env = zeros(1, N);
lower_env = zeros(1, N);

win_size = 20;  % 滑动窗口大小

for i = 1:N
    start_idx = max(1, i - win_size);
    end_idx = min(N, i + win_size);
    window_data = x(start_idx:end_idx);
    
    upper_env(i) = max(window_data);
    lower_env(i) = min(window_data);
end
end

2.6 应用示例 (hht_example.m)

%% HHT应用示例：故障诊断
clear; clc; close all;

% 生成轴承故障信号（模拟）
fs = 12000;          % 采样频率
t = 0:1/fs:1;       % 1秒数据
N = length(t);

% 故障特征频率
fault_freq = 100;   % 故障频率 (Hz)
mod_freq = 5;       % 调制频率 (Hz)

% 故障信号模型
fault_signal = sin(2*pi*fault_freq*t) .* (1 + 0.5*sin(2*pi*mod_freq*t));
noise = 0.1 * randn(size(t));
x = fault_signal + noise;

% HHT分析
[imf, H, f, t_hht] = fast_hht(x, fs, 'MaxIMF', 5);

% 可视化
figure('Position', [100, 100, 1200, 800]);

subplot(3, 2, 1);
plot(t, x, 'b-', 'LineWidth', 1.2);
xlabel('时间 (s)'); ylabel('幅值');
title('轴承故障信号');
grid on;

subplot(3, 2, 2);
imagesc(t_hht, f, abs(H));
set(gca, 'YDir', 'normal');
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('希尔伯特时频谱');
colorbar;

% 边际谱
marginal_spec = mean(abs(H), 2);
subplot(3, 2, 3);
plot(f, marginal_spec, 'r-', 'LineWidth', 1.5);
xlabel('频率 (Hz)'); ylabel('幅值');
title('边际谱');
grid on;

% IMF分量
for i = 1:min(3, size(imf, 1))
    subplot(3, 2, 3+i);
    plot(t, imf(i, :), 'k-', 'LineWidth', 1.2);
    xlabel('时间 (s)'); ylabel('幅值');
    title(sprintf('IMF%d', i));
    grid on;
end

% 瞬时频率分析
subplot(3, 2, 6);
inst_freq = instantaneous_frequency(imf(1, :), fs);
plot(t, inst_freq, 'b-', 'LineWidth', 1.2);
xlabel('时间 (s)'); ylabel('瞬时频率 (Hz)');
title('IMF1瞬时频率');
ylim([0, 200]);
grid on;

fprintf('故障特征频率: %.2f Hz\n', fault_freq);
fprintf('检测到的峰值频率: %.2f Hz\n', f(find(marginal_spec == max(marginal_spec), 1)));

三、优化

3.1 计算加速

% 使用GPU加速（需要Parallel Computing Toolbox）
if gpuDeviceCount > 0
    x_gpu = gpuArray(x);
    % HHT计算...
    H = gather(H);
end

% 并行计算IMF
parfor i = 1:num_imf
    % 每个IMF独立计算
end

3.2 内存优化

% 使用稀疏矩阵存储时频谱
H_sparse = sparse(freq_indices, time_indices, amplitude_values);

参考代码 希尔伯特黄变换的MATLAB算法程序 www.youwenfan.com/contentcnu/60547.html

四、常见问题解决

问题	原因	解决方案
端点效应严重	信号两端不连续	添加镜像延拓或波形匹配延拓
模态混叠	频率相近的IMF混合	改进停止准则或使用EEMD
计算速度慢	EMD迭代次数过多	优化包络提取算法
虚假IMF	停止准则不严格	调整停止阈值