Localization of eye Saccadic signatures in Electrooculograms using sparse representations with data driven dictionaries

Accepted Manuscript

Localization of Eye Saccadic Signatures in Electrooculograms using Sparse Representations with Data driven Dictionaries Suvodip Chakraborty, Anirban Dasgupta, Aurobinda Routray PII: DOI: Reference:

S0167-8655(17)30407-5 10.1016/j.patrec.2017.11.001 PATREC 6993

To appear in:

Pattern Recognition Letters

Received date: Revised date: Accepted date:

30 November 2016 30 October 2017 1 November 2017

Please cite this article as: Suvodip Chakraborty, Anirban Dasgupta, Aurobinda Routray, Localization of Eye Saccadic Signatures in Electrooculograms using Sparse Representations with Data driven Dictionaries, Pattern Recognition Letters (2017), doi: 10.1016/j.patrec.2017.11.001

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT 1 Highlights • Two methods for localizing saccadic signatures in EOG signals have been proposed. • The first method uses sparse representations in nonanalytical dictionaries to determine the presence of saccadic signatures. • The second method uses DTW to obtain the same.

AC

CE

PT

ED

M

AN US

• The findings of this work are that the first method is computationally more efficient, while the second one is more accurate.

CR IP T

• The proposed methods have been compared with the existing method on the Physiosig Dataset.

ACCEPTED MANUSCRIPT 2

Pattern Recognition Letters journal homepage: www.elsevier.com

Suvodip Chakrabortya,∗∗, Anirban Dasguptab , Aurobinda Routrayb a Advanced

CR IP T

Localization of Eye Saccadic Signatures in Electrooculograms using Sparse Representations with Data driven Dictionaries

Technology Development Center, Indian Institute of Technology Kharagpur, Kharagpur, India - 721302. of Electrical Engineering, Indian Institute of Technology Kharagpur, Kharagpur, India - 721302.

b Department

ABSTRACT

M

AN US

In this paper, we propose two methods for localizing saccadic eye movement signatures from Electrooculograms (EOG). The first approach uses a sparse representation of data-driven dictionaries of saccadic movements. In this approach, we match the EOG subsequence with the dictionary element using distance metrics to identify the saccades. The second approach is to compare a saccadic signature template with the EOG subsequence using Dynamic Time Warping (DTW). We find that the proposed methods have advantages over one another in context specific applications. The first method is significantly faster with considerable accuracy, while the second approach is more accurate. c Elsevier B.V. All rights reserved Copyright c 2017 Elsevier Ltd. All rights reserved.

1. Introduction

CE

PT

ED

Biomedical signals contain certain repeating patterns which often convey information related to some specific physiological phenomenon [1]. Detection and localization of such patterns in biomedical signals have developed great interest among the researchers over the last few decades [2]. Analysis of repeating patterns has been carried out using adaptive signal representation techniques in earlier works [3], [4]. Such representations have provided appreciable results in the case of audio signals [5]. However, such methods show limitations for biomedical signals due to the following reasons:

AC

• low signal-to-noise ratios (SNR) in most biomedical signals • presence of additive noise which may include interference from other biomedical signals • the patterns may be stretched or shrunk in time Biomedical signal analysis is conventionally carried out using signal decomposition methods such as the Short Time Fourier

∗∗ Corresponding

author. Tel.: +91 9126573570 e-mail: [email protected] (Suvodip Chakraborty), [email protected] (Anirban Dasgupta), [email protected]. (Aurobinda Routray)

Transform and Wavelet Transform [6]. The limitations of such representations are that they have fixed or analytical dictionaries, and hence are not suitable for representing any generic pattern. Such cases are usually dealt by feature selection based methods [7], which may lead to additional computational overhead depending on the feature vector size. Pursuit-type decompositions [8] which use analytical, redundant dictionaries, are designed to yield sparse representations. Such decomposition works appreciably for target signals which are in the same class as of the dictionary elements [9]. The pursuit search, however, performs poorly in the presence of realistic noise and clutter and has a significant computational burden. We can see that automatic isolation of a particular pattern from any generic biomedical signal is a difficult task. This is because the patterns are dependent on the nature of the biomedical signal. In this paper, we propose a case study on identifying saccadic signatures on Electrooculograms (EOG). The proposed methods will work for any biomedical signal for any kind of pattern if the training pattern is selected accordingly. EOG signals signify the standing potential difference existing between the cornea and the retina [10]. The magnitude of these signals indicates how far the eyes move from a reference position. Hence, EOG is useful in analyzing different kinds of eye movements such as saccades, gaze fixation, vergence, smooth pursuits, blinks, etc. [11]. A saccade is a fast simultaneous movement of both the eyes in the same direction [12]. A smooth

ACCEPTED MANUSCRIPT 3

1.1. Prior Art

AN US

A few earlier attempts has been made to automate the localization of saccades in EOG. The work by Behrens et al.[14] isolates saccades using thresholding technique, where the derivative of the EOG amplitude signal is thresholded and zerocrossings are subsequently identified. Jung et al.[15] used extended Independent Component Analysis to isolate the saccades from EOG. Kumar et al. [16] have used Discrete Wavelet Transform (DWT) in order to detect saccades based on the DWT coefficients. These methods, however, fail to distinguish saccadic and blink signatures isolated from EOG in certain cases. In such a situation where exclusively saccades are required, these methods may not be too accurate.

1.4. Signature detection in signals In this section we try to define the generic problem of signature detection in a signal. Fig. 1 shows that three similar patterns occur in a stochastic and non-stationary signal. The signals of interest are highlighted in figure 1. Given a signal x(n) of length N and a non-analytic data driven pattern m(n) of length L where N >> L, the problem is to find how close a match m(n) is to a window, xP (n) of length P, which is a subsequence of the signal x(n). Here, P = L is the case when the sample template length matches with the template present in the signal. In case, the signal is stretched or shrunk as compared to the template, the similarity of the signal to the template is evaluated for P > L or P < L respectively. To estimate the similarity between the signal x(n) and the nonanalytic data driven pattern m(n), one may use a distance metric D(x, m) as a dissimilarity score. The distance function D(x, m) is a measure of the divergence between the pattern m(n) and the sliding signal window xP (n) of the signal. The discovery of saccades in the EOG time-series first requires the selection of a smoothed sample saccade referred to as primitive template. This is obtained by observing the most commonly occurring saccade. Let the EOG signal be modeled as x(n) = [x(1) x(2) . . . x(N)] while the saccadic signature be m(n) = [m(1) m(2) . . . m(L)]. xP (n) is to be matched with a template saccade to find the probable saccades in the whole EOG signal. The sub-sequence is given by

CR IP T

pursuit movement occurs when someone carefully observes a target stimulus as it moves. Fixational eye movements are small oscillations of the eyeball when the person extracts visual information from a particular scene [10]. This implies that a saccadic movement occurs during the transition between two gaze fixations [13]. Saccadic eye movements are of prime importance in the fields of medical diagnosis as well as understanding cognition and emotion and, therefore, has been a field of interest to researchers.

M

The existing methods perform poorly for automated saccade isolation owing to the following challenges: • EOG data contains noise

ED

• Movements of body parts introduce artifacts to EOG data • Blinks and saccades carry similar signature

PT

• Saccadic movements may be bidirectional, which cause reversal of a unidirectional pattern

CE

• Temporal variations in signature of saccades in EOG

xP (n) = [m(p) m(p + 1) . . . m(p + P − 1)]

(1)

Here, xP (n) starts at index p and extends up to p + P − 1. Signal Amplitude

1.2. Research Issues

10

0

-10 0

500

1000

1500

Time axis

Fig. 1. An example of repeating patterns in time-series

1.3. The Motivation

AC

Saccadic movements are of significant importance when it comes to observing human behavior through eye movements [17]. It plays a vital role in the field of Human-computer interaction [18], biometrics [19], rehabilitation [20], [21] emotion classification [22] etc. The medical applications of eye saccades are the diagnosis of Parkinson’s syndrome [23], spinocerebellar ataxias, Niemann-Pick disease [24] etc. The peak saccadic velocity, saccadic latency and saccadic duration are some of the saccadic parameters which are validated indicators of visual attention [25]. Hence, it is essential to accurately identify and localize the saccadic signatures in EOG. However, the existing methods for automatic saccade isolation encounter a lot of issues as mentioned in the previous subsection. These issues create a need for a new method to automatically classify saccadic signatures accurately.

1.5. Objectives The objective of this paper is to develop a fast and accurate method to detect and localize the eye saccadic signatures from the EOG signal. 1.6. Contributions This work proposes two approaches for detecting and localizing the saccadic signatures from a noisy EOG signal. Each of the approaches has a merit over the other, and the suitable method may be employed based on the situation. Our first approach is to build a discriminative dictionary for eye saccadic signatures directly from data without relying on analytical pattern or wavelet. We match a sub-window of the EOG with the dictionary and present a comparison among different distance metrics.

ACCEPTED MANUSCRIPT

700 600 0

20

800 700 600 0

40

500

0 40

15

700 600 0

20

Shifted version 1 of upper saccade

1000

500

0

60

0

EOG sample points

10

20

30

0 0

20

40

60

500

0 0

Shifted version 2 of upper saccade

500

0

EOG sample points

10

20

30

80

50

100

EOG sample points

1000

0

60

Shifted version 1 of upper saccade

40

EOG Amplitude in uV

500

EOG Amplitude in uV

EOG Amplitude in uV

1000

40

1000

EOG sample points Shifted version 2 of upper saccade

20

EOG sample points

EOG Amplitude in uV

1000

EOG Amplitude in uV

EOG Amplitude in uV

Shifted version 1 of upper saccade

20

10

Table 1. EMD Algorithm for Denoising

800

EOG sample points

EOG sample points

0

5

Upscaled by 2

900

40

Shifted version 2 of upper saccade

1000

500

Input: EOG data Output: Sample patterns 1. Find all the extrema of the signal x s (n) 2. Interpolate the maxima of the x s (n) by using a natural cubic spine to obtain the envelope emax (n) 3. Similarly, interpolate the minima of the x s (n) by using a natural cubic spine to obtain the envelope emin (n) max (n) 4. The mean curve is then obtained as m(n) = emin (n)+e 2 5. The detail is obtained as d(n) = x s (n) − m(n) 6. The residual m(n) is iterated from step 1 until stopping criterion is reached.

CR IP T

800

Downscaled by 2

900

EOG Amplitude in uV

Primitive Saccade Motif

900

EOG Amplitude in uV

EOG Amplitude in uV

4

2.2. Denoising 0 0

EOG sample points

50

100

EOG sample points

Fig. 2. Sparse representation from a sample motif saccade

AN US

Our second approach is to use a non-linear time warping function that can take care of the time normalization effect and the fluctuations in the time axis. The prime contributions of this work are as follows:

The detrended EOG signal contains interference such as sensor noise, power-line interference, etc. We denoise the resulting signal using Empirical Mode Decomposition (EMD) as proposed in [27]. Our previous research [28] suggests that EMD performs better as compared to other denoising methods for EOG. EMD splits up the EOG signal into different components known as the Intrinsic Mode Functions (IMF). Table 1 contains the denoising algorithm using EMD. The decomposition equation can be written as

• sparse representation of the saccadic signature templates to form the non-analytical dictionary.

M

• optimization of the dictionary by observing if redundant patterns are present in the dictionary.

ED

• alternate technique for the localization of saccades using DTW which is more accurate.

PT

AC

CE

The EOG signal contains baseline wander, sensor noise and blink artifacts. We have used preprocessing for proper conditioning of the signal as briefly discussed in the following subsections. 2.1. Removal of Baseline Wander Baseline wander is the effect where the baseline of the EOG signal wanders up and down rather than being straight [26]. This effect is highly undesirable as it clamps the amplitude to some other level. We assume a linear trend in the acquired EOG signal and then detrend it. For the EOG sub-sequence xeog (n) of length Nw , the resulting detrended signal, x s (n) is obtained as x s (n) = xeog (n) −

N−1 1 X xeog (n − k) Nw k=0

N I MF X

fl (n) + r(n)

(3)

l=1

NI MF gives the total number of IMF’s in the noisy detrended EOG signal, x s (n) and r(n) gives the residual signal, while fl (n) corresponds to the lth IMF. In this case, we have obtained NI MF = 11 using the stopping criteria as mentioned in [10]. 3. Method 1 - Dictionary-based Detection and Localization

Both the approaches are described in details in section 3

2. Preprocessing

x s (n) =

(2)

In this proposed method, we create a dictionary of saccadic signatures. This process is carried out by selecting a sample saccade template manually by observing the prerecorded EOG signal. This saccade template is translated and scaled at a total of 10 different scales and 10 different shifts, thereby making a dictionary of 100 atoms. One such representation is shown in Fig. 3. We compare different distance metrics and finally select the best metric to match the template saccade with the sub-sequence. Given a signal x(n) of support of size S in an infinite dimensional space, the central problem is to compute a good approxe as a linear superposition of N basic elements picked imation m up in a huge collection of signals referred to as a dictionary. Biomedical signals in their unprocessed form usually have a highly complex underlying structure which makes it difficult to explicitly define the link between a class of signals and a dictionary. However, in our case, the patterns that we are specifically looking for are identical in nature but have different temporal scales. This condition makes it feasible to construct such a dictionary element from the signal itself. We call this unit element, an atom. Learning fundamental patterns can be used as an alternate to

ACCEPTED MANUSCRIPT 5

M = {mk } ;nk=1

(4)

of real valued,finite,discrete representations of M such that a dictionary G can be formed containing the effect of all translations on the generating function M. Let there be another operator θa such that,θa translates an signal. The translation operator when applied to Mk thus generates a dictionary G represented by, nn o o G = Θ p Mk , K = 1, .......n (5) 3.1. Sparse Representation For a dictionary G, let the atoms be represented as gi . Typically a dictionary learning problem can be represented as G is restricted to a close convex set G such that, o ∆ n G = G ∈ Rm∗k (6)

G = uG [G − δt 5G L(G)]

(7)

(8)

PT

ED

M

where δt is the gradient operator, and uG represents the projector to refine the dictionary set G represented by Equation 5. Now the problem remains in defining the dictionary elements gi which can be represented as a set of translated parent structures of interest M(L). The dictionary can be represented as, g1 = α11 m1 + α12 m2 + . . . α1n mn g2 = α21 m1 + α12 m2 + . . . α2n mn . . . gn = αn1 m1 + αn2 m2 + . . . αnn mn

αi = [α11 . . . α1n ]

(10)

A = [α1 . . . αn ]

(11)

AC

CE

  g1  g   2   .  G =    .   .    gn

∆

min E(G, A) =

G∈G,A

αˆ i = [αi1 αi2 . . . αiL ]

(13)

To utilize the update, we introduce another parameter I(αˆ i ) to further introduce sparsity to the dictionary learning framework P as minG∈G,A E(D, A) + µ ni=1 I(αˆ i ), where µ is the balance parameter. The acceptance parameter I(αˆ i ) may be defined as ( 0, if αˆ i = 0 I(αˆ i ) = (14) 1 otherwise If the sum of indices of all functions of the row of the dictionary G are non-zero, then the acceptance parameter is non-zero. We assume that no zero vectors are present in the dictionary. We use the objective function to introduce adaptation in the matrix of αi which represents G in terms of the matrix G = [A][M]

(15)

Here [A] is the sparse representation while [M] is the original pattern. To further facilitate the process and reduce data redundancy in data, we introduce the Euclidean distance parameter between the atom with non-zero lower-bound. This not only decreases the lower bound but also satisfies the conditions for the signal subsequence to be isolated.

AN US

In order to compute G, we use classical first-order projected stochastic gradient descent algorithm [29], in which G is updated iteratively. In each iteration,

the pattern M, the αˆ i as a sparse representation on the dictionary G.

CR IP T

the use of analytical bases. These data-driven patterns are shift invariant in nature. Here our aim is to construct a collection of patterns M such that

(9)

The use of the Euclidean distance metric as an optimal discriminator may be justified as follows. Let us assume two dictionary atoms as g∗m and g∗r with acceptance parameters I(αˆ m ) = 1 and I(αˆ n ) = 1. The parameters must satisfy the following conditions. nµλ2 kdm∗ − dr∗ k22 ≥ (16) kφ2

Here, µ is the balance parameter. The parameters φ and k may be obtained as L X 1 φ= kmi k22 (17) 2 i=1 λ k =1+ √ nφ

n

1X 1 [ kGαi − gi k22 + λkαi k1 ] n i=1 2

3.2. Proof of the Optimum discriminator

(12)

here λ is the regularization parameter, with M = [m1 , m2 . . . mL ] is the training dataset for finding the structure of interest with sparse co-efficients A = [α1 . . . αn ] over the dictionary G in Rm . The data fitting is represented in the dictionary as kDαi − mi k22 and kαi k1 is the sparsity regularization. Now the proposed structure modulates the data by representing

(18)

For the compactness of the dictionary, we can show that n n n X X 1 X 1 [ ( kgαi − mi k22 ) + λ ||αi ||1 + µn I(αˆ i )] n i=1 2 i=1 i=1 (19) Now we have

E(G, A) =

E(X1 + X2 + X3 ) = E(X1 ) + E(X2 ) + E(X3 )

(20)

This helps us to rewrite (19) as E1 (G, A) =

n 1 X 1 [ ( kgαi − mi k22 )] n i=1 2

(21)

ACCEPTED MANUSCRIPT

0 0

50

100

EOG sample points Sample Atoms

150 100 50 0 0

50

100

EOG sample points

100

50

0 0

50

100

EOG sample points Sample Atoms

300 200 100 0 0

100

200

0

50

100

EOG sample points Sample Atoms

200 100 0 -100 0

50

100

EOG sample points Sample Atoms

300 200 100

EOG sample points

0 0

50

100

EOG sample points

200 100 0 -100 0

50

100

EOG sample points Sample Atoms

200 100 0 -100 0

50

100

EOG sample points Sample Atoms

300 200 100 0 0

50

0 0

50

100

EOG sample points Sample Atoms

-150 -200 -250 -300 0

50

100

150

EOG sample points Sample Atoms

200 100 0 -100 0

50

100

150

EOG sample points Sample Atoms

-100

EOG Amplitude in uV

EOG Amplitude in uV

100

50

EOG Amplitude in uV

0

EOG Amplitude in uV

EOG Amplitude in uV

100

100

50

EOG sample points Sample Atoms

100

EOG Amplitude in uV

50

50

EOG sample points Sample Atoms

200

0

Sample Atoms

150

EOG Amplitude in uV

100

0

300

0

EOG Amplitude in uV

150

0

600

50

EOG Amplitude in uV

100

100

400

100

EOG Amplitude in uV

50

EOG sample points Sample Atoms

200

200

EOG sample points Sample Atoms

Sample Atoms

150

Sample Atoms

300 200 100 0 0

100

-200

-300

EOG sample points

0

50

100

150

100

200

EOG sample points Sample Atoms

-150 -200 -250 -300 0

50

100

150

EOG sample points Sample Atoms

200 100 0 -100 0

50

100

EOG sample points Sample Atoms

-100

CR IP T

0

0

EOG Amplitude in uV

0

300

-400

EOG Amplitude in uV

100

200

-200

EOG Amplitude in uV

200

100

EOG sample points Sample Atoms

0

EOG Amplitude in uV

300

0

EOG Amplitude in uV

100

EOG Amplitude in uV

50

EOG sample points Sample Atoms

0 -100

Sample Atoms

200

EOG Amplitude in uV

EOG Amplitude in uV EOG Amplitude in uV

0

100

EOG Amplitude in uV

0 -100

Sample Atoms

200

EOG Amplitude in uV

100

EOG Amplitude in uV

Sample Atoms

200

EOG Amplitude in uV

EOG Amplitude in uV

6

-200

-300

0

EOG sample points

50

100

150

EOG sample points

Fig. 3. A test case showing the localization of saccades using the dictionary based method

400

AN US

EOG Amplitude in uV

600

200 0 -200 -400 1000

2000

3000 4000 5000 EOG sample points

6000

7000

M

0

EOG Saccades

n

1 X ||αi ||1 ] [λ n i=1

PT

E2 (G, A) =

ED

Fig. 4. A test case showing the localization of saccades using the dictionary based method

E3 (G, A) = µ

n X

I(αˆ i )

(22)

(23)

Here, ||.||1 is an element-wise operator. Now, λX ∗ E2 (G− , A− ) = [||α ||1 + ||α∗m + α∗r ||1 ]] n iΩ i Here, Ω is the set of α’s without αm and αr .

CE

i=1

AC

In the earlier equations, E1 (G, A) is the data fitting term, E2 (G, A) is the sparsity inducing term, while E3 (G, A) is the acceptance parameter term. Let us assume that the optimal solution be represented as G− , A− . Considering the Euclidean distance metric as a measure of dissimilarity, we can construct another representation of the dictionary as G+ , A+ . Using triangular inequality 1 E1 (G− , A− ) ≥ E(G+ , A+ ) − ||αˆ i ||22 ||dm∗ − dr∗ ||22 − 2n n 1X ∗ p ||α ||1 E(G− , A− )||dm∗ − dr∗ ||22 n i=1 i E2 (G+ , A+ ) =

λ n

n X i=1

||α∗i ||1

λ E2 (G− , A− ) = E2 (G+ , A+ ) + ||α∗m + α∗r ||1 ]] n

(27)

As the triangular inequality gets satisfied by the norms ||α∗m ||1 + ||α∗r ||1 ≥ ||α∗m + α∗r ||1

(28)

For the third case, E3 (G− , A− ) may be represented as n X i=1

(24)

(26)

I(αi ) =

X iΩ

I(αˆ i ) + I(αˆ m ) + I(αˆ r ) − I(αˆ m + αˆ r )

(29)

Thus E3 (G− , A− ) = E3 (G+ , A+ ) + µ[I(αˆ m ) + I(αˆ r ) − I(αˆ m + αˆ r )] (30) Since, gm and gr are already in the dictionary, hence it is can concluded that I(αˆ m ) = I(αˆ r ) = 1. This proves that

(25)

I(αˆ m ) + I(αˆ r ) − I(αˆ m + αˆ r ) ≥ 1

(31)

ACCEPTED MANUSCRIPT 7 Now, to estimate the upper bound of each of the functions E1 (G− , A− ), E2 (G− , A− ), E3 (G− , A− ) and combining them we have

Mahalanobis Distance (MD). The Mahalanobis Distance (MD) [33] between two vectors with covariance matrix C can be represented as,

n

nµλ2 kdm∗ − dr∗ k22 = kφ2 Here, φ=

L X 1 i=1

2

(32)

(33)

kmi k22

(34)

p

(P(i) − Q(i))T C−1 (P(i) − Q(i))

(38)

In general the equation of Mahalanobis distance uses the full covariance matrix C, which includes the covariances between the vector components. It is helpful to point out that the Euclidean distance can be rewritten as a dot product operation. 4. Method 2 - Using Dynamic Time Warping (DTW) The first approach created a dictionary which was found to be considerably fast during the detection phase. However, the scales and time shifts are discrete which indicates missing of some saccadic signatures which may not bear much similarity among the dictionary elements. One solution to this issue is to warp the time axis and find the saccade. In such a case, we do not need to create a dictionary, as a single atom if warped suitably could find the match. Dynamic Time Warping (DTW) is such an approach to find an optimal match between two given signals which may alter in time or speed [34], [35]. The signals are compared in a nonlinear manner to align with each other. In this method, the sample pattern and the signal window are made to match using a rectangular grid. The entries of the i and j indices of the grid represents the distance between the corresponding signal amplitudes. Usually a Euclidean distance is used as the distance metric. The optimal match is the path through the grid which minimizes the total distance between them. This method has a computational complexity of O(LP) in computing the DTW grid matrix. We reduce this overhead by applying some interesting properties of the optimal DTW path [36]. These properties for a DTW grid matrix D(i, j) are as follows:

AN US

3.3. Matching with the Dictionary Once the dictionary is formed, the motifs are matched based on distance measure to obtain the divergence of the subsequence from the atom. First, we outline some standard distance metrics. The metrics may be classified as Euclidean and nonEuclidean ones. Some popular Euclidean metrics compared in our work as Kullback Leibler Divergence (KLD), Jesnen Shanon Divergence(JSD), while some non-Euclidean ones are Cayley Klein Hilbert Distance, Mahalanobis distance. Here, the template is represented by P(i) and the signal window is represented by Q(i).

D(s) =

CR IP T

1 X λ ||mi ||22 = E(G− , 0) ≥ E(G− , A− ) ≥ ||αˆ i ||22 2n i=1 n

D(s) =

n X

ED

M

Kullback Leibler Divergence (KLD). Kullback Leibler Divergence (KLD) [30] measures the discrepancy between two distributions Q(i) and P(i). From the aspect of information theory, KLD is the extra data that is necessary to transfer a random process to an improperly encoded process. This is a non-symmetric distribution and we have also used variants of this distribution throughout this paper. Q(i) log (Q(i)/P(i))

i=1

(35)

PT

Here, D(S ) represents the distance between Q(i) and P(i). In this method, we take into account the information gain that is achieved if P(i) is used instead of Q(i) , and D(s) provides us with a measure of the same.

CE

Jensen Shanon Divergence (JSD). Jensen Shanon divergence (JSD) [31] is a symmetric and smoothed version of the KLD. n

1X P(i)(log(P(i)) − log(H(i)))+ 2 i=0

AC

D(s) =

n

1X Q(i)(log(Q(i)) − log(H(i))) 2 i=0

Here, H(i) =

P(i)+Q(i) 2

(36)

Cayley Klein Hilbert Distance (CLHD). Cayley Klein Hilbert Distance (CLHD) [32] is a non-Euclidean metric which measures the distance between points inside a conic defined in terms of logarithm and projective cross-ratio function. (1 − Q(i)2 ) (1 − Q(i)2 − P(i)2 )2

• The DTW path can advance only by a unit step in both the i and j indices. • The DTW path begins at the bottom-left and ends at the top-right of the DTW grid. • An optimal path is unlikely to drift very far from the diagonal joining the beginning and ending of the path. • The DTW path should not be too steep or too shallow and is governed by the slope constraint condition.

for i = 1, 2......n

D(s) =

• Both the i (row) and j (column) indices either stay the same or monotonically increase.

(37)

The Slope constraint condition is usually expressed as a ratio n m where m is the number of steps in the x direction and n in the y direction. This condition implies that, after m steps in x direction, the DTW must make a step in y direction and viceversa. By applying the constraints, we have restricted the computation of unnecessary non-optimal paths. Hence the number of paths to be considered is reduced significantly.

ACCEPTED MANUSCRIPT 8 Table 2. Comparison with Existing Methods Time (s) 1.6512 2.4524 3.6347 5.3682 19.7316

% Accuracy 81.4 90.8 92.4 95.4 96.7

% False alarm 17.5 9.4 4.3 4.1 4.2

6. Conclusion

Table 3. Comparison of the Distance Metrics for Method 1

Distance metric KLD JSD CLHD MD

Accuracy (%) 92.6 92.3 92.7 93.8

Time of detection (ms) 5.003 5.059 6.2094 6.7691

For a P-length window of the EOG signal Q(i) and a primitive saccade atom, P( j) of length L, we form the DTW matrix D(i, j) using the Euclidean distances between P( j) and Q(i)) as, q D(i, j) = (P( j) − Q(i))2 ; 0 ≤ i ≤ L, 0 ≤ j ≤ P (39)

M

AN US

Using the symmetricity property of D(i, j), we only require to compute the upper triangular matrix. The warping path W is finally selected incorporating the constraints as described earlier. There are many warping paths that satisfy the above conditions. Hence we define warping cost S as v u t K X 1 Wi (40) S = min K i=1

This paper proposes two methods to identify eye saccadic signatures in EOG. The first approach relies on the creation of a dictionary of template saccades derived from the saccadic signatures in the training EOG signals. The sample atom is subsequently scaled and shifted in time to obtain a sparse representation. The concept is similar to the creation of a wavelet family from a mother wavelet. The sparse representations finally forms a non-analytic dictionary. The sub-sequences are matched with the dictionary atoms to detect and localize the saccade in the time axis. A series of distance functions have been compared and finally we select the Mahalanobis distance, which is found to provide highest accuracy. The second approach uses a single saccadic signature to match with the EOG subsequence by dynamically warping the time axis. The first method is found to be faster, while the second method has a superior accuracy. The process of matching a particular sub-sequence being independent processes can run in parallel on a GPU and may help achieve faster rates in future for real-time applications. The proposed techniques may be useful for other biomedical signals as well for selected repetitive patterns.

CR IP T

Method Amplitude Thresholding Independent Component Analysis Discrete Wavelet Transform Method 1 - Dictionary based detection Method 2 - Dynamic Time Warping (DTW)

Sample detection using the the first method is shown in Fig. 4. The proposed methods are tested on a 2.60 GHz, Core i5 processor. The processing time has been obtained as an average of ten trials. It is observed that both the methods have superiority in terms of accuracy over existing methods.

ED

The division by K is done to compensate the effect of varying path lengths. This optimal path is the one which minimizes the above cost function. 5. Results

AC

CE

PT

We have proposed two approaches for localizing saccades from EOG data. The first method involves the creation of a database of non-analytic patterns and comparison of each database element to a windowed EOG time series as shown in Fig. 3. The comparison among the distance metrics is given in Table 3. It can be seen that accuracy and computation time varies with the choice of distance metrics, and MD gave an accuracy of 93.8% with a small increase in computation time. The second method involves the application of DTW to isolate saccades from EOG time series. The first approach is found to be faster as dictionary creation (Eqn. 15), and optimization is carried out beforehand which reduces the time required for matching. This process is to be performed only once on the dataset. But, it is evident from Table 2 that DTW despite having a higher computational cost is more efficient in terms of accuracy. This method may be implemented where the computational cost is not an issue. A comparison with traditional methods can be seen in Table 2. The proposed methods have been tested on the EOG data of the Physiosig Database [10]. The standard database is available at https : //sites.google.com/view/physiosig − database/home.

References [1] Michael Unser and Akram Aldroubi. A review of wavelets in biomedical applications. Proceedings of the IEEE, 84(4):626–638, 1996. [2] Ulrich Parlitz, Sebastian Berg, Stefan Luther, Alexander Schirdewan, J¨urgen Kurths, and Niels Wessel. Classifying cardiac biosignals using ordinal pattern statistics and symbolic dynamics. Computers in biology and medicine, 42(3):319–327, 2012. [3] Richard G Baraniuk and Douglas L Jones. Shear madness: New orthonormal bases and frames using chirp functions. IEEE Transactions on Signal Processing, 41(12):3543–3549, 1993. [4] R´emi Gribonval. Fast matching pursuit with a multiscale dictionary of gaussian chirps. IEEE Transactions on signal processing, 49(5):994– 1001, 2001. [5] Khaled N Hamdy, Murtaza Ali, and AH Tewfi. Low bit rate high quality audio coding with combined harmonic and wavelet representations. In Acoustics, Speech, and Signal Processing, 1996. ICASSP-96. Conference Proceedings., 1996 IEEE International Conference on, volume 2, pages 1045–1048. IEEE, 1996. [6] U Wiklund, M Akay, and U Niklasson. Short-term analysis of heart-rate variability of adapted wavelet transforms. IEEE engineering in medicine and biology magazine, 16(5):113–118, 1997. [7] Yvan Saeys, I˜naki Inza, and Pedro Larra˜naga. A review of feature selection techniques in bioinformatics. bioinformatics, 23(19):2507–2517, 2007. [8] R´emi Gribonval. Sparse decomposition of stereo signals with matching pursuit and application to blind separation of more than two sources from a stereo mixture. In Acoustics, Speech, and Signal Processing (ICASSP), 2002 IEEE International Conference on, volume 3, pages III–3057. IEEE, 2002. [9] Scott Shaobing Chen, David L Donoho, and Michael A Saunders. Atomic decomposition by basis pursuit. SIAM review, 43(1):129–159, 2001. [10] Anirban Dasgupta, Suvodip Chakraborty, and Aurobinda Routray. A twostage framework for denoising electrooculography signals. Biomedical Signal Processing and Control, 31:231–237, 2017.

ACCEPTED MANUSCRIPT 9

[25] [26] [27]

[28]

[29] [30]

[31]

CR IP T

[24]

John T O’brien, and Ian G McKeith. Saccadic eye movement changes in parkinson’s disease dementia and dementia with lewy bodies. Brain, 128(6):1267–1276, 2005. Cecilia Bonnet and Hanuˇska. Horizontal and vertical eye movement metrics: what is important? Clinical Neurophysiology, 124(11):2216–2229, 2013. Burkhart Fischer and H Weber. Express saccades and visual attention. Behavioral and Brain Sciences, 16(03):553–567, 1993. Manuel Blanco-Velasco, Binwei Weng, and Kenneth E Barner. Ecg signal denoising and baseline wander correction based on the empirical mode decomposition. Computers in biology and medicine, 38(1):1–13, 2008. M Sanjeeva Reddy, A Sammaiah, B Narsimha, and K Subba Rao. Analysis of eog signals using empirical mode decomposition for eye blink detection. In Multimedia and Signal Processing (CMSP), 2011 International Conference on, volume 2, pages 293–297. IEEE, 2011. Anirban Dasgupta, Suvodip Chakrborty, Aritra Chaudhuri, and Aurobinda Routray. Evaluation of denoising techniques for eog signals based on snr estimation. In Systems in Medicine and Biology (ICSMB), 2016 International Conference on, pages 152–156. IEEE, 2016. L´eon Bottou. Stochastic gradient descent tricks. In Neural networks: Tricks of the trade, pages 421–436. Springer, 2012. Wenchao Zhang, Shiguang Shan, Xilin Chen, and Wen Gao. Local gabor binary patterns based on kullback–leibler divergence for partially occluded face recognition. IEEE signal processing letters, 14(11):875–878, 2007. Bent Fuglede and Flemming Topsoe. Jensen-shannon divergence and hilbert space embedding. In IEEE International Symposium on Information Theory, pages 31–31, 2004. AF Beardon. The klein, hilbert and poincar´e metrics of a domain. Journal of computational and applied mathematics, 105(1):155–162, 1999. Roy De Maesschalck, Delphine Jouan-Rimbaud, and D´esir´e L Massart. The mahalanobis distance. Chemometrics and intelligent laboratory systems, 50(1):1–18, 2000. Cory Myers, Lawrence Rabiner, and Aaron Rosenberg. Performance tradeoffs in dynamic time warping algorithms for isolated word recognition. IEEE Transactions on Acoustics, Speech, and Signal Processing, 28(6):623–635, 1980. John Thomas, Jing Jin, Justin Dauwels, Sydney S Cash, and M Brandon Westover. Clustering of interictal spikes by dynamic time warping and affinity propagation. In Acoustics, Speech and Signal Processing (ICASSP), 2016 IEEE International Conference on, pages 749–753. IEEE, 2016. Bin Huang and W Kinsner. Ecg frame classification using dynamic time warping. In Electrical and Computer Engineering, 2002. IEEE CCECE 2002. Canadian Conference on, volume 2, pages 1105–1110. IEEE, 2002.

AC

CE

PT

ED

M

AN US

[11] Anirban Dasgupta, Suvodip Chakraborty, Pritam Mondal, and Aurobinda Routray. Identification of eye saccadic signatures in electrooculography data using time-series motifs. In India Conference (INDICON), 2016 IEEE Annual, pages 1–5. IEEE, 2016. [12] Aritra Chaudhuri, Anirban Dasgupta, and Aurobinda Routray. Video & eog based investigation of pure saccades in human subjects. In Intelligent Human Computer Interaction (IHCI), 2012 4th International Conference on, pages 1–6. IEEE, 2012. [13] Suvodip Chakraborty, Aritra Chaudhuri, Anirban Dasgupta, and Aurobinda Routray. Determining ocular gaze point fixation captured during triggered full saccades using eog. In Systems in Medicine and Biology (ICSMB), 2016 International Conference on, pages 41–44. IEEE, 2016. [14] Frank Behrens and Lutz-R Weiss. An automated and modified technique for testing the retinal function (arden test) by use of the electrooculogram (eog) for clinical and research use. Documenta ophthalmologica, 96(4):283–292, 1998. [15] Tzyy-Ping Jung, Colin Humphries, Te-Won Lee, Scott Makeig, Martin J McKeown, Vicente Iragui, Terrence J Sejnowski, et al. Extended ica removes artifacts from electroencephalographic recordings. Advances in neural information processing systems, pages 894–900, 1998. [16] P Senthil Kumar, R Arumuganathan, K Sivakumar, and C Vimal. Removal of ocular artifacts in the eeg through wavelet transform without using an eog reference channel. Int. J. Open Problems Compt. Math, 1(3):188–200, 2008. [17] Mary Hayhoe and Dana Ballard. Eye movements in natural behavior. Trends in cognitive sciences, 9(4):188–194, 2005. [18] Anwesha Banerjee, Amit Konar, R Janarthana, and DN Tibarewala. Electro-oculogram based classification of eye movement direction. In Advances in Computing and Information Technology, pages 151–159. Springer, 2013. [19] Sherif N Abbas and M Abo-Zahhad. Eye blinking eog signals as biometrics. In Biometric Security and Privacy, pages 121–140. Springer, 2017. [20] Anwesha Banerjee, Sumantra Chakraborty, Pratyusha Das, Shounak Datta, Amit Konar, DN Tibarewala, and Ramadoss Janarthanan. Single channel electrooculogram (eog) based interface for mobility aid. In Intelligent Human Computer Interaction (IHCI), 2012 4th International Conference on, pages 1–6. IEEE, 2012. [21] Anwesha Banerjee, Pratyusha Das, Shounak Datta, Amit Konar, R Janarthanan, and DN Tibarewala. Real time electro-oculogram driven rehabilitation aid. In Proceedings of International Conference on Advances in Computing, pages 435–440. Springer, 2013. [22] Joao Perdiz, Gabriel Pires, and Urbano J Nunes. Emotional state detection based on emg and eog biosignals: A short survey. In Bioengineering (ENBENG), 2017 IEEE 5th Portuguese Meeting on, pages 1–4. IEEE, 2017. [23] Urs P Mosimann, Ren´e M M¨uri, David J Burn, Jacques Felblinger,

[32] [33]

[34]

[35]

[36]