Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

New EMD routine #1

Merged
merged 3 commits into from

3 participants

@Cadair

Hi,

Myself and my friend have had need of a fully functional EMD routine, therefore we expanded your initial code, a lot. We implemented a few different approaches to the cubic spline problem, but the mirror approach defiantly seemed the best.

I have cleaned up this function as much as I can, it should be pretty PEP8 etc.

I think our next move will be to write a Hilbert spectrum routine, then I will probably recreate the pyhht.py file.

Any suggestions etc ?

Stuart

@Cadair

Dropped the extrap method because it made the code ugly and is no where near as good as the mirror method.

@jaidevd jaidevd merged commit 76f2590 into jaidevd:master
@jaidevd
Owner
@jaidevd
Owner

Also, I don't know if I merged the request correctly. Don't even know what that means. Hope it has worked out fine and I didn't mess anything up.

@Cadair

I have a good plotting routine if you want it. I have pushed it to my fork, there is a example data generation test routine you can use.

It appears to be pulled ok!! We should decide where we want to go with this, it would be excellent if we could get it pulled into SciPy eventually.

@jaidevd
Owner
@raininja

hey, is anyone watching this thread? I want to use some of this code for a project! I have some tyro questions :dart:

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Apr 18, 2012
  1. @Cadair

    New EMD routine, comes with 3 different types of spline fitting, one …

    Cadair authored
    …standard, one that uses a Mirror and one that extrapolates
Commits on Apr 19, 2012
  1. @Cadair
  2. @Cadair

    Updated README

    Cadair authored
This page is out of date. Refresh to see the latest.
Showing with 192 additions and 135 deletions.
  1. +191 −0 EMD.py
  2. +1 −1  README
  3. +0 −134 scipy.tex
View
191 EMD.py
@@ -0,0 +1,191 @@
+"""
+To do:
+ Tests / Examples (same in some literature)
+ Instanous freequncuies
+ Hilbert-Huang Transform
+"""
+import numpy as np
+from scipy import interpolate
+
+__all__ = 'EMD'
+
+def emd(data, extrapolation='mirror', nimfs=12, shifting_distance=0.2):
+ """
+ Perform a Empirical Mode Decomposition on a data set.
+
+ This function will return an array of all the Imperical Mode Functions as
+ defined in [1]_, which can be used for further Hilbert Spectral Analysis.
+
+ The EMD uses a spline interpolation function to approcimate the upper and
+ lower envelopes of the signal, this routine implements a extrapolation
+ routine as described in [2]_ as well as the standard spline routine.
+ The extrapolation method removes the artifacts introduced by the spline fit
+ at the ends of the data set, by making the dataset a continuious circle.
+
+ Parameters
+ ----------
+ data : array_like
+ Signal Data
+ extrapolation : str, optional
+ Sets the extrapolation method for edge effects.
+ Options: None
+ 'mirror'
+ Default: 'mirror'
+ nimfs : int, optional
+ Sets the maximum number of IMFs to be found
+ Default : 12
+ shifiting_distance : float, optional
+ Sets the minimum variance between IMF iterations.
+ Default : 0.2
+
+ Returns
+ -------
+ IMFs : ndarray
+ An array of shape (len(data),N) where N is the number of found IMFs
+
+ Notes
+ -----
+
+ References
+ ----------
+ .. [1] Huang H. et al. 1998 'The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis.'
+ Procedings of the Royal Society 454, 903-995
+
+ .. [2] Zhao J., Huang D. 2001 'Mirror extending and circular spline function for empirical mode decomposition method'
+ Journal of Zhejiang University (Science) V.2, No.3,P247-252
+
+ .. [3] Rato R.T., Ortigueira M.D., Batista A.G 2008 'On the HHT, its problems, and some solutions.'
+ Mechanical Systems and Signal Processing 22 1374-1394
+
+
+ """
+
+ #Set up signals array and IMFs array based on type of extrapolation
+ # No extrapolation and 'extend' use signals array which is len(data)
+ # Mirror extrapolation (Zhao 2001) uses a signal array len(2*data)
+ if not(extrapolation):
+ base = len(data)
+ signals = np.zeros([base, 2])
+ nimfs = range(nimfs) # Max number of IMFs
+ IMFs = np.zeros([base, len(nimfs)])
+ ncomp = 0
+ residual = data
+ signals[:, 0] = data
+ #DON'T do spline fitting with periodic bounds
+ inter_per = 0
+
+ elif extrapolation == 'mirror':
+ #Set up base
+ base = len(data)
+ nimfs = range(nimfs) # Max number of IMFs
+ IMFs = np.zeros([base, len(nimfs)])
+ ncomp = 0
+ residual = data
+ #Signals is 2*base
+ signals = np.zeros([base*2, 2])
+ #Mirror Dataset
+ signals[0:base / 2, 0] = data[::-1][base / 2:]
+ signals[base / 2:base + base / 2, 0] = data
+ signals[base + base / 2:base * 2, 0] = data[::-1][0:base / 2]
+ # Redfine base as len(signals) for IMFs
+ base = len(signals)
+ data_length = len(data) # Data length is used in recovering input data
+ #DO spline fitting with periodic bounds
+ inter_per = 1
+
+ else:
+ raise Exception(
+ "Please Specifiy extrapolation keyword as None or 'mirror'")
+
+ for j in nimfs:
+# Extract at most nimfs IMFs no more IMFs to be found when Finish is True
+ k = 0
+ sd = 1.
+ finish = False
+
+ while sd > shifting_distance and not(finish):
+ min_env = np.zeros(base)
+ max_env = min_env.copy()
+
+ min_env = np.logical_and(
+ np.r_[True, signals[1:,0] > signals[:-1,0]],
+ np.r_[signals[:-1,0] > signals[1:,0], True])
+ max_env = np.logical_and(
+ np.r_[True, signals[1:,0] < signals[:-1,0]],
+ np.r_[signals[:-1,0] < signals[1:,0], True])
+ max_env[0] = max_env[-1] = False
+ min_env = min_env.nonzero()[0]
+ max_env = max_env.nonzero()[0]
+
+ #Cubic Spline by default
+ order_max = 3
+ order_min = 3
+
+ if len(min_env) < 2 or len(max_env) < 2:
+ #If this IMF has become a straight line
+ finish = True
+ else:
+ if len(min_env) < 4:
+ order_min = 1 #Do linear interpolation if not enough points
+
+ if len(max_env) < 4:
+ order_max = 1 #Do linear interpolation if not enough points
+
+#==============================================================================
+# Mirror Method requires per flag = 1 No extrapolation requires per flag = 0
+# This is set in intial setup at top of function.
+#==============================================================================
+ t = interpolate.splrep(min_env, signals[min_env,0],
+ k=order_min, per=inter_per)
+ top = interpolate.splev(
+ np.arange(len(signals[:,0])), t)
+
+ b = interpolate.splrep(max_env, signals[max_env,0],
+ k=order_max, per=inter_per)
+ bot = interpolate.splev(
+ np.arange(len(signals[:,0])), b)
+
+ #Calculate the Mean and remove from the data set.
+ mean = (top + bot)/2
+ signals[:,1] = signals[:,0] - mean
+
+ #Calculate the shifting distance which is a measure of
+ #simulartity to previous IMF
+ if k > 0:
+ sd = (np.sum((np.abs(signals[:,0] - signals[:,1])**2))
+ / (np.sum(signals[:,0]**2)))
+
+ #Set new iteration as previous and loop
+ signals = signals[:,::-1]
+ k += 1
+
+ if finish:
+ #If IMF is a straight line we are done here.
+ IMFs[:,j]= residual
+ ncomp += 1
+ break
+
+ if not(extrapolation):
+ IMFs[:,j] = signals[:,0]
+ residual = residual - IMFs[:,j]#For j==0 residual is initially data
+ signals[:,0] = residual
+ ncomp += 1
+
+ elif extrapolation == 'mirror':
+ IMFs[:,j] = signals[data_length / 2:data_length
+ + data_length / 2,0]
+ residual = residual - IMFs[:,j]#For j==0 residual is initially data
+
+ #Mirror case requires IMF subtraction from data range then
+ # re-mirroring for each IMF
+ signals[0:data_length / 2,0] = residual[::-1][data_length / 2:]
+ signals[data_length / 2:data_length + data_length / 2,0] = residual
+ signals[data_length
+ + data_length / 2:,0] = residual[::-1][0:data_length / 2]
+ ncomp += 1
+
+ else:
+ raise Exception(
+ "Please Specifiy extrapolation keyword as None or 'mirror'")
+
+ return IMFs[:,0:ncomp]
View
2  README
@@ -1 +1 @@
-DUMMY README
+A Hilbert Haung Transform package, should perform EMD as well as generate a Hilbert Spectra.
View
134 scipy.tex
@@ -1,134 +0,0 @@
-% This is the source code for my presentation on pyhht at SciPy India 2011 held in Mumbai in Dec 2011
-
-\documentclass[xcolor=dvipsnames]{beamer}
-\usetheme{Warsaw}
-\usecolortheme[named=Blue]{structure}
-\usepackage{graphics}
-
-\title{pyhht: A Python Toolbox for the Hilbert-Huang Transform}
-\author{Jaidev Deshpande \\ deshpande.jaidev@gmail.com}
-\institute{VIIT, Pune}
-\date{December 5, 2011}
-
-
-\begin{document}
-
-\begin{frame}
-\titlepage
-\end{frame}
-
-
-
-\begin{frame}{Introduction - Why this talk?}
-\begin{itemize}
-\item Introducing HHT
-\item An additional view of nonlinear and nonstationary phenomena
-\item HHT is almost entirely algorithmic
-\item It needs more math
-\item I need an audience
-\end{itemize}
-\end{frame}
-
-%3
-\begin{frame}{Motivation for the HHT}
-\begin{itemize}
-\item What's the Hilbert Transform good for?
-\begin{enumerate}
-\item Analytical signals
-\item Instantanous frequency
-\end{enumerate}
-\item Meaningful representations of data
-\item A check against harmonics
-\item Feature extraction
-\end{itemize}
-\end{frame}
-
-%4
-\begin{frame}{The Hilbert-Huang Transform}
-\begin{itemize}
-\item HHT = Empirical Mode Decomposition + Hilbert Spectral Analysis
-\item EMD: breaks a wideband signal down into piecewise narrowband signals, called IMFs
-\item Intrinsic Mode Functions:
-\begin{enumerate}
-\item Well behaved Hilbert Transforms
-\item Piecewise stationary
-\item Almost purely AM/FM
-\end{enumerate}
-\end{itemize}
-\end{frame}
-
-%5
-\begin{frame}{Intrinsic Mode Functions}
-\begin{definition}
-A function is called an Intrinsic Mode Function if
-
-\begin{enumerate}
-\item the number of zero crossings and local extrema in the function differ at most by unity (takes care of localized oscillations)
-\item the local mean of the enevelopes described by the local maxima and the local minima is zero at all times (required for meaningful instantaneous frequencies)
-\end{enumerate}
-\end{definition}
-
-\end{frame}
-
-%6
-\begin{frame}{Empirical Mode Decomposition}
-\begin{enumerate}
-\item Find all local extrema in the signal
-\item Join the maxima and minima with separate cubic splines, creating an upper and a lower envelope
-\item Calculate the mean of the envelopes
-\item Subtract mean from original signal
-\item Repeat steps 1-4 until the result is an IMF
-\item Subtract this IMF from the original signal
-\item Repeat steps 1-6 till there are no more IMFs left in the signal
-\end{enumerate}
-\end{frame}
-%
-%%7
-\begin{frame}{Hilbert Spectral Analysis}
-Analytical signals from the IMFs:
-\begin{equation}
-z(t)=x(t)+i*y(t)
-\end{equation}
-Instantaneous phase can be calculated as:
-\begin{equation}
-\theta = \tan ^{-1} {\frac{y(t)}{x(t)}}
-\end{equation}
-Instantaneous frequency, therefore:
-\begin{equation}
-\omega _{i}=\frac{d \theta}{dt}
-\end{equation}
-The original array is the real part of:
-\begin{equation}
-X(t)=\sum _{j=1} ^{N} x_{j}(t) \exp (i\int \omega_{j} (t).dt)
-\end{equation}
-\end{frame}
-%
-%%8
-\begin{frame}{Comparison with Fourier and Wavelet}
-%\includegraphics[scale=0.3]{Screenshot}
-\end{frame}
-%
-%%9
-\begin{frame}{Heuristics for HHT}
-\begin{itemize}
-\item Stoppage criteria for sifting
-\item Screening tools for IMFs based on statistical moments, information theoretic measures, etc
-\item Alternatives to sifting like SVD, ICA
-\item Better interpolation schemes
-\end{itemize}
-\end{frame}
-
-\begin{frame}{What Would be Awesome}
-\begin{itemize}
-\item Speeding it up
-\item A quantitative analysis of the EMD as a feature extraction, data compression method
-\item Use of wavelets as a handle on the HHT
-\item Establishing it as a veritable generalization of the Fourier method
-\end{itemize}
-\end{frame}
-%
-\begin{frame}
-Thank You
-\end{frame}
-
-\end{document}
Something went wrong with that request. Please try again.