✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab仿真内容点击👇

智能优化算法       神经网络预测       雷达通信       无线传感器        电力系统

信号处理              图像处理               路径规划       元胞自动机        无人机 

基于Matlab实现稳健的经验模式分解 (REMD)

clc; clear; close all;

fs = 10000; % sampling frequency

N = 30000; % data amount

t = (1:N)/fs; % time vector

x = (2+cos(2*pi*0.5*t)).*cos(2*pi*5*t+15*t.^2)+...

     cos(2*pi*2*t);

[imf,ort,fvs,iterNum] = emd_sssc(x,fs,'display',1);

function [imf,ort,fvs,iterNum] = emd_sssc(x,fs,varargin)

%

% ERPHM Code Log

% 2017-06-20 Created by Zhiliang Liu, Dandan Peng and Yaqiang Jin.

% 2019-01-19 Modified by Zhiliang Liu and Dandan Peng.

%

% This funcion performs the empirical mode decomposition(EMD) with a soft sifting stopping criterion (SSSC) on the input signal, and

% return intrinsic mode function (IMF). If you have any questions, please contact us via [email protected]

%

% SYNTAX:

% imf = emd_sssc(x);

% [imf,ort,fvs,iterNum] = emd_sssc(x,fs);

% [imf,ort,fvs,iterNum] = emd_sssc(x,fs,options);

% [imf,ort,fvs,iterNum] = emd_sssc(x,fs,'option_name1',option_value1,....);

%

% INPUTS:

% [1] x: a row signal vecter.

% [2] option: a struct contains options' name(option_name) and

%         corresponding values(option_value).

% [2.1] 'display'

%       plot imf or not,1 plot,0 do not.

%       default: 0

% [2.2] 'max_iter'

%       max number of iterations in one imf sift procedure.

%       default: 30

% [2.3] 'max_imfs'

%    max number of imf obtained in emd procedure.

%    default: 10

% [2.4] 'ext_ratio'

%       end extension length to original data length ratio.

%       can be any positive real number from (0,1]

%    default: 0.2

%

% OUTPUTS

% [1] imf

%       Intrinsic Mode Function (imf) Matrix of which each row is an imf and the last

%       row is the residual signal.

% [2] iterNum

%       iteration number for each imf.

% [3] fvs

%       funciton values.

% [4] ort

%       index of orthogonality.

%

% REFERENCE

% [1] Dandan Peng, Zhiliang Liu, Yaqiang Jin, Yong Qin. Improved EMD with a Soft

%     Sifting Stopping Criterion and Its Application to Fault Diagnosis of Rotating Machinery.

%     Journal of Mechanical Engineering. Accepted on Jan. 01, 2019.

% [2] Zhiliang Liu, Yaqiang Jin, Ming J. Zuo, and Zhipeng Feng. Time-frequency

%     representation based on robust local mean decomposition for multi-component

%     AM-FM signal analysis. Mechanical Systems and Signal Processing. 95: 468-487, 2017.

%

% EXAMPLE:

[x, display, ssc, max_iter, max_imfs, smooth_mode,ext_ratio, x_energy, imf, iterNum, fvs]= initial(x,varargin{:});

% Initialize main loop

i = 0;

xs = x; % copy x to xs for sifting process, reserve original input as x.

nx = length(x);

while   i < max_imfs && ~stopemd(xs, x_energy) % outer loop for imf selection

    i = i+1;

    % initialize variables used in imf sifting loop

    s_j = zeros(max_iter,nx);

    % imf sifting iteration loop

    j = 0;

    stop_flag_sifting_process = 0;

    s = xs;

    while  j < max_iter && ~stop_flag_sifting_process %  inner loop for sifting process

        j = j+1;

        if j==1

            [m_j, n_extr] = emd_mean(s, smooth_mode, ext_ratio);

        end

        % force to stop iteration if number of extrema of s is smaller than 3.

        if n_extr < 3

            break;

        end

        s = s-m_j; % remove mean.

        s_j(j, :) = s;

        [m_j, n_extr] = emd_mean(s, smooth_mode,ext_ratio);

        [stop_flag_sifting_process,fvs(i,:)] = is_sifting_process_stop(m_j, s,j, fvs(i,:), ssc);

    end

    switch ssc

        case {'liu'}

            [~, opt0] = min(fvs(i,1:j)); % ***Critical Step***

            opt_IterNum = min(j, opt0); % in case iteration stop for n_extr 2 % input options seperately.

    try

        in_optns = struct(varargin{:});

    catch

        error('wrong argmument syntax')

    end

else

    error('arguments error: maybe not enough or wrong syntax')

end

names = fieldnames(in_optns);% get input options' name and value

for k = names'

    if ~any(strcmpi(char(k), optn_fields))

        % find any wrong argument in syntax.

        error(['bad option field name: ',char(k)])

    end

    if ~isempty(eval(['in_optns.',char(k)]))

        % alter default option values with input, and empty input keep default.

        eval(['optns.',lower(char(k)),' = in_optns.',char(k)])

    end

end

display = optns.display;

ssc = optns.ssc;

max_iter = optns.max_iter;

max_imfs = optns.max_imfs;

smooth_mode = optns.smooth_mode;

ext_ratio = optns.ext_ratio;

% initialize x(input signal), x_energy, IMF, iterNum and fvs 

x = x(:)'; % make x a row vector.

nx = length(x);

x_energy = sum(x.^2); % energy = square summation.

imf = zeros(max_imfs,nx);

iterNum = zeros(1,max_imfs);

fvs = zeros(max_imfs,max_iter);

end

% Check whether there are enough (3) extrema to continue the decomposition

function stop = stopemd(xs, x_energy)

[indmin,indmax] = extr(xs);

peak = length(indmin) + length(indmax);

ratio = sum(xs.^2)/x_energy;

stop = peak < 3 | ratio < 0.001;

end

% Compute mean function of x in EMD

function [m, n_extr] = emd_mean(x, smooth_mode, ext_ratio)

% find extremum indices

[indmin, indmax, ~] = extr(x);

% total amount of extrema

n_extr = length(indmin)+length(indmax);

if n_extr < 3

    m = [];

    return

end

% extend original data to refrain end effect

[ext_indmin,ext_indmax,ext_x,cut_index] = extend(x, indmin, indmax, ext_ratio);

% compute local mean 

switch smooth_mode

    case 'spline'

        l = length(ext_x);

        ext_indmax(ext_indmax==1) = [];

        ext_indmax(ext_indmax==l) = [];

        ext_indmin(ext_indmin==1) = [];

        ext_indmin(ext_indmin==l) = [];

        max_ext_x = interp1([1,ext_indmax,l], ext_x([1,ext_indmax,l]),...

            1:l, 'spline'); % pchip

        min_ext_x = interp1([1,ext_indmin,l], ext_x([1,ext_indmin,1]),...

            1:l, 'spline');

        ext_m = 0.5*(max_ext_x+min_ext_x);

        m = ext_m(cut_index(1):cut_index(2));

    otherwise

        error(['mean computation method must be one of the following: ma,',...

            ' spline',' maspline'])

end

end

% sifting stopping criterion

function [stop_flag_sifting_process,fv_i] = is_sifting_process_stop(m, s, j, fv_i, ssc)

df = m;

[indmin, indmax, indzer] = extr(s);

lm = length(indmin);

lM = length(indmax);

nem = lm + lM;

nzm = length(indzer);

switch ssc

    case 'liu' % local optimal iteration.

        fv_i(j) =rms(df)+ abs(kurtosis(df)-3);

        if j >= 3 && abs(nzm-nem)= fv_i(j-1)) && (fv_i(j-1) >= fv_i(j-2)))

                stop_flag_sifting_process = 1;

                return;

            end

        end

end

stop_flag_sifting_process = 0;

end

% Plot IMFs

function emdplot(x,imf,fs)

t = (1:size(imf,2))/fs;

imfn = size(imf,1);

figure

subplot(imfn+1,1,1);

plot(t,x);

title('Input Signal');

for imfi = 1:imfn

    subplot(imfn+1,1,imfi+1);

    plot(t,imf(imfi,:));

    if imfi < imfn

        title(['IMF',num2str(imfi)] );

    else

        title('Residual');

    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,imf)

% ort = IO(x,pfs) computes the index of orthogonality

%

% inputs : - x   : analyzed signal

%          - pfs  : production function

n = size(imf,1);

s = 0;

for i = 1:n

    for j =1:n

        if i~=j

            s = s + abs(sum(imf(i,:).*conj(imf(j,:)))/sum(x.^2));

        end

    end

end

ort = 0.5*s;

end

% Extracts the indices of extrema

% ** Copied from emd toolbox by G.Rilling and P.Flandrin

% ** http://perso.ens-lyon.fr/patrick.flandrin/emd.html

function [indmin, indmax, indzer] = extr(x)

m = length(x);

if nargout > 2

    x1 = x(1:m-1);

    x2 = x(2:m);

    indzer = find(x1.*x2 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

% 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, ext_ratio)

if ext_ratio == 0 % do not extend x

    ext_indmin = indmin;

    ext_indmax = indmax;

    ext_x = x;

    cut_index = [1,length(x)];

    return

end

nbsym = ceil(ext_ratio*length(indmax)); % number of extrema in extending end

xlen = length(x);

t = 1:xlen;

% boundary conditions for interpolations :

% left end extend

if indmax(1) < indmin(1) % first extremum is max

    if x(1) > x(indmin(1)) % first point > first min extremum

        lmax = fliplr(indmax(2:min(end,nbsym+1)));

        lmin = fliplr(indmin(1:min(end,nbsym)));

        lsym = indmax(1);

    else                  % first point < first min extremum

        lmax = fliplr(indmax(1:min(end,nbsym)));

        lmin = [fliplr(indmin(1:min(end,nbsym-1))),1];

        lsym = 1;

    end

else                     % first extremum is maxmum

    if x(1) < x(indmax(1)) % first point < first maxmum

        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

if indmax(end) < indmin(end) % last extremum is minimum

    if x(end) < x(indmax(end)) % last point < last maxmum

        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 maxmum

        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 maxmum

    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

tlmin = 2*t(lsym)-t(lmin);

tlmax = 2*t(lsym)-t(lmax);

trmin = 2*t(rsym)-t(rmin);

trmax = 2*t(rsym)-t(rmax);

% in case symmetrized parts do not extend enough

if tlmin(1) > t(1) || tlmax(1) > t(1)

    if lsym == indmax(1)

        lmax = fliplr(indmax(1:min(end,nbsym)));

    else

        lmin = fliplr(indmin(1:min(end,nbsym)));

    end

    if lsym == 1

        error('bug')

    end

    lsym = 1;

    %     tlmin = 2*t(lsym)-t(lmin);

    %     tlmax = 2*t(lsym)-t(lmax);

end

if trmin(end) < t(xlen) || trmax(end) < t(xlen)

    if rsym == indmax(end)

        rmax = fliplr(indmax(max(end-nbsym+1,1):end));

    else

        rmin = fliplr(indmin(max(end-nbsym+1,1):end));

    end

    if rsym == xlen

        error('bug')

    end

    rsym = xlen;

    %     trmin = 2*t(rsym)-t(rmin);

    %     trmax = 2*t(rsym)-t(rmax);

end

l_end = max(max(lmax, lmin));

r_end = min(min(rmax, rmin));

new_lmax = l_end+1-lmax;

new_lmin = l_end+1-lmin;

new_rmax = rsym-rmax;

new_rmin = rsym-rmin;

lx_length = l_end-lsym;

lx = fliplr(x(lsym+1:l_end));

rx = fliplr(x(r_end:rsym-1));

ext_x = [lx, x(lsym:rsym), rx];

ext_indmin = [new_lmin,indmin+lx_length-lsym+1,new_rmin+lx_length-lsym+1+...

    rsym];

ext_indmax = [new_lmax,indmax+lx_length-lsym+1,new_rmax+lx_length-lsym+1+...

    rsym];

% Index for cutting extension of x

cut_index = [lx_length-lsym+2, length(x)+lx_length-lsym+1];

end

[1] Zhiliang Liu*, Dandan Peng, Ming J. Zuo, Jiansuo Xia, and Yong Qin. Improved Hilbert-Huang transform with soft sifting stopping criterion and its application to fault diagnosis of wheelset bearings. ISA Transactions. 125: 426-444, 2022​

❤️部分理论引用网络文献，若有侵权联系博主删除

❤️ 关注我领取海量matlab电子书和数学建模资料