emdcommnet.m

emd code with detail commnet

  pause
end
end

%---------------------------------------------------------------------------------------------------
% 演示分解过程（FIX和FIX_H停止准则）
function display_emd_fixe(t,m,mp,r,envmin,envmax,envmoy,nbit,k,display_sifting)
subplot(3,1,1)
plot(t,mp);hold on;
plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r');
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' before sifting']);
set(gca,'XTick',[])
hold  off
subplot(3,1,2)
plot(t,m)
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);
set(gca,'XTick',[])
subplot(3,1,3);
plot(t,r-m)
title('residue');
if display_sifting == 2
  pause(0.01)
else
  pause
end
end

%---------------------------------------------------------------------------------------
% 处理边界条件（镜像法）
function [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym)
% 实数情况下，x = z

lx = length(x);

% 判断极值点个数
if (length(indmin) + length(indmax) < 3)
  error('not enough extrema')
end

% 插值的边界条件
if indmax(1) < indmin(1)	% 第一个极值点是极大值
  if x(1) > x(indmin(1))	% 以第一个极大值为对称中心
    lmax = fliplr(indmax(2:min(end,nbsym+1)));
    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 x(1) < x(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 x(end) < x(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 x(end) > x(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);
tlmax = 2*t(lsym)-t(lmax);
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);
    
% 如果对称的部分没有足够的极值点
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 

% 延拓点上的取值       
zlmax = z(lmax); 
zlmin = z(lmin);
zrmax = z(rmax); 
zrmin = z(rmin);
     
% 完成延拓
tmin = [tlmin t(indmin) trmin];
tmax = [tlmax t(indmax) trmax];
zmin = [zlmin z(indmin) zrmin];
zmax = [zlmax z(indmax) zrmax];

end
    
%---------------------------------------------------------------------------------------------------
% 极值点和过零点位置提取
function [indmin, indmax, indzer] = extr(x,t)

if(nargin==1)
  t = 1:length(x);
end

m = length(x);

if nargout > 2
  x1 = x(1:m-1);
  x2 = x(2:m);
  indzer = find(x1.*x2<0);	% 寻找信号符号发生变化的位置

  if any(x == 0)	% 考虑信号采样点恰好为0的位置
    iz = find( x==0 );  % 信号采样点恰好为0的位置
    indz = [];
    if any(diff(iz)==1) % 出现连0的情况
      zer = x == 0;	% x=0处为1，其它地方为0
      dz = diff([0 zer 0]);	% 寻找0与非0的过渡点
      debz = find(dz == 1);	% 0值起点
      finz = find(dz == -1)-1;  % 0值终点
      indz = round((debz+finz)/2);	% 选择中间点作为过零点
    else
      indz = iz;	% 若没有连0的情况，该点本身就是过零点
    end
    indzer = sort([indzer indz]);	% 全体过零点排序
  end
end

% 提取极值点
d = diff(x);
n = length(d);
d1 = d(1:n-1);
d2 = d(2:n);
indmin = find(d1.*d2<0 & d1<0)+1;	% 最小值
indmax = find(d1.*d2<0 & d1>0)+1;	% 最大值


% 当连续多个采样值相同时，把最中间的一个值作为极值点，处理方式与连0类似
if any(d==0)

  imax = [];
  imin = [];

  bad = (d==0);
  dd = diff([0 bad 0]);
  debs = find(dd == 1);
  fins = find(dd == -1);
  if debs(1) == 1	% 连续值出现在序列开头
    if length(debs) > 1
      debs = debs(2:end);
      fins = fins(2:end);
    else
      debs = [];
      fins = [];
    end
  end
  if length(debs) > 0
    if fins(end) == m	% 连续值出现在序列末尾
      if length(debs) > 1
        debs = debs(1:(end-1));
        fins = fins(1:(end-1));

      else
        debs = [];
        fins = [];
      end
    end
  end
  lc = length(debs);
  if lc > 0
    for k = 1:lc
      if d(debs(k)-1) > 0	% 取中间值
        if d(fins(k)) < 0
          imax = [imax round((fins(k)+debs(k))/2)];
        end
      else
        if d(fins(k)) > 0
          imin = [imin round((fins(k)+debs(k))/2)];
        end
      end
    end
  end

  if length(imax) > 0
    indmax = sort([indmax imax]);
  end

  if length(imin) > 0
    indmin = sort([indmin imin]);
  end

end
end

%---------------------------------------------------------------------------------------------------

function ort = io(x,imf)
% ort = IO(x,imf) 计算正交指数
%
% 输入 : - x    : 分析信号
%        - imf  : IMF信号

n = size(imf,1);

s = 0;
% 根据公式计算
for i = 1:n
  for j = 1:n
    if i ~= j
      s = s + abs(sum(imf(i,:).*conj(imf(j,:)))/sum(x.^2));
    end
  end
end

ort = 0.5*s;
end

%---------------------------------------------------------------------------------------------------
% 函数参数解析
function [x,t,sd,sd2,tol,MODE_COMPLEX,ndirs,display_sifting,sdt,sd2t,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin)

x = varargin{1};
if nargin == 2
  if isstruct(varargin{2})
    inopts = varargin{2};
  else
    error('when using 2 arguments the first one is the analyzed signal X and the second one is a struct object describing the options')
  end
elseif nargin > 2
  try
    inopts = struct(varargin{2:end});
  catch
    error('bad argument syntax')
  end
end

% 默认停止条件
defstop = [0.05,0.5,0.05];

opt_fields = {'t','stop','display','maxiterations','fix','maxmodes','interp','fix_h','mask','ndirs','complex_version'};
% 时间序列，停止参数，是否演示，最大迭代次数，每一轮迭代次数，IMF个数，插值方法，每一轮迭代次数（带条件），mask信号，方向数，是否采用复数模式

defopts.stop = defstop;
defopts.display = 0;
defopts.t = 1:max(size(x));
defopts.maxiterations = 2000;
defopts.fix = 0;
defopts.maxmodes = 0;
defopts.interp = 'spline';
defopts.fix_h = 0;
defopts.mask = 0;
defopts.ndirs = 4;
defopts.complex_version = 2;

opts = defopts;

if(nargin==1)
  inopts = defopts;
elseif nargin == 0
  error('not enough arguments')
end

names = fieldnames(inopts);
for nom = names'
  if ~any(strcmpi(char(nom), opt_fields))
    error(['bad option field name: ',char(nom)])
  end
  if ~isempty(eval(['inopts.',char(nom)])) 
    eval(['opts.',lower(char(nom)),' = inopts.',char(nom),';'])
  end
end

t = opts.t;
stop = opts.stop;
display_sifting = opts.display;
MAXITERATIONS = opts.maxiterations;
FIXE = opts.fix;
MAXMODES = opts.maxmodes;
INTERP = opts.interp;
FIXE_H = opts.fix_h;
mask = opts.mask;
ndirs = opts.ndirs;
complex_version = opts.complex_version;

if ~isvector(x)
  error('X must have only one row or one column')
end

if size(x,1) > 1
  x = x.';
end

if ~isvector(t)
  error('option field T must have only one row or one column')
end

if ~isreal(t)
  error('time instants T must be a real vector')
end

if size(t,1) > 1
  t = t';
end

if (length(t)~=length(x))
  error('X and option field T must have the same length')
end

if ~isvector(stop) || length(stop) > 3
  error('option field STOP must have only one row or one column of max three elements')
end

if ~all(isfinite(x))
  error('data elements must be finite')
end

if size(stop,1) > 1
  stop = stop';
end

L = length(stop);
if L < 3
  stop(3) = defstop(3);
end

if L < 2
  stop(2) = defstop(2);
end

if ~ischar(INTERP) || ~any(strcmpi(INTERP,{'linear','cubic','spline'}))
  error('INTERP field must be ''linear'', ''cubic'', ''pchip'' or ''spline''')
end

% 使用mask信号时的特殊处理
if any(mask)
  if ~isvector(mask) || length(mask) ~= length(x)
    error('masking signal must have the same dimension as the analyzed signal X')
  end

  if size(mask,1) > 1
    mask = mask.';
  end
  opts.mask = 0;
  imf1 = emd(x+mask, opts);
  imf2 = emd(x-mask, opts);
  if size(imf1,1) ~= size(imf2,1)
    warning('emd:warning',['the two sets of IMFs have different sizes: ',int2str(size(imf1,1)),' and ',int2str(size(imf2,1)),' IMFs.'])
  end
  S1 = size(imf1,1);
  S2 = size(imf2,1);
  if S1 ~= S2	% 如果两个信号分解得到的IMF个数不一致，调整顺序
    if S1 < S2
      tmp = imf1;
      imf1 = imf2;
      imf2 = tmp;
    end
    imf2(max(S1,S2),1) = 0;	% 将短的那个补零，达到长度一致
  end
  imf = (imf1+imf2)/2;

end


sd = stop(1);
sd2 = stop(2);
tol = stop(3);

lx = length(x);

sdt = sd*ones(1,lx);
sd2t = sd2*ones(1,lx);

if FIXE
  MAXITERATIONS = FIXE;
  if FIXE_H
    error('cannot use both ''FIX'' and ''FIX_H'' modes')
  end
end

MODE_COMPLEX = ~isreal(x)*complex_version;
if MODE_COMPLEX && complex_version ~= 1 && complex_version ~= 2
  error('COMPLEX_VERSION parameter must equal 1 or 2')
end


% 极值点和过零点的个数
ner = lx;
nzr = lx;

r = x;

if ~any(mask) % 如果使用了mask信号，此时imf就已经计算得到了
  imf = [];
end
k = 1;

% 提取每个模式时迭代的次数
nbit = 0;

% 总体迭代次数
NbIt = 0;
end
%---------------------------------------------------------------------------------------------------