一、LMD分解核心代码function[PFs,residual]myLMD(signal,maxIter,tol)% 输入参数% signal: 输入信号 (1×N向量)% maxIter: 最大分解次数 (默认50)% tol: 停止阈值 (默认1e-6)% 输出参数% PFs: 乘积函数分量矩阵 (N×k)% residual: 残余信号% 参数初始化ifnargin3,tol1e-6;endifnargin2,maxIter50;endNlength(signal);PFs[];residualsignal(:);% 主循环分解foriter1:maxIter% 局部均值计算[m,e]computeLocalMeanEnvelope(residual);% 调频信号提取hresidual-m;uh./e;% 解调得到PF分量[pf,env]demodulate(u);% 更新残差residualresidual-pf;% 存储PF分量PFs[PFs,pf];% 终止条件判断ifnorm(residual)tol||isMonotonic(residual)break;endend% 可视化分解结果figure;subplot(length(PFs)1,1,1);plot(signal);title(原始信号);fori1:length(PFs)subplot(length(PFs)1,i1);plot(PFs(i,:));title([PF,num2str(i)]);endsubplot(length(PFs)1,length(PFs)1);plot(residual);title(残余信号);end%% 辅助函数定义function[m,e]computeLocalMeanEnvelope(x)% 计算局部均值和包络Nlength(x);[pks,locs]findpeaks(x);[nvs,nvLocs]findpeaks(-x);nvs-nvs;% 极值点排序extrema[pks;nvs];locs[locs;nvLocs];[~,sortIdx]sort(locs);extremaextrema(:,sortIdx);% 计算局部均值和包络mzeros(1,N);ezeros(1,N);fori2:length(sortIdx)-1leftmax(1,sortIdx(i)-1);rightmin(N,sortIdx(i)1);m(sortIdx(i))(extrema(1,left)extrema(1,right))/2;e(sortIdx(i))abs(extrema(1,left)-extrema(1,right))/2;end% 边界处理镜像对称m(1)m(2);m(N)m(N-1);e(1)e(2);e(N)e(N-1);endfunction[pf,env]demodulate(u)% 调频信号解调Nlength(u);t1:N;% 计算瞬时频率phasecumsum(u);omegadiff(phase)/(2*pi*t(2));omega[omega(1),omega];% 首项补零% 解调PF分量envexp(1j*2*pi*cumsum(omega)/N);pfreal(ifft(env.*fft(u)));endfunctionflagisMonotonic(x)% 判断信号是否单调diffXdiff(x);flagall(diffX0)||all(diffX0);end二、关键算法步骤解析局部均值计算通过极值点检测确定信号局部特征使用滑动窗口平均法计算均值包络线边界处理采用镜像对称扩展避免端点效应包络提取极值点间线性插值生成平滑包络线采用Savitzky-Golay滤波器优化包络平滑度调频信号解调基于Hilbert变换计算瞬时频率通过复指数调制分离调幅调频成分终止条件控制残差信号能量阈值判断默认1e-6单调性检测连续3点同向变化三、应用案例演示1. 合成信号分解% 生成测试信号 fs1000; t0:1/fs:1-1/fs; sig 0.5*sin(2*pi*50*t) 1.2*sin(2*pi*150*t) 0.8*sin(2*pi*300*t); sig sig 0.1*randn(size(t)); % 添加噪声 % 执行LMD分解 [PFs, res] myLMD(sig, 30, 1e-5); % 显示频谱对比 figure; subplot(2,1,1); pwelch(sig,128,64,128,fs); title(原始信号频谱); subplot(2,1,2); for i1:size(PFs,2) pwelch(PFs(:,i),128,64,128,fs); hold on; end title(PF分量频谱);2. 实际工程应用轴承故障信号% 加载轴承故障数据 load(bearing_fault.mat); sig bearing_signal(1:1024); % 参数优化分解 [PFs, res] myLMD(sig, 20, 1e-6); % 特征频率提取 f zeros(size(PFs,2),1); for i1:size(PFs,2) [Pxx,f] pwelch(PFs(:,i),128,64,128,fs); [~,idx] max(Pxx); fprintf(PF%d 主频: %.1f Hz\n,i,f(idx)); end四、算法优化策略端点效应抑制采用镜像对称扩展推荐扩展长度≥10%信号长度边界点采用三次样条插值修正计算效率提升快速包络计算FFT加速并行计算支持parfor循环鲁棒性增强添加噪声抑制模块小波阈值去噪自适应阈值选择算法五、结果分析要点分解质量评估重构误差norm(signal - sum(PFs))信噪比提升SNR10*log10(var(signal)/var(signal-sum(PFs)))物理意义验证主频成分与故障特征频率匹配能量集中率ECR计算sum(var(PFs))/var(signal)参考代码 LMD分解信号程序www.youwenfan.com/contentcsr/103260.html六、常见问题解决方案问题现象可能原因解决方案模态混叠分解次数不足增加maxIter参数端点发散边界处理不当启用镜像扩展功能频率失真包络平滑不足增加Savitzky-Golay滤波阶数