% EMD.M
%
% G. Rilling, July 2002
%
% computes EMD (Empirical Mode Decomposition) according to:
%
% N. E. Huang et al., "The empirical mode decomposition and the
% Hilbert spectrum for non-linear and non stationary time series analysis,"
% Proc. Royal Soc. London A, Vol. 454, pp. 903-995, 1998
%
% with variations reported in:
%
% G. Rilling, P. Flandrin and P. Gon鏰lv鑣
% "On Empirical Mode Decomposition and its algorithms"
% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
% NSIP-03, Grado (I), June 2003
%
% stopping criterion for sifting :
% at each point : mean amplitude < threshold2*envelope amplitude
% &
% mean of boolean array ((mean amplitude)/(envelope amplitude) > threshold) < tolerance
% &
% |#zeros-#extrema|<=1
%
% inputs: - x : analysed signal (line vector)分析信号（线矢量）
% - t (optional) : sampling times (line vector) (default : 1:length(x))采样时间（行向量）（默认值：1：长度（x））
% - stop (optional) : threshold, threshold2 and tolerance (optional)
% for sifting stopping criterion
% default : [0.05,0.5,0.05] Rilling停止条件：阈值，阈值2和容差（可选）筛选停止标准；默认：[0.05,0.5,0.05]
% - tst (optional) : if equals to 1 shows sifting steps with pause
% if equals to 2 no pause如果等于1，显示暂停的筛选步骤；如果等于2没有暂停
%
% outputs: - imf : intrinsic mode functions (last line = residual)固有模式函数（最后一行=残差）
% - ort : index of orthogonality正交指数
% - nbits : number of iterations for each mode每个模式的迭代次数
%
% calls: - extr (finds extrema and zero-crossings)
% - io : computes the index of orthogonality
function [imf,ort,nbits] = emd(x,t,stop,tst);
% default for stopping默认停止标准
defstop = [0.05,0.5,0.05];
if(nargin==1) %nargin为“number of input arguments”的缩写，用来判断输入变量个数的函数。输入变量只有x
t = 1:length(x);%length(x)为x的总数据量
stop = defstop;
tst = 0;
end
if(nargin==2)
stop = defstop;
tst = 0;
end
if (nargin==3)
tst=0;
end
S = size(x);%返回行数和列数 矢量[ ]
if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
error('x must have only one row or one column')%x必须只有一行或一列
end
if S(1) > 1
x = x';
end
S = size(t);
if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
error('t must have only one row or one column')%t必须只有一行或一列
end
if S(1) > 1
t = t';
end
if (length(t)~=length(x))%不相等
error('x and t must have the same length')%x和t长度必须相等
end
S = size(stop);
if ((S(1) > 1) & (S(2) > 1)) | (S(1) > 3) | (S(2) > 3) | (length(S) > 2)
error('stop must have only one row or one column of max three elements')%stop必须只有一行或一列，且最多三个元素
end
if S(1) > 1
stop = stop';
S = size(stop);
end
if S(2) < 3
stop(3)=defstop(3);%默认停止标准
end
if S(2) < 2
stop(2)=defstop(2);%默认停止标准
end
sd = stop(1);
sd2 = stop(2);
tol = stop(3);
if tst%如果tst不等于0，执行下述语句
%看这个程序，tst是判断语句，如果 tst=1，那么每一次计算IMF信号的过程都以fig的形式显示给用户看，就是如何一步步得到这个IMF分量的过程用图的形式表现出来。
%如果tst=2，那么直接跳过这个过程
figure
end
% maximum number of iterations 最大迭代次数=2000
MAXITERATIONS=2000;
% maximum number of symmetrized points for interpolations 插值的最大对称点数=2
NBSYM = 2;
lx = length(x);
sdt(lx) = 0;
sdt = sdt+sd;%sd为迭代部分的阈值上限
sd2t(lx) = 0;
sd2t = sd2t+sd2;%sd2为剩余部分的阈值上限
% maximum number of extrema and zero-crossings in residual 极值和残差的过零点的最大数量
ner = lx;%极值最大数量
nzr = lx;%过零点最大数量
r = x;
imf = [];
k = 1;
% iterations counter for extraction of 1 mode 单轮迭代次数
nbit=0;
% total iterations counter 总迭代计数器
NbIt=0;
while ner > 2%只要极值数大于2继续执行循环
% current mode
m = r;%原始信号
% mode at previous iteration 上一次迭代
mp = m;
sx = sd+1;%即为评价函数的值
% tests if enough extrema to proceed 测试是否有足够的极值进行
test = 0;
[indmin,indmax,indzer] = extr(m);%得出极小值、极大值、过零点的索引位置
lm=length(indmin);
lM=length(indmax);
nem=lm + lM;
nzm=length(indzer);
j=1;
% sifting loop 筛选循环
while ( mean(sx > sd) > tol | any(sx > sd2) | (abs(nzm-nem)>1)) & (test == 0) & nbit<MAXITERATIONS
%G.Rilling停止准则，极值点与零点数量相差超过1，单轮迭代次数小于2000
if(nbit>MAXITERATIONS/5 & mod(nbit,floor(MAXITERATIONS/10))==0)%mod取余，floor取比它小的整数 nbit大于400，且是200的倍数
disp(['mode ',int2str(k),' nombre d iterations : ',int2str(nbit)])%输出字符串：mode k nombre d iterations :nbit
disp(['stop parameter mean value : ',num2str(s)])
end
% % boundary conditions for interpolations : 插值的边界条件
% %往两边延拓
if indmax(1) < indmin(1)
if m(1) > m(indmin(1))
lmax = fliplr(indmax(2:min(end,NBSYM+1)));%fliplr
lmin = fliplr(indmin(1:min(end,NBSYM)));
lsym = indmax(1);
else
lmax = fliplr(indmax(1:min(end,NBSYM)));
lmin = [fliplr(indmin(1:min(end,NBSYM-1))),1];
lsym = 1;
end
else
if m(1) < m(indmax(1))
lmax = fliplr(indmax(1:min(end,NBSYM)));
lmin = fliplr(indmin(2:min(end,NBSYM+1)));
lsym = indmin(1);
else
lmax = [fliplr(indmax(1:min(end,NBSYM-1))),1];
lmin = fliplr(indmin(1:min(end,NBSYM)));
lsym = 1;
end
end
if indmax(end) < indmin(end)
if m(end) < m(indmax(end))
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
rmin = fliplr(indmin(max(end-NBSYM,1):end-1));
rsym = indmin(end);
else
rmax = [lx,fliplr(indmax(max(end-NBSYM+2,1):end))];
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
rsym = lx;
end
else
if m(end) > m(indmin(end))
rmax = fliplr(indmax(max(end-NBSYM,1):end-1));
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
rsym = indmax(end);
else
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
rmin = [lx,fliplr(indmin(max(end-NBSYM+2,1):end))];
rsym = lx;
end
end
tlmin = 2*t(lsym)-t(lmin);%带有l的参数与左边起点的延拓相关 极小值点索引位置
tlmax = 2*t(lsym)-t(lmax);%极大值点索引位置
trmin = 2*t(rsym)-t(rmin);%带有r的参数与右边终点的延拓相关
trmax = 2*t(rsym)-t(rmax);
% in case symmetrized parts do not extend enough 如果对称部分不够长，则将对称点移至左边起点处
if tlmin(1) > t(1) | tlmax(1) > t(1)
if lsym == indmax(1)
lmax = fliplr(indmax(1:min(end,NBSYM)));
else
lmin = fliplr(indmin(1:min(end,NBSYM)));
end
if lsym == 1
error('bug')
end
lsym = 1;
tlmin = 2*t(lsym)-t(lmin);
tlmax = 2*t(lsym)-t(lmax);
end
if trmin(end) < t(lx) | trmax(end) < t(lx)%如果对称部分不够长，则将对称点移至右边终点处
if rsym == indmax(end)
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
else
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
end
if rsym == lx
error('bug')
end
rsym = lx;
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);
end
mlmax =m(lmax); %左边延拓极大值点的纵坐标值
mlmin =m(lmin);%左边延拓极小值点的纵坐标值
mrmax =m(rmax); %右边延拓极大值点的纵坐标值
mrmin =m(rmin);%右边延拓极小值点的纵坐标值
% definition of envelopes from interpolation 插值包络的定义
envmax = interp1([tlmax t(indmax) trmax],[mlmax m(indmax) mrmax],t,'spline');%极大值点的包络曲线，spline意为三次样条插值
envmin = interp1([tlmin t(indmin) trmin],[mlmin m(indmin) mrmin],t,'spline');%极小值点的包络曲线
% envmax = interp1(t(indmax),m(indmax),t,'spline');%极大值点的包络曲线，spline意为三次样条插值
% envmin = interp1(t(indmin),m(indmin),t,'spline');%极小值点的包络曲线
envmoy = (envmax + envmin)/2;%均值信号
m = m - envmoy;%伪模态信号
[indmin,indmax,indzer] = extr(m);
lm=length(indmin
