Esmaeil Seraj
email: e.seraj@cse.shirazu.ac.ir

Reza Sameni
email: rsameni@shirazu.ac.ir

Robust Electroencephalogram Phase Estimation with Applications in Brain-Computer Interface Systems

Keywords: Electroencephalogram Phase, Electroencephalogram Frequency, Analytic Signal Representation, Electroencephalogram Hilbert Transform, Brain-Computer Interface

In this study, a robust method is developed for frequency-specific electroencephalogram (EEG) phase extraction using the analytic representation of the EEG. Based on recent theoretical findings in this area, it is shown that some of the phase variations-previously associated to the brain response-are systematic sideeffects of the methods used for EEG phase calculation, especially during low analytical amplitude segments of the EEG.

With this insight, the proposed method generates randomized ensembles of the EEG phase using minor perturbations in the zero-pole loci of narrow-band filters, followed by phase estimation using the signal's analytical form and ensemble averaging over the randomized ensembles to obtain a robust EEG phase and frequency. This Monte Carlo estimation method is shown to be very robust to noise and minor changes of the filter parameters and reduces the effect of fake EEG phase jumps, which do not have a cerebral origin.

As proof of concept, the proposed method is used for extracting EEG phase features for a brain computer interface (BCI) application. The results show significant improvement in classification rates using rather simple phase-related features and a standard K-nearest neighbors and random forest classifiers, over a standard BCI dataset. The average performance was improved between 4-7% (in absence of additive noise) and 8-12% (in presence of additive noise). The significance of these improvements was statistically confirmed by a paired sample t-test, with 0.01 and 0.03 p-values, respectively. The proposed method for EEG phase calculation is very generic and may be applied to other EEG phase-based studies.

Introduction

The phase analysis of electroencephalogram (EEG) signals has found great interest in the past decades. It has been shown to be a notable compliment for-and in cases more informative than-the EEG spectral amplitude. The EEG phase has been associated to neural propagation and synchronized firing of neuronal populations (neural assembly) during different mental and motor tasks. It is believed that a neural assembly corresponds to a group of functionally interconnected neurons within the brain, which interact with a reentry mechanism [START_REF] Gray | Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties[END_REF], to produce responses to specific intentions such as motor tasks. A reentry mechanism in terms requires an assembly to synchronize their firing rates [START_REF] Marshall | Statistical Analysis of EEG Phase Shift Events[END_REF]. These findings have led to an increased interest in utilizing features of cerebral phase signals, such as their synchronization.

In [START_REF] Ktonas | Instantaneous envelope and phase extraction from real signals: theory, implementation, and an application to EEG analysis[END_REF], various methods were presented for extracting the instantaneous EEG phase and amplitude of real EEG. The relationship between the different methods were studied using the theory of analytic signals. In [START_REF] Mormann | Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients[END_REF], the relation between phase synchronization in EEG signals and brain activity in patients with temporal lobe epilepsy was investigated. The authors reported a strong correlation between the shifts in phase synchrony and pathological activity. In [START_REF] Rudrauf | Frequency flows and the time-frequency dynamics of multivariate phase synchronization in brain signals[END_REF], the concept of frequency flows analysis (FFA) was introduced, as a new approach for studying the dynamics of phase synchrony in brain signals. The application of EEG phase and phase synchrony in brain recordings for brain-computer interface (BCI) systems was studied in [START_REF] Brunner | Phase relationships between different subdural electrode recordings in man[END_REF]. It was reported that additional information could be obtained by utilizing phase-related quantities for measuring brain synchrony. In [START_REF] Naruse | Statistical method for detecting phase shifts in alpha rhythm from human electroencephalogram data[END_REF] and [START_REF] Sauseng | Are event-related potential components generated by phase resetting of brain oscillations? A critical discussion[END_REF], two other important phase related quantities, namely phase shift and phase resetting and their relationship with event-related potentials were studied. 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], a phase synchronization method based on empirical mode decomposition (EMD) was proposed and evaluated for BCI applications. The phase-locking value (PLV) has been previously utilized to associate the EEG phase information in BCI systems [START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF][START_REF] Hamner | Phase-based features for motor imagery brain-computer interfaces[END_REF][START_REF] He | Combination of amplitude and phase features under a uniform framework with EMD in EEG-based Brain-Computer Interface[END_REF][START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] Carreiras | Phase-locking factor in a motor imagery brain-computer interface[END_REF]. More recently, several studies have utilized other EEG-based features for BCI applications, including steady-state visually evoked potential (SSVEP) [START_REF] Min | Decoding of top-down cognitive processing for SSVEP-controlled BMI[END_REF], event-related synchronization (ERS) and de-synchronization (ERD) [START_REF] Blankertz | The Berlin Brain-Computer Interface: Progress Beyond Communication and Control[END_REF][START_REF] Sturm | Interpretable Deep Neural Networks for Single-Trial EEG Classification[END_REF]. Herein, we focus on rather simple EEG phase features to show how a robust phase extraction mechanism can significantly improve the performance of BCI systems, which are based on phase synchronization and phase desynchronization [START_REF] Hamner | Phase-based features for motor imagery brain-computer interfaces[END_REF][START_REF] Carreiras | Phase-locking factor in a motor imagery brain-computer interface[END_REF].

In previous research, various methods such as wavelet transforms and analytic signal representations have been used for EEG phase extraction and PLV measurement [START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF][START_REF] Rosenblum | Phase synchronization of chaotic oscillators[END_REF]. In a recent study [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF], the authors presented a statistical framework for EEG phase analysis. Using an additive data model between the so-called background (spontaneous) and foreground EEG, probability density functions and other statistical properties of the instantaneous EEG envelope, phase and frequency were derived. It was analytically and numerically shown that in low analytical signal envelopes, the EEG phase is highly susceptible to the background EEG and noise. It was shown that although EEG phase variations convey important information regarding the EEG, some instantaneous phase jumps are merely systematic side effects of the processing stages used for EEG phase extraction, without any cerebral origin. A Monte Carlo approach was also proposed in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF] to detect and smooth the phase contents, during the time instants in which the estimated EEG phase is unreliable.

In this study, using our recent findings reported in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF], a new EEG phase extraction procedure is presented for extracting reliable phase sequences from the EEG during a BCI experiment. It is shown that the discrimination between true and fake EEG phase variations can significantly improve classification rates in BCI applications, even using rather basic features and conventional classifiers.

In Section 2, some preliminary backgrounds and limitations of classical EEG phase extraction techniques are reviewed. In Section 3, the proposed method and modifications to the conventional procedure of instantaneous EEG phase extraction are presented. As proof of concept, the proposed method is used for feature extraction in a BCI application in Section 4, showing significant improvement in the classification rates, followed by a discussion and future perspectives regarding the proposed scheme.

Background

The conventional phase estimation procedure

The conventional procedure for extracting the instantaneous phase sequence of a signal consists of two main stages: 1) narrow-band filtering, and 2) estimating the phase of the narrow-band signal [START_REF] Boashash | Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals[END_REF][START_REF] Picinbono | On instantaneous amplitude and phase of signals[END_REF].

For a unique and canonical definition, the instantaneous phase is extracted from very narrow frequency band signals [START_REF] Chavez | Towards a proper estimation of phase synchronization from time series[END_REF]. Moreover, the input signal's phase contents should not be affected by the filtering procedure. For the first stage, almost all previous studies on EEG phase extraction have used finite impulse response (FIR) filters to make the signal narrow-band in its frequency spectrum [START_REF] Mormann | Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients[END_REF][START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF][START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF][START_REF] Fell | The role of phase synchronization in memory processes[END_REF][START_REF] Freeman | Origin, structure, and role of background EEG activity. Part 2. Analytic phase[END_REF]. However, in order to have a reliable phase sequence, there are important considerations regarding this procedure, including: the filter's band-width, its phase response and the convolution process.

The second stage requires choosing a phase estimation method to extract the phase sequence from the narrow-band signal. The most common method for phase estimation is based on the analytic signal representation of the narrow-band signal [START_REF] Chavez | Towards a proper estimation of phase synchronization from time series[END_REF]. As shown in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF], the calculation of the instantaneous phase from the analytical representations becomes challenging and highly susceptible to noise in low analytical signal amplitudes; resulting in fake jumps and spikes in the extracted phase signal. In the following sections, this issue is further studied and partially solved by applying perturbations in the phase extraction procedure followed by ensemble averaging, as proposed in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF].

Unambiguous phase estimation conditions

For the signal x(t), its analytical form is defined as follows [START_REF] Picinbono | On instantaneous amplitude and phase of signals[END_REF]:

z x (t) = x(t) + jH{x(t)} (1)
where H{•} denotes the Hilbert transform. Using the analytical form, the instantaneous envelope (IE), the instantaneous phase (IP), and the instantaneous frequency (IF) are uniquely defined as follows:

IE: a x (t) = |z x (t)| = x(t) 2 + H{x(t)} 2 (2)
IP: φ x (t) = arctan H{x(t)} x(t) (3)
IF: f x (t) = 1 2π [ dφ x (t) dt mod 2π] (4)
For discrete-time implementations, ( 4) is replaced by a finite-difference approximation of the derivative operator (e.g., the first-order difference in the simplest case).

Unless the signal has a narrow-band spectral support, the pair (a x (t), φ x (t)) does not convey significant information regarding the instantaneous phase [START_REF] Boashash | Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals[END_REF]. This is due to the fact that only for narrow-band signals, the relative variations of the amplitude a x (t) are rather slow (and negligible) as compared with the variations of the phase φ x (t) [START_REF] Picinbono | On instantaneous amplitude and phase of signals[END_REF], i.e., dφ(t) dt

1 a(t) da(t) dt (5)
For cerebral signals, it is known that the EEG has a wide frequency range (0 Hz-150 Hz in the extreme case), which makes narrow bandpass filtering (of the order of Hertz) an essential prerequisite for extracting a meaningful instantaneous phase sequence.

Linear-phase filtering

Previous studies on EEG phase extraction have commonly employed linear-phase FIR filters to make the signal narrow-band in its frequency spectrum [START_REF] Mormann | Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients[END_REF][START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF][START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF][START_REF] Fell | The role of phase synchronization in memory processes[END_REF][START_REF] Freeman | Origin, structure, and role of background EEG activity. Part 2. Analytic phase[END_REF]. The advantage of linear phase filters is their constant group delay, which avoids phase distortions in the filtered signal. Nevertheless, in most FIR filter design techniques, the order of the filter proportionally increases with the inverse of their transition bandwidth, which means that narrow-band FIR filters have very long impulse responses and inputoutput delays [START_REF] Oppenheim | Discrete-time signal processing, ser. Prentice-Hall signal processing series[END_REF]. Moreover, highly narrow band filters are difficult to design and susceptible to design parameters.

To avoid these issues, previous studies have kept a trade-off between the order of the FIR filter and its bandwidth (BW). The BW is commonly chosen relatively wide, e.g. between 4 to 12 Hz, to have a low-order and practically realizable filter [START_REF] Mormann | Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients[END_REF][START_REF] Freeman | Origin, structure, and role of background EEG activity. Part 2. Analytic phase[END_REF][START_REF] Mormann | Epileptic seizures are preceded by a decrease in synchronization[END_REF]. However, as discussed before, using a bandwidth in this range, the envelope-phase pair {a(t), φ(t)} obtained from the Hilbert transform fails to correctly (and uniquely) define the instantaneous envelope. Therefore, the extracted phase becomes unreliable. To overcome this issue, very narrow band filters with reasonably low orders are required in practice to be less sensitive to noise and filter parameter variations. In Section 3, it is shown that for offline applications, zero-phase infinite impulse response (IIR) filters can be used as an alternative solution.

Low-amplitude analytical signal

The instantaneous phase sequence derived from the analytic representation of a signal is prone to fake jumps (without cerebral source) in low-amplitude analytical signal (LAAS) time instants [START_REF] Picinbono | On instantaneous amplitude and phase of signals[END_REF][START_REF] Chavez | Towards a proper estimation of phase synchronization from time series[END_REF][START_REF] Freeman | Origin, structure, and role of background EEG activity. Part 2. Analytic phase[END_REF]. As depicted in Fig. 1, the instantaneous frequency tends to have big jumps at LAAS epochs. The problem was rigorously studied in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF]. The main reason underlying this phenomenon could be linked to the arctan(•) operator used for phase calculation. According to (3), LAAS leads into very small numerator and denominator values. Consequently, any minor change in the real or imaginary parts of the analytic form (due to noise or background EEG fluctuations), results in a significant change in the estimated phase. This fact is illustrated in Fig. 2 for a sample EEG signal. It is seen that phase values corresponding to lower amplitudes tend to have bigger variations. More rigorously, it was shown in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF] that during LASS epochs the instantaneous phase tends to a uniform distribution over [-π, π] and in discrete-time implementations, the instantaneous frequency becomes uniform over the entire Nyquist band [0, f s ]. The findings of [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF] are in accordance with previous research [START_REF] Freeman | Origin, structure, and role of background EEG activity. Part 2. Analytic phase[END_REF], which reported that LAAS occurs more frequently in low power time-frequency regions, such as the high frequency bands of the EEG, for which the EEG power decays more rapidly with increasing frequency. To illustrate this point, in Fig. 3 the instantaneous phase differences of frequency components in the range of DC to 50 Hz have been extracted alongside their corresponding instantaneous envelopes. Accordingly, the first 5 s of the results are due to the FIR filter's transient response; resulting in very low-magnitude instantaneous amplitudes during this period. As it can be seen, the corresponding regions of the phase difference plot contains many phase jumps and spikes. For the rest of the signal, in lower frequencies, where the analytic form has higher amplitudes, the phase sequences are less contaminated with jumps and spikes. However, as the frequency increases and the power in EEG signal decays, the analytical signal envelope decreases and the rate of phase jumps increases once more.

Based on these findings, in the next section a robust method is proposed for the estimation of the instantaneous phase using perturbation of filter parameters and a Monte Carlo estimation process.

Method

The overall proposed scheme is summarized in Algorithm 1. In this algorithm, the bandpass filter is specified with a set of parameters: bandwidth (BW), transition-band width (TB), maximum pass-band ripple (PR) and minimum stop-band attenuation (SA) that we consider as the filter's design parameter set λ = {BW, TB, PR, SA}, which in term result in the zero-pole sets Z = {z 1 , z 2 , ..., z n } and P = {p 1 , p 2 , ..., p m }, depending on the filter design technique. The variables δ p i and δ z i are uniformly distributed random variables with very minor standard deviations (10 -4 in the sequel), used for making random perturbations in the zero-pole loci. The details of each stage of the algorithm are elaborated in the following subsections.

Step 1: Narrow-band zero-phase smoothing

The challenges in narrow-band FIR filter design were noted in Section 2.3. To overcome these issues, we propose to use forward-backward zero-phase IIR filters. Although the filter is performed offline in a non-causal manner (which is not a limiting issue for offline applications), the major advantage is that the order of a narrow band IIR filters is much lower than its FIR counterpart and by applying it in a forward-backward manner, the nonlinear phase response of the filter is compensated and zero-phase distortion-which is necessary for EEG phase analysis-is guaranteed.

Algorithm 1 Robust Instantaneous EEG Phase Extraction

Require: Input EEG signal x(t) Require: Bandpass filter prototype h BP (t) with design parameters λ = {BW, TB, PR, SA} and zero-pole sets Z and P , respectively. Require: Number of Monte Carlo iterations (N )

1:

for all i = 1 • • • N do 2:
Perturb the poles and zeros with minor random deviations δ p i and δ z i , while keeping the zero-pole pairs conjugate symmetric and preserving the poles inside the unit circle:

P i ← P ± δ p i , Z i ← Z ± δ z i 3:
Construct the bandpass filter h Zero-phase forward-backward filter the input signal

x i (t) ← h (i) BP (t) * x(t) 5:
Form the analytic representation of x i (t):

z i (t) ← x i (t) + jH{x i (t)} 6:
Calculate the instantaneous phase:

φ i (t) ← arctan( Im{z i (t)} Re{z i (t)} ) 7:
Unwrap the estimated phase sequences Φ i (t) ← unwrap{φ i (t)} 8: end for 9: Ensemble average over all iterations:

Φ(t) ← 1 N N i=1 Φ i (t)
Various types of IIR filters such as Chebyshev types 1 and 2, Butterworth and Elliptic filters were studied to determine the best filter for this application, and the Elliptic filter was chosen due to its steeper roll-off characteristics (as compared with Butterworth or Chebyshev filters) and its equi-ripple behavior in both the passband and stopband. In general, by allowing ripples in both passband and stopbands, Elliptic filters meet given performance specifications with the lowest order as compared with their counterparts [START_REF] Lutovac | Filter design for signal processing using MATLAB and Mathematica[END_REF]. In order to preserve the filter's frequency response over all frequency bands, instead of designing various bandpass filters in each band, a fixed narrow band lowpass filter prototype was designed and applied to the signal using a frequency domain shifting scheme illustrated in Fig. 4. Further details regarding this scheme is presented in Section 3.2. The prototype Elliptic IIR lowpass filter used in this study has the following characteristics: 0.3 Hz pass-band frequency, 0.5 Hz stopband cutoff frequency, 0.1 dB maximum pass-band ripple, and 70 dB minimum stopband attenuation, designed at a sampling frequency of 160 Hz (the sampling rate of the sample EEG). The order of this prototype filter was 6, which is computationally far more effective than any FIR filter with the same specifications. This filter was performed in a forward backward manner, which doubles its pass-band ripple and stop band attenuation in dB.

The procedure of zero-phase forward-backward smoothing (FBS) is shown in Fig. 5. This process can be implemented using the filtfilt function in major signal processing languages such as Matlab, Octave, or R. Accordingly, FBS uses the time-reversal property of the Fourier transform to perform zero-phase smoothing by processing the input signal in both the forward and reverse directions [START_REF] Oppenheim | Discrete-time signal processing[END_REF]. Considering H(jω) as the frequency response of the forward path digital filter, the overall response of FBS is

H eff (jω) = |H(jω)| 2 (6)
which is real-valued. Therefore, regardless of the nonlinear phase-response of the IIR filter, FBS has a zero-phase (and zero-group delay) response, which preserves the input signal's phase.

Step 2: Phase calculation

The next step is to compute the phase sequence. For this, we use the analytic representation of the filtered signal. In order to reduce the processing complexity and avoid the direct calculation of the Hilbert transform, this stage can be merged with the bandpass filtering as follows: as illustrated in Fig. 4, the proposed bandpass filtering scheme uses a lowpass filter prototype. To filter the signal x(t) around the center frequency ω 0 , x(t) is shifted in the frequency by multiplying the pure phase signal exp(jω 0 t), to obtain a complex valued signal x f (t). Next, the real and imaginary parts of x f (t) are given to the lowpass prototype H LP (jω) to obtain the narrow-band analytical signal y f (t). Finally, y f (t) is shifted back to the center frequency ω 0 , by multiplying the phase signal exp(-jω 0 t). This procedure provides x a (t), which is the narrow band analytical form of the original signal x(t) around the center frequency ω 0 . After computing the analytic form of the filtered EEG, the instantaneous phase is calculated as follows:

φ(t) = arctan Im(x a (t)) Re(x a (t)) (7)
For discrete-time signals, the instantaneous frequency can be approximated by the first order difference of the instantaneous phase:

f (t) ≈ f s φ(t) -φ(t -∆) 2π ( 8
)
where ∆ is the sampling time and f s = 1/∆ is the sampling frequency.

Step 3: Zero-pole perturbation of the prototype filter

In the proposed scheme, the previous two stages (narrow band filtering and phase calculation) are repeated several times, each time with very minor changes in the filter design parameters λ (and zero-pole positions). Here, the idea is to generate random ensembles of the signal's analytical form and EEG phase, using infinitesimal perturbations in the filter parameters. Apparently, clinically relevant EEG phase information should not be susceptible to minor filter design variations at the order of, e.g., 0.01 Hz. However, during LAAS epochs, even minor deviations in the filter parameters can significantly change the phase estimates, resulting into fake phase jumps.

The zero-pole plot of the utilized prototype IIR filter and the loci of its perturbed zeros and poles are illustrated in Fig. 6(c). It is known that all Elliptic IIR filter zeros are located on the unit-circle. Thus, in order to prevent any major change in the filter characteristics, these zeros are perturbed randomly only on the unit-circle. For this, the filter's zeros are taken into polar coordinates (ρ, θ), and the random perturbations are applied only to θ. Another important consideration is that any perturbation in the pole loci of a causal filter should not move its poles out of the unit-circle (to preserve its stability). Moreover, the conjugate symmetry of the zeros and poles should be preserved to guarantee the realness of the impulse response. These properties are achieved by applying random complex-valued perturbations δ p i (to the poles) and δ z i (to the zeros) to each pair of conjugate symmetric pair of poles/zeros. The complex values added to each conjugate symmetric pair should have equal real parts, and imaginary parts with equal magnitudes and opposite signs, to preserve the conjugate symmetry of the perturbed poles/zeros. The perturbed pole loci is finally checked to be inside the unit circle (Fig. 6(c)), to guarantee the filter's stability. The magnitude and phase response of 100 lowpass prototype filters obtained by the proposed zero-pole perturbation procedure are shown in Fig. 6(a) and Fig. 6(b), where it is seen that the zero-pole perturbations have had rather minor impact on the filters' response (irrelevant to most cerebral studies). However, it is later shown that even these minor changes can significantly change the EEG phase, especially during LAAS.

Step 4: Phase unwrapping

Calculating the phase using the four quadrant arc-tangent causes phase-wrapping [START_REF] Itoh | Analysis of the phase unwrapping algorithm[END_REF]. The amplitude of the phase sequence can take any value and even exceed the range [-π, π]. In cases where the phase exceeds this range, it is wrapped so that it stays within the principal range [START_REF] Itoh | Analysis of the phase unwrapping algorithm[END_REF][START_REF] Ghiglia | Two-dimensional phase unwrapping: theory, algorithms, and software[END_REF]. In such cases, the wrapped phase sequence will contain some phase-jumps greater than ±π. Therefore, for EEG phase analysis (either from the phase itself or from its time difference), an unwrapping procedure is required to obtain the phase sequence in its original continuous form. An unwrapped phase sequence typically diverges from zero over time (similar to a random-walk process). For better illustration of the phase fluctuations, one may either subtract the constant linear phase signal ω 0 t from the instantaneous phase to obtain its temporal fluctuations, or alternatively use the time difference of the phase signal, which is proportional to the instantaneous frequency.

Step 5: Ensemble averaging

The final stage of the algorithm is to average over the randomized ensembles of the phase sequences obtained from different filter responses, to obtain an average phase estimate. For illustration, Fig. 7 shows 100 ensembles of the instantaneous frequency of a sample EEG segment obtained by the aforementioned randomization scheme (zeropole perturbation) for a center frequency of 8 Hz, and the average of the randomized ensembles. This figure illustrates the importance of the proposed scheme and the significance of the zero-pole perturbation. Accordingly, without this procedure (by simply calculating the phase from a single filtered signal as in conventional methods), the obtained phase sequence is unreliable, since each set of filter parameters (even with minor differences) would lead to significantly different results; especially during low analytical signal envelope epochs in which the standard deviation of the instantaneous frequency is significantly high (Fig. 7).

Example

In this section, the robustness of the proposed phase/frequency calculation method is verified versus conventional methods (without random ensemble generation and averaging) for a sample EEG signal. The sample signal, represented in Fig. 9(a), consists of 24 s of an ongoing EEG acquired with a sampling frequency of f s =173.61 Hz, recorded during a BCI study. The complete description of the data is presented in [START_REF] Andrzejak | Indications of nonlinear deterministic and finite-dimensional structures in time series of brain electrical activity: Dependence on recording region and brain state[END_REF]. In the following, the phase robustness of this sample data is studied from three aspects: (1) impact of filter parameter variations, (2) robustness to non-stationary background noise and (3) low-amplitude analytic signal and phase jumps. The power spectral density (PSD) of this sample data is shown in Fig. 8. Accordingly, f 0 =8 Hz (the alpha-band peak) is chosen as the center frequency in the following evaluations.

Filter Parameter Variations:

In order to assess the sensitivity of the proposed method to filter parameter variations, the following three scenarios are considered for digital lowpass Elliptic IIR filter design. All filters are designed for a sampling frequency of f s =173.61 Hz (the sampling rate of the sample EEG), and a 0 dB DC gain. The instantaneous frequency obtained from (8) for these three scenarios are shown in Fig. 9(b), using both the conventional and proposed methods. It can be seen that the conventional method is very sensitive to variations in filter parameters and the results have significantly changed with minor changes in the bandpass filter design parameters. On the other hand, the proposed method has been more stable to such parameter variations and the filter design variation effects are almost removed during the frequency response perturbations and ensemble averaging.

3.6.2.

Low-amplitude Analytic Signal and Phase Jumps: Fig. 10 shows the instantaneous phase difference and the instantaneous analytical signal envelope calculated for the sample EEG in Fig. 9(a) using the conventional method, for frequency components in the range of DC to 30 Hz, using the three filtering schemes stated in Section 3.6.1. The same procedure is performed using the proposed method and the results are shown in Fig. 11 for comparison.

The comparison of Fig. 10 and Fig. 11, clearly shows the effects of LAAS on the EEG phase jumps using the classical and proposed methods. It is seen that the phase sequences estimated by the conventional method are prone to fake jumps at points where the corresponding analytic signal have lower amplitudes; while this issue has been significantly improved using the proposed method.

3.6.3. Noise Susceptibility: For evaluating the robustness of the proposed phase estimation method to noise, the results of phase estimations using both the proposed and conventional methods are tested in presence of additive non-stationary noise. Non-stationary noise with a time-variant variance and EEG-like spectra is used to model background non-stationary EEG activity. For this purpose, a time-varying autoregressive (TVAR) model of order 20 is trained over sliding windows of length 4 s using a sample 20 s EEG. The resulting TVAR model is next fed by white Gaussian noise, to generate non-stationary noises of the same length with a spectra similar to a real EEG and added to the EEG understudy at different signal-to-noise ratios (SNR). For this scenario, an Elliptic IIR filter with f 0 =8 Hz, BW=0.5 Hz, TB =0.5 Hz, δ p =0.1 dB and δ s =70 dB is used as the bandpass filter.

The corresponding results are shown in Fig. 12 for both the conventional and presented methods in SNR = 15 dB. It can be seen that while the conventional phase estimation procedure (gray line) becomes unreliable during noisy epochs, the proposed method (thick black line) yields a robust estimation of phase information despite the non-stationary background noise (which is statistically independent from the underlying mental task). Therefore, the proposed method significantly improves the reliability of the extracted instantaneous phase/frequency estimates. A rigorous discussion on the effect of noise and SNR level on the probability of correct and false phase detections has been presented in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF]. Fig. 13 shows the instantaneous phase difference and the instantaneous analytical signal envelope calculated for the sample EEG in Fig. 9(a), using the conventional method for frequency components in the range of DC to 30 Hz, at three SNR: infinity (no noise), SNR= 10 dB, and SNR= 0 dB. In Fig. 14, the same procedure is repeated using the proposed method for comparison; where it is seen that the proposed method has been considerably less susceptible to additive noise.

Case Study: Phase and Frequency Features for a Brain Computer Interface Application

In order to show the significance of the proposed method, its performance is evaluated for a visual evoked potential (VEP)-based BCI system, as proof of concept. The stateof-the-art classification procedure used in previous VEP-based BCI studies is employed for this purpose [START_REF] He | Combination of amplitude and phase features under a uniform framework with EMD in EEG-based Brain-Computer Interface[END_REF][START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] Townsend | Using phase information to reveal the nature of event-related desynchronization[END_REF], and we focus on the impact of EEG phase features using conventional versus the proposed schemes. It should be noted that this case study is only presented as a typical application for reliable and robust EEG phase/frequency extraction and we do not seek improvements in classification rates of state-of-the-art BCI systems. Therefore, only phase-related features and standard classifiers are used.

Dataset

The dataset used for this study is adopted from the Neuroelectric and Myoelectric Databases, which are online available on Physionet [START_REF] Goldberger | Physiobank, physiotoolkit, and physionet components of a new research resource for complex physiologic signals[END_REF]. This dataset includes one and two-minute recordings of 109 volunteers, performing a series of motor and motorimagery tasks. Each record contains sixty four channels of EEG recorded using the BCI2000 System, during a set of annotated mental tasks [START_REF] Schalk | BCI2000: a general-purpose brain-computer interface (BCI) system[END_REF]. The complete description of the dataset is available at [START_REF] Goldberger | Physiobank, physiotoolkit, and physionet components of a new research resource for complex physiologic signals[END_REF]. Each subject has performed a series of mental tasks: two one-minute baseline runs, with open and closed eyes and three two-minute runs of the four following tasks:

(i) A target appears on the left or right side of a screen in front of the subject. The subject opens and closes the corresponding fist until the target disappears. Then the subject relaxes.

(ii) A target appears on the left or right side of the screen. The subject imagines opening and closing the corresponding fist until the target disappears. Then the subject relaxes. (iii) A target appears on the top or bottom of the screen. The subject opens and closes both fists (if the target is on top) or both feet (if the target is on the bottom) until the target disappears. Then the subject relaxes.

(iv) A target appears on the top or bottom of the screen. The subject imagines opening and closing both fists (if the target is on top) or both feet (if the target is on the bottom) until the target disappears. Then the subject relaxes.

Fig. 15 shows the placement of the electrodes used for recording EEG signals in this dataset. Since the primary cortical regions involved in the task of motor imagery are the supplementary motor area (SMA) and the primary motor cortex area (M1), electrodes FCz, C3, and C4 are chosen for this study [START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] Wei | Amplitude and phase coupling measures for feature extraction in an EEG-based brain-computer interface[END_REF][START_REF] Wang | Phase synchrony measurement in motor cortex for classifying single-trial EEG during motor imagery[END_REF]. The annotations provided in the dataset consist of three classes for identifying rest versus left/up or right/down side activities: 1) T 0 corresponding to rest condition, 2) T 1 corresponding to motion (real or imagined) onset of either the left fist or both fists, and 3) T 2 corresponding to motion (real or imagined) onset of either the right fist or both feet. We therefore have a three-class classification problem.

The targets appeared on the screen every four seconds, resulting in thirty 4 s annotated EEG segments for each two minute records (per subject), each corresponding to a mental task trial.

Feature Extraction

A survey of previous studies on VEP-based BCI systems reveals that EEG phaserelated features are among the most discriminative and informative features for BCI applications [START_REF] He | Combination of amplitude and phase features under a uniform framework with EMD in EEG-based Brain-Computer Interface[END_REF][START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] Pan | Enhancing the classification accuracy of steadystate visual evoked potential-based brain-computer interfaces using phase constrained canonical correlation analysis[END_REF]. During the feature extraction phase, a broad range of features are commonly extracted from the frequency band of interest and passed to the feature selection and classification stages. However, in this study, a single phase-related feature, namely the PLV, is used to evaluate the robustness and feasibility of the proposed instantaneous phase estimation procedure.

PLV is an index for quantifying how constant the phase difference between two signals is. In order to calculate the PLV between two signals (or channels) x(t) and y(t), the following steps are applied [START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF][START_REF] Rosenblum | Phase synchronization of chaotic oscillators[END_REF]:

• Channelize the signals using narrow-band filters centered at f

• Calculate the instantaneous frequency-specific phase values φ x (t, f ) and φ y (t, f ).

• Calculate the instantaneous phase-difference between x(t) and y(t) and quantify the local stability of this phase-difference over time:

PLV xy (f ) = 1 T T t=1 exp (j[φ x (t, f ) -φ y (t, f )]) ( 9
)
where T is the signal length and the summation is taken over all temporal samples of the instantaneous phases.

PLV varies between 0 and 1, corresponding to completely non-synchronized signals and complete synchronization, respectively [START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF][START_REF] Rosenblum | Phase synchronization of chaotic oscillators[END_REF]. The notion of phase-locking using PLV has been previously studied in BCI applications [START_REF] Hamner | Phase-based features for motor imagery brain-computer interfaces[END_REF][START_REF] Carreiras | Phase-locking factor in a motor imagery brain-computer interface[END_REF]. For illustration, the envelope and phase sequences (phase derivatives), and the corresponding PLV index extracted from two 10 s segments of EEG acquired from channels C3 and C1 of the first subject in the dataset are shown in Figs. 16(a) and 16(b). During this segment, the subject is performing a right-hand task (starting between t = 55 and t =56 s) for four seconds and then relaxes. The PLV index obtained by averaging over 250 ms sliding windows (controlled by the window length T , in ( 9) ) is shown in Fig. 16(c), using phase values obtained from the conventional and proposed methods. Accordingly, the PLV computed using the conventional phase estimation method is reporting some de-synchronization around t =50.2 s, 50.6 s, 52.4 s, 57.7 s and 59 s, which have coincided with the envelope drops at the same instants (Fig. 16(a)). This shows that phase variations during low analytical signal epochs are unreliable in conventional methods, which is not the case for the PLV reported by the proposed method.

In order to evaluate the overall performance, the EEG phase and the corresponding PLV features were extracted using the conventional and proposed procedures, from 105 out of 109 subjects (four subjects were excluded due to data deficiency). The PLV was calculated for a single frequency band f = 10 Hz with an effective bandwidth of 1 Hz, which was identified as the dominant alpha-band peak by visual inspection of the EEG spectra. The index was calculated between the three possible combinations of the selected electrodes, i.e., FCz-C3, FCz-C4 and C3-C4 (as shown in Fig. 15), resulting in PLV feature vectors of length three. The feature vectors were computed from each of the thirty 4 s annotated temporal windows over all two-minutes records of each subject.

Classification

The PLV feature vectors calculated by both the conventional and proposed methods together with the annotations provided in the database were used for training and testing the classifiers. K-Nearest Neighbors (KNN) with K=30 (the number of nearest neighbors used in the classification), and Random Forest (RF) with number of trees equal to 10, were used as classifiers in a leave-one-out cross-validation process, in which the feature-set of one subject is considered as test data and the rest of the feature sets are used for training the classifiers.

Results

The comparison has been made both in absence and presence of an additive white Gaussian noise with SNR=5 dB, to investigate the robustness of the conventional and proposed procedures to noise. The average results of the noiseless and noisy cases over all subjects are reported in Fig. 17. Accordingly, the proposed method for extracting PLV features has improved the mean RF classification rates 7%, as compared with conventional methods of EEG phase extraction (in absence of additive noise) and 12% (in presence of additive noise). Using the KNN classifier, the mean classification rates have been improved 4% (in absence of additive noise) and 8% (in presence of additive noise). The significance of these results was tested by a paired sample t-test between the conventional and proposed methods, under the null hypothesis of equal average performance of PLV features obtained by both methods (conventional versus proposed).

The null hypothesis was rejected with p-values p = 0.01 and p = 0.03, for the noiseless and noisy scenarios, respectively, using 6442 degrees of freedom in both cases. The per-subject three-class accuracies using PLV as feature are shown and compared for the proposed and conventional methods in Fig. 18; where the outperformance of the proposed method is seen, both in terms of higher classification accuracies and smaller inter-subject variances. The higher performance can be In order to make the results reproducible, all source codes related to this study are online available in the open-source electrophysiological toolbox (OSET) [START_REF] Sameni | The Open-Source Electrophysiological Toolbox (OSET)[END_REF].

Discussion

The EEG phase is a rich source of information for various fields of brain studies. Conventional methods for calculating the instantaneous phase and frequency of EEG signals are unreliable in presence of spontaneous background EEG activity (and noise) and during low analytical signal envelopes [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF]. Therefore, robust methods for phase calculation are required.

Herein, a robust phase estimation procedure was proposed to overcome these issues. The proposed method has additional steps as compared to conventional methods: 1) zero-pole perturbation of the utilized bandpass filters, 2) offline forward-backward zerophase filtering, and 3) ensemble averaging between the perturbed phase estimates to obtain robust estimates.

The zero-pole perturbation decreases the effects of LAAS and provides more reliable instantaneous phase sequences. It was shown that the proposed method, significantly reduces the sensitivity to noise and variations in filter parameters, which have been commonly neglected in previous studies. However, as recently reported in [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF], the phase ambiguities due to LAAS are unavoidable and should be considered as an intrinsic limitation for phase estimation.

For BCI applications, the results presented in Section 4 demonstrate the significance and robustness of the proposed phase extraction procedure. The results show that phaserelated features obtained through this method not only outperform conventional phase features, but also are more robust to noise for BCI applications.

The proposed method is also useful for robust instantaneous EEG frequency estimation. For instance, previous studies have shown that the attention of an individual to a particular visual stimulus with a specific flashing rate can be detected by inspecting the peaks of the EEG frequency spectra at the corresponding flashing rate [START_REF] Middendorf | Brain-computer interfaces based on the steady-state visual-evoked response[END_REF]. Therefore, accurate estimation of the instantaneous frequency of the EEG over time is of particular importance in designing SSVEP-based BCI systems. Although the hereby utilized dataset was not based on SSVEP, the proposed method can be used as a potential instantaneous frequency estimator in SSVEP-based BCI studies.

The scope of the current study is not limited to BCI applications. In fact, without using the hereby proposed scheme, the effects of LAAS and narrow-band filtering lead to unreliable and ambiguous phase sequences, which result in wrong interpretations of phase-related quantities in any similar application. This highlights the necessity of the mentioned additional steps in phase estimation, to improve the reliability of the estimated instantaneous phase sequence of an EEG for different applications.

The theoretical findings of [START_REF] Sameni | A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis[END_REF] and the hereby reported results further highlight the importance of the analytical signal envelope in EEG phase-related studies. In future studies, the combination of phase and analytical signal envelopes can be used for improving performance in BCI and other applications, including sleep stage classification and pathological cases such as epileptic EEG. We also expect the hereby developed scheme to have considerable applications in brain connectivity and cognitive studies based on the EEG phase.

Figure 1 .

 1 Figure 1. A raw 60 s EEG signal used for phase extraction (top panel), the analytical signal's instantaneous envelope (IE) (middle panel), and the corresponding instantaneous frequency (IF) centered at 8 Hz obtained from the phase sequence (bottom panel).

Figure 2 .

 2 Figure 2. A 14 s segment of a raw EEG signal used for extracting the instantaneous phase sequence (top panel); the analytical signal envelope (middle panel); the polar representation of two segments of the analytical form, at time instants with (b) high and (c) low analytical signal envelope (bottom panel). The phase jump between successive samples (which is proportional to the instantaneous frequency) are considerably greater in lower analytic envelopes.

Figure 3 .

 3 Figure 3. Phase difference of the frequency components from DC to 50 Hz (left panel), the corresponding instantaneous analytical signal envelope (right panel), for a 25 s EEG segment.

  with perturbed zero-pole pairs Z i and P i 4:

Figure 5 .

 5 Figure 5. Block-diagram of the forward-backward filtering process.

Figure 6 .

 6 Figure 6. The (a) magnitude and (b) phase of the frequency responses of 100 lowpass prototype filters obtained by zero-pole perturbation; (c) zero-pole plot of the prototype Elliptic IIR filter. The contours specify the region of random perturbations of zeros and poles.

Figure 7 .

 7 Figure 7. (top panel) The gray shades show ensembles of the estimated instantaneous using 100 filters with minor zero-pole perturbations and the final averaged instantaneous frequency in solid black line; (bottom panel) the standard deviation of the instantaneous frequency for 100 ensembles at each time instant.

Figure 8 .

 8 Figure 8. The PSD of the sample EEG used for evaluation in Section 3.6. The peak at f 0 =8 Hz is chosen as the center frequency (in the alpha-band) in all evaluations.

(i) f 0 = 8

 8 Hz (center frequency), BW=0.5 Hz (pass-band bandwidth), TB=0.5 Hz (transient-band from either side of the center frequency), δ p =0.1 dB (maximum pass-band ripple), and δ s =-70 dB (minimum stop-band attenuation) (ii) f 0 = 8 Hz, BW=0.2 Hz, TB=0.2 Hz, δ p =0.1 dB, and δ s =-70 dB (iii) f 0 = 8 Hz, BW=1.0 Hz, TB=1.0 Hz, δ p =0.1 dB, and δ s =-70 dB

Figure 9 .

 9 Figure 9. (a) Sample ongoing EEG signal to assess the robustness of the proposed phase estimation method. (b) Instantaneous frequency around f 0 =8 Hz calculated by the conventional (gray) and proposed method (black) for three filter parameter sets described in Section 3.6.1 (from top to bottom)

Figure 10 .

 10 Figure 10. The instantaneous phase difference (top row) and the instantaneous amplitudes (bottom row) calculated from the EEG signal shown in Fig. 9(a) using the conventional phase extraction procedure, for frequency components in the range of DC to 30 Hz, (from left to right) with three different sets of filter parameters 1) BW=0.5 Hz, TB=0.5 Hz, 2) BW=0.2 Hz, TB=0.2 Hz and 3) BW=1.0 Hz, TB=1.0 Hz (detailed in Section 3.6.1).

Figure 11 .Figure 12 .

 1112 Figure 11. The instantaneous phase difference (top row) and the instantaneous amplitudes (bottom row) calculated from the EEG signal shown in Fig. 9(a) using the proposed phase extraction procedure, for frequency components in the range of DC to 30 Hz, (from left to right) with three different sets of filter parameters 1) BW=0.5 Hz, TB=0.5 Hz, 2) BW=0.2 Hz, TB=0.2 Hz and 3) BW=1.0 Hz, TB=1.0 Hz (detailed in Section 3.6.1).

Figure 13 .

 13 Figure 13. The instantaneous phase difference (top row) and instantaneous amplitudes (bottom row), calculated through the conventional phase extraction procedure for frequency components in the range of DC to 30 Hz for the sample EEG shown in Fig. 9(a), in three different cases: 1) no noise, 2) non-stationary noise with SNR=10 dB, and 3) non-stationary noise with SNR=0 dB from left to right, respectively.

Figure 14 .

 14 Figure 14. The instantaneous phase difference (top row) and instantaneous amplitudes (bottom row), calculated through the proposed phase extraction procedure for frequency components in the range of DC to 30 Hz for the sample EEG shown in Fig. 9(a), in three different cases: 1) no noise, 2) non-stationary noise with SNR=10 dB, and 3) non-stationary noise with SNR=0 dB from left to right, respectively.

Figure 15 .

 15 Figure 15. The location of the three (out of sixty four) leads used for the BCI classification case study. (a), (b) and (c) represent C4, C3 and FCz, respectively.

Figure 16 .

 16 Figure 16. The envelope (top panel) and phase (middle panel) sequences for 10 s EEG segments acquired from channels C3 and C1 of an EEG signal recorded during a righthand task from subject number one of the dataset; (bottom panel) the time-courses (dynamics) of the PLV values between C3-C1 (related to the task).

  The procedure of bandpass filtering and analytical form calculation using a lowpass prototype filter and frequency domain shifting of the input signal.

	x(t)	X	H LP (jω)		X	y(t)
		exp(jω 0 t)			exp(-jω 0 t)
	Figure 4. Time-Reversal x(t) z(t) H(jω)	w(t)	H(jω)	v(t)	Time-Reversal	y(t)

Table 1 .

 1 Per-subject three-class accuracies for the entire dataset (105 subjects), obtained using PLV as feature and KNN and RF as classifiers, for the conventional and proposed EEG phase extraction methods. Average and standard deviation of the results of Fig.18, for the proposed and conventional methods over all subjects robustness of the proposed EEG phase extraction method. The average and standard deviation of the per-subject results shown in Fig.18are reported in Table1, which show that KNN has slightly outperformed RF in all cases.

		80%
		75%
		70%
	Accuracy	60% 65%
			Conventional Method -RF
		55%	Conventional Method -KNN
			Proposed Method -RF
		50%	Proposed Method -KNN
			10	20	30	40	50	60	70	80	90	100
							subject index
			Figure 18. Method			Mean ± Standard Deviation
			Conventional Method -RF	65.45 ± 4.23
			Conventional Method -KNN 66.06 ± 4.03
			Proposed Method -RF	70.09 ± 2.56
			Proposed Method -KNN	70.90 ± 1.45
	associated to the

Acknowledgment

This research has been supported by the Cognitive Sciences and Technologies Council of Iran (COGC), under the grant number 2250.