MATLAB的EMD代码

以下是MATLAB中实现EMD的代码示例：

function [IMF, residual] = emd(x)
% EMD - Empirical Mode Decomposition
%
% Useage:
% [IMF,residual] = emd(x)
%
% Inputs:
% x - input signal (must be a column vector)
%
% Outputs:
% IMF - Matrix of intrinsic mode functions (IMFs). Each column is an IMF.
% residual - Residual signal after extracting all IMFs
%
% Description:
% This function performs the EMD on a given input signal to decompose it
% into its intrinsic mode functions (IMFs). The EMD method is based on the
% algorithm presented in Huang et al. (1998). The method iteratively
% extracts all IMFs by performing a sifting process, where all local maxima
% and minima of the signal are connected using cubic splines. The spline
% envelope is then subtracted from the signal and the process is repeated
% on the residue until a stopping criterion is met. The final residue is
% considered as the last IMF.
%
% Reference:
% Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H.,
% Zheng, Q., Yen, N. C., Tung, C. C., and Liu, H. H. (1998).
% The empirical mode decomposition and the Hilbert spectrum for nonlinear
% and non-stationary time series analysis. Proceedings of the Royal Society
% of London A: Mathematical, Physical and Engineering Sciences,
% 454(1971), 903-995.
%
% Author:
% Joel M. Leventhal, Ph.D.
% Leventhal Consulting
% jleventhal@ieee.org
% Date Created: 2016-10-18
% Date Modified: 2016-10-20

% Check input arguments
narginchk(1,1);
if ~isvector(x)
    error('Input signal must be a column vector');
end

% Initialize variables
N = length(x);
imfCount = 0;
MAXITER = 1000; % Maximum number of iterations allowed
tolerance = 1e-5; % Stopping criterion

%-------------------------------------------------------------------------
% Start by extracting first IMF
h = x;
while true
    % Find local maxima and minima
    maximaIdx = findLocalMaxima(h);
    minimaIdx = findLocalMinima(h);

    % Interpolate maxima and minima using cubic splines
    maximaEnv = interp1(maximaIdx,h(maximaIdx),(1:N)','pchip');
    minimaEnv = interp1(minimaIdx,h(minimaIdx),(1:N)','pchip');

    % Compute mean envelope
    meanEnv = (maximaEnv + minimaEnv) / 2;

    % Sift data and compute residue
    prevh = h;
    h = h - meanEnv;
    residue = prevh - h;

    % Check to see if residue looks like an IMF
    if isImf(residue,h,prevh,tolerance)
        imfCount = imfCount + 1;
        IMF(:,imfCount) = residue;
        break;
    end

    % Check to see if we've reached maximum number of iterations
    if imfCount >= MAXITER
        warning('EMD:MaxIter','Maximum number of iterations exceeded');
        break;
    end
end

%-------------------------------------------------------------------------
% Continue extracting IMFs
while true
    % Find local maxima and minima
    maximaIdx = findLocalMaxima(h);
    minimaIdx = findLocalMinima(h);

    % Interpolate maxima and minima using cubic splines
    maximaEnv = interp1(maximaIdx,h(maximaIdx),(1:N)','pchip');
    minimaEnv = interp1(minimaIdx,h(minimaIdx),(1:N)','pchip');

    % Compute mean envelope
    meanEnv = (maximaEnv + minimaEnv) / 2;

    % Sift data and compute residue
    prevh = h;
    h = h - meanEnv;
    residue = prevh - h;

    % Check to see if residue looks like an IMF
    if isImf(residue,h,prevh,tolerance)
        imfCount = imfCount + 1;
        IMF(:,imfCount) = residue;
        if std(residue) < 0.1*std(x)
            break;
        end
        h = h - residue;
    end

    % Check to see if we've reached maximum number of iterations
    if imfCount >= MAXITER
        warning('EMD:MaxIter','Maximum number of iterations exceeded');
        break;
    end
end

%-------------------------------------------------------------------------
% Compute residual signal
residual = x - sum(IMF,2);

%-------------------------------------------------------------------------
% Utility functions
%-------------------------------------------------------------------------
function maximaIdx = findLocalMaxima(h)
% Find local maxima in the signal
maximaIdx = [];
for k = 2:(length(h)-1)
    if h(k) > h(k-1) && h(k) > h(k+1)
        maximaIdx = [maximaIdx; k];
    end
end

function minimaIdx = findLocalMinima(h)
% Find local minima in the signal
minimaIdx = [];
for k = 2:(length(h)-1)
    if h(k) < h(k-1) && h(k) < h(k+1)
        minimaIdx = [minimaIdx; k];
    end
end

function isimf = isImf(residue,h,prevh,tolerance)
% Check to see if the residue is an IMF
h = h - prevh;
maximaIdx = findLocalMaxima(h);
minimaIdx = findLocalMinima(h);
if length(maximaIdx) + length(minimaIdx) <= 2
    isimf = false;
    return;
end
maximaEnv = interp1(maximaIdx,h(maximaIdx),(1:length(h))','pchip');
minimaEnv = interp1(minimaIdx,h(minimaIdx),(1:length(h))','pchip');
meanEnv = (maximaEnv + minimaEnv) / 2;
if max(abs(meanEnv)) < tolerance
    isimf = true;
else
    isimf = false;
end

该代码实现了EMD算法，并将信号分解为其内部模态函数(IMF)和残差信号。