National Chiao Tung University - Computer Science and Engineering

Report
National Chiao Tung University
MVEMD vs. MDEMD +
Applications in EEG & Gait Analyses
John K. Zao
Computer Science Dept. & Brain Research Center
National Chiao Tung University, Taiwan
2013/08/29
2013/8/29
MEMD Improvement & Apps
1
National Chiao Tung University
Agenda
 EMD vs. MVEMD vs. MDEMD
 MVEMD with PCA
 Application in Gait & EEG Analysis
 On-line & Light-weight Enhancements
2013/8/29
MEMD Improvement & Apps
2
National Chiao Tung University
Empirical Mode Decomposition (EMD)

Proposed by Dr. Norden E. Huang (1998)

Useful for non-linear non-stationary signal analysis

Decompose signals into Intrinsic Mode Functions
(IMFs) using sifting processing

原始訊號
找到局部極大、極小值,
且利用立方雲線找上下包
絡線
、
算出包絡線均值
k=k+1
IMFs capture oscillations at different speeds
No
k=0
n=n+1
是否為IMF
分量?
Yes
是否為單
調函數?
No
Yes
3
Empirical Mode Decomposition
National Chiao Tung University
Methodology : Original Signal
Source: NCU Lecture Slides 4
Empirical Mode Decomposition
National Chiao Tung University
Methodology : Original & m1 Signal
Source: NCU Lecture Slides 5
Empirical Mode Decomposition:
National Chiao Tung University
Methodology : Original & h1 Signal
Source: NCU Lecture Slides 6
National Chiao Tung University
M(V)EMD vs. MDEMD
 Multivariate EMD (MVEMD)
 Multidimensional EMD (MDEMD)
 Treats data from each channel as
the coordinate of a time-varying vector
in a vector space
 ( )
2013/8/29
⋯  ( )
 Treats data from each channel as
the value of a time-varying scalar
over a parameter space
  ⋯ ( )
MEMD Improvement & Apps
7
National Chiao Tung University
Multivariate Empirical Mode Decomposition (MVEMD)
 Decompose the trajectory of a vector into rotations at different speeds
 Find the envelop of trajectory
 Find the “center” of envelop
 Obtain the rotating component
by removing the trajectory of
the center
Questions:
 How to find the envelop?
 How o find i s “cen er”?
2013/8/29
MEMD Improvement & Apps
Source: BEMD & MEMD paper
8
National Chiao Tung University
Sifting based on Omnidirectional Projection
 Find the envelop of the trajectory by identifying the extrema of its projection in “evenly
spread” directions
 Evenly spread direction vectors in n-dimensional space can be found by placing
evenly distributed points on n-sphere using quasi-Monte Carlo methods based on
Hammersley sequences. Beware of he “curse of dimensionali y”!
 Extrema of the projection of the trajectory can be found using two methods:
a)
Find the centroids of the extrema  more sensitive to sampling errors
b) Find the mid-points of projection coordinates  more robust against sampling errors
 Algorithm (b) corresponds to 1D shifting along each projection directions
 Projections in evenly spread directions are used to reduce estimation errors of local mean
since trajectory orientation is unknown.
 Is it really needed?!
2013/8/29
MEMD Improvement & Apps
9
National Chiao Tung University
Multidimensional Empirical Mode Decomposition (MDEMD)
 Decompose the profile of a scalar field into n-dimensional oscillations
 Identify extrema of the profile
 Problems created by saddle points, ridges and valleys
 Create n-dimensional spline surfaces over the extrema
 No simple way to construct n-dimensional spline surfaces
 Several methods for 2D spline fitting
 Radial Based Function
 Thin Plate Interpretation
 Delaunay Triangulation
 By Slicing
 Non-Uniform Rational B-Spline
National Chiao Tung University
MDEMD based on EEMD & Min-Scale Combination
0
50
100
150
200
250
0
50
100
150
200
250
0
50
100
150
200
0
250
50
100
150
200
250
20
c1 ( x1 , y)
80
…
140
f ( x1, y)
EEMD
c2 ( x1, y)
..
.
200
2D
Image
c J ( s1 , t )
10
50
110
190
10
50
110
10
190
50
110
200
150
c1 ( x2 , y)
g1 ( x, y)
…
100
50
0
f ( x, y)
f (x2 , y)
EEMD
200
..
.
g2 ( x, y)
cJ (x2 , y)
..
.
…
100
50
0
…
g J ( x, y)
…
..
.
150
…
..
.
c2 (x2 , y)
200
c1(xM , y)
150
…
100
f (xM , y)
EEMD
c2 (xM , y)
50
..
.
cJ (xM , y)
0
Final 2DDecompositions:
2D-IMF1
2D-IMF2
2D-IMFn
2DResidual
190
National Chiao Tung University
MVEMD with PCA Preprocessing


Signal Re-orientation according to its Principal components
Signal Whitening according to its eigenvalues

−/   
Where ,  are eigenvalues and eigenvectors of covariance matrix  

Purposes:

Eliminate the effects of signal orientation and uneven power distribution
[Ques.] Can we simplify MVEMD algori hm when i ’s applied o whi ened
principal components? I think so.
2013/8/29
MEMD Improvement & Apps
12
National Chiao Tung University
PCA + MVEMD
 Separate 6D signals to two sets of 3D signals to do PCA (3D PCA)
 Recombine two sets of 3D principal components to do MEMD (6D
MEMD) and get same numbers IIMFs
Ax
Linear
Acceleration
Ay
3D PCA
Az
PCA1
PCA1 IMFs
PCA2
PCA2 IMFs
PCA3
PCA3 IMFs
6D MEMD
Gx
Angular
Velocity
Gy
Gz
2013/8/29
3D PCA
PCA1
PCA1 IMFs
PCA2
PCA2 IMFs
PCA3
PCA3 IMFs
MEMD Improvement & Apps
13
National Chiao Tung University
Principal Component Analysis (PCA)
 After analyzing, we can get
 eigenvectors
X
PCA1
0.4
0.3
0.2
0.2
 eigenvalues
0.1
0
0
 Use orthogonal transformation
-0.2
-0.1
-0.4
-0.2
1
 Reduce signal space dimensions
2
3
4
5
6
1
2
3
Y
4
5
6
4
5
6
4
5
6
PCA2
0.4
0.3
0.2
0.2
0.1
0
0
-0.2
-0.1
-0.4
-0.2
1
2
3
4
5
6
1
2
3
Z
PCA3
0.4
0.3
0.2
0.2
0.1
0
0
-0.2
-0.1
-0.4
-0.2
1
2
3
原資料
2013/8/29
MEMD Improvement & Apps
4
5
6
1
2
3
分析後
14
National Chiao Tung University
3D PCA
 Linear accelerations and angular velocities must be separated
 Do the whitening processing
 The unit-variance property of the whitened principal components
enhances the ability of MEMD
0.2
0.2
0.1
0.1
0
0
-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4
-0.4
-0.5
-0.5
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
18
20
(b)Principle Components
(a)Original signals
(a) is original signal,(b) is principal components
2013/8/29
MEMD Improvement & Apps
15
National Chiao Tung University
6D MVEMD
 Recombine two sets of 3D
0.02
0
-0.02
principal components
 Separate the each sets input
signals into a set of IMFs that
distinct frequency bands
 Each input signals will get
the same number of IMFs
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
0.05
0
-0.05
0.04
0.02
0
-0.02
-0.04
0.02
0
-0.02
0.02
0
-0.02
0
-3
x 10
2
0
-2
-4
0
-4
x 10
5
0
-5
0
-4
x 10
5
0
-5
-10
0
-4
x 10
4
2
0
-2
-4
0
-4
x 10
10
5
0
-5
0
2013/8/29
MEMD Improvement & Apps
16
National Chiao Tung University
Selection of PCA IMFs
IMF1
IMF2
IMF3
IMF4
IMF5
IMF6
IMF7
IMF8
IMF9
Residue
PCA1
0.3879
2.8707
2.6956
1.3987
2.0968
0.0340
0.0012
0.0031
0.0013
0.0069
PCA2
0.0717
0.2733
0.3455
0.7473
2.9635
0.2350
0.0469
0.3645
0.1729
0.6617
PCA3
0.1397
0.1252
0.0725
0.1417
0.0328
1.0051
0.0205
0.0302
0.0117
0.0944
3.5000
3.0000
2.5000
2.0000
PCA1
1.5000
PCA2
PCA3
1.0000
0.5000
0.0000
IMF1
2013/8/29
IMF2
IMF3
IMF4
IMF5
IMF6
IMF7
MEMD Improvement & Apps
IMF8
IMF9
IMF10
17
National Chiao Tung University
Construction of Characteristic Waveforms
 Derived from PCA IMFs of linear accelerations
 Gait cycle IMFs are selected first
 Remove gait cycles and trend IMFs
 Do the Gaussian distribution curve
fitting
 Impact IMFs are constructed from
IMFs fall into the main lobe of Gaussian distribution
2013/8/29
MEMD Improvement & Apps
18
National Chiao Tung University
Gaiting Characteristic Waveforms
Original sampled waveforms of 3D Linear Accelerations
Sampling rate: 50 samples/second
Waveforms of Dominant IMFs and “Shock Waves” extracted using PCA + MVEMD
es”
2013/8/29
MEMD Improvement & Apps
19
National Chiao Tung University
Feature Extractions
 Amplitude Modulation components
- signal’s time-varying amplitude
 Frequency Modulation components
- signal’s time-varying frequency
0.4
0.3
0.2
0.1
 Peak points
- when cause the stepping impacts
0
-0.1
 Phase Offset
-0.2
- whether the 3 axes are phase-locked
-0.3
 Trend
-0.4
5
5.5
6
6.5
7
7.5
- the changing direction of whole signal
2013/8/29
MEMD Improvement & Apps
20
National Chiao Tung University
Amplitude Modulation Components (AM)
 Find local extrema
 Perform cubic-spine interpolation through extrema
 Change of amplitudes reflects changes of step sizes
AM
AM
0.15
0.05
PCA1
PCA2
PCA3
0.1
0.04
0.05
0.03
0
0.02
-0.05
0.01
-0.1
0
0
10
15
20
25
MEMD Improvement & Apps
2013/8/29
5
5
5.5
6
6.5
7
7.5
8
21
8.5
9
9.5
10
National Chiao Tung University
Frequency Modulation components (FM)
 Calculate instantaneous frequency using Generalized Zero Crossing (GZC)
Observation
 Changes of frequency
reflect changes in gaiting speed
FM
5.5
PCA1
PCA2
PCA3
5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0
5
10
15
20
25
Time (s)
2013/8/29
MEMD Improvement & Apps
22
National Chiao Tung University
Phase Offset
 Deduced from time offsets between IMF zero-crossing points
Phase
2
PCA1 & PCA2
PCA1 & PCA3
PCA2 & PCA3
1.5
1
0.5
0
-0.5
-1
0
5
10
15
20
25
Time (s)
2013/8/29
MEMD Improvement & Apps
23
National Chiao Tung University
Impact Points
 Calculate instantaneous periods and use them as sliding windows
 Find the local maxima within the sliding windows
Observation
 Every impact point indicates
an impact of the feet with
the ground
0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0
2013/8/29
2
4
MEMD Improvement & Apps
6
8
10
12
14
16
18
20
24
National Chiao Tung University
Trend
 The last IMF corresponds to the trend signal
 Plot the trend signals into 3D space
Observation
 The trend of 3D linear acceleration corresponds to the general motion
directions of the human subject
-3
x 10
5
4
3
2
1
0
-1
-2
-3
-4
-5
-5
-4
-3
5
-2
4
-1
3
0
2
1
1
-3
x 10
0
2
-1
3
-2
-3
4
5
-4
-5
-3
x 10
2013/8/29
MEMD Improvement & Apps
25
National Chiao Tung University
SSVEP Stimulation
Color
RED
Frequency Luminanc
(Hz)
e (cd/m2)
32
Duty
Cycle (%)
153
20
50 sec recording
5~15 sec
Segment
(f10)
35~45 sec
Segment
(s10)
※MEEMD & MVEMD Analyses with 2 10-sec segments
National Chiao Tung University
Signal Processing
SSVEP Signal
Band Pass Filtering
1Hz ~ 100Hz
Down Sampling
MEEMD
Analysis
Select 6 Channels
(Fz, Fcz, Cz, Pz, Poz, Oz)
MVEMD Analysis
1000Hz → 500Hz
Noisy Channel & Epoch Removal
Select 6 Components
ICA
PCA
Stop condition: 1E-8
Channel Signal
Reconstruction
Bad ICA Component
Removal
Select 6 Good ICA Components
MVEMD Analysis
Channel Signal
Reconstruction
National Chiao Tung University
PCA Component Retrieval

EEGLAB function “runpca”

[pc,eigvec,sv] = runpca(EEG.data)

Select first 6 components from ‘pc’
P
C
A
_≒
3
2
H
z
National Chiao Tung University
0.01
f10中
20.8~22.8秒
0
-0.01
20.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
0.01
約32Hz波型圖
由上到下為
PCA1,
PCA2,
PCA3,
PCA4,
PCA5,
PCA6
0
-0.01
20.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
P
C
A
_≒
6
H
z
National Chiao Tung University
0.01
0
約16Hz波型圖
由上到下為
PCA1,
PCA2,
PCA3,
PCA4,
PCA5,
PCA6
-0.01
20.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
National Chiao Tung University
PCA__>64Hz
0.01
0
>64Hz波型圖
由上到下為
PCA1,
PCA2,
PCA3,
PCA4,
PCA5,
PCA6
-0.01
20.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
0.01
0
-0.01
20.8
National Chiao Tung University
PCA__Residue
0.02
0
Residue波型圖
由上到下為
PCA1,
PCA2,
PCA3,
PCA4,
PCA5,
PCA6
-0.02
20.8
20.9
21
21.1
21.2
21.3
21.4
21.5
21.6
21.7
21.8
20.9
21
21.1
21.2
21.3
21.4
21.5
21.6
21.7
21.8
20.9
21
21.1
21.2
21.3
21.4
21.5
21.6
21.7
21.8
20.9
21
21.1
21.2
21.3
21.4
21.5
21.6
21.7
21.8
20.9
21
21.1
21.2
21.3
21.4
21.5
21.6
21.7
21.8
20.9
21
21.1
21.2
21.3
21.4
21.5
21.6
21.7
21.8
0.02
0
-0.02
20.8
0.02
0
-0.02
20.8
0.02
0
-0.02
20.8
0.02
0
-0.02
20.8
0.02
0
-0.02
20.8
P
C
A
_≒
3
2
H
z
National Chiao Tung University
f10中
20.8~22.8
秒
約32Hz等高線圖
由上到下為
PCA1,
PCA2,
PCA3,
PCA4,
PCA5,
PCA6
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
P
C
A
_≒
6
H
z
National Chiao Tung University
約16Hz等高線圖
由上到下為
PCA1,
PCA2,
PCA3,
PCA4,
PCA5,
PCA6
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
National Chiao Tung University
PCA__>64Hz
>64Hz等高線圖
由上到下為
PCA1,
PCA2,
PCA3,
PCA4,
PCA5,
PCA6
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
National Chiao Tung University
PCA__Residue
Residue等高線圖
由上到下為
PCA1,
PCA2,
PCA3,
PCA4,
PCA5,
PCA6
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
National Chiao Tung University
Channel Signal Reconstruction

ICA and Bad Component Removal

EEGLAB -> Edit -> Select Data -> Data Range (Fz, FCz, Cz, Pz, POz, Oz)
National Chiao Tung University
LWH _ 32R – Fz、FCz、Cz、Pz、POz、Oz
Fz
Fz
Fz
Fz
FCz
FCz
FCz
FCz
Cz
Cz
Cz
Cz
Pz
Pz
Pz
Pz
POz
POz
POz
POz
Oz
Oz
Oz
Oz
0
100
200
300
原DATA
400
500
0
100
200
300
400
500
0
100
200
300
400
500
Fz
Fz
Fz
Fz
FCz
FCz
FCz
FCz
Cz
Cz
Cz
Cz
Pz
Pz
Pz
Pz
POz
POz
POz
POz
Oz
Oz
Oz
Oz
0
100
200
C1
300
400
500
0
100
200
300
400
500
0
100
200
300
400
500
Fz
Fz
Fz
Fz
FCz
FCz
FCz
FCz
Cz
Cz
Cz
Cz
Pz
Pz
Pz
Pz
POz
POz
POz
POz
Oz
Oz
Oz
Oz
0
100
200
300
400
500
0
100
200
300
400
500
0
100
200
300
400
500
0
100
200
300
400
500
0
100
200
300
400
500
0
100
200
300
400
500
0
100
200
300
400
500
C2
Fz
Fz
Fz
Fz
FCz
FCz
FCz
FCz
Cz
Cz
Cz
Cz
Pz
Pz
Pz
Pz
POz
POz
POz
POz
Oz
Oz
Oz
Oz
0
100
200
300
C3
400
500
0
100
200
300
400
500
0
100
200
300
400
500
C
a
e
l_M
E
D
≒
3
2
H
z
National Chiao Tung University
10
f10中
20.8~22.8
秒
0
-10
20.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
10
約32Hz波型圖
0
由上到下為
Fz,
-10
20.8
FCz,
10
Cz,
Pz,
0
POz,
-10
20.8
Oz
10
0
-10
20.8
10
0
-10
20.8
10
0
-10
20.8
C
a
e
l_M
E
D
≒
6
H
z
National Chiao Tung University
10
0
約16Hz波型圖
由上到下為
Fz,
FCz,
Cz,
Pz,
POz,
Oz
-10
20.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
10
0
-10
20.8
10
0
-10
20.8
10
0
-10
20.8
10
0
-10
20.8
10
0
-10
20.8
National Chiao Tung University
Channel__MEEMD__>64Hz
10
0
>64Hz波型圖
由上到下為
Fz,
FCz,
Cz,
Pz,
POz,
Oz
-10
20.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
10
0
-10
20.8
10
0
-10
20.8
10
0
-10
20.8
10
0
-10
20.8
10
0
-10
20.8
National Chiao Tung University
Channel__MEEMD__Residue
20
0
-20
Residue波型圖
由上到下為
Fz,
FCz,
Cz,
Pz,
POz,
Oz
20.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
20
0
-20
20.8
20
0
-20
20.8
20
0
-20
20.8
20
0
-20
20.8
20
0
-20
20.8
C
a
e
l_M
E
D
≒
3
2
H
z
National Chiao Tung University
f10中
20.8~22.8
秒
約32Hz等高線圖
由上到下為
Fz,
FCz,
Cz,
Pz,
POz,
Oz
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
C
ae
l_M
E
D
≒
6
H
z
National Chiao Tung University
約16Hz等高線圖
由上到下為
Fz,
FCz,
Cz,
Pz,
POz,
Oz
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
National Chiao Tung University
Channel__MEEMD__>64Hz
>64Hz等高線圖
由上到下為
Fz,
FCz,
Cz,
Pz,
POz,
Oz
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8
National Chiao Tung University
Channel__MEEMD__Residue
Residue等高線圖
由上到下為
Fz,
FCz,
Cz,
Pz,
POz,
Oz
21
21.2
21.4
21.6
21.8
22
22.2
22.4
22.6
22.8

similar documents