technique for EEG to correct the Ocular Artifact Mahesh Khadtare, - - PowerPoint PPT Presentation

Accelerated SWT based de-noising technique for EEG to correct the Ocular Artifact

Mahesh Khadtare, Pragati Dharmale

3/19/2015 2

Plan of Session

I. Introduction II. Overview III. Stationary Wavelet Transform (SWT) IV. Algorithmic & implementation of EEG denoise V. Results VI. Discussions

Introduction

• Every time we think, move, feel or remember something, our neurons are at
• work. That work is carried out by small electric signals that zip from neuron to

neuron as fast as 250 mph [source: Walker]

• Using data recorded from the brain, the BCI processes it, interprets the

intention of the user, and acts on it, Figure shows Awake state and Asleep state [source: Wikipedia].

• Eye movements and blinks cause a severe problem for EEG measurements

[Source: EEG Lead placements-EEGLab]

3/19/2015 3

Introduction

Inspiration

• Brain Computer Interface(BCI)

3/19/2015 4

Neurologists

Introduction The references used for research are Evenly spread

5

1880 1900 1920 1940 1960 1980 2000 2020 Hans Berger S.A. Hillyard and R. Galambos

• R. Verleger, T. Gasser, and J.…

J.C. Woestenburg, M.N.…

• R. J. Croft and R. J. Barry

T.P. Jung, C. Humphries

• A. Cichocki and S. Vorobyov
• T. P. Jung, S. Makeig, M.…
• M. Potter, N. Gadhok, and W.…

T-P. Jung, C. Humphries, M.… Jean-François Cardoso

• M. Rahalova, P. Sykacek, M.…

Jung T-P, Humphries C, Lee… Jung T-P, Makeig S, Lee T-…

• Y. Li, A. Cichocki, and S.…

Ruijiang Li, Weifeng Liu,… Guger C, Ramoser H and… Ramoser H, Muller-Gerking J…

• A. SchlÄogl, P. Anderer, S.J.…

Charles W Anderson, James… Zachary A. Keirn , Jorge I.… Stone, M.

• I. Daubechies and W. Sweldens
• O. A. Rosso, M. T. Martin, A.…

Nazareth P. Castellanos,… David A. Peterson, James N.… BasËar E, Demiralp T,… Bartnik EA, Blinowska KJ,… Akay M, Akay YM, Cheng P,… Ademoglu A, Micheli-… Puthusserypady S, Zhou Z Peter Driessen Payam Refaeilzadeh, Lei… R.Amod, R. Sumit , Z. Supriya

• I. Smita, K. Asmita, P. Shradhha

YEAR AUTHER NAME

YEAR

Most are recent years

3/19/2015

Overview

• Almost all existing approaches to ocular artifact (OA) detection and

removal use one or more electro oculogram (EOG) signals either directly or indirectly as a reference.

• Some feature of the output of the algorithm (e.g., signal-to-noise ratio,
• r SNR) is compared to the original artifact-free EEG.
• For real EEG, the artifact-free (“true”) EEG is not known, so the

performance of the algorithm on real data is usually reported ,often based on visual inspection of the resulting waveforms.

• So the challenge is removal of EOG artifact without or less EEG

data distortion

6 3/19/2015

Overview

• What is Wavelet

– A small wave

• Wavelet Transforms

– Convert a signal into a series of wavelets – Provide a way for analyzing waveforms, bounded in both frequency and duration – Allow signals to be stored more efficiently than by Fourier transform – Be able to better approximate real-world signals – Well-suited for approximating data with sharp discontinuities

7 3/19/2015

Overview Fourier Vs Wavelet Transform

Fourier Transform Wavelet Transform Frequency Information only, Time/space information is lost Joint Time and Frequency Information Single Basis Function Many Basis Functions Computational Cost High Low computational costs Analysis Structures:

• FS( periodic functions only)
• FT and DFT

Numerous Analysis structures:

 CWT,  DWT(2 Band and M Band),  WP  SWT,  Frame structure

3/19/2015 8

Standard DWT

• Classical DWT is not shift invariant: This means

that DWT of a translated version of a signal x is not the same as the DWT of the original signal.

• Shift-invariance is important in many applications

such as:

– Change Detection – Denoising – Pattern Recognition

• In DWT, the signal is convolved and decimated (the even

indices are kept.)

• The decimation can be carried out by choosing the odd

indices.

3/19/2015

Standard DWT (2 Stage)

3/19/2015

High Pass Low Pass Low pass High pass

Decimated Wavelet Transform-Analysis Stage

High Freq Details Low freq Approx Signal

2

2

SWT

• Apply high and low pass filters to the data

at each level

• Do not decimate
• Modify the filters at each level, by padding

them with zeroes

• Computationally more complex

3/19/2015

SWT (2 Stage)

3/19/2015

High Pass Low Pass Low pass High pass

High Freq Details Low freq Approx Signal

Different Implementations

• A Trous Algorithm: Upsample the filter

coefficients by inserting zeros

• Beylkin’s algorithm: Shift invariance, shifts by
• ne will yield the same result by any odd shift.

Similarly, shift by zeroAll even shifts.

– Shift by 1 and 0 and compute the DWT, repeat the same procedure at each stage – Not a unique inverse: Invert each transform and average the results

• Undecimated Algorithm: Apply the lowpass and

highpass filters without any decimation.

3/19/2015

Traditional parallel algorithm

14

EEG 1 EEG 2

EEG N

PE 1

EEG 1 EEG 2

EEG N

PE 2 PE N

Dependencies

• Scheme depends upon the PE’s Configuration
• PE’s Interconnectedness
• Message passing Protocols
• Data transmission bandwidth

Traditional distributed processing algorithm does not include any parallel processing algorithm for Signal processing

3/19/2015

Proposed parallel algorithm

15

EEG 1

PE 1 PE 2 PE N

EEG 1 EEG N

PE 1 PE 2 PE N

EEG N

3/19/2015

16

EEG 1

G-PE 1 G-PE 2 G-PE N

P-EEG 1 EEG N

G-PE 1 G-PE 2 G-PE N

P-EEG N

Proposed parallel algorithm

3/19/2015

Mathematical Model

Wavelet Coefficients

17

• Detail Coefficient:
• Approximate Coefficient:



 

n

n k h n x k ylow ] 2 [ ]. [ ] [

3/19/2015



 

n

n k g n x k yhigh ] 2 [ ]. [ ] [

Algorithm Flow

Signal domain basis functions Transformed domain, Coeffs

Transformation Analysis

Reconstruction

Signals Recon Signal

basis functions

Modified Coefficients

Information extraction

Wavelet Decomposition

• The maximum frequency of EEG data sample is 125 Hz.
• The original data can be decomposed in up to 10 detail levels (D1 –

D10) and a last approximation (Apr).

• The frequency limits of each scale are approx. calculated dividing by 2

the sampling rate. In the case of 250 Hz, these are (with a 5 scales decomposition):

• D1: 64 – 128 Hz;
• D2: 32 – 64 Hz
• D3: 16 – 32 Hz;
• D4: 8 – 16 Hz;
• D5: 4 – 8 Hz;
• A5: 0 – 4 Hz.

Note that D2, D3, D4, D5, A5 approximately correspond to the EEG frequency bands: Gamma,Beta,Alpha,Theta,Delta respectively.

19

GAMMA ALPHA BETA THETA DELTA

3/19/2015

General denoising procedure

• The involves three steps. The basic version of the

procedure follows the steps described below:

• Decompose: Select a wavelet, select a level N. Compute

the wavelet decomposition of the signal at level N.

• Zeroing or Thresholding detail coefficients: For each

level from 1 to N, either make detail coefficient zero or select a threshold value and apply to the detail coefficients.

• Reconstruct: Compute wavelet reconstruction using the
• riginal approximation coefficients of level N and the

modified detail coefficients of levels from 1 to N.

20 3/19/2015

Threshold method

• In wavelet denoising process, the threshold method is
• ne of the main methods.
• Threshold function has important relationship with the

continuity and accuracy of the reconstructed signals, and has a significant impact on wavelet denoising.

• At present, there are two choices of methods of which

are hard threshold and soft threshold.

• But the hard threshold method makes the wavelet

coefficients discontinuous in the threshold value position and leads to the oscillations of the reconstructed signals; while with the soft threshold method wavelet coefficient has improved continuity.

3/19/2015 21

Threshold method

• The key point of wavelet threshold denoising is the

selection of threshold and how to choose the threshold function when dealing with wavelet coefficient after decomposition.

• There are four commonly used methods to select

threshold of which are Donoho-Johnstone methods: 1.Fixed-form (default)

• 2. Heursure
• 3. Rigsure

4.Minimax

3/19/2015 22

EEG Data

• EEG Database –
• Physionet EEG database

– 325 Recording with 7 Channel with different task

• Baseline, Rotation, Letter composing, Counting,

Multiplication

• Duration: 8 min (Sampling Frequency (fs): 256Hz)
• Typically : 840MB data size per recording
• Realtime EEG database

– 60min Recording with 4 Channels – Eye blinking, all above activity

Source: http://sccn.ucsd.edu/~arno/fam2data/publicly_available_EEG_data.html

24

Signal Data

Sensor

SWT Forward Transform Thresholding SWT Inverse EEG Signal Processing

The rate of Execution Time*

60% 61% 30 %

*1 Thread, Conventional Corner Turn

Distribution Of Execution Time

3/19/2015

Recons- tructed Signal

Demo – Sample unit test

Algorithm Implementation

3/19/2015 25

Read All Channel EEG Data for Signal Processing Rearrange all channel Data ( Block formation and groups) Rearrange Image Data ( Block formation and Flipping ) DWT or SWT forward DWT or SWT Inverse

Denoise Thresholds

G P U C P U

Compute Thresholds

Code Snippet

CPU code

• Lifting Step Wavelet
• Wavelet levels and

direction CUDA Kernel

• Row Kernel –
• Column Kernel

3/19/2015 26

__global__ void columnsKernel( double *d_Dst, double *d_Src, double *d_Kernel, int RowSize, int ColSize, int kernelRadius )

__global__, void rowsKernel( double *d_Dst double *d_Src, double *d_Kernel, int RowSize, int ColSize, int kernelRadius )

lifting_step_ti(double[,,] x, double[] h, int dir, int dist) perform_wavelet_transf(double[, ,]x, int Jmin, int dir)

Test Results

3/19/2015 27

• 30000
• 20000
• 10000

10000 20000 30000 40000 1 1627 3253 4879 6505 8131 9757 11383 13009 14635 16261 17887 19513 21139 22765 24391 26017 27643 29269 30895 32521 34147 35773 37399 39025 40651 42277 43903 45529 47155 48781 50407 52033 53659 55285 56911 58537 60163 61789 63415 65041 66667 68293 69919 71545 73171 74797

Amp(microV) Time (s)

Input EEG

Input EEG

• 2000
• 1000

1000 2000 1 1960 3919 5878 7837 9796 11755 13714 15673 17632 19591 21550 23509 25468 27427 29386 31345 33304 35263 37222 39181 41140 43099 45058 47017 48976 50935 52894 54853 56812 58771 60730 62689 64648 66607 68566 70525 72484 74443

Amp(microV) Time(s)

Denoise O/p MATLAB/C# Ch1/ CUDA

Denoise O/p MATLAB

EEG Denoising Timing

CPU / GPU Core Details Timing(sec) Speedup Intel Core i3 – MATLAB 2 45.337 1x Intel Core i3 – C# 2 12.121 3x GeForce GT 525M 96 1.031 44x Data size

512 x 8

Results – EEG Denoise

3/19/2015

Signal to Noise Ratio (SNR)

) log( 10 Pnoise Psignal SNR 

3/19/2015 29

From SNR, We decided accuracy of the wavelet.

 It is defined as the ratio of signal power to the noise power

corrupting the signal. A ratio higher than 1:1 indicates more signal than noise.

Pnoise Psignal SNR 

SIGNAL POWER NOISE POWER.

 SNR in DB is calculated as-

) log( 20 Anoise Asignal SNR 

Comparison Five Wavelet Types

TASK ALL COEFF. Max SNR WAVELET TYPE CD1,CD2 Max SNR WAVELET TYPE BASELINE 119.6598 COIF1 13.2003 DMEY COUNTING 9.4811 DMEY 18.35 SYM6 LETTER COMPOSING 7.17 SYM6 14.29 SYM6 MULTIPLICATI ON 9.37 DMEY 12.76 DMEY ROTATION 7.8818 DB4 19.41 DMEY

30

3/19/2015 31

3/19/2015 32

References

• I. Daubechies, “The wavelet transform, time-frequency localization and signal

analysis”, IEEE Trans. Inform. Theory, vol. 36, pp. 961,1005, statiscal learning theory, Cambridge press MIT

• J. S. Sahambi, S. N. Tandon, and R. K. P. Bhatt, “Using wavelet transforms for ecg

characterization”, IEEE Eng. in Medicine and Biology, pp. 77, 83, Jan/Feb 1997.

• Mukherjee S., “Classifying microarray data using support vector machines,” in

Understanding And Using Microarray Analysis Techniques: A Practical Guide. Boston Kluwer Academic Publishers, 2003.