Esmaeil Seraj
email: esmaeil.seraj09@gmail.com

Cerebral Synchrony Assessment Tutorial: A General Review on Cerebral Signals' Synchronization Estimation Concepts and Methods 1 by

Keywords: Synchronization, Connectivity, Cerebral Connectivity, Phase synchrony, Phase Estimation, PSD, Correlation, Coherence, Magnitude Squared Coherence, MSC, Higher Order Spectra, HOS, Bispectrum, Multiple MSC, Thomson Multitapers, Imaginary MSC, Phase Lag Index, PLI, WPLI, Parametric Spectral Estimation, Phase Locking Value, PLV, PLS, Mean Phase Coherence, MPC, Wavelet Coherence, WC, WPLV, Empirical Mode Decomposition, EMD, Phase Analysis

The human brain is ultimately responsible for all thoughts and movements that the body produces. This allows humans to successfully interact with their environment. If the brain is not functioning properly many abilities of human can be damaged. The goal of cerebral signal analysis is to learn about brain function.

The idea that distinct areas of the brain are responsible for specific tasks, the functional segregation, is a key aspect of brain function. Functional integration is an important feature of brain function, it is the concordance of multiple segregated brain areas to produce a unified response. There is an amplified feedback mechanism in the brain called reentry which requires specific timing relations. This specific timing requires neurons within an assembly to synchronize their firing rates. This has led to increased interest and use of phase variables, particularly their synchronization, to measure connectivity in cerebral signals. Herein, we propose a comprehensive review on concepts and methods previously presented for assessing cerebral synchrony, with focus on phase synchronization, as a tool for brain connectivity evaluation.

Chapter 1 Introduction

Cerebral signals are among the easiest, most useful and efficient tools for studying brain and its different cognitive states. According to their concept (i.e. embedded information) and also the recording procedure, these signals are called in many different names such as EEG (electroencephalogram), EMG (magnetoencephalogram), ERP (event related potential) or etc.. Analyzing these signals needs some powerful and discriminative features extracted from them. Amplitude is one of these powerful features. Amplitude of Biological signals has been shown that is very informative and is used as main analysis factor in former studies.

In some applications, because of the environment being too much noisy, or presence of other biological contaminants, the amplitude information of captured signals become contaminated and inadequate. As an example of such applications we can refer to ERP signals studies in BCI purposes that the spontaneous EEG and other biological signals such as EOG play the contaminant role and make the amplitude of ERP signals inadequate in information. As an other instance, investigating fetal brain signals needs recordings from maternal abdominal. Also in this application there are some biological contaminants such as fetal and maternal cardiac activity signals (fECG, mECG). These signals are stronger, i.e. 10-100 fold, in amplitude than fetal brain signals. Thus, the amplitude solely can not yield adequate information. In such cases we need to look for new features to use and capture information from new aspects that are less contaminated and noisy.

Phase of this biological signals is an informative useful alternative feature in these cases that has been of interest in past few years. In many recent researches it has been shown that phase has valuable information and can be used as a complement for amplitude information. Thus, using phase of brain signals in cases that amplitude informations are contaminated or inadequate, or the phase itself has adjunct information over amplitude (i.e. brain cognitive response investigations), would be very helpful.

Phase of cerebral signals has been extracted and used through various procedures in past studies such as Fourier and Wavelet transforms in frequency domain, Analytic form of signals and also using multi-channel signals and the angles between signal subspaces. Each of these methods are applied on data sets related to a specific application such as brain computer interface purposes, cognitive researches, event related potentials studies and etc. Since for different applications the result of all these various methods will vary, the important confronting problem here is that there is not a comprehensive and reliable comparison between these methods. So we decided to first specify these different methods and extract a well-defined conceptual framework out of them to prevent from being perplexed in using phase analysis methods. After that we will implement some of these methods on specific data sets in specific applications to compare the results and specify the properties of each method in details. We hope that the results of this study will be useful in applications such as ERP signals analysis (i.e. BCI systems studies) and also brain cognitive responses investigations.

To this goal, in Chapter 1, the brain signals' phase analysis methods from two viewpoints are presented: [START_REF] Da | Biological aspects of eeg and magnetoencephalogram generation[END_REF] the early analysis and (2) the modern analysis. Different related quantities are defined and the statistics are presented and some benefits and limitations confronting in both the early and the modern viewpoints are discussed. Finally at the end of this section the related applications are introduced. Phase analysis is of an extended interest in this applications because of their particular characteristics.

In Chapter 2, phase analysis methods are presented and investigated afterward in different fields such as measuring connectivity in distinct brain areas, phase modulations and demodulations and high order statistics. The different approaches in each specific case such as Coherence, Phase Synchronization and Desynchronization and etc. are presented in this Section. The statistics and the motivation of each method are proposed and discussed for each part.

Cerebral Signals: Early Analysis

The human brain is ultimately responsible for all thought and movement that the body produces. This allows humans to successfully interact with their environment. If the brain is not functioning properly many abilities of human can be damaged. The goal of cerebral signal analysis is to learn about brain function.

Cerebral signals are susceptible to be contaminated by various noise sources such as muscular activities (i.e. blinking, jaw clenching or frowning), artifact generated by the electrical equipments which are near recording areas and also other biological signals (i.e. EMG or EOG). Artifact detection and removal is a necessary step for many types of cerebral signal analysis [START_REF] Da | Biological aspects of eeg and magnetoencephalogram generation[END_REF].

Early cerebral signal recording systems were very simple in comparison to the their modern species. The analogue EEG signals recorded at that time were analyzed visually and the qualitative features of them were been investigated usually. Features of the EEG signal were classified as either paroxysmal (transient, bursting) or on-going (background, spontaneous) activity. A simple feature of EEG recordings was measured by counting the number of oscillations in one second.

Because of recordings being too much noisy the quantification and interpretation of the amplitude of a desired frequency component were too difficult, since the peak-to-peak amplitude of the signals changed over time. One approximation used was to create an envelope for the signal by joining together all the peaks/troughs and taking an average of the resulting lines [START_REF] Walter | Visual imagery and electroencephalographic changes[END_REF]. Another used measure was the duration of time that a recording would spend oscillating in a particular frequency band [START_REF] Glass | A percent time alpha computer for the electroencephalograph[END_REF].

Many complex waveform that exist in EEG recording such as the alpha rhythm which was the first EEG waveform reported by Hans Berger [START_REF] Berger | Uber das elektrenkephalogramm des menschen[END_REF] (also called as the Berger waveform) and the K-complex which is a transient waveform that occurs during sleep were discovered [START_REF] Davis | Electrical reactions of the human brain to auditory stimulation during sleep[END_REF](many more examples can be found in [START_REF] Chatrian | Report of the committee on terminology[END_REF]). Also topographic maps of the scalp potentials were used to visually identify spatial patterns in EEG signals, however creating these maps were difficult from analogue recordings.

Often visual analysis was not sufficient to discriminate between components and furthermore, more complex events were very cumbersome with analogue recordings. Analogue-to-digital converters (ADC) are used to create digital EEG recordings. Digital EEG recordings opened new possibilities for analyzing the spectral content of brain signals. Creating topographic maps was made more practical with digital recordings and systems. Using the digital EEGs, computers used interpolation to create detailed contour plots of the scalp potential. However, there are also some ambiguities remaining when interpolation is used, due to projecting a 3D volume onto a 2D surface [START_REF] Binnie | EEG contour mapping I: Theory and practice[END_REF].

Based on the idea of Fourier analysis, using the Fourier transform leads to measuring the power spectral density(PSD) and so the mean squared coherence (MSC) (which employs the PSD) and also extracting significant frequency specific components such as instantaneous phase and amplitude.

1.1.0.1 Power Spectral Density Let x(t) be a stochastic stationary process, The Fourier transform of x(t) is shown by X(f ) and is determined as below:

X(f ) = ∞ -∞
x(t)e -j2πf t dt (1.1)

Then the power spectral density P SD x (f ) of x(t) is calculated as follow:

P SD x (f ) = ∞ i=-∞ |X i (f )| 2 (1.2)
The PSD can also be calculated between two distinct signals and it is called as cross-spectrum function. Cross-spectrum is a complex function that its amplitude and phase show the correlation between two signals. For instance, the cross-spectrum of two sinusoidal signals of the same frequency has a sharp peak in that frequency. The phase of the cross-spectrum at that frequency equals the phase difference between the records. If two signals share a sinusoidal activity, but each also contains other, unshared activity, their cross-spectrum has a peak at the shared frequency only.

Correlation Function

Another measure availed through digital signal recordings is the correlation function. Auto-correlation function R x (τ ) of x(t) is measured as below:

R x (τ ) = x(τ )x(-τ ) = ∞ -∞ x(t)x(t + τ )dt (1.3)
The Fourier transform of the Auto-correlation function does indeed equal to the P SD of our stochastic stationary signal x(t).

P SD x (f ) = +∞ -∞ R x (τ )e j2πf τ dτ (1.4)
The proof is quite straightforward, by mean of Fourier transform of the correlation function R x (τ ) we have:

f {R x (τ )} = ∞ -∞ R(τ )e -j2πf τ dτ = ∞ -∞ ∞ -∞ x(t)x(t + τ )e -j2πf τ dtdτ = ∞ -∞ x(t) ∞ -∞ x(t + τ )e -j2πf τ dτ dt = X(f ) ∞ -∞ x(t)e j2πf t dt = X(f )X * (f ) = |X(f )| 2 1.1.0.3 Phase & Amplitude Extraction
Two other quantities availed through digitalized early cerebral signal analysis are the instantaneous phase and amplitude. There are two conceptually distinct set of techniques for the purpose of capturing the amplitude and phase of a signal: (1) linear techniques and (2) nonlinear techniques.

In one hand, the linear techniques assume a constant amplitude and phase within the estimation window [START_REF] Bingham | Modern techniques of power spectrum estimation[END_REF]. As a common example, Fourier transform represents the stationary process x(t) in the frequency space as a phasor X(f ):

X(f ) = ∞ -∞ x(t)e -j2πf t dt = A(f ) cos θ(f ) + jA(f ) sin θ(f ) = A(f )e jθ(f )
where A(f ) and θ(f ) are the amplitude and phase of signal x(t) within the estimation window respectively.

In the other hand, nonlinear techniques try to measure the time-dependent instantaneous amplitude and phase which are illustrating the moment-to-moment change in signal. Instantaneous amplitude and phase are usually obtained using either the Hilbert or Wavelet transforms which has been illustrated that both produce similar results [START_REF] Van Quyen | Comparison of hilbert transform and wavelet methods for the analysis of neural synchrony[END_REF]. There are other recently proposed methods for phase extraction purpose which are represented here.

1.1.0.3.1 Hilbert Transform The Hilbert transform of a real valued process x(t) is defined as:

x h (t) = 1 π P V ∞ -∞ x(τ ) t -τ dτ (1.5)
where x h (t) = H{x(t)} represents the Hilbert transform of process x(t) and the integral is taken using the Cauchy principal value due to the potential singularity at t = τ . The Hilbert transform is the convolution of the signal x(t) with the function 1 πt which means the PSD of x h (t) is equal to the product of the PSD of x(t) and the PSD of 1 πt . Besides, the Fourier transform of the function 1 πt is equal to:

F { 1 πt } = -isgn(f )
which states that the Hilbert transform is simply a π 2 shift in the phase of the original signal. Thus using the analytic extension of a signal as presented below does not change the PSD of that signal [START_REF] Gabor | Theory of communication[END_REF]:

Z(t) = x(t) + jx h (t) (1.6)
Now, the instantaneous amplitude and phase of x(t) can be defined as the magnitude and argument of the analytic extension as follow:

Z(t) = A(t)e jφ(t) (1.7)
where A(t) is the instantaneous amplitude and φ(t) is the instantaneous phase which can be calculated as below:

φ(t) = arctan x h (t) x(t) (1.8)
To have a physically meaningful interpretation of the instantaneous phase and amplitude, they have to be extracted from a narrow-band signal or in a specific frequency of interest, so a filtering pre-step before applying Hilbert transform is needed [START_REF] Chavez | Towards a proper estimation of phase sychronization from time series[END_REF].

1.1.0.3.2 Wavelet Transform As mentioned before, another useful procedure for measuring the instantaneous amplitude and phase of a signal is through time-frequency transformations such as Wavelet transform. By mean of the standard Morlet-Wavelet, the Wavelet coefficients of signal x(t) can be determined as follow [START_REF] Mallat | A wavelet tour of signal processing[END_REF]:

W x (t, f ) = ∞ -∞ x(u)ψ * t,f (u)du (1.9)
where the ψ t,f (u) is a modulated Gaussian function (with center of t and variance of σ) by a sinusoidal wave with frequency of f :

ψ t,f (u) = f e j2πf (u-t) e -(u-t) 2 2σ 2
Then two important parameters obtained from Wavelet transform are spectrogram and Instantaneous Phase [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. We have:

Spectrogram = |W x (t, f )| 2 (1.10) and exp(jφ x (t, f )) = W x (t, f ) |W x (t, f )| (1.11)
which in latter the φ x (t, f ) represents the instantaneous phase.

1.1.0.3.3 Phase Extraction Based on Complex Energy Density Function "Rihaczek derived the signal energy distribution in time and frequency by application of the complex signal notation" [START_REF] Aviyente | A Phase Synchrony Measure for Quantifying Dynamic Functional Integration in the Brain[END_REF].

Considering two complex signals x(t) and y(t) which are at similar frequency, the total complex energy can be defined as:

CE(t, f ) = ∞ -∞ x(t)y * (t)dt (1.12)
Rihaczek utilized this idea to compute the interaction energy of a signal within some frequency band. This leads to the Complex Energy Density Function (CEDF), stated as below [START_REF] Rihaczek | Signal energy distribution in time and frequency[END_REF]:

CEDF (t, f ) = 1 √ 2π x(t)X * (f ) exp(-j2πf t) (1.13)
CEDF measures the complex energy of a signal around time t and frequency f . As stated in [START_REF] Aviyente | A Phase Synchrony Measure for Quantifying Dynamic Functional Integration in the Brain[END_REF], CEDF provides a better understanding of phase-modulated signals. Now, the time-varying phase for signal x(t) in time-frequency plane is defined as below: [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis[END_REF], and is based on the assumption that any signal consists of different simple intrinsic independent modes of oscillation. Each of these modes will have the same number of extrema and zero-crossings and there is only one extremum between successive zero-crossings [START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF]. Here we present the statistics given in [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis[END_REF] for EMD algorithm with respect to the signal x(t):

φ x (t, f ) = arg CEDF (t, f ) |CEDF (t, f )| (1.14)
• Find all the local extrema. Then connect all the local maxima together and also all the local minima by a cubic spline as the upper and lower envelope respectively, (e max (t)ande min (t)).

• Extract the detail as the difference between the signal and the mean of the upper and lower envelope values :

d(t) = x(t) -emax(t)+e min (t) 2 .
• If d(t) is not an IMF, treat it as original signal and repeat first two steps until d(t) becomes an IMF, then d(t) = I 1 (t) is the first IMF.

• After getting the first IMF, remove it from the original signal and obtain the residual r 1 :

x(t) -I 1 (t) = r 1 (t).

• Treat the r 1 (t) as the original signal and repeat the above steps to obtain the second IMF, I 2 (t).

• Repeat the process n times. The decomposition process can be stopped when no more IMF components can be extracted from the last residual.

Thus, the EMD algorithm decomposes the signal x(t) into n empirical modes as:

x(t) = n i=1 I i (t) + r n (t) (1.15)
The IMFs I 1 (t), ..., I n (t), include different frequency bands ranging from high to low [START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF]. Now, after decomposing signal x(t) into its relative IMFs, for the purpose of phase extraction, Hilbert transform can be applied to the IMFs and produce the instantaneous phase sequences.

The main difference between this method and other phase extraction method is that, here, phase sequences are extracted from IMFs and there is no need for band-pass filtering the signal anymore. Also EMD based methods have some other merits which are discussed in next chapter while introducing EMD based phase analysis method for assessing synchronization in brain.

Coherence

The idea that distinct areas of the brain are responsible for specific tasks, the functional segregation, is a key aspect of brain function. Functional integration is another important feature of brain function, it is the concordance of multiple segregated brain areas to produce a unified response [START_REF] Felleman | Distributed hierarchical processing in the primate cerebral cortex[END_REF], [START_REF] Finkel | Integration of distributed cortical systems by reentry: A computer simulation of interactive functionally segregated visual areas[END_REF].

To measure functional integration in brain, different procedures were introduced to quantify functional connectivity, a statistical dependence between different cerebral signals captured from different brain regions. The first famous measure of functional connectivity was coherence.

Coherence determines the correlation between the signals at specific frequencies [START_REF] Walter | Spectral analysis for electroencephalograms: Mathematical determination of neurophysiological relationships from records of limited duration[END_REF], [START_REF] Marshall | Statistical Analysis of EEG Phase Shift Events[END_REF]. Coherence in a simple definition is a measure of synchronization and correlation between two random process, stochastic, wide-sense stationary (WSS ) signals that is computable through the coherence function. Coherence is widely used to study dependency and relationship between different brain regions particularly during a specific cognitive task or while experiencing a specific stimuli [START_REF] Carter | Estimation of the Magnitude-Squared Coherence Function Via Overlapped Fast Fourier Transform Processing[END_REF].

In general Coherence function is a complex valued number with both an amplitude and phase that are used to measure phase synchronization in signals. The main advantage of coherence over other correlation calculating methods (i.e. in [START_REF] Do | Spectral Analysis for Electroencephalograms: Mathematical Determination of Neuro physiological Relationships from Records of Limited Duration[END_REF]) is that in fact coherence gives the correlation between two signals as a function of frequency that provides the possibility of studying spatial correlation in different frequency bands [START_REF] Carter | Estimation of the Magnitude-Squared Coherence Function Via Overlapped Fast Fourier Transform Processing[END_REF].

According to the different characteristics of signals (i.e. stochastic or deterministic) and also different procedures presented for determining PSDs in literature, several approaches for MSC calculation exist. Here we represent the most widely used approaches and in next Chapter after reintroducing the conventional MSC, The rest of the approaches will be discussed.

1.1.0.4.1 MSC Using FFT In [START_REF] Carter | Estimation of the Magnitude-Squared Coherence Function Via Overlapped Fast Fourier Transform Processing[END_REF] a method for measuring coherence function is presented. The method is based on weighted windowing of the Fourier transform of signals. Let x(t) and y(t) be two random process, zero mean, wide-sense stationary and ergodic signals with length l. A summary of the method in [START_REF] Carter | Estimation of the Magnitude-Squared Coherence Function Via Overlapped Fast Fourier Transform Processing[END_REF] is given below:

• Divide the signals into N equal parts with T samples (different part can be overlapped or disjoint)

• Multiply the samples of each ensemble in a weighting function (i.e. cosine)

• Take FFT of each weighted ensemble

• Measure power spectral densities of signals (PSD)

• Finally, measure the coherence function (or as said "Mean Squared Coherence (MSC)")

If X i (f ) and Y i (f ) be the Fourier transform of i-th ensemble, the PSDs are calculated as follow:

P SD xx (f ) = 1 NT N i=1 |X i (f )| 2 (1.16) P SD yy (f ) = 1 NT N i=1 |Y i (f )| 2
(1.17)

P SD xy (f ) = 1 NT N i=1 X i (f )Y * i (f ) (1.18)
With the ergodicity of signals being neglected we have:

P SD xx (f ) = E{|X i (f )| 2 } (1.19) P SD yy (f ) = E{|Y i (f )| 2 } (1.20) P SD xy (f ) = E{X i (f )Y * i (f )} (1.21)
where E{.} illustrates the mathematical expectation. Using these PSDs we compute the mean squared coherence MSC as:

|MSC(f )| 2 = |P SD xy (f )| 2 P SD xx (f )P SD yy (f ) = 1 N N i=1 |(P SD xy (f )) i | 2 P SD xx (f )P SD yy (f ) (1.22)
Since P SD xy (f ) is complex, the Cartesian or Polar coordinates can be used for averaging. We discussed This issue in Appendix A.

1.1.0.4.2 MSC Using Correlation Functions Another method for obtaining complex coherence is presented in [START_REF] Carter | Coherence and Time Delay Estimation[END_REF] which computes the PSDs using auto and cross correlations. In this method after dividing the signals into N separate parts, first the auto/cross correlations R xx , R yy and R xy are calculated and then the PSDs are measured by taking Fourier transform of these correlation functions:

R xx (τ ) = E{x(t)x(t + τ )} ⇒ P SD xx = +∞ -∞ R xx (τ )e j2πf τ dτ (1.23) R yy (τ ) = E{y(t)y(t + τ )} ⇒ P SD yy = +∞ -∞ R yy (τ )e j2πf τ dτ (1.24) R xy (τ ) = E{x(t)y(t + τ )} ⇒ P SD xy = +∞ -∞ R xy (τ )e j2πf τ dτ (1.25)
The equations (1.23), (1.24) and (1.25) can also be illustrated using inner product by the complex sinusoidal mother function as follow: This method also takes all the assumptions considered in former method. Here, as well, if the ergodicity assumption is neglected thus the expectation term in correlation calculation could be changed by time averaging. Finally the MSC is obtained as below:

P SD xx =< R xx ,
MSC(f ) = P SD xy (f ) P SD xx (f )P SD yy (f ) = F {R xy(f ) } F {R xx (f )}F {R yy (f )} (1.29)
where F {.} denotes the Fourier transform. MSC varies between 0 to 1 (0 < MSC < 1). The more two signals are related and synchronized, the grater the MSC value will become. In fact a MSC value close to 1 means the activity between corresponding signals follows a linear transformation in that frequency of interest and a MSC near zero shows two non-related signals in that frequency.

As discussed before, cross-spectrum is a complex function that its amplitude and phase show the correlation between two signals. For instance, the cross-spectrum of two sinusoidal signals of the same frequency has a sharp peak in that frequency. The phase of the cross-spectrum at that frequency equals the phase difference between the records. If two signals share a sinusoidal activity, but each also contains other, unshared activity, their cross-spectrum has a peak at the shared frequency only. The MSC has the same interpretation to the cross-spectrum and concluded from the discussion above, it is said that the phase of MSC is meaningful and equals the phase difference between two corresponding signals, if the amplitude in that frequency is considerable. The large value of coherence amplitude in a phase of 0 degree indicates the large positive correlation and in a phase of 180 degree it's the vise versa and the the rest of the phase values are interpreted similarly due to the amplitude in that frequency [START_REF] Do | Spectral Analysis for Electroencephalograms: Mathematical Determination of Neuro physiological Relationships from Records of Limited Duration[END_REF].

The discriminant factor between methods presented in [START_REF] Carter | Estimation of the Magnitude-Squared Coherence Function Via Overlapped Fast Fourier Transform Processing[END_REF] and [START_REF] Carter | Coherence and Time Delay Estimation[END_REF] is the procedure of measuring power spectral densities of signals. As reported in [START_REF] Carter | Coherence and Time Delay Estimation[END_REF], the PSD can be estimated using non-parametric methods such as Blackman-Tukey method, weighted overlapped segment averaging (WOSA), also known as Welch's method, Lag-reshaping method or parametric ones, as discussed later in next Chapter.

Cerebral Signal: Modern Analysis-Cerebral Connectivity

Brain connectivity depicts patterns of links in the brain. Most of the brain functions, for instance the ones involved in learning, memory, behavior adaptation to stimuli, emotions, as well as pathological processes in some brain disorders like epilepsy, autism or schizophrenia, are based on interactions between neuronal assemblies distributed within and across distinct cerebral regions [START_REF] Hassan | Graph-based analysis of brain connectivity during spelling task[END_REF]. For instance, it has been shown that segregated areas may activate in response to a particular cognitive task. Neural assemblies are the basic concept of current scientific model for how cerebral functional integration is achieved. Integration of cerebral areas can be measured by assessing brain connectivity [START_REF] Hassan | Graph-based analysis of brain connectivity during spelling task[END_REF], [START_REF] Wendling | From EEG signals to brain connectivity: A model-based evaluation of interdependence measures[END_REF], [START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF]. Brain connectivity can be separated into three different concepts: (1) structural (or neuroanatomical ) connectivity, which denotes anatomical links, (2) functional connectivity, which rely on statistical dependencies between signals from different brain areas and (3) effective connectivity, which introduces causal interactions among these signals [START_REF] Hassan | Graph-based analysis of brain connectivity during spelling task[END_REF], [START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF], [START_REF] Sporns | Networks of the Brain[END_REF]. There are reasons discussed in [START_REF] Van Ooyen | Competition in the development of nerve connections: a review of models[END_REF] which make structural connectivity difficult to define intrinsically. functional connectivity is defined as the temporal correlation among the activity of different neural assemblies whereas the direct or indirect effect that one neural system exerts over another is defined as effective connectivity [START_REF] Fingelkurts | Functional connectivity in the brain is it an elusive concept?[END_REF], [START_REF] Horwitz | The elusive concept of brain connectivity[END_REF].

Functional and effective connectivity techniques are intensely dependent on measuring the correspondence of neural signals over time. Thus cerebral signal recording techniques such as EEG and MEG with good temporal resolution, are optimal for calculating such connectivity. Most frequently used measures for connectivity are: correlation, coherence, mutual information, transfer entropy, generalized synchronization [START_REF] Pereda | Causal Influence: Nonlinear multivariate analysis of neurophysical signals[END_REF], continuity measure [START_REF] Netoff | Detecting coupling in the presence of noise and nonlinearity[END_REF], synchronization likelihood [START_REF] Stam | Synchronization likelihood: an unbiased measure of generalized synchronization in multivariate data sets[END_REF] and phase synchronization [START_REF] Rosenblum | Phase Synchronization of Chaotic Oscillators[END_REF].

There is an amplified feedback mechanism called reentry which requires specific timing relations. This specific timing requires neurons within an assembly to synchronize their firing rates. This has led to increased interest and use of phase variables, particularly their synchronization, to measure connectivity in cerebral signals [START_REF] Marshall | Statistical Analysis of EEG Phase Shift Events[END_REF]. Moreover, our focus in this work is based on signal phase and its analysis methods. Thus, we will introduce the phase synchrony as a subset of modern methods for cerebral signal analysis and discuss various approaches presented to achieve this quantity in next section.

The next section is dedicated to introducing synchronization and its different aspects. Each case is introduced and its characteristics are represented briefly. Finally the phase synchronization as one of the major investigated fields in this thesis and its various approaches are presented.

Synchronization

Connectivity measures are widely used to determine the level of synchronization between distinct brain regions using cerebral recordings. As discussed before, since structural connectivity is difficult to define inherently, only the two remaining connectivity measures (functional connectivity and effective connectivity) are used for determining synchronization in brain.

As discussed in [START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF], methods to measure effective connectivity can be subdivided into two main categories: (1) model-based and (2) data-driven techniques. The former case assumes theoretical models that describe how brain areas interact and influence each other, whereas the latter case does not assume any specific underlying model or prior knowledge concerning spatial or temporal relationships. Granger-causality (GC), directed transfer function (DTF), directed coherence (DC) and its extension partial directed coherence (PDC) are the most common data-driven connectivity techniques. It is notable that DTF, DC and PDC are developed out of the GC method. The PDC technique will be presented in next chapter as an extension of coherence methods.

The functional connectivity methods are subdivided into three main categories: (1) linear, (2) nonlinear and (3) information-based techniques [START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF]. Each of these categories and the corresponding most common approaches are introduced in next subsections.

Linear Synchronization Methods

Cross Correlation

Cross-correlation and MSC are the most commonly used linear synchronization methods. As we discussed before (see Sections 1.2.1.2 and 1.2.1.4), assuming two discrete time processes x n and y n where n = 1, ..., N, then the cross-correlation function R xy is defined as:

R xy (τ ) = 1 N -τ N -τ n=1 ( x n -µ x σ x )( y n+τ -µ y σ y ) (1.30)
where µ and σ denote mean and variance, respectively, while τ is the time lag.

MSC

Another linear technique for synchronization measurements is the MSC or simply coherence. MSC is defined precisely in Section (1.1.0.4). To remind what the basic idea is, the main formula to calculate MSC is given below:

|MSC(f )| 2 = 1 N N i=1 |(P SD xy (f )) i | 2 P SD xx (f )P SD yy (f ) (1.31)
which is averaged over trials.

Although MSC is a successful measure of functional integration, it has several limitations [START_REF] Marshall | Statistical Analysis of EEG Phase Shift Events[END_REF]:

• Stationarity: Because of using FFT, existence of a PSD depends on the stationarity of a process. As reported in [START_REF] Cohen | Stationarity of the human electroencephalogram[END_REF], there is evidence that EEG recordings longer than a particular period of time are non-stationary.

• Linearity: Coherence is based on linear correlation in the process. Thus, if information is being transferred in a nonlinear manner, it may not be detected.

• Direction: When there is synchronization between two signals, coherence does not clarify which signal is driving the flow of information.

Nonlinear Synchronization Methods

Many crucial neural processes have nonlinear characteristics. Due to the lack of linear synchronization methods to detect this nonlinear connectivity, basic motivation for developing nonlinear methods were not put on outperforming linear methods but rather providing complementary information under certain assumptions [START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF].

Within the field of neuroscience, synchronization is commonly illustrated in the concept of phase synchrony. There are various approaches to obtain this quantity based on different characterization procedures. The main ideas are described briefly in next subsection and more comprehensively in next chapter. Another nonlinear connectivity measure termed as Generalized synchronization is also introduced very briefly.

Phase Synchronization

Phase synchronization (PS) is defined as temporal adjustment of rhythms of a pair of signals in a frequency band of interest, whereas the amplitudes can remain uncorrelated [START_REF] Rudrauf | Frequency flows and the time-frequency dynamics of multivariate phase synchronization in brain signals[END_REF]. PS also can be determined through successive frequency bands to provide a complete measure of synchrony between two signals. Phase locking-value (PLV) is the most commonly used measure of PS. In order to detect periods of phase synchrony (i.e. PLV) in cerebral signals, there are four main steps summarized as follow [START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF], [START_REF] Lachaux | Studying single-trials of phase-synchronous activity in the brain[END_REF], [START_REF] Tass | Detection of n:m phase locking from noisy data: application to magnetoencephalography[END_REF], [START_REF] Palus | Detecting phase synchronization in noisy systems[END_REF], [START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF]:

• The first step is to calculate a sequence consisting instantaneous frequency specific phase values within a temporal window.

• The second step simply is to determine the instantaneous phase-differences (in frequency f and time t) from the corresponding phase sequences captured in first step.

• Quantifying the local stability of this phase-differences across trials of cerebral signals.

• Determining the degree of statistical significance of each quantity obtained in step three.

The first step can be performed through one of the two common signal phase extraction methods: (1) Hilbert Transform or (2) Wavelet transform (see Section 1.2.1.3). The second step actually depends on the third step. The quantification in step three can be performed through various statistical dependence parameters such as mean phase difference, circular variance, standard deviation, Shannon entropy, or mutual information [START_REF] Rudrauf | Frequency flows and the time-frequency dynamics of multivariate phase synchronization in brain signals[END_REF]. Depending on which of these procedures to be selected, the process in second step might be required or not. If the phase differences have to be found, there are different ways to this goal depending on the method used for phase extraction. This might differ between simply subtracting the corresponding time-dependent phases captured through Hilbert Transform or using description below due to phase extraction through Wavelet Transform:

exp(j(φ y (t, f ) -φ x (t, f ))) = W x (t, f )W * y (t, f ) |W x (t, f )||W y (t, f )| (1.32)
where the term φ y (t, f ) -φ x (t, f ) represents the phase-difference between processes y and x. The main conventional approaches to assess the synchrony between neural signals are distinct in their technique to quantify the captured instantaneous phase-differences. Finally, for the last step commonly surrogate data are utilized. There are different techniques presented for generating surrogate data in literature such as shifting or scrambling the original series [START_REF] Lachaux | Studying single-trials of phase-synchronous activity in the brain[END_REF], [START_REF] Palus | Detecting phase synchronization in noisy systems[END_REF]. The degree of statistical significance of the PLVs is then determined by comparing them to values obtained from surrogate data.

Generally speaking, for two neural signals x(t) and y(t) and their corresponding phases φ x (t) and φ y (t), phase synchrony in its most general form is defined as:

nφ x (t) -mφ y (t) = const. (1.33)
where n and m are integers indicating the ratios of possible frequency locking [START_REF] Rosenblum | Phase Synchronization of Chaotic Oscillators[END_REF]. In what follows we assume n = m = 1 for simplicity. Because of problems such as Volume Conduction (introduced later in current chapter) or background spontaneous brain signals the true synchrony is buried in a considerable background noise [START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF]. In this case, as in [START_REF] Rosenblum | Phase Synchronization of Chaotic Oscillators[END_REF], this condition of phase locking can be replaced by the weaker condition as stated below:

|nφ x (t) -mφ y (t)| < const. (1.34)
As discussed in [START_REF] Rudrauf | Frequency flows and the time-frequency dynamics of multivariate phase synchronization in brain signals[END_REF], such approaches to phase synchrony impose two main limitations as presented below:

• As mentioned, only synchrony between pairs of signals can be directly studied. This extremely increases the computational cost. For instance, when there is

n signals, d = n 2 -n 2
pairs of signal exist.

• As illustrated in [START_REF] Rudrauf | Frequency flows and the time-frequency dynamics of multivariate phase synchronization in brain signals[END_REF], synchronizations with frequency non-stationarity can be observed in brain signals (frequency of synchronization alters continuously through time, although the phasedifference remains stable). Thus, defining phase synchronization in this manner will prevent such phase synchronizations from being detected.

• It has been demonstrated in [START_REF] Breakspear | Nonlinear phase desynchronization in human electroencephalographic data[END_REF] that synchrony-measurement methods based on this definition may be insensitive to very short periods of phase synchronization (see also [START_REF] Freeman | Fine temporal resolution of analytic phase reveals episodic synchronization by state transitions in gamma EEGs[END_REF]).

The relation between coherence and phase synchronization is investigated in [START_REF] Mezeiova | Comparison of coherence and phase synchronization of the human sleep electroencephalogram[END_REF] and according to the reported results there is a very close relationship between these concepts. In many studies it has been reported that a complete phase synchronization is manifested by highly coherent phases and correlated amplitudes. In fact, both the phase synchronization and the coherence analysis are looking for periods of phase synchrony in various frequency bands but there is some differences between their sensitivity to detecting them. With focusing on Signals phase, the coherence analysis has some objections (see Section 1.3.1.2), however, in [START_REF] Mezeiova | Comparison of coherence and phase synchronization of the human sleep electroencephalogram[END_REF] it is reported that phase synchronization analysis is not superior than coherence considerably.

Generalized synchronization

After the successful application of PS in EEG analysis, the concept of Generalized synchronization (GS) was introduced [START_REF] Arnhold | A robust method for detecting interdependences: application to intracranially recorded EEG[END_REF]. GS approaches were proposed to investigate the dependencies between nonlinear signals without any knowledge about the governing equations [START_REF] Wendling | From EEG signals to brain connectivity: A model-based evaluation of interdependence measures[END_REF]. GS represents how much neighborhoods of one chaotic attractor can be mapped onto the other [START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF]. Although this mapping is prone to stationarity deficits, it is considered to be a robust way of evaluating the GS [START_REF] Arnhold | A robust method for detecting interdependences: application to intracranially recorded EEG[END_REF]. Generating such attractors needs delay vectors to be constructed from time series using a procedure known as time-delay embedding [START_REF] Takens | Detecting strange attractors in turbulence, Dynamical systems and Turbulence[END_REF].

Information-based Techniques

In probability theory and information theory, Mutual Information (MI) measures the mutual dependencies between two random-variable processes by quantifying the amount of information gained about one signal from measuring the other. If two signals have a statistical dependence, observing one signal gains information about the other and reduces the entropy because of the knowledge of other signal. The benefit of using information-based techniques is that these techniques are sensitive to both linear and nonlinear statistical dependencies between signals [START_REF] Marshall | Statistical Analysis of EEG Phase Shift Events[END_REF], [START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF].

MI analysis can be used to assess functional connectivity between cerebral signals. For this purpose a measure of entropy is needed [START_REF] Vastano | Information transport in spatiotemporal system[END_REF]. Defining the entropy for a signal x(t) as:

H(x) = - x P (x) log(P (x)) (1.35)
we have the mutual information between signals x(t) and y(t) determined as:

I(x, y) = y x
P (x, y) log( P (x, y) P (x)P (y) )dxdy (1.36) where P (x, y) is the joint probability density function of x(t) and y(t), and P (x) and P (y) are the marginal probability density functions of x(t) and y(t) respectively [START_REF] Vastano | Information transport in spatiotemporal system[END_REF]. MI has been used for evaluating functional connectivity in EEG recordings (see [START_REF] Xu | Information transmission in human cerebral cortex[END_REF] and [START_REF] Mars | Spread of epileptic seizure activity in humans[END_REF]). Also MI and entropy based features have been used previously for CAP detection in sleep EEG recordings [START_REF] Karimzadeh | Presenting efficient features for automatic CAP detection in sleep EEG signals[END_REF].

Cerebral Signal: Modern Analysis-Higher order Spectra

The previous introduced statistical tools utilize first and second order statistics to extract information from random signals. However, the main problem confronts in the presence of nonlinearity in systems where the first and second order statistics are unable to adequately analyze many signals [START_REF] Sanaullah | A Review of Higher Order Statistics and Spectra in Communication Systems[END_REF]. Presenting High-order Spectra (Statistics) (HOS) is mainly motivated to overcome this problem and plays an important role in digital signal processing. HOS methods are very useful in problems where non-Gaussian, non-minimum phase, phase coupling or nonlinear behavior and robustness to additive noise are important [START_REF] Sanaullah | A Review of Higher Order Statistics and Spectra in Communication Systems[END_REF]. They are also used in detection and classification in communication and pattern recognition applications. "In general, there are three motivations behind the use of HOS in signal processing which can be considered as its property over lower order spectra: (1) suppressing Gaussian noise of unknown mean and variance; (2) detecting and characterize nonlinearities in the data and (3) reconstructing and preserving the phase as well as the magnitude response of signals or systems" [START_REF] Petropulu | Higher-Order Spectral Analysis[END_REF]. Also the applications of HOS on biomedical signals are based on these properties. Particularly in this contribution we focus on the latter property of HOS. In fact, measuring Power spectrum causes phase relations and signals phase informations to be suppressed, while HOS, based on the ability of cumulant spectra to preserve the Fourier-phase of signals, contains such informations [START_REF] Nikias | Bispectrum Estimation -A Digital Signal Processing Framework[END_REF], [START_REF] Nikias | Higher-order spectra analysis: A nonlinear signal processing framework[END_REF] and [START_REF] Nikias | Signal processing with higher-order spectra[END_REF].

Higher order spectra are functions of two or more component frequencies while the power spectrum and lower order spectra are functions of a single frequency [START_REF] Sanaullah | A Review of Higher Order Statistics and Spectra in Communication Systems[END_REF]. The bispectrum and bicoherence and their related extensions (i.e. auto/cross bispectrum and bicoherence) are the third order spectra and also the most used HOS quantities as they are easiest computationally. The other less common HOS quantities than the 3-rd orders are the (auto/cross) trispectrum and tricoherence. First, we need to review approaches for calculating moments and cumulants of order n as high order spectra quantities' computation is based on them.

Moments and Cumulants of order n

For a stationary discrete time random process X(k), the moments of order n are given by:

m n (τ 1 , τ 2 , ..., τ n-1 ) = E{X(k)X(k + τ 1 )...X(k + τ n-1 )} (1.37)
where E{.} denotes expectation.

The n -th order cumulants are functions of the moments of order up to n, for example the first order cumulants can be obtained as below:

c 1 = m 1 = E{X(k)} (1.38)
which is the mean. Consequently, the second order cumulants are given by:

c 2 (τ 1 ) = m 2 (τ 1 ) -(m 1 ) 2 (1.39)
which is the covariance. Also for the third order cumulants we have:

c 3 (τ 1 , τ 2 ) = m 3 (τ 1 , τ 2 ) -(m 1 )[m 2 (τ 1 ) + m 2 (τ 2 ) + m 2 (τ 2 -τ 1 )] + 2(m 1 ) 3 (1.40)
where m 3 (τ 1 , τ 2 ) in the third order moment sequence. Following this manner, the n -th order cumulants can be obtained. The general relations between cumulants and moments are given in [START_REF] Nikias | Higher-order spectra analysis: A nonlinear signal processing framework[END_REF].

1.4.2 Bispectrum 1.4.2.1 Auto-bispectrum
The bispectrum of a signal x(t) can be defined in two ways. First is the triple product of the DFT's [START_REF] Nikias | Higher-order spectra analysis: A nonlinear signal processing framework[END_REF]:

B xxx (f 1 , f 2 ) = lim T →∞ 1 T E{X T (f 1 )X T (f 2 )X * T (f 1 + f 2 )} (1.41)
The bispectrum is complex. It means that it contains Fourier magnitude and phase information [START_REF] Fackrell | DETECTING PHASE COUPLING IN SPEECH SIG-NALS[END_REF]. With the above definition we have:

B xxx (f 1 , f 2 ) = |B x xx(f 1 , f 2 )| exp(jφ(f 1 , f 2 )) (1.42)
where the magnitude of bispectrum obtains as:

|B xxx (f 1 , f 2 )| = |X(f 1 )||X(f 2 )||X * (f 1 + f 2 )| (1.43)
and the phase of bispectrum (or biphase [START_REF] Fackrell | DETECTING PHASE COUPLING IN SPEECH SIG-NALS[END_REF]) is:

φ(f 1 , f 2 ) = φ(f 1 ) + φ(f 2 ) -φ(f 1 + f 2 ) (1.44)
Second definition of bispectrum is the double DFT of the 3rd order cumulant. Two methods for calculating bispectrum based on this definition in real data are presented later in current section.

The auto-bispectrum is a function of two frequency components, f 1 and f 2 which would have a small value if the biphase varies over the different realizations and conversely, would have a large value if this phase does not vary, which is an indication of quadratic coupling between f 1 and f 2 [START_REF] Chabalko | Identification of Transient Nonlinear Aeroelastic Phenomena[END_REF].

Cross-bispectrum

Following the definition of the auto-bispectrum, the cross-bispectrum is used to determine the quadratic coupling between frequency components in two different signals. Considering two signals with Fourier transforms as X(f ) and Y (f ), the cross-bispectrum is defined by [START_REF] Chabalko | Identification of Transient Nonlinear Aeroelastic Phenomena[END_REF]:

B xxy (f 1 , f 2 ) = lim T →∞ 1 T E{X T (f 1 )X T (f 2 )Y * T (f 1 + f 2 )} (1.45)
Through this definition, the coupling level between two frequency components in X(f ) namely f 1 and f 2 and their algebraic sum in Y (f ) is determined which is based on the phase relation between these components in the different realizations. A large value of the cross-bispectrum indicates high Quadratic Phase Coupling (QPC) between X(f ) and Y (f ). In fact QPC normally is obtained through the normalized extension of bispectrum known as bicoherence. Bicoherence as a tool for phase analysis is presented within next Chapter [START_REF] Chabalko | Identification of Transient Nonlinear Aeroelastic Phenomena[END_REF].

Methods for Calculating Bispectrum

There are two main methods for calculating bispectrum through its definition based on 3rd order cumulants: (1) indirect method and (2) direct method.

1.4.2.3.1 Indirect Method Let x(k), k = 1, ..., L be the available discrete signal. First, we have to segment the data into N parts each with M samples. Now x i (k), k = 1, ..., M is the i -th segment. Then, After subtracting the mean of each segment the 3rd order moment of each segment is computed as below:

m x i 3 (τ 1 , τ 2 ) = 1 M l 2 l=l 1 X i (l)X i (l + τ 1 )X i (l + τ 2 ) (1.46)
where l 1 = max(0, -τ 1 , -τ 2 ) and l 2 = min(M -1, M -2). Since we made each segment zero-mean, its third-order moments and cumulants are equal. Thus, now the average cumulants have to be calculated:

c x 3 (τ 1 , τ 2 ) = 1 N N i=1 m x i 3 (τ 1 , τ 2 ) (1.47)
Finally, the third order spectrum (bispectrum) can be estimated as [START_REF] Petropulu | Higher-Order Spectral Analysis[END_REF]:

B xxx (f 1 , f 2 ) = L τ 1 =-L L τ 2 =-L c x 3 (τ 1 , τ 2 ) exp(-j(f 1 τ 1 + f 2 τ 2 ))ω(τ 1 , τ 2 ) (1.48)
where L < M -1 and ω(τ 1 , τ 2 ) is a two-dimensional window, introduced to smooth out edge effects [START_REF] Petropulu | Higher-Order Spectral Analysis[END_REF]. A complete description of appropriate windows that can be used and their properties can be found in [START_REF] Nikias | Higher-order spectra analysis: A nonlinear signal processing framework[END_REF].

1.4.2.3.2 Direct Method Regarding to assumed signal x(k) in previous method, again in this method we have to segment the data and zero-mean each segment similar to previous method. After these steps, the Fourier transform of each segment based on M points have to be computed:

F x i (k) = M -1 l=0 x i (l) exp(-j 2π M lk) (1.49)
where k = 0, 1, ..., M -1. Now the third order spectrum (bispectrum) of each segment is obtained as:

C x i 3 (k 1 , k 2 ) = 1 M F x i (k 1 )F x i (k 2 )F * x i (k 1 + k 2 ) (1.50)
which is similar to the first definition of bispectrum. The C x i 3 (k 1 , k 2 ) need to be computed only in the triangular region 0 ≤ k 2 ≤ k 1 , K -1 + k 2 < M/2. Also a smoothing window for reducing the variance can be performed around each frequency. With these we have:

Ĉx i 3 (k 1 , k 2 ) = 1 W 2 (W/2)-1 l 1 =-W/2 (W/2)-1 l 2 =-W/2 C x i 3 (k 1 + l 1 , k 2 + l 2 ) (1.51)
where the smoothing window is of size (W × W ). Finally, the third order spectrum (bispectrum) is given as [START_REF] Petropulu | Higher-Order Spectral Analysis[END_REF]:

B xxx (f 1 , f 2 ) = 1 N N i=1 Ĉx i 3 (f 1 , f 2 ) (1.52)
Both the direct and the indirect methods produce asymptotically unbiased and consistent bispectrum estimates [START_REF] Petropulu | Higher-Order Spectral Analysis[END_REF].

Major Obstacles

There are some main obstacles confronting through modern cerebral signal analysis methods, particularly when measuring connectivity in brain. Each of them are introduced below and discussed through next chapter within each represented method.

Important Issues Regarding Instantaneous Phase Estimation: Phase Slipping Problem

Phase Slipping is a notorious problem common in phase estimations from cerebral signals which makes the phase sequence to contain fake (unrelated to brain activity) jumps during low foreground (narrow-band filtered signal) SNRs. This problem is rigorously studied in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF][START_REF] Seraj | Robust Electroencephalogram Phase Estimation with Applications in Brain-computer Interface Systems[END_REF]. In [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF], a Monte Carlo based statistical framework is proposed for EEG phase and frequency estimation to overcome the problems associated to previous deterministic approaches (introduced in section 1.1.0.3). Also, in [START_REF] Seraj | Robust Electroencephalogram Phase Estimation with Applications in Brain-computer Interface Systems[END_REF], despite the basic illustration of this problem, a straightforward approach for EEG phase estimation based on findings in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF] is proposed and is evaluated in a real BCI problem. To help other researches use the findings of these studies, a comprehensive cerebral signal phase analysis toolbox is generated by the authors of [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF] and [START_REF] Seraj | Robust Electroencephalogram Phase Estimation with Applications in Brain-computer Interface Systems[END_REF]. Please refer to [START_REF] Seraj | Cerebral Signal Instantaneous Parameters Estimation MATLAB Toolbox -User Guide Version 2[END_REF] for a detailed description of the provided toolbox

Volume Conduction

One of the most important obstacles confronting in investigating correlation and synchronization between different brain regions using EEG/MEG recordings is a phenomenon known as Volume Conduction (VC). The activity of a particular brain area will be observed in some adjacent electrodes because of conduction in brain volume which produces spurious synchrony. This substantial problem provides false results and shows incorrect related brain regions that can not be neglected. The problem of VC is especially large for scalp EEG and MEG data, because of their low spatial resolution. The difference between estimated connectivity directly from neural sources and the corresponding estimates from scalp recordings is an effect of VC [START_REF] Marshall | Statistical Analysis of EEG Phase Shift Events[END_REF].

Several recent methods are presented to overcome this substantial problem and as some of them reported, they obtained successful results. Many of these methods are represented in Section 2.

Presence of a Common Reference

EEG is a bipolar signal recorded from the scalp while MEG is a unipolar record. In a unipolar record, the recording channels are reference free whereas in a bipolar record, channels are constructed by subtracting two unipolar signals captured from two distinct scalp electrodes (i.e. a reference and an arbitrary electrode).

As reported in [START_REF] Guevara | Phase Synchronization Measurements Using Electroencephalographic Recordings: What Can We Really Say About Neuronal Synchrony?[END_REF], the confronting problem here is that bipolar recordings impose a distortion in observed synchrony values, such that in general, they destroy the intended physical interpretation of phase synchrony. It is notable to emphasize that the problem of a common reference is an obstacle when using EEG (and not MEG) recording due to their intrinsic bipolarity. The limitations due to this problem in coherence measurements [START_REF] Fein | Common reference coherence data are confounded by power and phase effects[END_REF], [START_REF] Zaveri | The effect of a scalp reference signal on coherence measurements of intracranial electroencephalograms[END_REF], and phase synchrony [START_REF] Guevara | Phase Synchronization Measurements Using Electroencephalographic Recordings: What Can We Really Say About Neuronal Synchrony?[END_REF], are widely investigated previously. It is reported that any contamination in the reference channel will significantly affect the coherence measurements [START_REF] Fein | Common reference coherence data are confounded by power and phase effects[END_REF], [START_REF] Zaveri | The effect of a scalp reference signal on coherence measurements of intracranial electroencephalograms[END_REF]. As mentioned before, one of the superb benefits of using PS over conventional coherence is that PS separates the amplitude information from the phase. Nevertheless, problem of common reference imposes limitations to PS too. It has been shown that if two signals are synchronized, their differences from a third signal are not necessarily synchronized. Based on this, the choice of the reference electrode can considerably affect the synchrony values obtained, up to the point that they almost span the entire interval [0 1] (see [START_REF] Guevara | Phase Synchronization Measurements Using Electroencephalographic Recordings: What Can We Really Say About Neuronal Synchrony?[END_REF]).

There are a number of different ways to choose the reference electrode such as typical EEG reference montages, which all voltage differences refer to one specific channel, bipolar montage, which reference electrode for every channel is different and normally spatially close to it, average reference and the Laplacian reference which is a mathematical approximation to a reference-free signal. All these techniques are susceptible to common reference problem but it is preferred to use the Laplacian reference method [START_REF] Guevara | Phase Synchronization Measurements Using Electroencephalographic Recordings: What Can We Really Say About Neuronal Synchrony?[END_REF].

Also as reported in [START_REF] Guevara | Phase Synchronization Measurements Using Electroencephalographic Recordings: What Can We Really Say About Neuronal Synchrony?[END_REF], the significance-level in PS measurements will not help and again high synchrony is detected when the reference electrode amplitude is relatively high. As a conclusion, utilizing MEG data in PS measurements and more generally in brain connectivity studies will produce much more reliable results.

Noise Sources and Spurious Synchrony

A noise source or an artifact in cerebral signal is referred to any unwanted contribution to the field potential. This might vary due to different applications. For example, when studying ERPs the spontaneous background EEG/MEG is considered as an artifact while in many other applications the background EEG/MEG is the subject of study. Also there are some common noise sources in EEG recordings that require a preprocessing step where the artifacts are detected and removed. Artifacts generated by the electrical equipment which are near the recording area, contaminations caused bay muscle activity (i.e. blinking, jaw clenching or frowning) and also other contaminant biological signals (i.e. EOG) are such artifacts. Neglecting these kind of noise sources which are uncorrelated to the signals under study in phase synchronization and brain connectivity measurements might impose incorrect PS values and cause Spurious Synchrony.

Phase Enslaving

"Narrow-band filtering tends to transfer phase behavior consistent with the target frequency (or frequency band) from the time region where this holds to adjacent regions where corresponding oscillations are lacking. Conversely, if the signal exhibits distinct oscillations of some particular frequency f , these may strike through to the phase course even if f is incompatible with the target frequency. Such tends to happen particularly in time-frequency regions where power is low. Since power in EEG signals decays rapidly with increasing frequency, phase enslaving may thus especially affect higher frequency (i.e. gamma) oscillations" [START_REF] Ehm | Detecting Frequency-specific Phase Synchrony in Brain Electrical Oscillations[END_REF].

Applications

Phase analysis has a broad usage in many fields. It has been used widely in communication investigations (i.e. signal phase modulation/demodulation) and also in biomedical researches, particularly in field of neuroscience (i.e. brain functional integration and cerebral signal synchrony measurements).

The most frequent applications in brain connectivity measurements are typically the Brain-Computer Interface (BCI), Event Related Potential (ERP) and the Cognitive studies. Recent applications include investigation of cognitive processes such as visual perception, mental rotation [START_REF] Bhattacharya | EEG gamma-band phase synchronization between posterior and frontal cortex during mental rotation in humans[END_REF], or pathological states such as epileptic seizures [START_REF] Mormann | Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients[END_REF], schizophrenia [START_REF] Bob | EEG phase synchronization in patients with paranoid schizophrenia[END_REF], [START_REF] Spencer | Abnormal neural synchrony in schizophrenia[END_REF], attention deficit (hyperactivity) disorder [START_REF] Tcheslavski | Phase synchrony and coherence analyses of EEG as tools to discriminate between children with and without attention deficit disorders[END_REF] or migraine [START_REF] Angelini | Steady-state visual evoked potentials and phase synchronization in migraine patients[END_REF].

Summary

In this chapter an overview of frequent cerebral signals, their history and their early and modern analysis methods, particularly in cerebral connectivity measurements were presented.

Most of the early quantities gained by digitalizing the EEG records were introduced. Quantities such as PSD, Correlation Function, signals' Instantaneous Phase and Amplitude and also the concept of Coherence. In modern analysis the main problem discussed was estimating functional integration between neural regions from spatio-temporal patterns in cerebral signals. To this purpose, the concept of Synchronization was introduced. It had been discussed that synchronization can be classified into three main categories: (1) linear, (2) nonlinear and (3) information-based methods. Conventional approaches in these categories were also represented. Finally, the major obstacles confronting in this new fast growing field were discussed and many possible and commonly used applications to this concept were introduced.

The next Chapter is dedicated to delicately representing many introduced Phase Analysis methods in various distinct fields which are applicable to brain studies.

Chapter 2 Phase Analysis Methods

Phase information is generally considered to be "purer"and less contaminated or more informative than the amplitude of cerebral recordings, which is more influenced by the impedance of the skull or various artifacts such as eyes or face muscles movements [START_REF] Galka | Topics in Nonlinear Time Series Analysis[END_REF]. Biological signal phase is an informative useful alternative feature in aforementioned cases that has been of interest in many different research fields through past few years. Thus, using phase of brain signals in cases that amplitude informations are contaminated or inadequate, or in cases that the phase itself has valuable adjunct information over amplitude (i.e. brain cognitive response investigations), would be very helpful.

As mentioned before, measuring brain connectivity is of a particular interest in neuroscience. Because of the reentry mechanism, phase variables and particularly their synchronization for the purpose of measuring connectivity in cerebral signals became more interesting [START_REF] Marshall | Statistical Analysis of EEG Phase Shift Events[END_REF].

Various methods within different fields and frameworks are presented in literature for the purpose of cerebral signal phase analysis. The main frameworks and basic concepts were introduced in previous Chapter. In this Chapter, we try to represent and elaborate these various methods within their corresponding field and framework.

Coherence

Conventional MSC

As introduced earlier in Section (1.1.0.4.1), the conventional method presented in [START_REF] Carter | Estimation of the Magnitude-Squared Coherence Function Via Overlapped Fast Fourier Transform Processing[END_REF] for measuring coherence is based on weighted windowing of the Fourier transform of signals. Let x(t) and y(t) be two random process, zero mean, wide-sense stationary and ergodic signals with length l. A summary of the method in [START_REF] Carter | Estimation of the Magnitude-Squared Coherence Function Via Overlapped Fast Fourier Transform Processing[END_REF] is given below:

• Divide the signals into N equal parts with T samples (different part can be overlapped or disjoint)

• Multiply the samples of each ensemble in a weighting function (i.e. cosine)

• Take FFT of each weighted ensemble

• Measure power spectral densities of signals (PSD)

• Finally, measure the coherence function (or as said "Mean Squared Coherence (MSC)")

If X i (f ) and Y i (f ) be the Fourier transform of i-th ensemble, the PSDs are calculated as follow:

P SD xx (f ) = 1 NT N i=1 |X i (f )| 2 (2.1) P SD yy (f ) = 1 NT N i=1 |Y i (f )| 2 (2.
2)

P SD xy (f ) = 1 NT N i=1 X i (f )Y * i (f ) (2.3)
where the signals are sampled at f s (f s > 2BW ) and the BW is the band width. N is the number of ensembles with 50% overlap, T is the number of samples in ensembles and a p -pointF F T is taken (p = T f s ).

With the ergodicity of signals being neglected we have:

P SD xx (f ) = E{|X i (f )| 2 } (2.4) P SD yy (f ) = E{|Y i (f )| 2 } (2.5) P SD xy (f ) = E{X i (f )Y * i (f )} (2.6)
where E{.} illustrates the mathematical expectation. Using these PSDs we compute the mean squared coherence MSC as:

|MSC(f )| 2 = |P SD xy (f )| 2 P SD xx (f )P SD yy (f ) = | N i=1 X i (f )Y * i (f )| 2 N i=1 |X i (f )| 2 N i=1 |Y i (f )| 2 (2.7)
and the complex coherence function is obtained as below:

MSC(f ) = P SD xy (f ) P SD xx (f )P SD yy (f ) (2.8)
It is notable that T have to be chosen large enough to reduce the bias and standard deviation of the measurements. Also the Fourier transform of the weight function have to be narrow. Thus a good resolution for PSD measurements is provided when N and T both are chosen large (the best resolution is obtained in 5o% overlap between ensembles).

The important question confronting here is that "what is the optimum way for averaging the MSC over ensembles? ". As reported in [START_REF] Brushe | On the Computation of an Averaged Coherence Function[END_REF], the procedure described above for computing averaged MSC might prevent some coherencies in case of non-stationary phase from being detected. Another method presented for averaging MSC over trials is described below:

|MSC(f )| 2 = 1 N N i=1 |(P SD xy (f )) i | 2 P SD xx (f )P SD yy (f ) (2.9)

MSC for Deterministic Evoked Signals

In some applications, for example detecting brain evoked potentials during a rhythmic periodic stimuli the MSC could be used to measure the correlation and synchronization between evoked potentials and background EEG signal. For Instance, in [START_REF] Mauricio | A Statistical Technique for Measuring Synchronism Between Cortical Regions in the EEG During Rhythmic Stimulation[END_REF] they assumed x(t) to be the evoked potential in response to a rhythmic periodic stimuli and y(t) as the background EEG. If the windowing process is done in the way that all the ensembles contain equal number of the periods of signal x(t) thus a modification as below is possible. Due to main equation presented for MSC in Section (1.2.1.4) we have:

|MSC(f )| 2 = | N i=1 Y i (f )X * i (f )| 2 N i=1 |X i (f )| 2 N i=1 |Y i (f )| 2 since x(t)
is deterministic, it can be put out from the summation as below:

|MSC(f )| 2 | x(t)periodic = | N i=1 Y i (f )| 2 N N i=1 |Y i (f )| 2 (2.10)
where 0 < MSC(f ) < 1, and it reaches to its maximum value (i.e. 1), when there is only the evoked potentials observing in the windows.

Multiple MSC

The extension of coherence concept and increasing in number of existing signals, particularly in case of multi-channel EEG recordings presents another aspect of coherence known as multiple coherence.

Multiple coherence shows the relation and synchronization between a signal and a group of other signals. This concept is used to increase the evoked potential detection rate for multi-channel EEG recordings in [START_REF] Mauricio | A Matrix Based Algorithm for Estimating Multiple Coherence of a Periodic Signal and Its Application to the Multichannel EEG During Sensory Stimulation[END_REF]. Similar to previous method, signal x(t) assumed to be the brain response of a rhythmic periodic stimuli (i.e. evoked potential) and so deterministic and y i (t)s(i = 1, 2, ..., L) as the different EEG signals recorded from different electrodes. The coherence function is computed as below:

|MSC(f )| 2 | x:y 1 ,..y L = P SD H yx (f )P SD -1 yy (f )P SD yx (f ) P SD xx (f ) (2.11)
where H indicates Hermitian of a matrix and the cross-spectrum vector and auto-spectrum matrix are as follow:

P SD yx (f ) =       P SD y 1 x (f ) . . . P SD y L x (f )      
(2.12)

P SD yy (f ) =     P SD y 1 y 1 (f ) . . P SD y 1 y L (f ) . . . . . . . . P SD y L y 1 (f ) . . P SD y L y L (f )     (2.13)
With x(t) being deterministic, same as previous method, its Fourier transform has an identical value in all ensembles and a simplification could be made:

P SD yx (f )| x(t)periodic =       N i=1 Y * 1i (f )X(f ) . . . N i=1 Y * Li (f )X(f )      
(2.14)

P SD yx (f )| x(t)periodic = X(f )       N i=1 Y * 1i (f ) . . . N i=1 Y * Li (f )       P SD yx (f )| x(t)periodic = X(f )V (f ) (2.15)
Putting this result into equation ( 25) results:

|MSC(f )| 2 | x(t)periodic = V H (f )P SD -1 yy (f )V (f ) N (2.16)
This procedure for measuring MSC is helpful calculating phase synchronization in multi-channel EEG recordings and particularly in application of detecting evoked potentials among a group of other channels of background records. As an example of recent studies using MSC for investigating cognitive sources of specific EEG based features within brain, we can mention [START_REF] Seraj | An Investigation on the Utility and Reliability of Electroencephalogram Phase Signal upon Interpreting Cognitive Responses in the Brain: A Critical Discussion[END_REF].

MSC Based on Thomson Multitapers

A novel method based on Thomson Multitapers [START_REF] Thomson | Spectrum estimation and harmonic analysis[END_REF], for determining coherence function is presented in [START_REF] Hansson-Sandsten | Cross-Spectrum And Coherence Function Estimation Using Time-Delayed Thomson Multitapers[END_REF]. In this method, the left and right singular vectors of the cross-covariance matrix corresponding to a known cross-spectrum model are evaluated as multitapers for estimation of cross-spectrum and coherence. A random process, zero mean, stationary signal x(t), is filtered through a causal linear filter h(n) and a stationary sequence y(t) is created and then the MSC function is measured using a multitaper spectrum estimator as:

|MSC(f ) mt | 2 = |P SD mt xy (f )| 2 P SD mt xx (f )P SD mt yy (f ) (2.17)
where -0.5 < f < 0.5 and:

P SD mt xy (f ) = 1 N N i=1 u H i F H (f )XY H F (f )v i (2.18) P SD mt xx (f ) = 1 N N i=1 u H i F H (f )XX H F (f )u i (2.19) P SD mt yy (f ) = 1 N N i=1 v H i F H (f )Y Y H F (f )v i (2.20)
where n samples of signals, x = [x(0), ..., x(n -1)] T and y = [y(0), ..., y(n -1)] T are used, and H is showing conjugate transpose and the left and right windows are

u i = [u i (0), ..., u i (n -1)] T and v i = [v i (0), ..., v i (n -1)
] T respectively and F = diag[1, e -j2πf , ..., e -j2π(n-1)f ] is the Fourier transform matrix.

According to [START_REF] Hansson-Sandsten | Cross-Spectrum And Coherence Function Estimation Using Time-Delayed Thomson Multitapers[END_REF], using this procedure for measuring MSC has some advantages over former methods (i.e. based on weighted overlapped segment averaging) in many applications.

Imaginary Part of Coherence

A new method based on imaginary part of the complex coherence function is presented in [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF] and it has been shown that the results of this method can not be caused by Volume Conduction (VC). This method assumes that there is no time lag between the potential recorded on the scalp and the source of this potential. This assumption is proved through an investigation in [START_REF] Stinstra | The volume conductor may act as a temporal filter on the ECG and EEG[END_REF]. During an accurate study, they do not observed any time lag or any shift in phase of the recorded potential compared with source of the corresponding signal for frequencies under 100Hz.

A very striking point discussed in [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF] is that they argued about using coherency for non-stationary signals and reported that the coherence analysis is feasible in non-stationary processes as strong as stationary ones.

The basic idea used in method presented by [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF] is using the real and imaginary part of the coherency instead of its amplitude and phase. Since the real and imaginary parts of coherency are just another representation of complex coherency function, no new quantities are being calculated and it is just a new perspective.

According to [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF], the (complex) coherency of the non-interacting areas of brain are necessarily real. Therefore the imaginary part of the coherency can be a good choice for investigation on interacting areas of brain. To see this, assume that the M possible active sources at any time are given as follow:

S = [S 1 , S 2 , ..., S M ] T and the L scalp electrodes as:

X = [X 1 , X 2 , ..., S L ] T now we have:       x 1 (t) . . . x L (t)       =     a 11 . . a 1M . . . . . . . . a L1 . . a LM           s 1 (t) . . . s L (t)       (2.21)
where a lm denote the degree to which source s m project to electrode x l . With the coefficient a lm assume to be constant, we have in frequency domain:

      X 1 (f ) . . . X L (f )       =     a 11 . . a 1M . . . . . . . . a L1 . . a LM           S 1 (f ) . . . S L (f )       (2.22)
Now assume that the signals recorded from two adjacent electrodes, x(t) and y(t) are superposed of K independent sources, then according to above equations we have:

X(f ) = K i=1 a i S i (f ) (2.23) Y (f ) = K j=1 b j S j (f ) (2.24)
Now based on the main assumption of this method, while computing the cross-spectrum of these signals we will have:

P SD xy (f ) = E{X(f )Y * (f )} = i j a i b j E{S i (f )S * j (f )} (2.25) ⇒ P SD xy (f ) = i a i b i E{S i (f )S * i (f )} = i a i b i |S i (f )| 2 (2.26)
which is real. Thus the normalization to auto-spectrum terms will also be real and so will be the coherence function. Based on this result and the assumption, if the VC causes no time lag (which has been shown in [START_REF] Stinstra | The volume conductor may act as a temporal filter on the ECG and EEG[END_REF]), the imaginary part of the coherency has no effect of VC artifact in it. This also has been shown in [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF] by adding non-interacting signal (i.e. noise) and experiencing a decrease in imaginary part of coherency.

In fact, the imaginary part of coherency shows the phase synchronization between time-lagged signals. This claim is coming from the fact that time-lag in time domain causes a phase shift in frequency domain and vise versa. Therefore, it could be said that the imaginary part of coherency specifically shows the synchronization between time-lagged signals.

According to [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF], however the amplitude and phase of the complex coherence function provide similar information to its real and imaginary part, but looking to the imaginary part instead of phase have some advantages as below:

• Non-interacting sources do not lead to small but rather to random phases. We cannot interpret a phase without having an estimate of its significance at the same time.

• One usually calculates coherency with respect to a baseline (a rest condition). Since in the individual coherencies the real parts are typically much larger than the imaginary parts, the phase flips by π depending on whether the real part of coherency is larger in the rest or active condition. The interesting structure is easily obscured by this rather meaningless effect.

• Phase is usually regarded as an additional information about time delay between two processes. However, volume conduction strongly affects the real part but does not create an imaginary part. Processes can appear to be synchronized with almost vanishing time delay while it is only the volume conducted copies of the signals which do not have a time delay.

Finally, stochastic given by [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF] to measure the imaginary part of coherency is presented in follow. It has been shown that for two signals x(t) and y(t), the imaginary part of power spectral density in determined as below:

Imag{X(f )Y * (f )} = K i=1 K j=1 A i A j (a i b j -a j b i )sin(∆θ ij ) (2.27)
where i and j indicate the independent sources and finally the imaginary part of coherency (ImagC) is obtained as:

ImagC xy (f ) = E{Imag{X(f )Y * (f )} E{A x } 2 E{A y } 2 = E{Imag{P SD xy (f )} E{A x } 2 E{A y } 2 (2.28)
which varies between 0 and 1 (similar to MSC(f )).

Phase Lag Index

The ImagC normalizes the components of cross-spectrum to the amplitudes of signals. Thus it is not true if one says it is only related to the imaginary components of cross-spectrum. Adding uncorrelated noise sources causes the amplitude of cross-spectrum to increase while the numerator of ImagC fraction stays unaffected. This problem is investigated in [START_REF] Stam | Phase Lag Index: Assessment of Functional Connectivity From Multi Channel EEG and MEG With Diminished Bias From Common Sources[END_REF] and by simulations it has been shown that the normalization of ImagC to the amplitudes of signals increases its sensitivity to additional uncorrelated noise sources and changes in phase of coherency. Due to this, another method for measuring phase synchronization is represented in [START_REF] Stam | Phase Lag Index: Assessment of Functional Connectivity From Multi Channel EEG and MEG With Diminished Bias From Common Sources[END_REF] to confront these problems. As a potential improvement on the ImagC, [START_REF] Stam | Phase Lag Index: Assessment of Functional Connectivity From Multi Channel EEG and MEG With Diminished Bias From Common Sources[END_REF] proposed the phase lag index (PLI). The PLI estimates, for a particular frequency, to what extent the phase leads and lags between signals from two sensors are non-equiprobable, irrespective of the magnitude of the phase leads and lags. In simulations, the PLI performed better than the ImagC in detecting true changes in phase synchronization, and was less sensitive to the addition of volume conducted noise sources [START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF]. The main purpose of PLI was to solve two major problems (1) volume conduction and (2) common reference electrode.

"The central idea is to discard phase differences that center around 0 mod π. One way to realize this is to define an asymmetry index for the distribution of phase differences, when the distribution is centered around a phase difference of zero. If no phase coupling exists between two time series, then this distribution is expected to be flat. Any deviation from this flat distribution indicates phase synchronization. Asymmetry of the phase difference distribution means that the likelihood that the phase difference ∆φ will be in the interval -π < ∆φ < 0 is different from the likelihood that it will be in the interval 0 < ∆φ < π. This asymmetry implies the presence of a consistent, nonzero phase difference (lag) between the two time series. The existence of such a phase difference or time lag, however, cannot be explained by the influences of volume conduction from a single strong source or an active reference, since these influences are effectively instantaneous. The distribution is expected to be symmetric when it is flat (no coupling), or if the median phase difference is equal to or centers around a value of 0 mod π" [START_REF] Stam | Phase Lag Index: Assessment of Functional Connectivity From Multi Channel EEG and MEG With Diminished Bias From Common Sources[END_REF].

An index of the asymmetry of the phase difference distribution can be obtained from a time series of phase differences ∆φ t k , k = 1, ..., N in the following way:

P LI = |E{sign[∆φ t k ]}| (2.29)
Phase lag index measures phase differences as they occur on the unit circle by first thresholding each ∆φ using the signum function and then averaging over successive data points [START_REF] Stam | Phase Lag Index: Assessment of Functional Connectivity From Multi Channel EEG and MEG With Diminished Bias From Common Sources[END_REF]. The PLI ranges between 0 and 1. A PLI of zero indicates either no coupling or coupling with a phase difference centered around 0 mod π. A PLI of 1 indicates perfect phase locking at a value of ∆φ different from 0 mod π. The stronger this nonzero phase locking is, the larger PLI will be. Note that PLI does no longer indicate, which of the two signals is leading in phase. Whenever needed, however, this information can be easily recovered, for instance, by omitting the absolute value in above equation [START_REF] Stam | Phase Lag Index: Assessment of Functional Connectivity From Multi Channel EEG and MEG With Diminished Bias From Common Sources[END_REF].

Weighted Phase Lag Index

The PLI is said to has some defects. Recently, in a comprehensive study in [START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF], it has been shown that PLI has some problems with VC, noise sources and detecting changes in phase synchronization due to its discontinuity. It is said that in PLI a small perturbation turns phase lags into phase leads and vise versa. Based on this, another quantity for measuring phase synchronization known as Weighted Phase Lag Index (WPLI) is presented in [START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF]. It is said that WPLI is an improvement over both the conventional PLI and the ImagC. According to [START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF], the advantage of WPLI over PLI is that it weights the observed phase lags and leads by the magnitude of the imaginary components of cross-spectrum that causes less sensitivity to additional uncorrelated noise sources and more power of detecting phase synchronization statistically.

As mentioned before (see Section 1.4), in case of phase synchronization measurements four problems are well known:

• The presence of a common reference signal

• Volume Conduction

• The presence of noise sources • Sample-size bias VC, and in case of EEG (but not MEG) data, the use of a common reference, can spuriously inflate phase synchronization indices. The problem of VS is especially large for scalp EEG and MEG data, because of their low spatial resolution.

Due to these problems, to increase the detection rate of true changes in phase synchronization and also to decrease the effects of additional uncorrelated noise sources, coherency phase changes and the common reference problem, the represented statistics for WPLI are as follow:

W P LI(f ) = |E{Imag{X(f )Y * (f )}}| E{|Imag{X(f )Y * (f )}|} = |E{Imag{P SD xy (f )}}| E{|Imag{P SD xy (f )}|} (2.30) W P LI(f ) = |E{|Imag{P SD xy (f )}|sign{Imag{P SD xy (f )}}| E{|Imag{P SD xy (f )}|} (2.31)
As it is clear, unlike the ImagC, the WPLI uses the components of imaginary part of the crossspectrum to normalize and also unlike the conventional PLI, sign{Imag{P SD xy (f )}} is weighted by |Imag{P SD xy (f )}|. The WPLI also varies between 0 to 1 and find its maximum value in:

sign{Imag{P SD xy (f )}} = ±1 when: |E{|Imag{P SD xy (f )}|sign{Imag{P SD xy (f )}}| = E{|Imag{P SD xy (f )}|}
The performance of WPLI is evaluated in confront with the mentioned phase synchronization measurement difficulties in a wide simulation experiment in [START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF] and it has been reported that the WPLI performed much better in comparison with both PLI and ImagC. As an instance of recent studies using PLI for BCI application we can mention [START_REF] Seraj | Improved Detection Rate in Motor Imagery Based BCI Systems Using Combination of Robust Analytic Phase and Envelope Features[END_REF].

MSC Based on Parametric Spectral Estimation Methods

Each of these represented method for analyzing of signal phase exhibit some sensitivity to both changes in background power and the instantaneous effects caused by volume conduction. As an improvement to all these methods, the traditional spectral estimation based on Fourier Transform of the raw EEG data, can be replaced by a parametric estimation derived using the Yule-Walker AR approach. After this as before, we have to look for periods of phase synchronization [START_REF] Gordon | Comparing parametric and nonparametric methods for detecting phase synchronization in EEG[END_REF].

Let p be the model order, v i (t) be the noise term. For two assumptive signals x(t) and y(t) the AR model could be derived as:

x(t) = p q=1 b 11 (q)x(t -q) + p q=1 b 12 (q)y(t -q) + v 1 (t)
(2.32)

y(t) = p q=1 b 21 (q)x(t -q) + p q=1 b 22 (q)y(t -q) + v 2 (t) (2.33)
which in frequency domain turns into:

X(f ) = p q=1 B 11 (q)X(f )e -j2πf q + p q=1 B 12 (q)Y (f )e -j2πf q + V 1 (f ) (2.34) Y (f ) = p q=1 B 21 (q)X(f )e -j2πf q + p q=1 B 22 (q)Y (f )e -j2πf q + V 2 (f ) (2.35)
To use this model in cross-spectrum measurement and then coherence function estimation using EEG/MEG signals we can extend this to a multi-channel case (i.e. multi-channel EEG). Let x(t) = [x 1 (t), ..., x L (t)] T be the scalp electrodes. Derivation of Multy-channel AR Model gives (B(0) = -I, F.T.{v(t)} = N(f ) [START_REF] Gordon | Comparing parametric and nonparametric methods for detecting phase synchronization in EEG[END_REF]:

x(t) = p q=1 b(q)x(t -q) + v(t) (2.36)
taking Fourier transform of this equation gives:

X(f ) = p q=1 B(q)X(f )e -j2πf q + N(f ) (2.37) - p q=1 B(q)X(f )e -j2πf q + X(f ) = N(f ) (I - p q=1 B(q)e -j2πf q )X(f ) = N(f ) (- p q=0 B(q)e -j2πf q )X(f ) = N(f ) B(f )X(f ) = N(f ) (2.38)
where B(f ) = -p q=0 B(q)e -j2πf q . If H(f ) = B-1 (f ) multiplying in equation above gives the cross-spectrum as:

X(f ) = H(f )N(f ) X(f )X * (f ) = H(f )N(f )N * (f )H * (f ) (2.39)
which is an estimation of cross-spectrum needed for establishing the coherence function. Normalization of the temporal averaged version of this cross-spectrum (i.e. expectation) to the amplitude of the estimated AR coefficients gives the coherence function. To estimate the coefficients of AR model, as presented in [START_REF] Tcheslavski | Phase synchrony and coherence analyses of EEG as tools to discriminate between children with and without attention deficit disorder[END_REF], the least squares modified Yule-Walker equation (LSMYWE ) estimator [START_REF] Kay | Modern Spectral Estimation. Theory and Application[END_REF], can be used:

b = -( RH RH ) -1 RH r (2.40) R =       rxx,
where rxx,n shows the i-th element of auto-correlation estimation for signal x n , p is the AR model order and M = 2p.

As reported in [START_REF] Gordon | Comparing parametric and nonparametric methods for detecting phase synchronization in EEG[END_REF], the parametric and non-parametric spectral estimation methods perform similarly in low noise power while in high noise power the parametric methods act much better in detecting periods of phase synchronization. Using this parametric model in first step of the presented phase synchronization detection methods will improve the results of them and make them less impressible to noise, background power changing and volume conduction effects. The main problem confronting here is dealing with number of channels and the optimum model order.

It is notable that using non-parametric methods have some important advantages too. For instance, parametric methods assume that signals are generated by a linear auto-regressive process but the non-parametric methods do not consider this assumption. Moreover, there are less parameters to be adjusted in non-parametric methods (i.e. the window length), but in parametric methods, after selecting the number of channels and adjusting the model order, the AR model coefficients have to be estimated and after that other parameters such as window length and etc. are set.

Phase Synchronization

As discussed before, in synchronization concept, coherence deals with several objections. The linearity, stationarity and also the lack of ability to separate the effects of amplitude and phase informations are the main abjections confronting. Thus, another concept called phase locking was introduced to overcome these obstacles.

Phase Locking Value

There are two major approaches to calculate this quantity. Both of these methods are based on the framework presented in Section (1.3.2.1). It means in first step they utilize one of the phase extraction methods and then after simply calculating the phase-differences, they quantify the local stability in this phase-differences and finally, in significance-level the significance of captured values are determined. The major difference between methods presented for Phase Locking Value (PLV) is in the manner they use to deal with the first and third steps. Due to different techniques utilized in these steps, we represent various statistics for PLV as follow.

Phase-Locking Statistics

This method, called Phase-Locking Statistics (PLS), measures the significance of the phase-covariance between two signals. PLS separates the phase and amplitude components and can be directly interpreted in the framework of neural integration [START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF].

To have a meaningful interpretation from instantaneous phase, it has to be extracted through a narrow filter which isolates the frequency of interest. In this method, First of all the signal is band-pass filtered in its frequency of interest and a short arbitrary interval around it (i.e. f ± 2). After that, its convolution with a complex Gabor-Wavelet centered at frequency f is calculated and phase of this convolution is extracted for all time-bins t, trials n = 1, ..., N, and for each of the pair of electrodes as discussed before in Section (1.2.1.3.2). As mentioned before, when using Wavelet Transform description below gives back the phase-differences which is calculated across trials, as in the second step [START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF], [START_REF] Lachaux | Studying single-trials of phase-synchronous activity in the brain[END_REF]:

exp(j(φ y (t, f ) -φ x (t, f ))) = W x (t, f )W * y (t, f ) |W x (t, f )||W y (t, f )| (2.43)
where the term φ y (t, f ) -φ x (t, f ) represents the phase-difference between processes y(t) and x(t).

As the third step, the stability of phase-differences across trials is quantified by a phase locking value:

P LV (t, f ) = | 1 N N i=1 exp(j(φ y,i (t, f ) -φ x,i (t, f )))| (2.44)
where N is the total number of trials and i indicates the index of trial. PLV varies between 0 and 1 which indicate the completely non-synchronized and synchronized signals respectively. PLV is said to has phase enslaving problem when the power is low.

As the final step in significance-level, the degree of statistical significance of each PLV value is determined by comparing it to the values obtained by the surrogate data which is a statistical test with the null hypothesis that states: "the value of original data is not significantly different from the specific class of signals" [START_REF] Mezeiova | Comparison of coherence and phase synchronization of the human sleep electroencephalogram[END_REF]. Efficient surrogate data are generated in a way that they simulate all signals' characteristics, but destroy only the one property suspected to be the source of outstanding values of measure computed [START_REF] Mezeiova | Comparison of coherence and phase synchronization of the human sleep electroencephalogram[END_REF]. For this purpose, as in [START_REF] Lachaux | Studying single-trials of phase-synchronous activity in the brain[END_REF], the surrogate data is generated by shifting the trials of one of the main signals (i.e. y(t)) in such a way that the phase-differences are no longer computed in simultaneously recorded trials. It means we have:

P LV surrogate (t, f ) = 1 K K j=1 | 1 N N i=1 exp(j(φ y,j(i) (t, f ) -φ x,i (t, f )))| (2.45)
The proportion of surrogate values greater than the original PLV for a tim t is called PLS. Normally a criterion of 5%(P LS < 5%) is used to characterize significant synchrony [START_REF] Lachaux | Studying single-trials of phase-synchronous activity in the brain[END_REF]. As direct applications of PLV in recent studies we can mention [START_REF] Seraj | Robust Electroencephalogram Phase Estimation with Applications in Brain-computer Interface Systems[END_REF][START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF][START_REF] Karimzadeh | A distributed classification procedure for automatic sleep stage scoring based on instantaneous electroencephalogram phase and envelope features[END_REF][START_REF] Seraj | An Investigation on the Utility and Reliability of Electroencephalogram Phase Signal upon Interpreting Cognitive Responses in the Brain: A Critical Discussion[END_REF] and [START_REF] Seraj | Improved Detection Rate in Motor Imagery Based BCI Systems Using Combination of Robust Analytic Phase and Envelope Features[END_REF]. Further implementation details and MATLAB codes related to these studies can be found in [START_REF] Seraj | Cerebral Signal Instantaneous Parameters Estimation MATLAB Toolbox -User Guide Version 2[END_REF] and [START_REF] Sameni | The Open-Source Electrophysiological Toolbox (OSET)[END_REF].

Single-trial Phase Locking Statistics

As an extension of PLS, Single-trial Phase Locking Statistics (S-PLS) is presented in [START_REF] Lachaux | Studying single-trials of phase-synchronous activity in the brain[END_REF] in order to estimate PLS in single-trials. S-PLS shares the whole idea with PLS except at a slight cost of temporal resolution, in this method the phase-difference variability is calculated through successive time-steps (not across trials). Based on this, the S-PLV for each individual trial is obtained as below:

S -P LV (t, f ) = | 1 δ t+ δ 2 t-δ 2 exp(j(φ y (τ, f ) -φ x (τ, f )))dτ | (2.46)
As PLV, the S-PLV varies between 0 and 1 with 1 indicating the strongest phase locking [START_REF] Lachaux | Studying single-trials of phase-synchronous activity in the brain[END_REF].

The significance-level in this method is performed similar to that of PLS except that here surrogate values are generated by estimating the maximum PLV between pairs of independent white noise signals. Finally, S-PLS is obtained through a similar criteria to that of PLS.

Phase Synchrony Based on Shannon Entropy

The instantaneous phase sequence which is used to estimate the PLV can also be captured through the Hilbert transform (statistics are given in Section 1.2.1.3.1). In a method presented in [START_REF] Tass | Detection of n:m phase locking from noisy data: application to magnetoencephalography[END_REF], as the first step they extract the instantaneous phase of signals by Hilbert transform. A merit of the this approach is that the phase can be easily obtained for an arbitrary broad-band signal. Nevertheless, instantaneous amplitude and phase have a clear physical interpretation only if the signals are narrowband [START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF]. Therefore, an extra filtering pre-step is required to separate the frequency band of interest. After the pre-step filtering and extracting the phase sequence from the frequency band of interest, Tass et al. [START_REF] Tass | Detection of n:m phase locking from noisy data: application to magnetoencephalography[END_REF], proposed that Shannon Entropy can be used to characterize the statistical strength of phase synchrony. For this purpose, "the deviation of the actual distribution of the phase-difference between recording signals from a uniform one must be quantified" [START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF]. The synchronization index proposed by Tass et al. in [START_REF] Tass | Detection of n:m phase locking from noisy data: application to magnetoencephalography[END_REF] can be illustrated as follow:

γ = H max -H H max (2.47)
where H is the entropy and defined by (see Section 1.3.3):

H = - N i=1 P i ln P i (2.48)
which N is the total number of trials, and H max = ln(N) is the maximal entropy, and P i the relative frequency of finding the phase-differences within the i -th trial [START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF].

As reported in [START_REF] Tass | Detection of n:m phase locking from noisy data: application to magnetoencephalography[END_REF], the optimal number of trials can be obtained as:

N = exp(0.626 + 0.4 ln(n -1)) (2.49)
where n is the number of samples. Finally, the degree of statistical significance of the phase-locking values was determined through generating surrogate values obtained from surrogate data. Surrogate data is constructed by scrambling the original series. This randomization destroys any temporal structure, if present in the original series [START_REF] Tass | Detection of n:m phase locking from noisy data: application to magnetoencephalography[END_REF], [START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF]. The γ is a normalized quantity (0 < γ < 1) which measures the strength of phase synchronization and the γ = 1 indicates the perfect synchrony.

Phase Synchrony Based on Mutual Information

Another phase-locking statistic was introduced in [START_REF] Palus | Detecting phase synchronization in noisy systems[END_REF]. This method is also similar to previous methods with the whole idea except in procedure presented for quantification the phase synchrony strength. The phase synchrony index in this method is as below:

ρ = I I max (2.50)
where I max=ln(N ) and N is total number of trials. I is the mutual information between phase sequences and can be determined as follow:

I(φ x , φ y ) = π -π π -π
P x,y (φ x , φ y ) log( P x,y (φ x , φ y ) P x (φ x )P y (φ y ) )dφ x dφ y (2.51) which P x (φ x ) and P x (φ y ) are the probability distribution of the phases φ x and φ y , and P x,y (φ x , φ y ) is their joint distribution. The Mutual Information (MI) in this method measures the dependence between phase sequences of signals x(t) and y(t).

It has been reported that MI produces surplus amount of information caused by an incorrect assumption that two processes are independent [START_REF] Shannon | The Mathematical Theory of Information[END_REF]. Also, as discussed in Section (1.3.3), MI can be rewritten as entropy as follow:

I(φ x , φ y ) = H(φ x ) + H(φ y ) -H(φ x , φ y ) (2.52)
Finally, the statistical significant of obtained phase synchrony values is tested through significantlevel similar to that of previous method. Normalized in this way, just like previous methods, ρ varies between 0 and 1 where latter indicates the complete synchrony [START_REF] Palus | Detecting phase synchronization in noisy systems[END_REF].

It seems that this method does not provide any further information in comparison to previous method. Although, more accurate results will be presented in next Chapter after implementing these methods using similar data sets and in a unique situation.

Mean Phase Coherence

Another statistical measure for phase synchrony called as Mean Phase Coherence (MPC) is presented in [START_REF] Mormann | Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients[END_REF]. This method utilize circular variance [START_REF] Mardia | Probability and Mathematical Statistics: Statistics of Directional Data[END_REF] to characterize the dependences between instantaneous phases of two signals. The method presented in [START_REF] Mormann | Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients[END_REF] shares the two first steps with previous methods where as first step the instantaneous phases are extracted through Hilbert transform and as second step the phase-differences are calculated by simply subtracting the phase sequences. Finally, as a measure for synchronization the MPC of an angular distribution is defined as:

MP C = | 1 N N -1 i=0 exp(jφ x,y (i∆t))| (2.53)
where φ x,y indicates the phase differences between phases φ x and φ y and 1 ∆t is the sampling rate of discrete series. Whit utilizing the Euler's formula, the above equation turns into:

MP C = ([ 1 N N -1 i=0 cos(φ x,y (i∆t))] 2 + [ 1 N N -1 i=0 sin(φ x,y (i∆t))] 2 ) 1 2
(2.54)

Clearly from this notation MP C is restricted to the interval [0, 1] [START_REF] Mormann | Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients[END_REF]. As same as all previous methods presented for phase synchrony, reliability of these obtained synchrony values is evaluated by the surrogate data approach. In [START_REF] Mezeiova | Comparison of coherence and phase synchronization of the human sleep electroencephalogram[END_REF] they offered to use Z-score statistics for the statistical test in significance-level. Z-score statistic calculates how different is the value of original data from the mean value of all surrogates in the units of standard deviation [START_REF] Mezeiova | Comparison of coherence and phase synchronization of the human sleep electroencephalogram[END_REF]:

Z = MP C org -µ(MP C surr ) σ(MP C surr ) (2.55)
where µ and σ are the mean value and the standard deviation of all surrogates respectively.

Wavelet Coherence

Another method presented in [START_REF] Lachaux | Studying single-trials of phase-synchronous activity in the brain[END_REF] called as Wavelet Coherence (WC) uses the wavelet coefficients of signals and the coherence function to provide a weighted version of PLV. The basic idea for presenting this method is that instantaneous phase may be meaningless where power is low and former methods are susceptible to a notorious problems called as phase enslaving [START_REF] Ehm | Detecting Frequency-specific Phase Synchrony in Brain Electrical Oscillations[END_REF]. AS in [START_REF] Lachaux | Studying single-trials of phase-synchronous activity in the brain[END_REF], WC for time delays t and frequency(bands) f can be obtained as below:

W C(t, f ) = | 1 N N i=1 W x i (t, f )W * y i (t, f )| 1 N N i=1 |W x i (t, f )| 2 1 N N i=1 |W y i (t, f )| 2 (2.56)
Since for both signals x(t) and y(t) we have:

exp(jφ(t, f )) = W (t, f ) |W (t, f )| (2.57)
and also from the original PLV equation:

P LV (t, f ) = | 1 N N i=1 exp(j(φ y,i (t, f ) -φ x,i (t, f )))| (2.58)
then we got:

W x i (t, f )W * y i (t, f ) = exp(j∆ i (t, f ))|W x i (t, f )||W y i (t, f )| (2.59)
then assuming:

λ i (t, f ) = |W x i (t, f )||W y i (t, f )| 1 N N i=1 |W x i (t, f )| 2 1 N N i=1 |W y i (t, f )| 2 (2.60)
thus WC can be rewritten as a weighted variant of PLV as below [START_REF] Ehm | Detecting Frequency-specific Phase Synchrony in Brain Electrical Oscillations[END_REF]:

W C(t, f ) = | 1 N N i=1
exp(j∆ i (t, f ))λ i (t, f )| (2.61) Due to multiplication by λ i (t, f ), wavelet amplitudes are low and this clearly reduces the potential consequences of phase enslaving. Nevertheless, WC has some defects such as it does not separate the effects of amplitude and phase in the interrelation between signals [START_REF] Ehm | Detecting Frequency-specific Phase Synchrony in Brain Electrical Oscillations[END_REF].

Weighted Phase Locking Value

The basic idea in this method presented in [START_REF] Ehm | Detecting Frequency-specific Phase Synchrony in Brain Electrical Oscillations[END_REF] is completely similar to WC method. The only different is between procedures for obtaining the weights. Description below known as Weighted Phase Locking Value (WPLV) is presented to measure the statistical strength of phase synchrony between signal pairs:

W P LV (t, f ) = | 1 N N i=1
exp(j∆ i (t, f ))ω x i (t, f )ω y i (t, f )| (2.62)

In contrast to WC, here weights are not proportional to amplitudes. As offered in [START_REF] Ehm | Detecting Frequency-specific Phase Synchrony in Brain Electrical Oscillations[END_REF], the weights differ from 1 -which gains the conventional PLV-only when the amplitudes are considered "not large enough"and that depends on the respective spectral characteristics. For calculating the weights assume that for signal x(t) we have (the procedure and statistics for y(t) is repeated):

P x (t, f ) = 1 f 2 -f 1 f 2 f 1 |W x (t, f )| 2 df
(2.63)

P x (f ) = 1 T T P x (t, f )dt (2.64)
where P x (t, f ) denotes the mean wavelet power in frequency band f at time t and P x (f ) is its time average across the relevant trial T . Now, assessment of weights can be summarized as below [START_REF] Ehm | Detecting Frequency-specific Phase Synchrony in Brain Electrical Oscillations[END_REF]:

• Large enough: ω x (t, f ) = 1, if P x (t, f ) > 3 2 P x (f ) • Too small: ω x (t, f ) = 0, if P x (t, f ) ≤ 3 10 P x (f )

• Between: ω x (t, f ) = log( Px (t,f ) Px(f ) )-log( 3 10 ) log [START_REF] Davis | Electrical reactions of the human brain to auditory stimulation during sleep[END_REF] For both the previous and current methods, since the wavelet transform has to be computed anyway, it is preferred to use it also for assessing the instantaneous phases. As discussed before in previous Chapter, the wavelet phases for frequency bands can be determined as the phase angle of the band-averaged complex wavelet coefficients as follow:

exp(jφ x (t, f )) = W x (t, f ) |W x (t, f )| (2.65)
where

W x (t, f ) = 1 f 2 -f 1 f 2 f 1 W x (t, f ) (2.66)
It is notable that phase enslaving is less serious problem in case of cerebral signals and it is normally serious when test signals with a sharply defined spectrum is used. Nevertheless, it is expected that WPLV yield more accurate results than the conventional PLV because of discussed problem.

Reduced Interference Rihaczek Distribution Phase Synchrony

This phase synchronization method is based on Cohen's Class of Distributions and the Complex Energy Density Function (CEDF) phase extraction method. Cohen's class of distributions are bilinear time-frequency distributions and can be illustrated as below [START_REF] Cohen | Time-Frequency Analysis: Theory and Applications[END_REF]:

CEDF (t, f ) = ∞ -∞ ∞ -∞
φ(θ, τ )A(θ, τ ) exp(-j(θt + 2πτ f ))dτ dθ (2.67) where φ(θ, τ ) is the kernel function and A(θ, τ ) is the ambiguity function and can be obtained as below:

A(θ, τ ) = ∞ -∞ x(u + τ 2 )x * (u - τ 2
) exp(jθu) (2.68)

"The major differences between Cohen's class of distributions compared to other time-frequency representations such as the wavelet transform are the nonlinearity of the distribution, energy preservation, and the uniform resolution over entire time-frequency plane" [START_REF] Aviyente | A Phase Synchrony Measure for Quantifying Dynamic Functional Integration in the Brain[END_REF]. This distributions describe energy of a signal over time and frequency simultaneously but they can not be used for describing phase information of a signal. Thus, we need complex time-frequency distributions to obtain both the energy and the phase information. For this purpose, as discussed earlier in Section (1.2.1.3.3), Rihaczek derived the signal energy distribution in time and frequency by application of the complex signal notation. This CEDF, known as Rihaczek distribution, encounters a major disadvantage due to producing crossterms in case of multi-component signals. These crossterms are at the same time and frequency as the original signals and will lead to biased energy and phase estimates [START_REF] Aviyente | A Phase Synchrony Measure for Quantifying Dynamic Functional Integration in the Brain[END_REF]. Assume z(t) = x(t) + y(t), then the Rihaczek distribution is as follows: where the last two terms are crossterms. Therefore, a filtering is needed to remove these crossterms out of calculations. For this purpose, this method proposes a method based on Cohen's class of distributions and uses its kernel function as the required filter. This is called Reduced Interference Rihaczek Distribution (RID) and utilizes the Choi-Williams kernel to filter out the crossterms. As in [START_REF] Aviyente | A Phase Synchrony Measure for Quantifying Dynamic Functional Integration in the Brain[END_REF], the Reduced Interference Complex Energy Density Function or simply (RI-CEDF) can be expressed as:

CEDF (t, f ) = 1 √ 2π (x(
RI -CEDF (t, f ) = exp(- (θτ ) 2 σ
) exp(j θτ 2 )A(θ, τ ) exp(-j(θt + 2πτ f ))dτ dθ (2.70)

where exp(j θτ 2 ) is the Rihaczek distribution kernel function and the integrals are from -∞ to -∞. After defining this RI-CEDF, as introduced previously, the phase in time-frequency plane can easily be defined as:

φ x (t, f ) = arg( RI -CEDF (t, f ) |RI -CEDF (t, f )| ) (2.71)
Now, as the second step in phase synchronization evaluation process, we need to estimate the phase-difference between two signals x(t) and y(t). The phase-difference between time-varying phases can be computed as follow: Finally in this method, as the final step to assess phase synchrony, the conventional PLV measure is used to quantify the phase synchrony strength between pairs of electrodes.

Phase Synchrony Based on Empirical Mode Decomposition

Another set of methods for evaluating phase synchrony utilize Empirical Mode Decomposition (EMD) and separate data into Intrinsic Mode Functions (IMFs) and then, extract the instantaneous phases of the IMFs and finally use a phase synchrony quantification procedure, for example in [START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF], [START_REF] Sweeny-Reed | A novel approach to the detection of synchronisation in EEG based on empirical mode decomposition[END_REF] and [START_REF] Looney | Measuring Phase Synchrony Using Complex Extensions of EMD[END_REF].

Extracting signals phase sequences through EMD algorithm was discussed in previous Chapter. First, signals are decomposed into their relative IMFs and then the Hilbert transform is applied to these IMFs to extract phases. As the next step, for evaluating phase synchrony, the phase-differences have to be calculated. Assuming signals x(t) and y(t) and their corresponding phases φ x (t) and φ y (t) which are extracted using EMD, the phase-difference is simply calculated as: δφ x,y (t) = φ y (t) -φ x (t) (2.73) Finally, as in [START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF] and [START_REF] Sweeny-Reed | A novel approach to the detection of synchronisation in EEG based on empirical mode decomposition[END_REF], the phase synchrony strength is evaluated by this phase-difference through measuring the conventional PLV.

Phase synchrony measurements using EMD have some advantages over other introduced previous methods [START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF]. First, it is said that the discriminative band of the phase synchrony phenomenon differs not only with subjects, but also with time. Nevertheless, EMD algorithm can decompose the nonstationary signal into its relative IMFs. Therefore, no subject-specific or frequency-specific factors need to be considered to get the discriminative band [START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF]. Second, band-pass filter is a required necessary step in previous methods which causes time-delay in output signal with respect to the input signal while EMD based methods prevent this time-delay.

As reported in [START_REF] Looney | Measuring Phase Synchrony Using Complex Extensions of EMD[END_REF], there are two main obstacles confronting while using this method for phase synchrony purposes: (1) it is computationally inefficient (2) misleading results can be obtained as components of different frequencies can temporarily appear synchronized. To overcome these problems, [START_REF] Looney | Measuring Phase Synchrony Using Complex Extensions of EMD[END_REF] proposes to use the Complex Extension of EMD. For more information refer to [START_REF] Looney | Measuring Phase Synchrony Using Complex Extensions of EMD[END_REF].

  φ x,y (t, f ) = arg( RI -CEDF x (t, f )RI -CEDF * y (t, f ) |RI -CEDF x (t, f )||RI -CEDF y (t, f )| ) (2.72)

Appendix A Which Coordinates? since P SD xy (f ) is complex, the Cartesian or polar coordinates can be used for averaging. In [START_REF] Brushe | On the Computation of an Averaged Coherence Function[END_REF] this problem is investigated. Using Cartesian coordinate we got:

where P SD c xy (f ) indicates the cross-spectrum averaged in Cartesian coordinates and its amplitude and phase are obtaining as follow:

Similarly, in Polar coordinates we have:

where P SD p xy (f ) indicates the cross-spectrum averaged in Polar coordinates and ρ its amplitude and φ its phase are obtaining as follow:

It is clear that equations (A.3) and (A.6) and also equations (A.4) and (A.7) are equal when all a i (f )s are equal and all b i (f )s are equal too. In all any other conditions the amplitude and phase averaged in Polar coordinates are larger than the amplitude and phase averaged in Cartesian coordinates. Due to EEG signal non-stationarity this condition will not occur and thus as reported in [START_REF] Brushe | On the Computation of an Averaged Coherence Function[END_REF], the optimum averaged coherence function performance can be obtained through presentation and averaging the power spectral density functions in Polar coordinates.