emd.m

emd经验模态分析 很好的信号分析包 希望大家支持

function [IMF] = emd(signal,options);

% Computes the IMFs using variants of EMD.
% Signal: Input signal to get decomposed.
% options: Structure having the proper parameters generated with the aid of
%         emdoptimset function. Use 'help emdoptimset', or just
%         'emdoptimset' in Matlab command line to see details.
% IMF: A Matrix containing the resulting IMFs in its lines.
% Examples:
%1) x=randn(1,1000);IMF=emd(x);
%2) options=emdoptimset('stop','Double_T','T',[0.05,0.5,0.05]);
%   x=randn(1,1000);IMF=emd(x,options);
%
% Ver. 1.0
% Yannis Kopsinis, 23/04/2008
% [e-mail]
%
% Some code (quite many actually) of the functions
% "edge_extrapolation" and "extremes" has been 'borrowed'
% from the first version of Rilling emd.m code. 
% His new emd versions with extra functionings can be found in 
% http://perso.ens-lyon.fr/patrick.flandrin/emd.html

t=1:length(signal);
if size(signal,2)==1, signal=signal.'; end
if nargin==1, options=emdoptimset; end


if strcmpi(options.Nodes,'DI_Extrema')
    if isempty(options.in_N) || isempty(options.in_Nskip) || isempty(options.in_Order) 
        error('The values of ''in_N'' and/or ''in_Nskip'' and/or ''in_Order'' are wrong.')
    end
end

if strcmpi(options.Stopping,'Single_T')
    if length(options.T)~=1
        error('The value of ''T'' is wrong.')
    end
end

if strcmpi(options.Stopping,'Double_T')
    if length(options.T)~=3
        error('The value of ''T'' is wrong.')
    end
end


NofExtr = inf; %nofextr: Number of Extrema
IMFnum = 1; %Current IMF number.


totalstopping=0;
while ~totalstopping & IMFnum<=options.IMFs
    tic
    
    iter = 0;
    h = signal; %Initial IMF value.
    Stoppingflag = 0; %Stoppinflag takes value 1 when the IMF criterion of options.stopping have been fulfilled
    zeromeaned = 0;

    while ~Stoppingflag & iter <= options.MaxN
        iter = iter + 1;

        options.iter = iter;
        siftdetails = sifting(h,t,options);
        
        if siftdetails.valid==1
            m = siftdetails.m;
            h = h - m;
        end
        
        [Stoppingflag,totalstopping,zeromeaned]=stoppingcriteria(siftdetails,h,options,zeromeaned,iter);

        % Display progress
        if options.Disp==1
            a=toc;
            if a>3
                fprintf('   sifting number: %d \n', iter)
                tic
            end
        end
    end

    if iter==1
        disp('')
    end
    
    if options.Disp==1
       fprintf('IMF number %d has been computed with %d sifting iterations\n',IMFnum,iter-1) 
    end
    
    IMF(IMFnum,:) = h;

    signal = signal - h;

    NofExtr = siftdetails.NofExtr;
    IMFnum = IMFnum + 1;
end



%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function siftdetails=sifting(signal,t,options);

[maxi, mini, zerocross] = extremes(signal,t,options);

NofExtr = length(maxi)+length(mini);
siftdetails.NofExtr=NofExtr;
siftdetails.valid=1;
if length(maxi) >= 2 & length(mini) >= 2
    N_extrap=floor(options.Order/2)+1;
    [maxp_ex, maxv_ex, minp_ex, minv_ex] = edge_extrapolation(maxi, mini, signal, t,N_extrap);

    %The total number of extrapolated maxima and minima should be larger than the spline
    %order (It takes effect for spline orders higher than 3.)
    if length(maxp_ex) <= options.Order | length(minp_ex) <= options.Order
        m=zeros(1,length(signal));
        siftdetails.valid=0;
        return
    end
else
    m=zeros(1,length(signal));
    siftdetails.valid=0;
    return
end


envmax = interpolation(maxp_ex,maxv_ex,t,options); % envelop of the maxima
envmin = interpolation(minp_ex,minv_ex,t,options); % envelop of the minima



siftdetails.m = 1/2*(envmax + envmin);
siftdetails.envmax=envmax;
siftdetails.envmin=envmin;
siftdetails.maxi=maxi;
siftdetails.mini=mini;
siftdetails.zerocross=zerocross;
siftdetails.NofExtr=NofExtr;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [maxima , minima, zerocross] = extremes(signal,t,options)

if nargin < 3
    kind = 'extrema';
    differentiation='diff_2p';
else
    kind = lower(options.Nodes);
    differentiation=options.Differentiation;
end

if strcmpi(lower(kind),'di_extrema')

    N_extrap=floor(options.in_Order/2)+1;

    if mod(options.iter,options.in_Nskip)~=0 & options.iter~=1
        kind='extrema';
    else


        h=t(2)-t(1);

        dsignal=diff_5p(signal,h); %first derivative of signal.
        options.Nodes='extrema'; %It is used in the "extremes" function for the internal sifting iterations
        options.Order=options.in_Order; %It is used in the "interpolation" function for the internal sifting iterations
        for nofsift=1:options.in_N
            
            [maxifd , minifd, zcrsfd] = extremes(dsignal,t,options);
            
            if length(maxifd)>=ceil(options.in_Order/2) & length(minifd)>=ceil(options.in_Order/2)
                [maxp_ex, maxv_ex, minp_ex, minv_ex] = edge_extrapolation(maxifd, minifd, dsignal, t, N_extrap);
            else
                break % for the rare case that some extrema in the dfastsignal get lost and lead to less than the required extrema for the specific interpolation order.
            end
            envmax = interpolation(maxp_ex,maxv_ex,t,options); % envelop of the max
            envmin = interpolation(minp_ex,minv_ex,t,options); % envelop of the min

            m = 1/2*(envmax + envmin);
            dsignal=dsignal-m;
        end
        if nofsift==1
            kind='extrema';
        else
            
            dif = sign(dsignal);

            % It considers the appearence of single zeros only.
            if any(dif==0)
                zdif=find(dif==0);      %in case of zero the sign of the previous
                if zdif(1)==1, zdif(1)=[]; end
                if zdif(end)==length(dif), zdif(end)=[]; end
                dif(zdif)=dif(zdif-1);  %point is adopted.
            end

            N = length(dif);

            sdif = dif(1:N-1).*dif(2:N);
            extremespos = find(sdif==-1);
            extremestype = dif(extremespos);

            maxima = extremespos(find(extremestype>0))+1;
            minima = extremespos(find(extremestype<0))+1;
            zerocross = sort([maxifd, minifd]);
            if signal(maxima(1))<signal(minima(1))
                tmp=maxima(1);maxima(1)=minima(1);minima(1)=tmp;
            end
            if signal(maxima(end))<signal(minima(end))
                tmp=maxima(end);maxima(end)=minima(end);minima(end)=tmp;
            end
        end
    end
end




if strcmpi(kind,'extrema')
    if strcmpi(differentiation,'diff_2p')==1
         dif = diff(signal); %approximate derivative
    else
        h=t(2)-t(1);
        dif=diff_5p(signal,h);
    end        
  
    dif = sign(dif);

    % It considers the appearence of single zeros only.
    if any(dif==0)
        zdif=find(dif==0);      %in case of zero the sign of the previous
        if zdif(1)==1, zdif(1)=[]; end
        if zdif(end)==length(dif), zdif(end)=[]; end
        dif(zdif)=dif(zdif-1);  %point is adopted.
    end

    N = length(dif);

    sdif = dif(1:N-1).*dif(2:N);
    extremespos = find(sdif==-1);
    extremestype = dif(extremespos);

    maxima = extremespos(find(extremestype>0))+1;
    minima = extremespos(find(extremestype<0))+1;

    zerocross = find(signal(1:N-1).*signal(2:N)<0);




end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function YI=interpolation(X,Y,XI,options)

h = diff(X);
if any(h == 0)
    [Kx,Lx]=find(h==0);
    X(Lx)=[];
    Y(Lx)=[];
end

if options.Order~=3
    order = options.Order + 1;
    kovw=order/2-1;
    knotseqmax=[X(1)*ones(1,order), X(3+kovw-1:end-1-kovw), X(end)*ones(1,order)];
    spout = spapi(knotseqmax,X,Y);
    YI = fnval(spout,XI);

elseif  options.Order==3
    YI = interp1(X,Y,XI,'spline');

end