希尔伯特黄变换（HHT）MATLAB程序
希尔伯特黄变换（HHT）MATLAB程序。HHT包含两个核心步骤：经验模态分解（EMD） 和希尔伯特谱分析（HSA），能够有效分析非线性、非平稳信号。
HHT 算法完整实现
以下是两个版本的实现，第一个是简洁的基础版，第二个是功能更完善的增强版。
版本1：基础简洁版HHT程序
%% 希尔伯特黄变换 (HHT) 基础实现
clear; clc; close all;
%% 1. 生成测试信号
fs = 1000; % 采样频率
t = 0:1/fs:1-1/fs; % 时间向量
% 组合信号：包含多个频率分量
signal = 2*sin(2*pi*5*t) + ... % 5Hz低频
0.5*sin(2*pi*50*t) + ... % 50Hz中频
0.2*sin(2*pi*120*t) + ... % 120Hz高频
0.1*randn(size(t)); % 噪声
%% 2. 经验模态分解 (EMD) - 获取本征模态函数IMF
[imf, residue] = emd_basic(signal);
num_imf = size(imf, 1);
%% 3. 希尔伯特谱分析 - 计算各IMF的瞬时频率和幅值
[hht_time, hht_freq, hht_amp, hilbert_spectrum] = hilbert_spectrum(imf, fs);
%% 4. 可视化结果
figure('Position', [100, 100, 1200, 800]);
% 4.1 原始信号
subplot(3, 2, 1);
plot(t, signal, 'b', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('幅值');
title('原始信号'); grid on;
% 4.2 IMF分量
subplot(3, 2, 2);
for i = 1:num_imf
plot(t, imf(i,:) + (num_imf-i)*0.5, 'LineWidth', 1);
hold on;
end
xlabel('时间 (s)'); ylabel('IMF分量');
title('经验模态分解结果'); grid on;
hold off;
% 4.3 希尔伯特谱（时频谱）
subplot(3, 2, [3, 4]);
imagesc(t, hht_freq, 20*log10(hht_amp+eps));
axis xy; colorbar;
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('希尔伯特谱 (Hilbert Spectrum)');
ylim([0, fs/2]);
% 4.4 边际谱（频率分布）
subplot(3, 2, 5);
marginal_spectrum = sum(hht_amp, 2);
plot(hht_freq, marginal_spectrum, 'r', 'LineWidth', 2);
xlabel('频率 (Hz)'); ylabel('幅值');
title('边际谱 (Marginal Spectrum)');
xlim([0, fs/2]); grid on;
% 4.5 各IMF瞬时频率
subplot(3, 2, 6);
colors = jet(num_imf);
for i = 1:num_imf
plot(t, squeeze(hilbert_spectrum(2,:,i)), 'Color', colors(i,:), 'LineWidth', 1.5);
hold on;
end
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('各IMF瞬时频率');
legend(arrayfun(@(x) sprintf('IMF%d', x), 1:num_imf, 'UniformOutput', false));
grid on; hold off;
sgtitle('希尔伯特黄变换 (HHT) 分析结果', 'FontSize', 14, 'FontWeight', 'bold');
%% 5. 基础EMD函数实现
function [imf, residue] = emd_basic(signal, max_imf, max_iter)
% signal: 输入信号
% max_imf: 最大IMF数量（默认10）
% max_iter: 每次筛分最大迭代次数（默认1000）
if nargin < 2, max_imf = 10; end
if nargin < 3, max_iter = 1000; end
signal = signal(:)';
residue = signal;
imf = [];
for imf_idx = 1:max_imf
h = residue;
for iter = 1:max_iter
% 找到极值点
[max_peaks, max_locs] = findpeaks(h);
[min_peaks, min_locs] = findpeaks(-h);
if length(max_peaks) < 3 || length(min_peaks) < 3
break;
end
% 插值得到上下包络
t = 1:length(h);
max_env = interp1(max_locs, max_peaks, t, 'spline', 'extrap');
min_env = interp1(min_locs, -min_peaks, t, 'spline', 'extrap');
% 计算均值包络
mean_env = (max_env + min_env) / 2;
% 更新h
h_prev = h;
h = h - mean_env;
% 停止准则：均值包络接近零
if sum(abs(mean_env)) / sum(abs(h_prev)) < 0.1
break;
end
end
% 保存IMF，更新残差
imf = [imf; h];
residue = residue - h;
% 停止准则：残差为单调函数或极值点太少
if length(findpeaks(residue)) < 2 || length(findpeaks(-residue)) < 2
break;
end
end
end
%% 6. 希尔伯特谱分析函数
function [time_vec, freq_vec, amp_mat, hilbert_data] = hilbert_spectrum(imf, fs)
% imf: IMF矩阵 (每行一个IMF)
% fs: 采样频率
% 返回: 时间向量、频率向量、幅值矩阵、希尔伯特数据
[num_imf, sig_len] = size(imf);
time_vec = (0:sig_len-1) / fs;
% 初始化
hilbert_data = zeros(3, sig_len, num_imf); % 幅值、频率、相位
amp_mat = zeros(round(fs/2)+1, sig_len);
% 频率向量
freq_vec = 0:0.5:fs/2;
for i = 1:num_imf
% 希尔伯特变换
analytic_signal = hilbert(imf(i,:));
% 计算瞬时幅值、相位、频率
inst_amplitude = abs(analytic_signal);
inst_phase = unwrap(angle(analytic_signal));
inst_frequency = diff(inst_phase) * fs / (2*pi);
inst_frequency = [inst_frequency(1), inst_frequency]; % 保持长度
% 存储数据
hilbert_data(1,:,i) = inst_amplitude;
hilbert_data(2,:,i) = inst_frequency;
hilbert_data(3,:,i) = inst_phase;
% 构建希尔伯特谱
for t = 1:sig_len
freq_idx = round(inst_frequency(t) * 2) + 1; % 0.5Hz分辨率
if freq_idx >= 1 && freq_idx <= length(freq_vec)
amp_mat(freq_idx, t) = amp_mat(freq_idx, t) + inst_amplitude(t);
end
end
end
end
版本2：增强版HHT程序（推荐用于实际分析）
%% 希尔伯特黄变换 (HHT) 增强实现 - 包含更多功能和优化
function hht_enhanced_demo()
clear; clc; close all;
%% 1. 生成更复杂的测试信号
fs = 2000; % 采样频率
t = 0:1/fs:2-1/fs; % 2秒时长
% 非线性调频信号 + 冲击 + 噪声
chirp_signal = chirp(t, 10, 2, 100); % 线性调频 10-100Hz
impulse_signal = zeros(size(t));
impulse_signal(500:520) = 2; % 局部冲击
impulse_signal(1500:1520) = 1.5;
signal = 1.5*sin(2*pi*20*t) + ... % 20Hz正弦
0.8*chirp_signal + ... % 调频成分
0.6*impulse_signal + ... % 冲击成分
0.15*randn(size(t)); % 高斯白噪声
%% 2. 使用改进的EMD分解
[imf, residue, info] = emd_enhanced(signal, 'MaxIMF', 8, 'MaxIterations', 500);
num_imf = size(imf, 1);
%% 3. 计算HHT谱（使用更高分辨率）
[hht_time, hht_freq, hht_amp, inst_freq, inst_amp] = ...
compute_hht_spectrum(imf, fs, 'FrequencyResolution', 1);
%% 4. 综合可视化分析
visualize_hht_results(t, signal, imf, hht_time, hht_freq, hht_amp, ...
inst_freq, inst_amp, fs);
%% 5. 信号重构与误差分析
reconstructed = sum(imf, 1) + residue;
reconstruction_error = rms(signal - reconstructed);
fprintf('信号重构均方根误差: %.6f\n', reconstruction_error);
%% 6. 计算各IMF的统计特征
compute_imf_statistics(imf, fs);
end
%% 改进的EMD实现（处理端点效应和模态混叠）
function [imf, residue, info] = emd_enhanced(signal, varargin)
% 解析输入参数
p = inputParser;
addParameter(p, 'MaxIMF', 10, @isnumeric);
addParameter(p, 'MaxIterations', 1000, @isnumeric);
addParameter(p, 'StopThreshold', 0.1, @isnumeric);
addParameter(p, 'BoundaryTreatment', 'mirror', @ischar);
parse(p, varargin{:});
signal = signal(:)';
N = length(signal);
residue = signal;
imf = [];
info.num_iterations = [];
info.stop_criteria = [];
% 镜像延拓处理端点效应
if strcmp(p.Results.BoundaryTreatment, 'mirror')
ext_len = round(N/10);
residue_ext = [fliplr(residue(1:ext_len)), residue, fliplr(residue(end-ext_len+1:end))];
else
residue_ext = residue;
end
for imf_idx = 1:p.Results.MaxIMF
h = residue_ext;
prev_h = h;
sift_num = 0;
for iter = 1:p.Results.MaxIterations
% 寻找极值点（使用更稳健的方法）
[max_vals, max_pos] = find_local_maxima(h);
[min_vals, min_pos] = find_local_maxima(-h);
if length(max_vals) < 3 || length(min_vals) < 3
info.stop_criteria{imf_idx} = 'Insufficient extrema';
break;
end
% 三次样条插值包络
t_ext = 1:length(h);
if length(max_pos) >= 4
max_env = spline(max_pos, max_vals, t_ext);
min_env = spline(min_pos, min_vals, t_ext);
else
max_env = interp1(max_pos, max_vals, t_ext, 'linear', 'extrap');
min_env = interp1(min_pos, min_vals, t_ext, 'linear', 'extrap');
end
mean_env = (max_env + min_env) / 2;
% 更新分量
prev_h = h;
h = h - mean_env;
% 改进的停止准则
sd = sum((prev_h - h).^2) / sum(prev_h.^2);
sift_num = sift_num + 1;
if sd < p.Results.StopThreshold || iter == p.Results.MaxIterations
info.num_iterations(imf_idx) = sift_num;
% 提取有效的IMF（去掉延拓部分）
if strcmp(p.Results.BoundaryTreatment, 'mirror')
valid_imf = h(ext_len+1:ext_len+N);
else
valid_imf = h;
end
imf = [imf; valid_imf];
residue = residue - valid_imf;
% 更新延拓后的残差
if strcmp(p.Results.BoundaryTreatment, 'mirror')
residue_ext = [fliplr(residue(1:ext_len)), residue, fliplr(residue(end-ext_len+1:end))];
else
residue_ext = residue;
end
break;
end
end
% 全局停止准则
if is_monotonic(residue) || length(find_local_maxima(residue)) < 2
info.stop_criteria{imf_idx} = 'Residue monotonic';
break;
end
end
info.total_imf = size(imf, 1);
end
%% 寻找局部极值点（改进版本）
function [vals, positions] = find_local_maxima(x)
x = x(:)';
N = length(x);
positions = [];
vals = [];
for i = 2:N-1
if x(i) > x(i-1) && x(i) > x(i+1)
positions = [positions, i];
vals = [vals, x(i)];
end
end
% 考虑端点
if x(1) > x(2)
positions = [1, positions];
vals = [x(1), vals];
end
if x(N) > x(N-1)
positions = [positions, N];
vals = [vals, x(N)];
end
end
%% 检查序列是否单调
function result = is_monotonic(x)
diff_x = diff(x);
result = all(diff_x >= 0) || all(diff_x <= 0);
end
%% 计算HHT谱（高分辨率）
function [time_vec, freq_vec, hht_matrix, inst_freq, inst_amp] = ...
compute_hht_spectrum(imf, fs, varargin)
p = inputParser;
addParameter(p, 'FrequencyResolution', 0.5, @isnumeric); % Hz
parse(p, varargin{:});
[num_imf, sig_len] = size(imf);
time_vec = (0:sig_len-1) / fs;
% 设置频率范围和分辨率
freq_res = p.Results.FrequencyResolution;
freq_vec = 0:freq_res:fs/2;
num_freq_bins = length(freq_vec);
% 初始化
hht_matrix = zeros(num_freq_bins, sig_len);
inst_freq = zeros(num_imf, sig_len);
inst_amp = zeros(num_imf, sig_len);
for i = 1:num_imf
% 希尔伯特变换
analytic = hilbert(imf(i,:));
amplitude = abs(analytic);
phase = unwrap(angle(analytic));
% 瞬时频率（使用中心差分提高精度）
inst_freq(i,2:end-1) = fs/(4*pi) * (phase(3:end) - phase(1:end-2));
inst_freq(i,1) = inst_freq(i,2);
inst_freq(i,end) = inst_freq(i,end-1);
inst_freq(i, inst_freq(i,:) < 0) = 0;
inst_freq(i, inst_freq(i,:) > fs/2) = fs/2;
inst_amp(i,:) = amplitude;
% 填充希尔伯特谱矩阵
for t = 1:sig_len
freq_bin = round(inst_freq(i,t) / freq_res) + 1;
if freq_bin >= 1 && freq_bin <= num_freq_bins
hht_matrix(freq_bin, t) = hht_matrix(freq_bin, t) + amplitude(t);
end
end
end
end
%% 综合可视化函数
function visualize_hht_results(t, signal, imf, hht_time, hht_freq, hht_amp, ...
inst_freq, inst_amp, fs)
num_imf = size(imf, 1);
figure('Position', [50, 50, 1400, 900]);
% 1. 原始信号与IMFs
subplot(4, 3, [1, 2]);
plot(t, signal, 'b', 'LineWidth', 1.8);
xlabel('时间 (s)'); ylabel('幅值');
title('原始信号', 'FontSize', 12, 'FontWeight', 'bold');
grid on; box on;
% 2. EMD分解结果
subplot(4, 3, 3);
offset = 0;
colors = lines(num_imf);
for i = 1:num_imf
plot(t, imf(i,:) + offset, 'Color', colors(i,:), 'LineWidth', 1.2);
hold on;
offset = offset + max(abs(imf(i,:))) * 1.5;
end
xlabel('时间 (s)'); ylabel('IMF分量');
title('经验模态分解 (EMD)', 'FontSize', 12, 'FontWeight', 'bold');
grid on; hold off;
% 3. 希尔伯特谱（等高线图）
subplot(4, 3, [4, 5, 6]);
contourf(hht_time, hht_freq, 20*log10(hht_amp+eps), 30, 'LineColor', 'none');
colorbar; colormap(jet);
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('希尔伯特谱 (Hilbert Spectrum)', 'FontSize', 12, 'FontWeight', 'bold');
ylim([0, fs/4]); % 只显示低频部分
% 4. 边际谱（累积频率分布）
subplot(4, 3, 7);
marginal = sum(hht_amp, 2);
plot(hht_freq, marginal, 'r', 'LineWidth', 2);
xlabel('频率 (Hz)'); ylabel('能量');
title('边际谱', 'FontSize', 12);
grid on; xlim([0, fs/4]);
% 5. 瞬时频率轨迹
subplot(4, 3, 8);
for i = 1:min(num_imf, 5) % 显示前5个IMF
plot(t, inst_freq(i,:), 'LineWidth', 1.5);
hold on;
end
xlabel('时间 (s)'); ylabel('瞬时频率 (Hz)');
title('瞬时频率变化', 'FontSize', 12);
legend(arrayfun(@(x) sprintf('IMF%d', x), 1:min(num_imf,5), 'UniformOutput', false));
grid on; hold off;
% 6. 能量分布饼图
subplot(4, 3, 9);
imf_energy = sum(imf.^2, 2);
pie(imf_energy, arrayfun(@(x) sprintf('IMF%d', x), 1:num_imf, 'UniformOutput', false));
title('各IMF能量占比', 'FontSize', 12);
% 7. 时频联合分布（3D视图）
subplot(4, 3, [10, 11, 12]);
[T, F] = meshgrid(hht_time, hht_freq(1:2:end));
Z = hht_amp(1:2:end, :);
surf(T, F, Z, 'EdgeColor', 'none');
view(45, 30);
xlabel('时间 (s)'); ylabel('频率 (Hz)'); zlabel('幅值');
title('时频分布 3D 视图', 'FontSize', 12);
colormap(jet); colorbar;
sgtitle('希尔伯特黄变换 (HHT) 综合分析报告', 'FontSize', 16, 'FontWeight', 'bold');
end
%% 计算IMF统计特征
function compute_imf_statistics(imf, fs)
num_imf = size(imf, 1);
fprintf('\n======= IMF统计特征 =======\n');
fprintf('%-8s %-12s %-12s %-12s %-12s\n', ...
'IMF', '均值', '方差', '主频(Hz)', '能量占比(%)');
fprintf('%s\n', repmat('-', 60, 1));
total_energy = sum(imf(:).^2);
for i = 1:num_imf
% 基本统计
mean_val = mean(imf(i,:));
var_val = var(imf(i,:));
% 频谱分析求主频
[Pxx, f] = periodogram(imf(i,:), [], length(imf(i,:)), fs);
[~, idx] = max(Pxx);
main_freq = f(idx);
% 能量占比
imf_energy = sum(imf(i,:).^2);
energy_ratio = 100 * imf_energy / total_energy;
fprintf('%-8d %-12.4f %-12.4f %-12.2f %-12.2f\n', ...
i, mean_val, var_val, main_freq, energy_ratio);
end
end
算法核心解析与使用
1. HHT 的两个核心步骤
| 步骤 | 原理 | 关键点 |
|---|---|---|
| EMD | 自适应地将信号分解为多个本征模态函数 (IMF) | 1. 识别信号极值点 2. 插值得到上下包络线 3. 迭代筛分直到满足IMF条件 |
| 希尔伯特谱分析 | 对每个IMF进行希尔伯特变换得到瞬时频率和幅值 | 1. 计算解析信号 2. 提取瞬时频率和幅值 3. 组合成时频分布 |
2. 程序使用说明
% 快速使用示例
fs = 1000;
t = 0:1/fs:1;
x = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t);
% 使用基础版
[imf, residue] = emd_basic(x);
[hht_time, hht_freq, hht_amp] = hilbert_spectrum(imf, fs);
% 使用增强版（推荐）
hht_enhanced_demo(); % 运行增强版演示
3. 参数调整建议
| 参数 | 推荐值 | 说明 |
|---|---|---|
| 最大IMF数 | 8-10 | 过多会导致过分解 |
| 筛分停止阈值 | 0.1-0.3 | 控制分解精度 |
| 频率分辨率 | 0.5-2 Hz | 影响希尔伯特谱精度 |
| 端点处理 | 'mirror' | 减少端点效应 |
4. 专业工具箱推荐
对于严肃的科研工作，建议使用成熟的EMD工具箱：
-
G. Rilling的EMD工具箱（最经典）：
% 下载地址：https://perso.ens-lyon.fr/patrick.flandrin/emd.html imf = emd(signal); % 使用工具箱函数 -
MATLAB官方信号处理工具箱（2018a+）：
[imf, residue, info] = emd(signal, 'Interpolation', 'pchip');
参考代码 希尔伯特黄变换(HHT)的 完整 MATLAB程序 www.3dddown.com/cna/64242.html
注意事项与常见问题
-
端点效应处理
- 基础版本存在端点效应，增强版提供了镜像延拓
- 实际应用中可考虑使用极值点延拓或神经网络预测
-
模态混叠问题
- 当信号频率成分过于接近时会发生模态混叠
- 解决方法：使用集合经验模态分解 (EEMD) 或互补集合经验模态分解 (CEEMD)
-
计算复杂度
- EMD的计算量较大，长信号处理较慢
- 可考虑分段处理或使用并行计算
-
结果验证
% 验证IMF的正交性 orthogonality_index = zeros(size(imf,1)); for i = 1:size(imf,1) for j = i+1:size(imf,1) orthogonality_index(i,j) = abs(sum(imf(i,:).*imf(j,:))) / ... (norm(imf(i,:))*norm(imf(j,:))); end end
应用场景扩展
HHT特别适合以下应用：
- 机械故障诊断：轴承、齿轮振动信号分析
- 生物医学信号：EEG、ECG的非平稳特征提取
- 地震信号分析：地震波时频特性研究
- 金融时间序列：股票价格波动模式分析
浙公网安备 33010602011771号