% EMD: Emprical mode decomposition(对信号进行EMD分解，并将各个imf绘出，以矩阵形式返回各个imf)
%
% imf = emd(x)
%
% x - input signal (must be a column vector or a row vector)
%
% This version will calculate all the imf's (longer)
%
% imf - Matrix of intrinsic mode functions (each as a row)
% with residual in last row.
%
% See: Huang et al, Royal Society Proceedings on Math, Physical,
% and Engineering Sciences, vol. 454, no. 1971, pp. 903-995,
% 8 March 1998
%
% Author: Ivan Magrin-Chagnolleau <[email protected]>
%
function imf = emd1(x);
tic;
num=size(x);
if num(1)==1
c=x;
else
c = x(:)'; % copy of the input signal (as a row vector)
end
N = length(x);
%-------------------------------------------------------------------------
% loop to decompose the input signal into successive IMF
imf = []; % Matrix which will contain the successive IMF, and the residue
while (1) % the stop criterion is tested at the end of the loop
%-------------------------------------------------------------------------
% inner loop to find each imf
h = c; % at the beginning of the sifting process, h is the signal
SD = 1; % Standard deviation(背离) which will be used to stop the sifting process
while SD > 0.3
% while the standard deviation is higher than 0.3 (typical value)
% find local max/min points
d = diff(h); % approximate derivative
maxmin = []; % to store the optima（最适宜的） (min and max without distinction so far)
for i=1:N-2%注：默认步长是1
if d(i)==0 % we are on a zero注：如果d(i)等于零，则是极值点的位置。
maxmin = [maxmin, i];
elseif sign(d(i))~=sign(d(i+1)) % we are straddling a zero so注：sign是用来判断一个值与零的关系，返回值有－1，0，1三种情况。
maxmin = [maxmin, i+1]; % define zero as at i+1 (not i)
end
end
if size(maxmin,2) < 2 % then it is the residue
break
end
% divide maxmin into maxes and mins注意此处区别极大值与极小值的方法。很简单。
if maxmin(1)>maxmin(2) % first one is a max not a min
maxes = maxmin(1:2:length(maxmin));
mins = maxmin(2:2:length(maxmin));
else % is the other way around
maxes = maxmin(2:2:length(maxmin));
mins = maxmin(1:2:length(maxmin));
end
% make endpoints both maxes and mins
maxes = [1 maxes N];
mins = [1 mins N];%此处是用maxes/mins加上两端点重新构成新的行向量。二者起点和终点相同。
%-------------------------------------------------------------------------
% spline interpolate（插值） to get max and min envelopes(包络); form imf此处用三次样条曲线拟合得到上下包络
maxenv = spline(maxes,h(maxes),1:N);
minenv = spline(mins, h(mins),1:N);%此处能直接用spline命令来得到上下包络吗？怎么调用此命令。
m = (maxenv + minenv)/2; % mean of max and min enveloppes简单的求上下包络的平均值
prevh = h; % copy of the previous value of h before modifying it后面求r(t)残差时还用的到原始信号。
h = h - m; % substract mean to h
% calculate standard deviation
eps = 0.0000001; % to avoid zero values
SD = sum ( ((prevh - h).^2) ./ (prevh.^2 + eps) );%此处好像是用方差之类的东西来得到判断条件。不是！
end
imf = [imf; h]; % store the extracted IMF in the matrix imf
% if size(maxmin,2)<2, then h is the residue
% stop criterion of the algo.
if size(maxmin,2) < 2
break
end
c = c - h; % substract the extracted IMF from the signal
end
a=1;
b=1;
c=size(imf);
while b<=c(1)
if a==1
k=2;
figure(a);
subplot(5,1,1)
plot(x);
title('原始信号');
while k<=5
if b>c(1)
break;
end
subplot(5,1,k)
plot(imf(b,:));
k=k+1;
b=b+1;
end
else
k=1;
while k<=5
if b>c(1)
break;
end
figure(a);
subplot(5,1,k)
plot(imf(b,:));
k=k+1;
b=b+1;
end
end
a=a+1;
end
toc;
return
评论0