时序分解 | Matlab实现RLMD鲁棒性局部均值分解-编程知识

时序分解 | Matlab实现RLMD鲁棒性局部均值分解

效果一览

在这里插入图片描述

基本介绍

Matlab实现RLMD鲁棒性局部均值分解,可直接替换 Matlab语言
1.算法新颖小众，用的人很少，包含分解图
2.直接替换数据即可用适合新手小白注释清晰～
3.附赠excel测试数据直接运行main一键出图～

鲁棒性局部均值分解（Robust Locally Mean Decomposition，简称RLMD）是一种信号处理方法，用于对复杂信号进行分解和分析。它是局部均值分解（Locally Mean Decomposition，简称LMD）的改进版本，旨在提高对噪声和干扰的鲁棒性。

局部均值分解是一种将信号分解为局部平稳信号和调制信号的技术。它通过对信号进行迭代滤波和振幅调整来实现分解。然而，传统的局部均值分解对于存在噪声和干扰的信号可能会产生较差的分解结果。

RLMD的主要目标是提高对噪声和干扰的抵抗力。它在局部均值分解的基础上引入了鲁棒性因子，通过自适应调整滤波器的参数来减小噪声和干扰对分解结果的影响。这种方法可以提高信号分解的准确性和稳定性，使得分解结果更加可靠。

程序设计

完整源码和数据获取方式资源处下载Matlab实现RLMD鲁棒性局部均值分解。

%%  清空环境变量
warning off             % 关闭报警信息
close all               % 关闭开启的图窗
clear                   % 清空变量
clc                     % 清空命令行
function [pfs, ams, fms, ort, fvs, iterNum] = lmd_public(x,varargin)
switch ma_iter_mode  case 'fixed' % stick stepcntr = 0; % count times of moving averagemax_c = ceil(smax/span)*15; % theoretic%         max_c = ceil(smax/(span-1));%         nm = length(x);k = (span+1)/2;kmax = nm - (span-1)/2;while (k < kmax) && (cntr < max_c) % find flat stepif x(k) == x(k+1);x = smooth(x, span);cntr = cntr+1;k = k-1;endk = k+1;end%         while ~isempty(find(diff(x)==0)) && (cntr < max_c)% %             find(diff(x)==0)%             x = smooth(x, span);%             cntr = cntr+1;%         endotherwiseerror('No specifications for ma_iter_mode.');end
end% Extend original data to refrain end effect
% ** Modified on emd by G.Rilling and P.Flandrin
% ** http://perso.ens-lyon.fr/patrick.flandrin/emd.html
function [ext_indmin, ext_indmax, ext_x, cut_index] = extend(x, indmin,...indmax, extd_r)
if extd_r == 0 % do not extend xext_indmin = indmin;ext_indmax = indmax;ext_x = x;cut_index = [1,length(x)];return
end
nbsym = ceil(extd_r*length(indmax)); % number of extrema in extending endif x(1) < x(indmax(1)) % first point < first maximumlmax = fliplr(indmax(1:min(end,nbsym)));lmin = fliplr(indmin(2:min(end,nbsym+1)));lsym = indmin(1);else                   % first point > first minimumlmax = [fliplr(indmax(1:min(end,nbsym-1))),1];lmin = fliplr(indmin(1:min(end,nbsym)));lsym = 1;end
end% right end extension
if indmax(end) < indmin(end) % last extremum is minimumif x(end) < x(indmax(end)) % last point < last maximumrmax = fliplr(indmax(max(end-nbsym+1,1):end));rmin = fliplr(indmin(max(end-nbsym,1):end-1));rsym = indmin(end);else                       % last point > last maximumrmax = [xlen, fliplr(indmax(max(end-nbsym+2,1):end))];rmin = fliplr(indmin(max(end-nbsym+1,1):end));rsym = xlen;end
else                         % last extremum is maximumif x(end) > x(indmin(end)) % last point > last minimumrmax = fliplr(indmax(max(end-nbsym,1):end-1));rmin = fliplr(indmin(max(end-nbsym+1,1):end));rsym = indmax(end);else                       % last point < last minimumrmax = fliplr(indmax(max(end-nbsym+1,1):end));rmin = [xlen, fliplr(indmin(max(end-nbsym+2,1):end))];rsym = xlen;end
end% when two or more successive points have the same value we consider only
% one extremum in the middle of the constant area (only works if the signal
% is uniformly sampled)if any(d==0)imax = [];imin = [];bad = (d==0);dd = diff([0 bad 0]);debs = find(dd == 1);fins = find(dd == -1);if debs(1) == 1if length(debs) > 1debs = debs(2:end);fins = fins(2:end);elsedebs = [];fins = [];endendif ~isempty(debs)if fins(end) == mif length(debs) > 1debs = debs(1:(end-1));fins = fins(1:(end-1));elsedebs = [];fins = [];endendendlc = length(debs);if lc > 0for k = 1:lcif d(debs(k)-1) > 0if d(fins(k)) < 0%           imax = [imax round((fins(k)+debs(k))/2)];endelseif d(fins(k)) > 0%           imin = [imin round((fins(k)+debs(k))/2)];endendendendif ~isempty(imax)indmax = sort([indmax imax]);endif ~isempty(imin)indmin = sort([indmin imin]);endend
end% Compute the index of orthogonality
% ** Copied from emd toolbox by G.Rilling and P.Flandrin
% ** http://perso.ens-lyon.fr/patrick.flandrin/emd.html
function ort = io(x,pfs)
% ort = IO(x,pfs) computes the index of orthogonality
%
% inputs : - x   : analyzed signal
%          - pfs  : production functionn = size(pfs,1);s = 0;for i = 1:nfor j =1:nif i~=js = s + abs(sum(pfs(i,:).*conj(pfs(j,:)))/sum(x.^2));endend
endort = 0.5*s;
end% Plot PF, Amplititude Signal and FM Signal
function lmdplot(pfs, ams, fms, smooth_mode)
t = 1:size(pfs,2);
pfn = size(pfs,1);
figure
for pfi = 1:pfnsubplot(pfn,1,pfi);plot(t,pfs(pfi,:));if pfi < pfntitle(['PF',num2str(pfi),' (',smooth_mode,')']);elsetitle(['Residual',' (',smooth_mode,')']);end
end
figure
for ai = 1:pfn-1subplot(pfn-1,1,ai);plot(t,ams(ai,:));title(['Amplititude Signal',num2str(ai),' (',smooth_mode,')']);
end
figure
for fsi = 1:pfn-1subplot(pfn-1,1,fsi);plot(t,fms(fsi,:));title(['FM Signal',num2str(fsi),' (',smooth_mode,')']);
endend