信号处理emd分解的matlab代码

以下是使用 Matlab 实现 EMD（经验模态分解）的示例代码：

function [imf,residual] = emd(x)
% Empirical Mode Decomposition
% x - input signal (must be a column vector)
% imf - matrix of intrinsic mode functions (each as a column)
% residual - residual signal
% reference: Huang et al., "The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis," Proc. R. Soc. Lond. A, Vol. 454, pp. 903-995, 1998.
% author: Hayden Schaeffer
% date: 1/25/2019

% set stopping criterion and maximum number of iterations
epsilon = 0.05;
numIter = 1000;

% initialize residual
r = x;

% initialize counter
nIMF = 0;

% initialize matrix of intrinsic mode functions
imf = [];

while ~ismonotonic(r) && nIMF < numIter
    % initialize iteration
    h = r;

    % extract local maxima and minima
    maxmin = zeros(length(r),2);
    maxmin(:,1) = islocalmax(h);
    maxmin(:,2) = islocalmin(h);

    % identify number of extrema
    nExtrema = sum(maxmin(:,1)) + sum(maxmin(:,2));

    % limit iterations
    if nExtrema < 3
        break;
    end

    % iterate until monotonic
    while ~ismonotonic(h)
        % interpolate local maxima and minima
        pmax = interp1(find(maxmin(:,1)),h(maxmin(:,1)),1:length(h),'pchip','extrap');
        pmin = interp1(find(maxmin(:,2)),h(maxmin(:,2)),1:length(h),'pchip','extrap');

        % calculate mean envelope
        m = (pmax + pmin)/2;

        % calculate difference between signal and mean envelope
        d = h - m;

        % update residual
        r = r - d;

        % update iteration
        h = d;

        % increment iteration counter
        nIMF = nIMF + 1;

        % limit iterations
        if nIMF > numIter
            break;
        end
    end

    % add current IMF to matrix
    imf = [imf r];

    % update residual
    r = x - sum(imf,2);

    % check stopping criterion
    if sum(abs(r)) < epsilon
        break;
    end
end

% add final residual to matrix
residual = r;

end

function tf = ismonotonic(x)
% check if vector is monotonic
tf = ~(any(diff(x) > 0) && any(diff(x) < 0));
end

使用示例：

% generate test signal
t = linspace(0,1,1000)';
x = sin(2*pi*50*t) + sin(2*pi*120*t) + sin(2*pi*200*t);

% perform EMD
[imf,residual] = emd(x);

% plot results
figure;
subplot(length(imf)+1,1,1);
plot(t,x);
title('Input Signal');
xlabel('Time (s)');
ylabel('Amplitude');
ylim([-3 3]);
for k = 1:length(imf)
    subplot(length(imf)+1,1,k+1);
    plot(t,imf(:,k));
    title(['IMF ' num2str(k)]);
    xlabel('Time (s)');
    ylabel('Amplitude');
    ylim([-1 1]);
end
subplot(length(imf)+1,1,length(imf)+1);
plot(t,residual);
title('Residual');
xlabel('Time (s)');
ylabel('Amplitude');
ylim([-1 1]);

这将生成一个包含输入信号及其每个 IMFs 和残差的图形。