imfpick.m

fastest algorithm to find EMD.

function [xx,index] = IMFpick(x,t,cumimf,tpow);
% [y,idx] = IMFpick(x,t,cumimf,tpow) is a trace picker.  The cumulative EMD 
% for all traces in 'x' from MAT_EMD is used to determine how many IMFs
% are required to produce the desired trace. 'tpow' scales the display
% by t^tpow
% 
% author: Bradley Matthew Battista
%   University of South Carolina
%   Department of Geological Sciences
%   701 Sumter Street, EWS 617
%   Columbia, SC. 29208
%
% COPYRIGHT: see the associated COPYRIGHT.txt file, and also
% http://software.seg.org/disclaimer2.txt
% This source code may be found online at:
% http://software.seg.org/2007/0003
%

% determine how many IMFs represent each trace
numIMF = size(cumimf,2)/size(x,2);

% preallocate memory for the output traces
xx = zeros(size(x));
index = zeros(1,size(x,2));

% use a user-provided time-gain (or not)
if ~isempty(tpow)
    scale = repmat(t.^tpow,1,numIMF);
else
    scale = 1;
end

% loop through the cumulative-sum IMFs
for n=1:size(x,2)
    cla;
    idx = (n-1)*numIMF+1;

    % plot the cumulative-sum IMFs for a given trace
    y = cumimf(:,idx:idx+numIMF-1).*scale;
    for m = 1:numIMF
        plot(y(:,m)/max(y(:,m))+m,t);hold on
    end
    
    % make t0 the top of the graph and set axis limits
    set(gca,'ydir','reverse',...
        'xlim',[0 numIMF+1],'ylim',[min(t) max(t)]);
    title(num2str(n));
        
    % allow user to click the desired output trace from the figure
    [buf1,buf2] = ginput(1);set(gcf,'Pointer','crosshair')
    
    % get the trace number and output it
    index(n) = round(buf1);
    xx(:,n) = cumimf(:,idx+round(buf1)-1);
end