Skip to content

Instantly share code, notes, and snippets.

@rushlinneu

rushlinneu/pyeeg

Created June 26, 2016 22:30
Show Gist options
  • Save rushlinneu/1e2c8eb559dfc219bbd992ca905d0695 to your computer and use it in GitHub Desktop.
Save rushlinneu/1e2c8eb559dfc219bbd992ca905d0695 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
pyeeg
"""Copyleft 2010 Forrest Sheng Bao http://fsbao.net
PyEEG, a Python module to extract EEG features, v 0.02_r2
Project homepage: http://pyeeg.org
**Data structure**
PyEEG only uses standard Python and numpy data structures,
so you need to import numpy before using it.
For numpy, please visit http://numpy.scipy.org
**Naming convention**
I follow "Style Guide for Python Code" to code my program
http://www.python.org/dev/peps/pep-0008/
Constants: UPPER_CASE_WITH_UNDERSCORES, e.g., SAMPLING_RATE, LENGTH_SIGNAL.
Function names: lower_case_with_underscores, e.g., spectrum_entropy.
Variables (global and local): CapitalizedWords or CapWords, e.g., Power.
If a variable name consists of one letter, I may use lower case, e.g., x, y.
Functions listed alphabetically
--------------------------------------------------
"""
from numpy.fft import fft
from numpy import zeros, floor, log10, log, mean, array, sqrt, vstack, cumsum, \
ones, log2, std
from numpy.linalg import svd, lstsq
import time
######################## Functions contributed by Xin Liu #################
def hurst(X):
""" Compute the Hurst exponent of X. If the output H=0.5,the behavior
of the time-series is similar to random walk. If H<0.5, the time-series
cover less "distance" than a random walk, vice verse.
Parameters
----------
X
list
a time series
Returns
-------
H
float
Hurst exponent
Examples
--------
>>> import pyeeg
>>> from numpy.random import randn
>>> a = randn(4096)
>>> pyeeg.hurst(a)
>>> 0.5057444
"""
N = len(X)
T = array([float(i) for i in xrange(1,N+1)])
Y = cumsum(X)
Ave_T = Y/T
S_T = zeros((N))
R_T = zeros((N))
for i in xrange(N):
S_T[i] = std(X[:i+1])
X_T = Y - T * Ave_T[i]
R_T[i] = max(X_T[:i + 1]) - min(X_T[:i + 1])
R_S = R_T / S_T
R_S = log(R_S)
n = log(T).reshape(N, 1)
H = lstsq(n[1:], R_S[1:])[0]
return H[0]
######################## Begin function definitions #######################
def embed_seq(X,Tau,D):
"""Build a set of embedding sequences from given time series X with lag Tau
and embedding dimension DE. Let X = [x(1), x(2), ... , x(N)], then for each
i such that 1 < i < N - (D - 1) * Tau, we build an embedding sequence,
Y(i) = [x(i), x(i + Tau), ... , x(i + (D - 1) * Tau)]. All embedding
sequence are placed in a matrix Y.
Parameters
----------
X
list
a time series
Tau
integer
the lag or delay when building embedding sequence
D
integer
the embedding dimension
Returns
-------
Y
2-D list
embedding matrix built
Examples
---------------
>>> import pyeeg
>>> a=range(0,9)
>>> pyeeg.embed_seq(a,1,4)
array([[ 0., 1., 2., 3.],
[ 1., 2., 3., 4.],
[ 2., 3., 4., 5.],
[ 3., 4., 5., 6.],
[ 4., 5., 6., 7.],
[ 5., 6., 7., 8.]])
>>> pyeeg.embed_seq(a,2,3)
array([[ 0., 2., 4.],
[ 1., 3., 5.],
[ 2., 4., 6.],
[ 3., 5., 7.],
[ 4., 6., 8.]])
>>> pyeeg.embed_seq(a,4,1)
array([[ 0.],
[ 1.],
[ 2.],
[ 3.],
[ 4.],
[ 5.],
[ 6.],
[ 7.],
[ 8.]])
"""
N =len(X)
if D * Tau > N:
print "Cannot build such a matrix, because D * Tau > N"
exit()
if Tau<1:
print "Tau has to be at least 1"
exit()
Y=zeros((N - (D - 1) * Tau, D))
for i in xrange(0, N - (D - 1) * Tau):
for j in xrange(0, D):
Y[i][j] = X[i + j * Tau]
return Y
def in_range(Template, Scroll, Distance):
"""Determines whether one vector is the the range of another vector.
The two vectors should have equal length.
Parameters
-----------------
Template
list
The template vector, one of two vectors being compared
Scroll
list
The scroll vector, one of the two vectors being compared
D
float
Two vectors match if their distance is less than D
Bit
Notes
-------
The distance between two vectors can be defined as Euclidean distance
according to some publications.
The two vector should of equal length
"""
for i in range(0, len(Template)):
if abs(Template[i] - Scroll[i]) > Distance:
return False
return True
""" Desperate code, but do not delete
def bit_in_range(Index):
if abs(Scroll[Index] - Template[Bit]) <= Distance :
print "Bit=", Bit, "Scroll[Index]", Scroll[Index], "Template[Bit]",\
Template[Bit], "abs(Scroll[Index] - Template[Bit])",\
abs(Scroll[Index] - Template[Bit])
return Index + 1 # move
Match_No_Tail = range(0, len(Scroll) - 1) # except the last one
# print Match_No_Tail
# first compare Template[:-2] and Scroll[:-2]
for Bit in xrange(0, len(Template) - 1): # every bit of Template is in range of Scroll
Match_No_Tail = filter(bit_in_range, Match_No_Tail)
print Match_No_Tail
# second and last, check whether Template[-1] is in range of Scroll and
# Scroll[-1] in range of Template
# 2.1 Check whether Template[-1] is in the range of Scroll
Bit = - 1
Match_All = filter(bit_in_range, Match_No_Tail)
# 2.2 Check whether Scroll[-1] is in the range of Template
# I just write a loop for this.
for i in Match_All:
if abs(Scroll[-1] - Template[i] ) <= Distance:
Match_All.remove(i)
return len(Match_All), len(Match_No_Tail)
"""
def bin_power(X,Band,Fs):
"""Compute power in each frequency bin specified by Band from FFT result of
X. By default, X is a real signal.
Note
-----
A real signal can be synthesized, thus not real.
Parameters
-----------
Band
list
boundary frequencies (in Hz) of bins. They can be unequal bins, e.g.
[0.5,4,7,12,30] which are delta, theta, alpha and beta respectively.
You can also use range() function of Python to generate equal bins and
pass the generated list to this function.
Each element of Band is a physical frequency and shall not exceed the
Nyquist frequency, i.e., half of sampling frequency.
X
list
a 1-D real time series.
Fs
integer
the sampling rate in physical frequency
Returns
-------
Power
list
spectral power in each frequency bin.
Power_ratio
list
spectral power in each frequency bin normalized by total power in ALL
frequency bins.
"""
C = fft(X)
C = abs(C)
Power =zeros(len(Band)-1);
for Freq_Index in xrange(0,len(Band)-1):
Freq = float(Band[Freq_Index]) ## Xin Liu
Next_Freq = float(Band[Freq_Index+1])
Power[Freq_Index] = sum(C[floor(Freq/Fs*len(X)):floor(Next_Freq/Fs*len(X))])
Power_Ratio = Power/sum(Power)
return Power, Power_Ratio
def first_order_diff(X):
""" Compute the first order difference of a time series.
For a time series X = [x(1), x(2), ... , x(N)], its first order
difference is:
Y = [x(2) - x(1) , x(3) - x(2), ..., x(N) - x(N-1)]
"""
D=[]
for i in xrange(1,len(X)):
D.append(X[i]-X[i-1])
return D
def pfd(X, D=None):
"""Compute Petrosian Fractal Dimension of a time series from either two
cases below:
1. X, the time series of type list (default)
2. D, the first order differential sequence of X (if D is provided,
recommended to speed up)
In case 1, D is computed by first_order_diff(X) function of pyeeg
To speed up, it is recommended to compute D before calling this function
because D may also be used by other functions whereas computing it here
again will slow down.
"""
if D is None: ## Xin Liu
D = first_order_diff(X)
N_delta= 0; #number of sign changes in derivative of the signal
for i in xrange(1,len(D)):
if D[i]*D[i-1]<0:
N_delta += 1
n = len(X)
return log10(n)/(log10(n)+log10(n/n+0.4*N_delta))
def hfd(X, Kmax):
""" Compute Hjorth Fractal Dimension of a time series X, kmax
is an HFD parameter
"""
L = [];
x = []
N = len(X)
for k in xrange(1,Kmax):
Lk = []
for m in xrange(0,k):
Lmk = 0
for i in xrange(1,int(floor((N-m)/k))):
Lmk += abs(X[m+i*k] - X[m+i*k-k])
Lmk = Lmk*(N - 1)/floor((N - m) / float(k)) / k
Lk.append(Lmk)
L.append(log(mean(Lk)))
x.append([log(float(1) / k), 1])
(p, r1, r2, s)=lstsq(x, L)
return p[0]
def hjorth(X, D = None):
""" Compute Hjorth mobility and complexity of a time series from either two
cases below:
1. X, the time series of type list (default)
2. D, a first order differential sequence of X (if D is provided,
recommended to speed up)
In case 1, D is computed by first_order_diff(X) function of pyeeg
Notes
-----
To speed up, it is recommended to compute D before calling this function
because D may also be used by other functions whereas computing it here
again will slow down.
Parameters
----------
X
list
a time series
D
list
first order differential sequence of a time series
Returns
-------
As indicated in return line
Hjorth mobility and complexity
"""
if D is None:
D = first_order_diff(X)
D.insert(0, X[0]) # pad the first difference
D = array(D)
n = len(X)
M2 = float(sum(D ** 2)) / n
TP = sum(array(X) ** 2)
M4 = 0;
for i in xrange(1, len(D)):
M4 += (D[i] - D[i - 1]) ** 2
M4 = M4 / n
return sqrt(M2 / TP), sqrt(float(M4) * TP / M2 / M2) # Hjorth Mobility and Complexity
def spectral_entropy(X, Band, Fs, Power_Ratio = None):
"""Compute spectral entropy of a time series from either two cases below:
1. X, the time series (default)
2. Power_Ratio, a list of normalized signal power in a set of frequency
bins defined in Band (if Power_Ratio is provided, recommended to speed up)
In case 1, Power_Ratio is computed by bin_power() function.
Notes
-----
To speed up, it is recommended to compute Power_Ratio before calling this
function because it may also be used by other functions whereas computing
it here again will slow down.
Parameters
----------
Band
list
boundary frequencies (in Hz) of bins. They can be unequal bins, e.g.
[0.5,4,7,12,30] which are delta, theta, alpha and beta respectively.
You can also use range() function of Python to generate equal bins and
pass the generated list to this function.
Each element of Band is a physical frequency and shall not exceed the
Nyquist frequency, i.e., half of sampling frequency.
X
list
a 1-D real time series.
Fs
integer
the sampling rate in physical frequency
Returns
-------
As indicated in return line
See Also
--------
bin_power: pyeeg function that computes spectral power in frequency bins
"""
if Power_Ratio is None:
Power, Power_Ratio = bin_power(X, Band, Fs)
Spectral_Entropy = 0
for i in xrange(0, len(Power_Ratio) - 1):
Spectral_Entropy += Power_Ratio[i] * log(Power_Ratio[i])
Spectral_Entropy /= log(len(Power_Ratio)) # to save time, minus one is omitted
return -1 * Spectral_Entropy
def svd_entropy(X, Tau, DE, W = None):
"""Compute SVD Entropy from either two cases below:
1. a time series X, with lag tau and embedding dimension dE (default)
2. a list, W, of normalized singular values of a matrix (if W is provided,
recommend to speed up.)
If W is None, the function will do as follows to prepare singular spectrum:
First, computer an embedding matrix from X, Tau and DE using pyeeg
function embed_seq():
M = embed_seq(X, Tau, DE)
Second, use scipy.linalg function svd to decompose the embedding matrix
M and obtain a list of singular values:
W = svd(M, compute_uv=0)
At last, normalize W:
W /= sum(W)
Notes
-------------
To speed up, it is recommended to compute W before calling this function
because W may also be used by other functions whereas computing it here
again will slow down.
"""
if W is None:
Y = EmbedSeq(X, tau, dE)
W = svd(Y, compute_uv = 0)
W /= sum(W) # normalize singular values
return -1*sum(W * log(W))
def fisher_info(X, Tau, DE, W = None):
""" Compute Fisher information of a time series from either two cases below:
1. X, a time series, with lag Tau and embedding dimension DE (default)
2. W, a list of normalized singular values, i.e., singular spectrum (if W is
provided, recommended to speed up.)
If W is None, the function will do as follows to prepare singular spectrum:
First, computer an embedding matrix from X, Tau and DE using pyeeg
function embed_seq():
M = embed_seq(X, Tau, DE)
Second, use scipy.linalg function svd to decompose the embedding matrix
M and obtain a list of singular values:
W = svd(M, compute_uv=0)
At last, normalize W:
W /= sum(W)
Parameters
----------
X
list
a time series. X will be used to build embedding matrix and compute
singular values if W or M is not provided.
Tau
integer
the lag or delay when building a embedding sequence. Tau will be used
to build embedding matrix and compute singular values if W or M is not
provided.
DE
integer
the embedding dimension to build an embedding matrix from a given
series. DE will be used to build embedding matrix and compute
singular values if W or M is not provided.
W
list or array
the set of singular values, i.e., the singular spectrum
Returns
-------
FI
integer
Fisher information
Notes
-----
To speed up, it is recommended to compute W before calling this function
because W may also be used by other functions whereas computing it here
again will slow down.
See Also
--------
embed_seq : embed a time series into a matrix
"""
if W is None:
M = embed_seq(X, Tau, DE)
W = svd(M, compute_uv = 0)
W /= sum(W)
FI = 0
for i in xrange(0, len(W) - 1): # from 1 to M
FI += ((W[i +1] - W[i]) ** 2) / (W[i])
return FI
def ap_entropy(X, M, R):
"""Computer approximate entropy (ApEN) of series X, specified by M and R.
Suppose given time series is X = [x(1), x(2), ... , x(N)]. We first build
embedding matrix Em, of dimension (N-M+1)-by-M, such that the i-th row of Em
is x(i),x(i+1), ... , x(i+M-1). Hence, the embedding lag and dimension are
1 and M-1 respectively. Such a matrix can be built by calling pyeeg function
as Em = embed_seq(X, 1, M). Then we build matrix Emp, whose only
difference with Em is that the length of each embedding sequence is M + 1
Denote the i-th and j-th row of Em as Em[i] and Em[j]. Their k-th elments
are Em[i][k] and Em[j][k] respectively. The distance between Em[i] and Em[j]
is defined as 1) the maximum difference of their corresponding scalar
components, thus, max(Em[i]-Em[j]), or 2) Euclidean distance. We say two 1-D
vectors Em[i] and Em[j] *match* in *tolerance* R, if the distance between them
is no greater than R, thus, max(Em[i]-Em[j]) <= R. Mostly, the value of R is
defined as 20% - 30% of standard deviation of X.
Pick Em[i] as a template, for all j such that 0 < j < N - M + 1, we can
check whether Em[j] matches with Em[i]. Denote the number of Em[j],
which is in the range of Em[i], as k[i], which is the i-th element of the
vector k. The probability that a random row in Em matches Em[i] is
\simga_1^{N-M+1} k[i] / (N - M + 1), thus sum(k)/ (N - M + 1),
denoted as Cm[i].
We repeat the same process on Emp and obtained Cmp[i], but here 0<i<N-M
since the length of each sequence in Emp is M + 1.
The probability that any two embedding sequences in Em match is then
sum(Cm)/ (N - M +1 ). We define Phi_m = sum(log(Cm)) / (N - M + 1) and
Phi_mp = sum(log(Cmp)) / (N - M ).
And the ApEn is defined as Phi_m - Phi_mp.
Notes
-----
#. Please be aware that self-match is also counted in ApEn.
#. This function now runs very slow. We are still trying to speed it up.
References
----------
Costa M, Goldberger AL, Peng CK, Multiscale entropy analysis of biolgical
signals, Physical Review E, 71:021906, 2005
See also
--------
samp_entropy: sample entropy of a time series
Notes
-----
Extremely slow implementation. Do NOT use if your dataset is not small.
"""
N = len(X)
Em = embed_seq(X, 1, M)
Emp = embed_seq(X, 1, M + 1) # try to only build Emp to save time
Cm, Cmp = zeros(N - M + 1), zeros(N - M)
# in case there is 0 after counting. Log(0) is undefined.
for i in xrange(0, N - M):
# print i
for j in xrange(i, N - M): # start from i, self-match counts in ApEn
# if max(abs(Em[i]-Em[j])) <= R:# compare N-M scalars in each subseq v 0.01b_r1
if in_range(Em[i], Em[j], R):
Cm[i] += 1 ### Xin Liu
Cm[j] += 1
if abs(Emp[i][-1] - Emp[j][-1]) <= R: # check last one
Cmp[i] += 1
Cmp[j] += 1
if in_range(Em[i], Em[N-M], R):
Cm[i] += 1
Cm[N-M] += 1
# try to count Cm[j] and Cmp[j] as well here
# if max(abs(Em[N-M]-Em[N-M])) <= R: # index from 0, so N-M+1 is N-M v 0.01b_r1
# if in_range(Em[i], Em[N - M], R): # for Cm, there is one more iteration than Cmp
# Cm[N - M] += 1 # cross-matches on Cm[N - M]
Cm[N - M] += 1 # Cm[N - M] self-matches
# import code;code.interact(local=locals())
Cm /= (N - M +1 )
Cmp /= ( N - M )
# import code;code.interact(local=locals())
Phi_m, Phi_mp = sum(log(Cm)), sum(log(Cmp))
Ap_En = (Phi_m - Phi_mp) / (N - M)
return Ap_En
def samp_entropy(X, M, R):
"""Computer sample entropy (SampEn) of series X, specified by M and R.
SampEn is very close to ApEn.
Suppose given time series is X = [x(1), x(2), ... , x(N)]. We first build
embedding matrix Em, of dimension (N-M+1)-by-M, such that the i-th row of Em
is x(i),x(i+1), ... , x(i+M-1). Hence, the embedding lag and dimension are
1 and M-1 respectively. Such a matrix can be built by calling pyeeg function
as Em = embed_seq(X, 1, M). Then we build matrix Emp, whose only
difference with Em is that the length of each embedding sequence is M + 1
Denote the i-th and j-th row of Em as Em[i] and Em[j]. Their k-th elments
are Em[i][k] and Em[j][k] respectively. The distance between Em[i] and Em[j]
is defined as 1) the maximum difference of their corresponding scalar
components, thus, max(Em[i]-Em[j]), or 2) Euclidean distance. We say two 1-D
vectors Em[i] and Em[j] *match* in *tolerance* R, if the distance between them
is no greater than R, thus, max(Em[i]-Em[j]) <= R. Mostly, the value of R is
defined as 20% - 30% of standard deviation of X.
Pick Em[i] as a template, for all j such that 0 < j < N - M , we can
check whether Em[j] matches with Em[i]. Denote the number of Em[j],
which is in the range of Em[i], as k[i], which is the i-th element of the
vector k.
We repeat the same process on Emp and obtained Cmp[i], 0 < i < N - M.
The SampEn is defined as log(sum(Cm)/sum(Cmp))
References
----------
Costa M, Goldberger AL, Peng C-K, Multiscale entropy analysis of biolgical
signals, Physical Review E, 71:021906, 2005
See also
--------
ap_entropy: approximate entropy of a time series
Notes
-----
Extremely slow computation. Do NOT use if your dataset is not small and you
are not patient enough.
"""
N = len(X)
Em = embed_seq(X, 1, M)
Emp = embed_seq(X, 1, M + 1)
Cm, Cmp = zeros(N - M - 1) + 1e-100, zeros(N - M - 1) + 1e-100
# in case there is 0 after counting. Log(0) is undefined.
for i in xrange(0, N - M):
for j in xrange(i + 1, N - M): # no self-match
# if max(abs(Em[i]-Em[j])) <= R: # v 0.01_b_r1
if in_range(Em[i], Em[j], R):
Cm[i] += 1
# if max(abs(Emp[i] - Emp[j])) <= R: # v 0.01_b_r1
if abs(Emp[i][-1] - Emp[j][-1]) <= R: # check last one
Cmp[i] += 1
Samp_En = log(sum(Cm)/sum(Cmp))
return Samp_En
def dfa(X, Ave = None, L = None):
"""Compute Detrended Fluctuation Analysis from a time series X and length of
boxes L.
The first step to compute DFA is to integrate the signal. Let original seres
be X= [x(1), x(2), ..., x(N)].
The integrated signal Y = [y(1), y(2), ..., y(N)] is otained as follows
y(k) = \sum_{i=1}^{k}{x(i)-Ave} where Ave is the mean of X.
The second step is to partition/slice/segment the integrated sequence Y into
boxes. At least two boxes are needed for computing DFA. Box sizes are
specified by the L argument of this function. By default, it is from 1/5 of
signal length to one (x-5)-th of the signal length, where x is the nearest
power of 2 from the length of the signal, i.e., 1/16, 1/32, 1/64, 1/128, ...
In each box, a linear least square fitting is employed on data in the box.
Denote the series on fitted line as Yn. Its k-th elements, yn(k),
corresponds to y(k).
For fitting in each box, there is a residue, the sum of squares of all
offsets, difference between actual points and points on fitted line.
F(n) denotes the square root of average total residue in all boxes when box
length is n, thus
Total_Residue = \sum_{k=1}^{N}{(y(k)-yn(k))}
F(n) = \sqrt(Total_Residue/N)
The computing to F(n) is carried out for every box length n. Therefore, a
relationship between n and F(n) can be obtained. In general, F(n) increases
when n increases.
Finally, the relationship between F(n) and n is analyzed. A least square
fitting is performed between log(F(n)) and log(n). The slope of the fitting
line is the DFA value, denoted as Alpha. To white noise, Alpha should be
0.5. Higher level of signal complexity is related to higher Alpha.
Parameters
----------
X:
1-D Python list or numpy array
a time series
Ave:
integer, optional
The average value of the time series
L:
1-D Python list of integers
A list of box size, integers in ascending order
Returns
-------
Alpha:
integer
the result of DFA analysis, thus the slope of fitting line of log(F(n))
vs. log(n). where n is the
Examples
--------
>>> import pyeeg
>>> from numpy.random import randn
>>> print pyeeg.dfa(randn(4096))
0.490035110345
Reference
---------
Peng C-K, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling
exponents and crossover phenomena in nonstationary heartbeat time series.
_Chaos_ 1995;5:82-87
Notes
-----
This value depends on the box sizes very much. When the input is a white
noise, this value should be 0.5. But, some choices on box sizes can lead to
the value lower or higher than 0.5, e.g. 0.38 or 0.58.
Based on many test, I set the box sizes from 1/5 of signal length to one
(x-5)-th of the signal length, where x is the nearest power of 2 from the
length of the signal, i.e., 1/16, 1/32, 1/64, 1/128, ...
You may generate a list of box sizes and pass in such a list as a parameter.
"""
X = array(X)
if Ave is None:
Ave = mean(X)
Y = cumsum(X)
Y -= Ave
if L is None:
L = floor(len(X)*1/(2**array(range(4,int(log2(len(X)))-4))))
F = zeros(len(L)) # F(n) of different given box length n
for i in xrange(0,len(L)):
n = int(L[i]) # for each box length L[i]
if n==0:
print "time series is too short while the box length is too big"
print "abort"
exit()
for j in xrange(0,len(X),n): # for each box
if j+n < len(X):
c = range(j,j+n)
c = vstack([c, ones(n)]).T # coordinates of time in the box
y = Y[j:j+n] # the value of data in the box
F[i] += lstsq(c,y)[1] # add residue in this box
F[i] /= ((len(X)/n)*n)
F = sqrt(F)
Alpha = lstsq(vstack([log(L), ones(len(L))]).T,log(F))[0][0]
return Alpha
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment