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

时序分解 | 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 step
        cntr = 0; % count times of moving average
        max_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 step
            if x(k) == x(k+1);
                x = smooth(x, span);
                cntr = cntr+1;
                k = k-1;
            end
            k = k+1;
        end
        %         while ~isempty(find(diff(x)==0)) && (cntr < max_c)
        % %             find(diff(x)==0)
        %             x = smooth(x, span);
        %             cntr = cntr+1;
        %         end
    otherwise
        error('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 x
    ext_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 end

    if x(1) < x(indmax(1)) % first point < first maximum
        lmax = fliplr(indmax(1:min(end,nbsym)));
        lmin = fliplr(indmin(2:min(end,nbsym+1)));
        lsym = indmin(1);
    else                   % first point > first minimum
        lmax = [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 minimum
    if x(end) < x(indmax(end)) % last point < last maximum
        rmax = fliplr(indmax(max(end-nbsym+1,1):end));
        rmin = fliplr(indmin(max(end-nbsym,1):end-1));
        rsym = indmin(end);
    else                       % last point > last maximum
        rmax = [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 maximum
    if x(end) > x(indmin(end)) % last point > last minimum
        rmax = fliplr(indmax(max(end-nbsym,1):end-1));
        rmin = fliplr(indmin(max(end-nbsym+1,1):end));
        rsym = indmax(end);
    else                       % last point < last minimum
        rmax = 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) == 1
        if length(debs) > 1
            debs = debs(2:end);
            fins = fins(2:end);
        else
            debs = [];
            fins = [];
        end
    end
    if ~isempty(debs)
        if fins(end) == m
            if length(debs) > 1
                debs = debs(1:(end-1));
                fins = fins(1:(end-1));
                
            else
                debs = [];
                fins = [];
            end
        end
    end
    lc = length(debs);
    if lc > 0
        for k = 1:lc
            if d(debs(k)-1) > 0
                if d(fins(k)) < 0
                    %           imax = [imax round((fins(k)+debs(k))/2)];
                end
            else
                if d(fins(k)) > 0
                    %           imin = [imin round((fins(k)+debs(k))/2)];
                end
            end
        end
    end
    
    if ~isempty(imax)
        indmax = sort([indmax imax]);
    end
    
    if ~isempty(imin)
        indmin = sort([indmin imin]);
    end
    
end
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 function

n = size(pfs,1);

s = 0;

for i = 1:n
    for j =1:n
        if i~=j
            s = s + abs(sum(pfs(i,:).*conj(pfs(j,:)))/sum(x.^2));
        end
    end
end

ort = 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:pfn
    subplot(pfn,1,pfi);
    plot(t,pfs(pfi,:));
    if pfi < pfn
        title(['PF',num2str(pfi),' (',smooth_mode,')']);
    else
        title(['Residual',' (',smooth_mode,')']);
    end
end
figure
for ai = 1:pfn-1
    subplot(pfn-1,1,ai);
    plot(t,ams(ai,:));
    title(['Amplititude Signal',num2str(ai),' (',smooth_mode,')']);
end
figure
for fsi = 1:pfn-1
    subplot(pfn-1,1,fsi);
    plot(t,fms(fsi,:));
    title(['FM Signal',num2str(fsi),' (',smooth_mode,')']);
end

end