example1.m

fastest algorithm to find EMD.

function imf = example1(t,x);
% imf = example1;
% imf = example1(t,x);
%
%   x = time-series to decompose
%   t = times for each sample in x
%
% author: Brad M. 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
%

% Create a time series if needed
if nargin ~= 2
    t = [0:.001:1]';
    x = sum(cos(2*pi*[2 20 100]'*t'))';
end
tn=length(t);

t1=min(t(1),t(tn));
t2=max(t(1),t(tn));

% Perform the EMD
imf = semd(x);

%%%%%%%%%%%%%%%%%%% Plot the EMD %%%%%%%%%%%%%%%%%%%
figure;
clf;

imfStart = 1;
imfStop = size(imf,2);
imfLast = imfStop+1;

for j=imfStart:imfLast;
    subplot(imfLast-imfStart+1,1,j+1-imfStart);
    if j==imfStart
        z=sum(imf(:,imfStart:imfStop),2);
        plot(t,z);
        yMax = max(z);
        yMin = min(z);
    else
        plot(t,imf(:,j-1));
        yMax = max(imf(:,j-1));
        yMin = min(imf(:,j-1));
    end
    if( abs(yMax-yMin) > 1e-10 )
        axis([t1 t2 yMin yMax]);
    else
        axis([t1 t2 -1 1]);
    end
   s1='c';
   s2=[s1 num2str(j-1)];
   ylabel(s2);
   if( j == 1 )
       title( 'Time Series and IMF Components' );
   end
   if( j < imfLast );
      set(gca,'xticklabel','');
   end;
end;
xlabel('Time');


%%%%%%%%%%%%%%%%%%% Plot the Hilbert Spectrum %%%%%%%%%%%%%%%%%%%

% Perform the Hilbert transform and 
% get instantaneous attributes
[omega,amp,mag,phase,w] = fhilbert(imf,t1,t2,[]);

% Plot the Hilbert Spectrum
[m,n] = size(mag);
h     = full(sum(mag,2)/n);

% Need time,frequency and amplitude to be column vectors
T = repmat(t(:),size(imf,2),1);
O = omega(:);
A = amp(:);

% Sum amp for identical freqs in different IMFs
stoa = sortrows([T O A],[2 1]);
T = stoa(:,1);
O = stoa(:,2);
A = stoa(:,3);
myeps = max(max(abs(T)),max(abs(O)))*eps^(1/3);
ind = [0; ((abs(diff(O)) < myeps) & (abs(diff(T)) < myeps)); 0];
if sum(ind) > 0
  disp(['Summing amplitudes at identical frequencies.']);
  fs = find(ind(1:end-1) == 0 & ind(2:end) == 1);
  fe = find(ind(1:end-1) == 1 & ind(2:end) == 0);
  for n = 1 : length(fs)
    % summing z values
    A(fe(n)) = sum(A(fs(n):fe(n)));
  end
  T = T(~ind(2:end));
  O = O(~ind(2:end));
  A = A(~ind(2:end));
end

figure;
ax1 = subplot(3,1,1);plot(t,sum(imf,2));
ylabel('Amplitude');set(gca,'XTickLabel','','xlim',[min(t) max(t)]);
grid on;box on;

ax2 = subplot(3,1,2);scatter(T,O,10,20*log10(A),'o','filled');
set(ax2,'ydir','normal','xlim',[min(t) max(t)]);

xlabel('Time');ylabel('Frequency');
%set(gca,'yscale','log');
grid on;box on

ax4 = colorbar('horiz');axes(ax4);
xlabel('Power (dB)')
grid on;box on;

ax3 = subplot(3,1,3);plot(h,w);
xlabel('Amp');set(gca,'YTickLabel','');
grid on;box on

set(ax1,'Position',[0.1300 0.5822 0.6770 0.3222]);
set(ax2,'Position',[0.1300 0.2100 0.6770 0.3222]);
set(ax3,'Position',[0.8480 0.2100 0.0770 0.3222]);
set(ax4,'Position',[0.1300 0.1000 0.6770 0.0100]);

% Plot the regenerated signal against the original signal
sig = real((sum(amp.*exp(i*phase),2)));%+imf(:,end);
figure;plot(t,[sum(imf,2) sig]);
xlabel('Time');
legend('\Sigma IMFs','\Sigma a_je^i^\Theta^t');