MATLAB 的ICEEMDAN分解代码实现

0、前言

        本文讲解ICEEMDAN分解方法，并分享代码。

1、ICEEMDAN实现

         下面为主函数部分：

ecg=data;%data为待分解的一个信号数据，请替换为自己数据就行
%% 参数设置
Nstd = 0.2;
NR = 1;
MaxIter = 5000;
%% ICEEMDAN
[modes]=iceemdan(ecg,Nstd,NR,MaxIter,1);%iceemdan
modes=modes';
t=1:length(ecg);
[a b]=size(modes);
figure;
subplot(a+1,1,1);
plot(t,ecg);% the ECG signal is in the first row of the subplot
ylabel('original')
set(gca,'xtick',[])
title('ICEEMDAN')
axis tight;

for i=2:a
    subplot(a+1,1,i);
    plot(t,modes(i-1,:));
    ylabel (['IMF ' num2str(i-1)]);
    set(gca,'xtick',[])
    xlim([1 length(ecg)])
end

subplot(a+1,1,a+1)
plot(t,modes(a,:))
ylabel(['IMF ' num2str(a)])
xlim([1 length(ecg)])
xlabel('样本点')

 子函数iceemdan的代码：

function [modes]=iceemdan(x,Nstd,NR,MaxIter,SNRFlag)
% The current is an improved version, introduced in:

%[1] Colominas MA, Schlotthauer G, Torres ME. "Improve complete ensemble EMD: A suitable tool for biomedical signal processing" 
%       Biomedical Signal Processing and Control vol. 14 pp. 19-29 (2014)

%The CEEMDAN algorithm was first introduced at ICASSP 2011, Prague, Czech Republic

%The authors will be thankful if the users of this code reference the work
%where the algorithm was first presented:

%[2] Torres ME, Colominas MA, Schlotthauer G, Flandrin P. "A Complete Ensemble Empirical Mode Decomposition with Adaptive Noise"
%       Proc. 36th Int. Conf. on Acoustics, Speech and Signa Processing ICASSP 2011 (May 22-27, Prague, Czech Republic)

%Author: Marcelo A. Colominas
%contact: macolominas@bioingenieria.edu.ar
%Last version: 25 feb 2015
desvio_x=std(x);
x=x/desvio_x;
[a,b]=size(x);
temp=zeros(b,1);
 modes=zeros(b,1);
 aux=zeros(a,b);
for i=1:NR
    white_noise{i}=randn(size(x));%creates the noise realizations
end;

for i=1:NR
    modes_white_noise{i}=emd(white_noise{i},'display',0);%calculates the modes of white gaussian noise
end;
% save interval modes_white_noise
for i=1:NR %calculates the first mode
    xi=x+Nstd*modes_white_noise{i}(:,1)'/std(modes_white_noise{i}(:,1));
    [temp, o, it]=emd(xi,'MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
    aux=aux+(xi-temp')/NR;% nnnnnnnnnnnnnnnnJub局部包络
end;

modes= (x-aux)'; %saves the first mode
medias = aux; %  r1
k=1;
aux=zeros(a,b);
es_imf = min(size(emd(medias(1,:),'SiftMaxIterations',MaxIter,'display',0)));

while es_imf>1 %calculates the rest of the modes
    for i=1:NR
        tamanio=size(modes_white_noise{i});
        if tamanio(2)>=k+1
            noise=modes_white_noise{i}(:,k+1);
            if SNRFlag == 2
                noise=noise/std(noise); %adjust the std of the noise
            end;
            noise=Nstd*noise;
            try
                [temp,o,it]=emd(medias(1,:)+std(medias(1,:))*noise','MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
            catch    
                temp=emd(medias(1,:)+std(medias(1,:))*noise','MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
            end;
        else
            try
                [temp, o, it]=emd(medias(1,:),'MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
            catch
                temp=emd(medias(1,:),'MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
            end;
        end;
        aux=aux+(medias(1,:)+std(medias(1,:))*noise'-temp')/NR;% r2 r3 r...
    end;
    modes=[modes (medias(1,:)-aux)'];
    medias = aux;
    aux=zeros(size(x));
    k=k+1;
    es_imf = min(size(emd(medias(1,:),'SiftMaxIterations',MaxIter,'display',0)));
end;
modes = [modes (medias(1,:))'];
modes=modes*desvio_x;

        上述代码可以直接运行，实现信号分解。