M File_Empirical Mode Decomposition

Hello fellas!

In data analysis techniques, a new method had discovered by N.E. Huang, by combining Hilbert - Huang Transformation (HHT), which is named as Empirical Mode Decomposition (EMD). This method decompose a signal into its Intrinsic Mode Function (IMF) form. This method applies to all types of data stationary or nonstationary and nonlinear. Unlike the Fourier Transformation, HHT is more of an algorithm (empirical approach) which can be applied to the data set.

By using the EMD method, any complicated data set can be decomposed into a finite or limited, and often a small number of components. These components are formed approaching the basis of the orthogonal of original signal, or can be called with the IMF. Without prejudice to the domain time, EMD is a method that is adaptive and having high efficiency. Because decomposition based on the characteristic scale of local time of the data, the EMD can be used to nonlinear and nonstationary processes.While to ensure the IMF it is said to be qualified as following:1. In the entire data set, the number of extrema are equal to the number of zero-crossing or has a difference of at least one.2. At any point, the mean or average of the envelope defined by the local maxima and the envelope defined by local minima is zero.

Here is the m-file for plotting HHT or IMF (1) and for processing the EMD (2) :
(1)
function plot_hht(x,Ts)
% Plot the HHT.
% plot_hht(x,Ts)
%
% :: Syntax
%    The array x is the input signal and Ts is the sampling period.
%    Example on use: [x,Fs] = wavread('Hum.wav');
%                    plot_hht(x(1:6000),1/Fs);
% Func : emd

% Get HHT.
Ts=1/5000;
imf = emd(x);
for k = 1:length(imf)
   b(k) = sum(imf{k}.*imf{k});
   th   = angle(hilbert(imf{k}));
   d{k} = diff(th)/Ts/(2*pi);
end
[u,v] = sort(-b);
b     = 1-b/max(b);

% Set IMF plots.
M = length(imf);
N = length(x);
c = linspace(0,(N-1)*Ts,N);
for k1 = 0:4:M-1
   figure
   for k2 = 1:min(4,M-k1), subplot(4,1,k2), plot(c,imf{k1+k2}); set(gca,'FontSize',8,'XLim',[0 c(end)]); end
   xlabel('Time');
end


(2)
function imf = emd(x)
% Empiricial Mode Decomposition (Hilbert-Huang Transform)
% imf = emd(x)
% Func : findpeaks

x   = transpose(x(:));
imf = []
figure(1)
hold on
sdmat=[];
while ~ismonotonic(x)
   x1 = x;
   sd = Inf;
  
   while (sd > 0.1) | ~isimf(x1)
      s1 = getspline(x1);
      s2 = -getspline(-x1);
      x2 = x1-(s1+s2)/2;
     
      sd = sum((x1-x2).^2)/sum(x1.^2)
      sdmat=[sdmat,sd];
      %figure (1)
      %plot(sd,'bO')
      x1 = x2;
   end
 
   imf{end+1} = x1;
   x          = x-x1;
end
sdmat
hold off
figure
plot (sdmat,'bo')
[f,xi] = ksdensity(sdmat);
figure
area(xi,f);
title('PDF')
imf{end+1} = x;

% FUNCTIONS

function u = ismonotonic(x)

u1 = length(findpeaks(x))*length(findpeaks(-x));
if u1 > 0, u = 0;
else,      u = 1; end

function u = isimf(x)

N  = length(x);
u1 = sum(x(1:N-1).*x(2:N) < 0);
u2 = length(findpeaks(x))+length(findpeaks(-x));
if abs(u1-u2) > 1, u = 0;
else,              u = 1; end

function s = getspline(x)

N = length(x);
p = findpeaks(x);
s = spline([0 p N+1],[0 x(p) 0],1:N);

Hope it will help you fellas!
Good luck

sources:

 R. Valles-Novo, J. Rangel-Magdaleno, J. Ramirez-Cortes, H. Peregrina-Barreto, R. Morales-Carporal, “Empirical Mode Decomposition Analysis for Broken-Bar Detection on Squirrel Cage Induction Motors”, IEEE Transactions on Instrumentation and Measurement, December, 2014.

N. E. Huang, Zheng Shen, Steven R. Long, M. C. Wu, H. H. Shih, Quanan Zheng, Nai-Chyuan Yen, C. C. Tung, H. H. Liu, “The Empirical Mode Decomposition and The Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis”, The Royal Society, Proc. R. Soc. Lond. A 454, 903-995, England, 1998.