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

基于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);
        end
        
        imf = [imf, h];
        residual = residual - h;
        
        if is_monotonic(residual) || length(findpeaks(residual)) < 2
            break;
        end
    end
end

function is_mono = is_monotonic(x)
    % 判断信号是否单调
    diff_x = diff(x);
    is_mono = all(diff_x >= -1e-6) || all(diff_x <= 1e-6);
end

function [H, f, t] = hilbert_spectrum(imf, Fs)
    % 希尔伯特谱计算
    [P, f] = pburg(imf, [], [], [], Fs);  % Burg功率谱估计
    t = (0:length(imf)-1)/Fs;
    H = abs(hilbert(imf));                % 希尔伯特变换
end

---

说明

1. 信号预处理

   • 使用Chebyshev II滤波器（cheby2）进行带通滤波，去除高频噪声干扰。
   • 滤波参数可根据实际信号调整（wp, ws, rp, rs）。

2. 经验模态分解（EMD）

   • 通过循环筛选过程分离IMF分量，停止条件为标准差（sd < 0.3）或信号单调性。
   • 关键步骤包括包络线拟合（三次样条插值）和残差更新。

3. 希尔伯特变换

   • 使用hilbert函数计算解析信号，提取瞬时幅值和频率。
   • Burg功率谱估计（pburg）用于时频谱分析，平衡分辨率和计算效率。

4. 可视化模块

   • 分别展示原始信号、IMF分量及希尔伯特时频谱。
   • 时频谱采用伪彩色图（imagesc），横轴为时间，纵轴为频率，颜色表示能量强度。

   参考代码 希尔伯特变换(HHT)的 完整 MATLAB程序

---

*应用场景示

1. 机械故障诊断

   • 分析轴承振动信号中的冲击特征，通过IMF分量定位故障频率。

2. 生物医学信号处理

   • 提取心电信号（ECG）中的P波、QRS波群等特征，用于心律失常检测。

3. 地球物理勘探

   • 处理地震信号中的复杂波形，识别地层结构变化。