一、核心代码框架
%% HHT完整实现（含EMD分解与希尔伯特谱分析）clear;clc;close all;%% 参数设置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)+0.1*randn(size(t));% 含噪信号%% 信号预处理（带通滤波）wp=[40160]/(Fs/2);% 通带边界ws=[30170]/(Fs/2);% 阻带边界rp=1;rs=40;% 滤波器参数[n,wn]=cheb2ord(wp,ws,rp,rs);[b,a]=cheby2(n,rp,wn);x_filt=filter(b,a,x);%% 经验模态分解（EMD）imf=emd(x_filt,'MaxNumIMF',5,'Display',1);% 分解为5个IMF分量%% 希尔伯特变换analytic=cell(size(imf));fori=1:size(imf,1)analytic{i}=hilbert(imf(i,:));% 构造解析信号end%% 瞬时频率与幅值计算[inst_amp,inst_phase]=deal(zeros(size(imf)));fori=1:size(imf,1)phase=angle(analytic{i});inst_phase(:,i)=unwrap(phase);% 相位解包裹inst_amp(:,i)=abs(analytic{i});% 瞬时幅值endinst_freq=zeros(size(imf));fori=1:size(imf,1)delta_t=diff(t);delta_phase=diff(inst_phase(:,i));inst_freq(:,i)=[inst_phase(1,i);delta_phase./delta_t];% 瞬时频率end%% 三维时频谱构建[~,f_idx]=min(abs(fft_freq(Fs,1/Fs)-f1:f2:end));% 频率索引t_res=0.01;% 时间分辨率Z=zeros(length(t_res),length(f_idx));fori=1:size(imf,1)[t_f,f_f]=tfridge(inst_amp(:,i),inst_freq(:,i),Fs);Z(:,f_idx)=Z(:,f_idx)+interp1(t_f,f_f,t,'linear',0);end%% 可视化figure;% IMF分量展示subplot(2,2,1);plot(t,x_filt,'k',t,imf','LineWidth',1.5);title('原始信号与IMF分量');xlabel('时间(s)');ylabel('幅值');legend('原始信号','IMF1','IMF2','IMF3','IMF4','IMF5');% 瞬时频率谱subplot(2,2,2);imagesc(t,1:size(imf,1),inst_freq);set(gca,'YDir','normal');title('瞬时频率分布');xlabel('时间(s)');ylabel('IMF序号');colorbar;% 三维时频谱subplot(2,2,3);surf(linspace(0,1,length(t)),linspace(1,size(imf,1),size(imf,1)),Z);shading interp;title('三维时频谱');xlabel('时间(s)');ylabel('IMF序号');zlabel('幅值');% Hilbert边际谱subplot(2,2,4);hs=hht(imf,Fs);imagesc(hs.time,hs.f,hs.power);title('Hilbert边际谱');xlabel('时间(s)');ylabel('频率(Hz)');colorbar;二、关键函数实现
1. EMD分解函数（增强版）
functionimfs=emd(signal,varargin)% 参数解析p=inputParser;addParameter(p,'MaxNumIMF',10,@(x)isscalar(x));addParameter(p,'SiftRelativeTol',0.1,@(x)isscalar(x));addParameter(p,'Display',0,@(x)islogical(x));parse(p,varargin{:});imfs={};residual=signal;num_imf=0;whiletrue num_imf=num_imf+1;imf_candidate=residual;% 极值点检测[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;% 筛选IMFimf_candidate=residual-mean_env;% IMF条件验证ifis_imf(imf_candidate,p.SiftRelativeTol)imfs{num_imf}=imf_candidate;residual=residual-imf_candidate;ifnorm(residual)<0.1*norm(signal)break;endelseresidual=imf_candidate;end% 显示进度ifp.Displayfprintf('IMF%d: 残差能量比=%.4f\n',num_imf,norm(residual)/norm(signal));endendendfunction[max_peaks,min_peaks]=find_extrema(x)% 极值点检测n=length(x);max_peaks=[];min_peaks=[];fori=2:n-1ifx(i)>x(i-1)&&x(i)>x(i+1)max_peaks=[max_peaks;i,x(i)];elseifx(i)<x(i-1)&&x(i)<x(i+1)min_peaks=[min_peaks;i,x(i)];endendendfunctionis_imf=is_imf(imf,tol)% IMF条件验证diff_imf=diff(imf);num_extrema=sum(diff_imf(1:end-1).*diff_imf(2:end)<=0);mean_diff=mean(abs(diff_imf));is_imf=(num_extrema-1)<=1&&mean_diff<tol;end三、核心参数优化
| 参数 | 推荐值 | 作用说明 |
|---|---|---|
MaxNumIMF | 5-10 | 控制分解深度与计算量 |
SiftRelativeTol | 0.1-0.3 | 分解精度控制 |
Display | 1 (True) | 显示分解过程 |
| 滤波器类型 | Chebyshev II | 高频噪声抑制 |
参考代码 希尔伯特黄变换(HHT)的 完整 MATLAB程序www.youwenfan.com/contentcsq/64242.html
四、工程应用扩展
1. 模态混叠抑制（EEMD改进）
functionimfs=eemd(signal,noise_level)% 集成经验模态分解imfs={};residual=signal;fori=1:10noise=noise_level*randn(size(signal));imf=emd(residual+noise);residual=residual-imf{end};imfs{end+1}=imf{end};endend2. 边界效应处理（镜像延拓）
functionx_pad=boundary_extension(x,n_pad)% 镜像延拓x_pad=[flipud(x(1:n_pad));x;flipud(x(end-n_pad+1:end))];end五、结果分析示例
IMF分量特征
IMF1：高频噪声（0-30Hz）
IMF2：主频50Hz正弦波
残余分量：低频趋势项（约0.1Hz）
三维时频谱解读
X轴：时间（0-1s）
Y轴：IMF序号
Z轴：幅值强度
颜色深度：能量密度
六、注意事项
数据预处理：建议先进行去趋势项处理（
detrend函数）频率校准：使用
resample函数统一采样率实时处理：分段处理时需重叠50%以上（参考Hanning窗）
七、参考文献
Huang N E, et al. Hilbert-Huang Transform and Its Applications. World Scientific, 2005.
王明阳, 等. 基于HHT的机械故障诊断方法. 振动工程学报, 2010.
MATLAB官方HHT工具箱文档（R2023a）
