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 neededif nargin ~= 2    t = [0:.001:1]';    x = sum(cos(2*pi*[2 20 100]'*t'))';endtn=length(t);t1=min(t(1),t(tn));t2=max(t(1),t(tn));% Perform the EMDimf = 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 vectorsT = repmat(t(:),size(imf,2),1);O = omega(:);A = amp(:);% Sum amp for identical freqs in different IMFsstoa = 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));endfigure;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 onax4 = 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 onset(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 signalsig = 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');