eemd分解matlab代码

以下是使用 Matlab 实现 EEMD 分解的代码：

function [IMF,residual] = emd(x)
% EMD - Perform Empirical Mode Decomposition
%
% [IMF,residual] = emd(x)
%
% Inputs:
%   x         - Input signal (must be a column vector)
%
% Outputs:
%   IMF       - Matrix of intrinsic mode functions (one IMF per row)
%   residual  - Residual signal (last IMF + residual = original signal)
%
% Example:
%
% t=linspace(0,1,2^14);
% x=sin(2*pi*50*t)+sin(2*pi*120*t);
% x=x+randn(size(t));
% [imf,residual]=emd(x);
% subplot(length(imf)+1,1,1);
% plot(t,x);
% for k=1:length(imf)
%     subplot(length(imf)+1,1,k+1);
%     plot(t,imf(k,:));
% end
% 
% Algorithm based on:
% Huang et al, "The empirical mode decomposition and the Hilbert spectrum
% for nonlinear and non-stationary time series analysis", Proc. Royal Soc.
% London A, Vol. 454, pp. 903-995, 1998.

%-------------------------------------------------------------------------
% Preprocess input
%-------------------------------------------------------------------------

% Force x to be a column vector
x = x(:);

%-------------------------------------------------------------------------
% Set parameters and initialize variables
%-------------------------------------------------------------------------

% Maximum number of iterations
nMax = 500;

% Sifting tolerance (stop criterion)
tol = 1e-5;

% Number of sifting iterations
nIMF = 0;

% Extract first IMF
x1 = x;
h = x;
while true
    
    nIMF = nIMF + 1;
    
    % Extract local maxima and minima
    [maxtab,mintab] = peakdet(x1,0.05);
    
    % Interpolate to get envelopes
    if isempty(maxtab) || isempty(mintab)
        break
    end
    max_env = interp1(maxtab(:,1),maxtab(:,2),1:length(x1));
    min_env = interp1(mintab(:,1),mintab(:,2),1:length(x1));
    
    % Calculate mean envelope
    mean_env = (max_env + min_env)/2;
    
    % Extract IMF
    IMF(nIMF,:) = x1 - mean_env;
    
    % Calculate residual
    x1 = x1 - IMF(nIMF,:);
    
    % Check for convergence
    if sum(abs(IMF(nIMF,:))) < tol || nIMF >= nMax
        break
    end
    
end

% Calculate residual
residual = x1;

end

function extrema = peakdet(v, delta)
%PEAKDET Detect peaks in a vector
%        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
%        maxima and minima ("peaks") in the vector V.
%        MAXTAB and MINTAB consists of two columns. Column 1
%        contains indices in V, column 2 the found values.
%
%        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
%        in MAXTAB and MINTAB are replaced with the corresponding
%        X-values.
%
%        A point is considered a maximum peak if it has the maximal
%        value, and was preceded (to the left) by a value lower by
%        DELTA. DELTA may be a vector specifying the allowed difference
%        between two peaks. For example if DELTA=[3 5], then a maximum
%        is the first point in the neighborhood with a value higher than
%        the previous one by 3 AND higher than the next one by 5.
%
%        Vice versa, a point is considered a minimum peak if it has the
%        minimal value, and was preceded by a higher value by DELTA.
%        (c) 2002,2004 Copyright (C) by Thomas C. O'Haver
%        http://www.mathworks.com/matlabcentral/fileexchange/12275-peakdet
%        Version 4, 1 June 2016


%-------------------------------------------------------------------------
% Parse and check input arguments
%-------------------------------------------------------------------------

% Check number of input arguments
narginchk(2,3);

% Check dimensions of input arguments
assert(isvector(v),'Input argument "v" must be a vector');
assert(isnumeric(delta) && isvector(delta),'Input argument "delta" must be a vector');

% Check values of input arguments
assert(all(delta > 0),'Values of input argument "delta" must be positive');

% Check for complex input
if ~isreal(v)
    warning('Input argument "v" is complex; imaginary part ignored');
    v = real(v);
end

% Force input to be a row vector
v = v(:)';

%-------------------------------------------------------------------------
% Set default values for optional input arguments
%-------------------------------------------------------------------------

x = 1:length(v);

%-------------------------------------------------------------------------
% Find peaks
%-------------------------------------------------------------------------

% Preallocate output
maxtab = [];
mintab = [];

% Loop over each specified delta
for d = delta
    
    % Find all maxima and their indices
    if d > 0
        % Find local maxima
        maxloc = (v(2:end-1) > v(1:end-2)) & (v(2:end-1) > v(3:end));
        % Add first and last point
        maxloc = [0, maxloc, 0];
        % Find indices
        maxind = find(maxloc);
        % Remove maxima below delta
        maxind(v(maxind) - v(maxind-1) < d) = [];
        maxind(v(maxind) - v(maxind+1) < d) = [];
        % Store peak values and indices
        maxtab = [maxtab; x(maxind)', v(maxind)'];
    end
    
    % Find all minima and their indices
    if d > 0
        % Find local minima
        minloc = (v(2:end-1) < v(1:end-2)) & (v(2:end-1) < v(3:end));
        % Add first and last point
        minloc = [0, minloc, 0];
        % Find indices
        minind = find(minloc);
        % Remove minima below delta
        minind(v(minind) - v(minind-1) > -d) = [];
        minind(v(minind) - v(minind+1) > -d) = [];
        % Store peak values and indices
        mintab = [mintab; x(minind)', v(minind)'];
    end
    
end

% Sort output
maxtab = sortrows(maxtab,1);
mintab = sortrows(mintab,1);

%-------------------------------------------------------------------------
% Replace indices with corresponding X-values
%-------------------------------------------------------------------------

if nargin == 3
    maxtab(:,1) = x(maxtab(:,1));
    mintab(:,1) = x(mintab(:,1));
end

%-------------------------------------------------------------------------
% Return output
%-------------------------------------------------------------------------

extrema = {maxtab, mintab};

end

这个代码实现了 EMD 的分解，并返回每个分量的 IMF 和剩余部分（residual）。在主函数 `emd` 中，使用了一个 while 循环来不断提取每个 IMF，直到满足收敛条件或达到最大迭代次数。在循环中，使用 `peakdet` 函数来寻找局部极大值和极小值，进而计算出 IMF。最终返回的 `IMF` 是一个矩阵，每一行代表一个 IMF。