改进版程序ceemd，抑制噪音不错

改进版程序ceemd，加入噪声，抑制噪声，有不错的效果。

function [imftresidtspecwaddmeancountimfnresidn] = ceemdnew（xthrshsthreshNalpha）
% Empiricial Mode Decomposition （Hilbert-Huang Transform）
% [imfresidspecwaphsmeancount] = emd（xthrshsthreshNalpha）
%   Inputs:
%         thrsh=thershold variance [0.2 to 0.3 usually]
 %        sthrh=threshold for first iteration [between 4/8]
%         N=number of trials for the ensemble emd
%         alpha=multiplier for white noise （defaults to .1*std（x） as  suggested by Huang）
%   Outputs:
%         imf= intrinsic mode functions （imf（:end））
%         resid=residue; spec is the spectrum （A*e^jwt） w is the
%         frequency a is the ampitlude phs is the instphase
x   = transpose（x（:））;
imf = [];
k=0;
xt=x;
if nargin<5
%     alpha=.1*std（xt）;
alpha=0.5*std（xt）;% 按照需要修改，huang的建议是0.1，但是0.1会造成端点噪声大。
end
cx=0;
s=0;
i=1;
lim=1;
imfprev=[];
while i    if i~=1
        Noise=alpha.*randn（size（xt））;
    else
        Noise=0;
    end
    for t=1:lim
        count=0;
        x=xt+（-1）^t*Noise;
        while ~ismonotonic（x）
            count=count+1;
            if i==1
                k=k+1;
                imf{count}=zeros（size（x））;
                meancount（k）=0;
            end
           x1 = x;
           VA = Inf;
           cn=0;
           while （VA > thrsh） | ~isimf（x1）
               cn=cn+1;
              s1 = getspline（x1）;
              s2 = -getspline（-x1）;
              x2 = x1-（s1+s2）/2;
              VA = var（（s1+s2）./（2*x1））;
              x1 = x2;
              if i==1 cx=cx+1; elseif cn>cx*2 break; end
              if i==1 && cx>numel（xt）/2
                  s=s+1;
                  if s>sthresh
                      x1=xt+alpha.*randn（size（xt））;
                      s=0;
                      cx=0;
                      disp（‘Restarting Sifting‘）;
                  end
              end
           end
%            if cn>cx*2 && i~=1
%                disp（‘Too long for ensemble member breaking count‘）
%                flag=1;
%                break;
%            else
%                flag=0;
%            end
           
           if count<=k 
               if i==1 
                   imf{count} = x1+imf{count}; 
               else
                   imf{tcount}=x1+imf{tcount}; 
               end;
           else
               break;
           end
           x=x-x1;
        end
        
        if i==1
            resid=repmat（x21）;
            imf=repmat（imf21）; 
%         elseif flag==0
        else
            resid（t:） = x+resid（t:）;
            if t==lim
                meancount=meancount+1;
            end
%             residprev=x;
%         else
%             if t==lim
%                 resid（1:）=resid（1:）-residprev;
%                 for l=1:numel（imf（1:））
%                     imf{1l}=imf{1l}-imfprev{l};
%                 end
%             end
%             i=i-1;
        end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        4623  2019-05-13 11:39  ceemd （1）.m