当前位置：首页 > news >正文

matlab实现希尔伯特变换(HHT)

news 2025/9/5 8:52:50

基于Hilbert-Huang变换（HHT）的完整MATLAB实现代码，包含经验模态分解（EMD）和希尔伯特谱分析，适用于非平稳信号处理。

HHT 代码

%% HHT变换主程序
% 输入参数
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);  % 合成信号（含噪声）%% 信号预处理（可选）
% 带通滤波去除高频噪声
wp = [40 160]/(Fs/2);  % 通带边界
ws = [30 170]/(Fs/2);  % 阻带边界
rp = 1; rs = 40;       % 通带/阻带波纹
[n,wn] = cheb2ord(wp,ws,rp,rs);  % 计算Chebyshev II滤波器阶数
[b,a] = cheby2(n,rp,wn);        % 设计滤波器
x_filtered = filter(b,a,x);     % 应用滤波器%% 经验模态分解（EMD）
imf = emd(x_filtered);  % 执行EMD分解
[~,num_imf] = size(imf);% 获取IMF数量%% 希尔伯特变换
[H, f, t_hilbert] = hilbert_spectrum(imf, Fs);  % 计算希尔伯特谱%% 可视化结果
figure;
subplot(3,1,1);
plot(t, x, 'b', t, x_filtered, 'r--');
title('原始信号与滤波后信号');
xlabel('时间(s)'); ylabel('幅值');
legend('原始信号', '滤波后信号');subplot(3,1,2);
plot(t, imf');
title('IMF分量分解结果');
xlabel('时间(s)'); ylabel('幅值');
legend(arrayfun(@(i) sprintf('IMF%d',i),1:num_imf,'UniformOutput',false));subplot(3,1,3);
imagesc(t_hilbert, f, abs(H));
title('希尔伯特时频谱');
xlabel('时间(s)'); ylabel('频率(Hz)');
colorbar;%% 自定义函数定义
function imf = emd(signal)% 经验模态分解（EMD）实现% 输入: signal - 原始信号% 输出: imf - 分解后的本征模态函数N = length(signal);imf = [];residual = signal;while true% 筛选过程h = residual;sd = Inf;while sd > 0.3 && ~is_monotonic(h)max_peaks = findpeaks(h);min_peaks = findpeaks(-h);min_peaks = -min_peaks(:,2);upper_env = interp1([0,N], [h(1), max_peaks, h(end)], 1:N, 'pchip');lower_env = interp1([0,N], [h(1), min_peaks, h(end)], 1:N, 'pchip');mean_env = (upper_env + lower_env)/2;h = h - mean_env;sd = sum((mean_env(2:end-1) - (upper_env(2:end-1) + lower_env(2:end-1))/2).^2) / ...sum((upper_env(2:end-1) + lower_env(2:end-1))/2.^2);endimf = [imf, h];residual = residual - h;if is_monotonic(residual) || length(findpeaks(residual)) < 2break;endend
endfunction is_mono = is_monotonic(x)% 判断信号是否单调diff_x = diff(x);is_mono = all(diff_x >= -1e-6) || all(diff_x <= 1e-6);
endfunction [H, f, t] = hilbert_spectrum(imf, Fs)% 希尔伯特谱计算[P, f] = pburg(imf, [], [], [], Fs);  % Burg功率谱估计t = (0:length(imf)-1)/Fs;H = abs(hilbert(imf));                % 希尔伯特变换
end

说明

信号预处理
- 使用Chebyshev II滤波器（cheby2）进行带通滤波，去除高频噪声干扰。
- 滤波参数可根据实际信号调整（wp, ws, rp, rs）。
经验模态分解（EMD）
- 通过循环筛选过程分离IMF分量，停止条件为标准差（sd < 0.3）或信号单调性。
- 关键步骤包括包络线拟合（三次样条插值）和残差更新。
希尔伯特变换
- 使用hilbert函数计算解析信号，提取瞬时幅值和频率。
- Burg功率谱估计（pburg）用于时频谱分析，平衡分辨率和计算效率。
可视化模块
- 分别展示原始信号、IMF分量及希尔伯特时频谱。
- 时频谱采用伪彩色图（imagesc），横轴为时间，纵轴为频率，颜色表示能量强度。
参考代码希尔伯特变换(HHT)的完整 MATLAB程序 www.youwenfan.com/contentcsf/80149.html