function allmode=eemd(Y,Nstd,NE)
% This is an EMD/EEMD program
%
% 输入：
% Y，Inputted data：1-d data only
% Nstd，附加噪声的标准差与Y的标准差之比：Nstd = (0.1 ~ 0.4)*std(Y).
% NE，EEMD中添加噪声的总体次数：NE = 10-50.
%
% 输出：
% N*（m+1）的矩阵，其中N是输入数据Y的长度，m=fix(log2(N))-1.
% 输出的矩阵中，第一列为原始输入数据，第2,3，...，m列是从高频到低频的IMF分量，m+1列是残余分量。
%
% 注意：
% 当Nstd设置为0，NE设置为1时，该EEMD程序就退化为EMD程序。
% 此代码限制筛选数=10，停止标准不能更改。
%
% References:
% Wu, Z., and N. E Huang (2008),
% Ensemble Empirical Mode Decomposition: a noise-assisted data analysis method.
% Advances in Adaptive Data Analysis. Vol.1, No.1. 1-41.
%
% code writer: Zhaohua Wu.
% footnote:S.C.Su 2009/03/04
%
% 这段代码中有三个循环耦合在一起。
% 1.读取数据，找出标准差，用标准差对所有数据进行分离
% 2.将TNM评估为IMF总数--eq1.
% TNM2=TNM+2,original data and residual included in TNM2
% assign 0 to TNM2 matrix
% 3.Do EEMD NE times-----------loop EEMD start
% 4.添加噪声
% 5.筛选前给出初始值
% 6.开始寻找IMF------IMF loop start
% 7.筛选10次以获得IMF------sift loop start and end
% 8.10次筛后 --we got IMF
% 9.从数据中减去IMF，让残差通过循环找到下一个IMF
% 6.在拥有了所有的imf之后-------------IMF loop end
% 9.经TNM-IMFs处理后，残差呈总体趋势
% 3.总结NE分解结果--------loop EEMD end
% 10.Devide EEMD summation by NE,std be multiply back to data
%% Association: no
% this function ususally used for doing 1-D EEMD with fixed
% stoppage criteria independently.
%
% Concerned function: extrema.m
% above mentioned m file must be put together
%function allmode=eemd(Y,Nstd,NE)
%part1.read data, find out standard deviation ,devide all data by std
xsize=length(Y);
dd=1:1:xsize;
Ystd=std(Y);
Y=Y/Ystd;
%part2.evaluate TNM as total IMF number,ssign 0 to N*TNM2 matrix
TNM=fix(log2(xsize))-5; % TNM=m
TNM2=TNM+2;
for kk=1:1:TNM2
for ii=1:1:xsize
allmode(ii,kk)=0.0;
end
end
%part3 Do EEMD -----EEMD loop start
for iii=1:1:NE %EEMD loop NE times EMD sum together
%part4 --Add noise to original data,we have X1
for i=1:xsize
temp=randn(1,1)*Nstd; % add a random noise to Y
X1(i)=Y(i)+temp;
end
%part4 --assign original data in the first column
for jj=1:1:xsize
mode(jj,1) = Y(jj); % assign Y to column 1of mode
end
%part5--give initial 0to xorigin and xend
xorigin = X1; %
xend = xorigin; %
%part6--start to find an IMF-----IMF loop start
nmode = 1;
while nmode <= TNM
xstart = xend; %last loop value assign to new iteration loop
%xstart -loop start data
iter = 1; %loop index initial value
%part7--sift 10 times to get IMF---sift loop start
while iter<=10
[spmax, spmin, flag]=extrema(xstart); %call function extrema
%the usage of spline ,please see part11.
upper= spline(spmax(:,1),spmax(:,2),dd); %upper spline bound of this sift
lower= spline(spmin(:,1),spmin(:,2),dd); %lower spline bound of this sift
mean_ul = (upper + lower)/2; %spline mean of upper and lower
xstart = xstart - mean_ul; %extract spline mean from Xstart
iter = iter +1;
end
%part8--subtract IMF from data ,then let the residual xend to start to find next IMF
xend = xend - xstart;
nmode=nmode+1;
%part9--after sift 10 times,that xstart is this time IMF
for jj=1:1:xsize
mode(jj,nmode) = xstart(jj);
end
end
%part10--after gotten all(TNM) IMFs ,the residual xend is over all trend
% put them in the last column
for jj=1:1:xsize
mode(jj,nmode+1)=xend(jj);
end
%after part 10 ,original + TNM IMFs+overall trend ---those are all in mode
allmode=allmode+mode;
end %part3 Do EEMD -----EEMD loop end
%part11--devide EEMD summation by NE,std be multiply back to data
allmode=allmode/NE;
allmode=allmode*Ystd;
评论2