EMD分解与希尔伯特变换能量谱分析

如何对信号进行经验模态分解(EMD)，然后对各个本征模态函数(IMF)进行希尔伯特变换，最终获得能量谱的完整MATLAB实现。

算法原理概述

EMD + Hilbert变换流程

原始信号 → EMD分解 → 多个IMF分量 + 残差 → 对每个IMF进行Hilbert变换 → 解析信号 → 计算瞬时频率和瞬时幅度 → 构建Hilbert谱 → 计算能量谱密度

MATLAB完整实现

1. 主程序：EMD-Hilbert能量谱分析

function[hht_spectrum,energy_density,imfs]=emd_hilbert_energy_analysis(signal,fs,varargin)% EMD分解与希尔伯特变换能量谱分析% 输入:% signal - 输入信号% fs - 采样频率% varargin - 可选参数% 输出:% hht_spectrum - Hilbert谱% energy_density - 能量谱密度% imfs - IMF分量%% 参数设置p=inputParser;addParameter(p,'num_imfs',10,@isnumeric);% 最大IMF数量addParameter(p,'max_iterations',1000,@isnumeric);% 最大迭代次数addParameter(p,'stop_threshold',0.05,@isnumeric);% 停止阈值addParameter(p,'plot_results',true,@islogical);% 是否绘图addParameter(p,'energy_band',[0,fs/2],@isnumeric);% 能量谱频带parse(p,varargin{:});params=p.Results;fprintf('开始EMD-Hilbert能量谱分析...\n');fprintf('信号长度: %d, 采样率: %.1f Hz\n',length(signal),fs);%% 步骤1: EMD分解tic;fprintf('正在进行EMD分解...\n');imfs=emd_decomposition(signal,...'max_imfs',params.num_imfs,...'max_iterations',params.max_iterations,...'stop_threshold',params.stop_threshold);decomposition_time=toc;fprintf('EMD分解完成! 得到 %d 个IMF分量, 耗时: %.2f 秒\n',size(imfs,1),decomposition_time);%% 步骤2: 希尔伯特变换与能量计算fprintf('正在进行希尔伯特变换...\n');[hht_spectrum,instantaneous_freq,instantaneous_amp]=...hilbert_spectrum_analysis(imfs,fs);%% 步骤3: 能量谱计算fprintf('计算能量谱密度...\n');energy_density=calculate_energy_spectrum(hht_spectrum,...instantaneous_freq,fs,params.energy_band);%% 步骤4: 结果显示ifparams.plot_resultsplot_emd_hilbert_results(signal,fs,imfs,hht_spectrum,...instantaneous_freq,energy_density);endfprintf('EMD-Hilbert能量谱分析完成!\n');end

2. EMD分解核心函数

functionimfs=emd_decomposition(signal,varargin)% EMD经验模态分解% 实现经典的EMD算法将信号分解为IMF分量% 参数解析p=inputParser;addParameter(p,'max_imfs',10,@isnumeric);addParameter(p,'max_iterations',1000,@isnumeric);addParameter(p,'stop_threshold',0.05,@isnumeric);parse(p,varargin{:});max_imfs=p.Results.max_imfs;max_iterations=p.Results.max_iterations;stop_threshold=p.Results.stop_threshold;signal=signal(:);% 确保列向量N=length(signal);residual=signal;imfs=[];forimf_idx=1:max_imfsfprintf(' 提取IMF %d...',imf_idx);h=residual;% 当前待提取的IMFiteration_count=0;whileiteration_count<max_iterations iteration_count=iteration_count+1;% 寻找局部极值点[max_locs,min_locs]=find_local_extrema(h);% 检查是否满足IMF条件iflength(max_locs)<2||length(min_locs)<2fprintf('极值点不足，停止分解\n');break;end% 计算上下包络线upper_env=compute_envelope(h,max_locs,'spline');lower_env=compute_envelope(h,min_locs,'spline');% 计算均值包络线mean_env=(upper_env+lower_env)/2;% 更新hh_prev=h;h=h-mean_env;% 检查停止条件stop_value=compute_stop_criteria(h_prev,h,mean_env);ifstop_value<stop_thresholdbreak;endend% 保存IMFimfs=[imfs;h'];residual=residual-h;fprintf('完成 (迭代次数: %d)\n',iteration_count);% 检查残差是否足够小ifsum(residual.^2)<1e-10*sum(signal.^2)||...length(find_local_extrema(residual))<2fprintf('残差已足够小，分解完成\n');break;endend% 添加残差作为最后一个分量if~isempty(residual)imfs=[imfs;residual'];endendfunction[max_locs,min_locs]=find_local_extrema(x)% 寻找局部极值点diff_x=diff(x);sign_changes=diff(sign(diff_x));% 极大值点 (二阶导数为负)max_locs=find(sign_changes<0)+1;% 极小值点 (二阶导数为正)min_locs=find(sign_changes>0)+1;% 确保端点被正确处理ifx(1)>x(2)max_locs=[1;max_locs];elsemin_locs=[1;min_locs];endifx(end)>x(end-1)max_locs=[max_locs;length(x)];elsemin_locs=[min_locs;length(x)];endendfunctionenvelope=compute_envelope(x,extr_locs,method)% 计算包络线t=1:length(x);ifstrcmpi(method,'spline')% 三次样条插值iflength(extr_locs)>=4envelope=spline(extr_locs,x(extr_locs),t);else% 极值点太少，使用线性插值envelope=interp1(extr_locs,x(extr_locs),t,'linear','extrap');endelse% 线性插值envelope=interp1(extr_locs,x(extr_locs),t,'linear','extrap');endendfunctionstop_value=compute_stop_criteria(h_prev,h,mean_env)% 计算停止准则% 使用基于标准差的方法numerator=sum((h_prev-h).^2);denominator=sum(h_prev.^2);ifdenominator==0stop_value=0;elsestop_value=numerator/denominator;endend

3. 希尔伯特变换与谱分析

function[hht_spectrum,instantaneous_freq,instantaneous_amp]=...hilbert_spectrum_analysis(imfs,fs)% 希尔伯特谱分析% 对每个IMF进行希尔伯特变换，计算瞬时频率和幅度[num_imfs,N]=size(imfs);t=(0:N-1)/fs;% 初始化输出变量instantaneous_freq=zeros(num_imfs,N);instantaneous_amp=zeros(num_imfs,N);hht_spectrum=zeros(num_imfs,N);fori=1:num_imfsfprintf(' 处理IMF %d 的希尔伯特变换...\n',i);current_imf=imfs(i,:);% 希尔伯特变换获得解析信号analytic_signal=hilbert(current_imf);% 计算瞬时幅度instantaneous_amp(i,:)=abs(analytic_signal);% 计算瞬时相位（解卷绕）instantaneous_phase=unwrap(angle(analytic_signal));% 计算瞬时频率 (导数方法)instantaneous_freq(i,:)=diff([instantaneous_phase(1),instantaneous_phase])*fs/(2*pi);% 限制频率范围为正值instantaneous_freq(i,instantaneous_freq(i,:)<0)=0;instantaneous_freq(i,instantaneous_freq(i,:)>fs/2)=fs/2;% 构建Hilbert谱 (幅度平方表示能量)hht_spectrum(i,:)=instantaneous_amp(i,:).^2;endend

4. 能量谱计算函数

functionenergy_density=calculate_energy_spectrum(hht_spectrum,...instantaneous_freq,fs,energy_band)% 计算能量谱密度[num_imfs,N]=size(hht_spectrum);% 频率分辨率设置freq_resolution=1;% Hzfreq_bins=energy_band(1):freq_resolution:energy_band(2);% 初始化能量谱矩阵energy_density=zeros(num_imfs,length(freq_bins)-1);fori=1:num_imfsfprintf(' 计算IMF %d 的能量谱...\n',i);fort=1:N current_freq=instantaneous_freq(i,t);current_energy=hht_spectrum(i,t);% 找到对应的频率区间freq_bin_idx=find(freq_bins<=current_freq,1,'last');if~isempty(freq_bin_idx)&&freq_bin_idx<length(freq_bins)energy_density(i,freq_bin_idx)=energy_density(i,freq_bin_idx)+current_energy;endend% 归一化能量谱ifsum(energy_density(i,:))>0energy_density(i,:)=energy_density(i,:)/sum(energy_density(i,:));endend% 计算总能量谱total_energy_density=sum(energy_density,1);ifsum(total_energy_density)>0total_energy_density=total_energy_density/sum(total_energy_density);endenergy_density=[energy_density;total_energy_density];end

5. 结果可视化函数

functionplot_emd_hilbert_results(signal,fs,imfs,hht_spectrum,...instantaneous_freq,energy_density)% 绘制EMD-Hilbert分析结果N=length(signal);t=(0:N-1)/fs;num_imfs=size(imfs,1);% 创建图形窗口figure('Position',[100,100,1400,1000]);%% 子图1: 原始信号subplot(4,3,1);plot(t,signal,'b','LineWidth',1.5);title('原始信号','FontSize',12,'FontWeight','bold');xlabel('时间 (s)');ylabel('幅值');grid on;%% 子图2: EMD分解结果subplot(4,3,[2,3]);plot_offset=0;colors=parula(num_imfs);fori=1:num_imfsplot(t,imfs(i,:)+plot_offset,'Color',colors(i,:),'LineWidth',1.2);hold on;text(t(end)+0.02,plot_offset,sprintf('IMF%d',i),...'FontSize',10,'HorizontalAlignment','left');ifi<num_imfs plot_offset=plot_offset+max(abs(imfs(i,:)))*1.5;endendtitle('EMD分解结果','FontSize',12,'FontWeight','bold');xlabel('时间 (s)');ylabel('IMF分量');grid on;xlim([t(1),t(end)]);%% 子图3: Hilbert谱 (时频谱)subplot(4,3,[4,5]);[T,F]=meshgrid(t,1:num_imfs);surf(T,F,hht_spectrum,'EdgeColor','none');view(2);colormap(jet);colorbar;title('Hilbert谱 (时频分布)','FontSize',12,'FontWeight','bold');xlabel('时间 (s)');ylabel('IMF序号');%% 子图4: 瞬时频率subplot(4,3,6);hold on;fori=1:num_imfsplot(t,instantaneous_freq(i,:),'Color',colors(i,:),'LineWidth',1);endtitle('瞬时频率','FontSize',12,'FontWeight','bold');xlabel('时间 (s)');ylabel('频率 (Hz)');grid on;legend(arrayfun(@(x)sprintf('IMF%d',x),1:num_imfs,'UniformOutput',false),...'Location','best','FontSize',8);%% 子图5: 各IMF能量谱subplot(4,3,[7,8]);freq_bins=1:size(energy_density,2);bar_data=energy_density(1:end-1,:)';h_bar=bar(freq_bins,bar_data,'stacked');title('各IMF分量能量谱分布','FontSize',12,'FontWeight','bold');xlabel('频率区间');ylabel('归一化能量');legend(arrayfun(@(x)sprintf('IMF%d',x),1:num_imfs,'UniformOutput',false),...'Location','best','FontSize',8);grid on;%% 子图6: 总能量谱subplot(4,3,[9,10]);total_energy=energy_density(end,:);stem(freq_bins,total_energy,'filled','LineWidth',1.5,'MarkerSize',4);title('总能量谱密度','FontSize',12,'FontWeight','bold');xlabel('频率区间');ylabel('能量密度');grid on;%% 子图7: 能量分布统计subplot(4,3,[11,12]);% 计算各IMF能量贡献imf_energy_contribution=zeros(1,num_imfs);fori=1:num_imfsimf_energy_contribution(i)=sum(hht_spectrum(i,:));endtotal_energy=sum(imf_energy_contribution);imf_energy_contribution=imf_energy_contribution/total_energy*100;pie(imf_energy_contribution,arrayfun(@(x)sprintf('IMF%d\n%.1f%%',x,...imf_energy_contribution(x)),1:num_imfs,'UniformOutput',false));title('各IMF能量贡献比例','FontSize',12,'FontWeight','bold');sgtitle('EMD分解与希尔伯特能量谱分析','FontSize',14,'FontWeight','bold');end

6. 测试用例与演示

%% 测试函数：模拟信号分析functiontest_emd_hilbert_analysis()% 测试EMD-Hilbert能量谱分析fprintf('=== EMD-Hilbert能量谱分析测试 ===\n\n');% 生成测试信号fs=1000;% 采样频率 1000 Hzt=0:1/fs:2;% 2秒信号% 多分量测试信号signal=2*sin(2*pi*50*t).*(1+0.5*sin(2*pi*5*t))+...% 调幅信号1.5*sin(2*pi*120*t)+...% 高频分量0.8*chirp(t,80,2,20)+...% 扫频信号0.3*randn(size(t));% 噪声fprintf('生成测试信号:\n');fprintf(' - 50Hz调幅信号\n');fprintf(' - 120Hz正弦信号\n');fprintf(' - 80-20Hz扫频信号\n');fprintf(' - 高斯白噪声\n\n');% 执行EMD-Hilbert分析[hht_spectrum,energy_density,imfs]=emd_hilbert_energy_analysis(...signal,fs,...'num_imfs',8,...'max_iterations',500,...'plot_results',true,...'energy_band',[0,200]);% 显示能量分析结果display_energy_analysis(imfs,hht_spectrum);endfunctiondisplay_energy_analysis(imfs,hht_spectrum)% 显示能量分析结果fprintf('\n=== 能量分析结果 ===\n\n');num_imfs=size(imfs,1);total_energy=sum(hht_spectrum(:));fprintf('%-8s %-12s %-12s %-15s\n',...'IMF','能量','能量百分比','累计百分比');fprintf('%s\n',repmat('-',1,50));cumulative_percentage=0;fori=1:num_imfs imf_energy=sum(hht_spectrum(i,:));percentage=imf_energy/total_energy*100;cumulative_percentage=cumulative_percentage+percentage;fprintf('IMF%-5d %-12.4e %-11.2f%% %-14.2f%%\n',...i,imf_energy,percentage,cumulative_percentage);endfprintf('%s\n',repmat('-',1,50));fprintf('总计 %-12.4e %-11.2f%%\n',total_energy,100.0);% 显示主要能量集中区域fprintf('\n主要能量集中分析:\n');energy_threshold=0.05;% 5%能量阈值significant_imfs=[];fori=1:num_imfs imf_energy=sum(hht_spectrum(i,:))/total_energy;ifimf_energy>energy_threshold significant_imfs=[significant_imfs,i];endendif~isempty(significant_imfs)fprintf('主要能量集中在: IMF%s\n',mat2str(significant_imfs));fprintf('这些IMF包含了 %.2f%% 的总能量\n',...sum(hht_spectrum(significant_imfs,:),'all')/total_energy*100);elsefprintf('能量分布较为均匀，无明显主导IMF\n');endend%% 实际应用示例：故障诊断信号分析functionbearing_fault_analysis()% 轴承故障信号分析示例fprintf('=== 轴承故障诊断信号分析 ===\n\n');% 模拟轴承故障信号 (简化模型)fs=12000;% 12 kHz采样频率t=0:1/fs:1;% 1秒信号% 轴承故障特征频率 (示例)bpfo=120;% 外圈故障频率 120 Hzbpfi=145;% 内圈故障频率 145 Hzbsf=65;% 滚动体故障频率 65 Hzftf=15;% 保持架故障频率 15 Hz% 生成故障信号fundamental=0.8*sin(2*pi*30*t);% 转频 30 Hz% 故障冲击序列fault_impacts=zeros(size(t));impact_interval=round(fs/bpfo);% 外圈故障冲击间隔impact_locations=1:impact_interval:length(t);fault_impacts(impact_locations)=1.5;% 共振响应resonance_freq=3000;% 3 kHz 共振频率resonance_signal=filter([1,-1.8*cos(2*pi*resonance_freq/fs),0.81],...[1,-1.2*cos(2*pi*resonance_freq/fs),0.25],fault_impacts);% 合成信号signal=fundamental+resonance_signal+0.1*randn(size(t));fprintf('轴承故障信号特征:\n');fprintf(' 转频: 30 Hz\n');fprintf(' 外圈故障频率: 120 Hz\n');fprintf(' 系统共振频率: 3000 Hz\n\n');% 执行EMD-Hilbert分析[hht_spectrum,energy_density,imfs]=emd_hilbert_energy_analysis(...signal,fs,...'num_imfs',6,...'energy_band',[0,5000],...'plot_results',true);% 故障特征提取extract_fault_features(imfs,hht_spectrum,fs,[bpfo,bpfi,bsf,ftf]);endfunctionextract_fault_features(imfs,hht_spectrum,fs,fault_freqs)% 提取故障特征fprintf('\n=== 故障特征分析 ===\n\n');% 计算各IMF的频谱num_imfs=size(imfs,1);fori=1:num_imfs% FFT分析imf_fft=fft(imfs(i,:));N=length(imf_fft);f=(0:N-1)*fs/N;% 寻找主要频率成分[~,dominant_freq_idx]=max(abs(imf_fft(1:floor(N/2))));dominant_freq=f(dominant_freq_idx);fprintf('IMF%d 主要频率成分: %.2f Hz\n',i,dominant_freq);% 检查是否接近已知故障频率forj=1:length(fault_freqs)ifabs(dominant_freq-fault_freqs(j))<10% 10 Hz容差switchjcase1fprintf(' → 检测到外圈故障特征!\n');case2fprintf(' → 检测到内圈故障特征!\n');case3fprintf(' → 检测到滚动体故障特征!\n');case4fprintf(' → 检测到保持架故障特征!\n');endendendendend

参考代码 emd分解之后再进行希尔伯特变换，获得能量谱www.3dddown.com/csa/77762.html

算法特点与优势

EMD-Hilbert方法的优势

1. 自适应分解- EMD根据信号特性自动分解
2. 局部特征提取- 能够捕捉非平稳信号的时变特性
3. 能量定位- 精确确定能量在时频域的分布
4. 物理意义明确- 每个IMF代表特定的振动模式

应用场景

• 机械故障诊断- 轴承、齿轮箱故障检测
• 生物信号处理- EEG、ECG信号分析
• 地震信号分析- 地层结构识别
• 语音信号处理- 语音特征提取