emd_intermit.m

一种新的时频分析方法的matlab源程序。

function imf=emd_intermit(x,s)
% imf=emd_intermit(x,s): EMD with sign intermittency

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

nx = max(size(x));
imf = zeros(0, nx);

nIMF = 0;

intermittencyActive = 1;

fig1 = gcf; % save current figure
fig2 = figure; % get a new figure

while any(x)
    if min(x) == max(x)
        nIMF = nIMF+1;
        imf(nIMF,:) = c;
        break;
    end
    c = emdstep(x,s);
    if intermittencyActive
        figure(fig1)
        subplot(3,1,1);
        plot(c);
        axis([0 nx-1 min(c) max(c)]);
        title('Sifted component');
        
        subplot(3,1,2);
        [dist,distidx] = getSignConsec(c);
        %hist(dist);
        plot(distidx, dist);
        if any(dist)
            axis([0 nx-1 0 max(dist)]);
            title('Zero-crossing spacing distribution');
            
            disp('');
            disp(' Interactive intermittency');
            disp(' -------------------------');
            disp('');
            disp(' Press Enter to accept the component');
            disp(' Select a positive number to set intermittency threshold.');
            disp(' ''0'' to cancel intermittency for this component');
            disp(' ''-1'' to cancel intermittency for this and all successive components');
            disp(' ''-2'' to abort the EMD and store the remaining signal as the last component');
            disp('');
            
            disp(['Spacing: mean=' num2str(mean(dist)) ', std=' num2str(std(dist))]);
            
            ic = c;
            reply = 1;
            while reply > 0
                figure(fig1)
                subplot(3,1,3);
                plot(ic);
                axis([0 nx-1 min(ic) max(ic)]);
                title('Candidate component');
                
                reply = input('Enter to accept or select a new parameter: ');
                if isempty(reply)
                    reply = -3;
                elseif reply == 0
                    ic = c;
                elseif reply == -1
                    ic = c;
                    intermittencyActive = 0;
                elseif reply == -2
                    ic = x;
                elseif reply > 0
                    ic = signIntermittency(c, reply);
                else
                    disp('Invalid negative reply.');
                    disp('');
                end
            end % intermittency loop
            
            c = ic; % update component
        end
    end % intermittency active
    
    nIMF = nIMF+1;
    imf(nIMF,:) = c;
    x = x - c;
    
    figure(fig2)
    if max(x) ~= min(x)
        plotcs([imf; x]);
    else
        plotcs(imf);
    end
end