Sergi Ventosa

Sylvain Le Roy

Irène Huard

Antonio Pica

Hérald Rabeson

Patrice Ricarte

Laurent Duval

Adaptive multiple subtraction with wavelet-based complex unary Wiener filters

Keywords: cwt, Morlet wavelet transform, unary adaptive multiple removal, Wiener filer, complex continuous wavelet frame 1

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Multiples correspond to unwanted coherent events related to wave field reflection bounces on given surfaces. We refer to the comprehensive survey in Verschuur [START_REF] Verschuur | Seismic multiple removal techniques: past, present and future[END_REF] for a detailed description. Their attenuation [START_REF] Verschuur | Adaptive surface related multiple elimination[END_REF] represents one of the greatest challenges in past and present seismic processing, since the January 1948 issue of Geophysics. Their processing mainly relies on two types of approaches:

(1) primary/multiple discrimination based on different velocities and reflections and (2) periodicity and predictability.

The first approach maps data to a transformed domain with minimum overlap between primaries and multiples. Transformed domain filters subsequently attenuate multiples or select primaries of interest. The filtered data is finally mapped back to data domain, using an appropriate inverse transform. Common tools include stacking combined with Normal Move Out (NMO), homomorphic filtering [START_REF] Buttkus | Homomorphic filtering -theory and practice[END_REF], f -k [START_REF] Wu | A case study of f -k demultiple on 2D offshore seismic data: The Leading Edge[END_REF], and τ -p representations [START_REF] Nuzzo | Improvement in GPR coherent noise attenuation using τ -p and wavelet transforms[END_REF], or alternative breeds -parabolic, hyperbolic -of the Radon transform [START_REF] Hampson | Inverse velocity stacking for multiple elimination[END_REF][START_REF] Trad | Latest views of the sparse Radon transform[END_REF][START_REF] Nowak | Amplitude preservation of Radon-based multiple-removal filters[END_REF], with pros & cons discussion in Kabir and Verschuur [START_REF] Kabir | Expert answers: Does parabolic Radon transform multiple removal hurt amplitudes for AVO analysis?[END_REF]. a sparsifying transform and its associated matched filtering. Its flexibility both resides in choices for a quasi-analytic wavelet and its scale discretization, to better adapt variations in data character and sampling. The algorithm alleviates some constraints of the curvelets (either real or complex) such as a relatively high memory footprint and its fixed dyadic frequency decomposition. However, as strictly 1D, it does not capitalize the information given by neighboring traces. Despite this, the closed-form of the matching filter may be beneficial in keeping low complexity, computational speed and possible parallelization. Being nevertheless rather competitive to 2D techniques, it may fit sparsely or non-uniformly sampled datasets.

Methodology

A classical 1D trace observation model is:

d[n] = p[n] + m[n] + w[n] , (1)
where

d[n], p[n], m[n]
and w[n] denote the recorded data, primary events, multiples and background noise, respectively, at discrete time index n.

Complex wavelet transform decomposition

We perform a time-scale decomposition of each data d[n] and multiple model x[n] traces with discrete wavelet frame approximations to continuous wavelet transform (CWT). We refer to Jorgensen and Song [START_REF] Jorgensen | Comparison of discrete and continuous wavelet transforms[END_REF] and references therein for comparison between discrete and continuous wavelet transforms. We choose the complex Morlet wavelet since it yields a simple interpretation of amplitude and phase delay in the transformed domain.

The mother Morlet wavelet is approximately analytic and writes:

ψ(t) = π -1/4 e -iω 0 t e -t 2 /2 , (2)
where ω 0 is the central frequency of the modulated Gaussian, t denoting the continuous time variable. A standard choice of π 2/ ln(2) for the central frequency halves side-lobe amplitude with respect to the main lobes. The associated discrete family of functions is defined as a sampling of the mother wavelet:

ψ v r,j [n] = 1 √ 2 j+v/V ψ nT -r2 j b 0 2 j+v/V , (3)
with T the sampling rate and V the number of voices per octave. Indices r, j ∈ Z and v ∈ [0, . . . , V -1] denote, respectively, discretized time, octave, and voice. Finally, b 0 stands for the sampling period at scale zero, j = v = 0. Its redundancy, of approximately 2V /b 0 , may be adapted to the expected computational efficiency. The time-scale representation of trace d[n], for instance, is given by the inner product:

d = d v r,j = d[n], ψ v r,j [n] = n d[n]ψ v r,j [n] , (4)
and written in bold, with indices dropped for more generality, accounting for the whole trace, or a windowed portion of it. Windowing in time-scale domain improves the localization in complex unary filters estimation. Expression ψ v r,j [n] denotes the complex conjugate of ψ v r,j [n].

Suppose d results from some time-scale correction or processing of d, here model matching and subtraction. Then the resulting filtered trace is synthesized back to the time domain as a sum of the dual frame components,

ψ v r,j [n], d[n] = r j,v dv r,j ψ v r,j [n] , (5)
weighted by the corrected signal dv r,j , for instance here with multiples subtracted. Illustrations of corresponding time-scale raw and processed data are provided in [START_REF] Ventosa | Complex wavelet adaptive multiple subtraction with unary filters[END_REF]. With sufficient redundancy, ψ[n] can be approximated by ψ[n] up to a constant multiplicative factor [see, e.g., 48].

The discriminative power of the wavelet frame helps reformulate the design of a long matching filter into combined global and local optimum complex unary filters, minimizing the error between multiple events and their matched model. The straightforward group-delay estimation allowed by Morlet wavelets, together with frame redundancy, helps avoiding reconstruction artifacts observed with non-redundant discrete wavelet processing [START_REF] Gilloire | Adaptive filtering in sub-bands with critical sampling: analysis, experiments, and application to acoustic echo cancellation[END_REF].

Unary filter estimation: amplitude and delay

When the group delay difference between model and multiple sequences is less than half a period at all scales, it can be corrected in the time-scale domain following an LSE approach. The optimum unary filter at a given wavelet scale, either in local or global windowed portions of the trace, is defined as the complex scalar a which, multiplied by the time-scale decomposed multiple model x, minimizes the mean square error ξ(a) with the decomposed input dataset d:

a opt = arg min a ξ(a) = arg min a d -ax 2 (6)
where the phase of a opt compensates small delays. This method is efficient, provided the above small group delay condition is satisfied [START_REF] Neelamani | Adaptive subtraction using complexvalued curvelet transforms[END_REF]. Otherwise, i.e., when the cross-correlation peak departs from the zero delay sample, faulty a opt estimation in equation 6 requires a redesign of the filter.

Two options can be conceived to add robustness against high fractional delays, while keeping with the unary filter approach: either to design filters with a high group delay, or to add a gross group delay correction first and then to apply the previous unary filter for fine group delay correction. In contrast with the Fourier transform, wavelets' fast decay strongly limits the maximum group delay correction on the first option. Hence, we focus exclusively on the second one. However, high group delay filters enable a relatively simple scale-wise unary filter design, that could reinforce matching in noisy scenarios, in applications where the specifications of group-delay range are smaller than the effective time support of the chosen wavelets.

We thus introduce a delay term into the above unary filter to yield higher delay robustness, while keeping an LSE criterion. This amounts to finding optimum a and delay l that minimize:

ξ(a, l) = r d v j [r] -a v j x v j [r -l] 2 = d -ax l 2 , (7)
where x l denotes a delay of x by l samples, sequences d v j [r] and x v j [r] are complex and the symbol denotes a locally weighted sum along odd L consecutive samples around index r. We opt to drop all subscripts for a opt for clarity. It is understood that this algorithm may be applied either locally or globally, depending on problem constraints.

For a given delay l, optimum a opt yields orthogonality between filtered signal y = d -a opt x l and adapted model ax l , which writes in inner product form as:

d -a opt x l , ax l = 0 ∀a = 0 . (8)
By developing and applying linearity on the first argument and conjugate linearity on the second one, we obtain:

a ( d, x l -a opt x l , x l ) = 0 ∀a = 0 (9)
and therefore

a opt [l] = d, x l x l 2 = r d v j [r]x v j [r -l] r x v j [r -l] 2 , (10)
i.e., the unary Wiener filter solution for complex signals.

Note that instead of performing a direct measurement of the fractional delay, or equivalently of the group delay, we estimate the phase delay at each scale or frequency component. When the signal-to-noise ratio (S/N) is high, all meaningful phase components are sufficiently well estimated. However, when the S/N decreases, moderate phase errors made at key scale components may induce huge errors on group delay. As a consequence, we have to ensure that the phase of a opt varies smoothly in scale for important components.

Several criteria, well adapted to the nature of the seismic signals, can be chosen in the selection of the optimum delay. The optimum integer delay, in LSE sense, is selected among all sub-optimal a opt [l] as the one minimizing ξ(a, l):

l opt = arg min l ξ(a opt [l]) . (11)
Alternatively, due to the importance of waveform over amplitude, we can define the optimum delay as the one that maximizes the normalized cross-correlation between the adapted model and the data, commonly called coherence [START_REF] Neidell | Semblance and other coherency measures for multichannel data[END_REF]39,[START_REF] Schimmel | Noise reduction and detection of weak, coherent signals through phase-weighted stacks[END_REF],

l opt = arg max l C[l] (12)
being

C[l] = Re d, a opt [l]x l d a opt [l]x l , (13)
where we focus on the real part of the normalized cross-correlation. The use of unary filters enables two important simplifications of the previous equations. Applying conjugate linearity in the second argument of the inner product and positive scalability of the norm gives:

l opt = arg max l Re a opt [l] |a opt [l]| d, x l d x l . (14)
Namely, the normalized cross-correlation between the adapted model and the data only depends on the phase of the optimum unary filter and the normalized cross-correlation between the original model and the data. Slighty better results can be obtained with an independent fine selection of the window functions used in the estimation of the unary filter and the normalized cross-correlation, equations 10 and 13. However, even more important to us, when these weighted sums are chosen equal, the normalized cross-correlation of the corrected sequence is equivalent to the cosine of the angle between d and x l . This can be shown noticing that a opt [l]/ |a opt [l]| can be written depending on d and x l from its definition in equation 10,

a opt [l] |a opt [l]| = d, x l | d, x l | . ( 15
)
As a consequence, the normalized cross-correlation is totally independent of the optimum unary filter in this particular case,

l opt = arg max l | d, x l | d x l (16)
and the coherence value measured is always real and positive. When no window is used, this criterion reduces to maximum cross-correlation between data and predicted multiples. This configuration enables an important reduction on the computational cost when the normalized cross-correlation estimator is used, as it enables a direct estimation of the optimum integer delay, without having to estimate the optimum unary filter for each candidate. Concisely, we first gauge an integer delay, equation 16, to rectify large time shifts; we then estimate the optimum unary filter, equation 10, to shrink amplitude and remaining delay differences.

Examples

Figure 1a represents synthetic data with two intersecting dipping events with wavelets of equal amplitude and frequency content; Figure 1b shows the multiple model. As expected, 1D algorithms (1D filters computed in 1D data windows) behave similarly and fail to properly cope with the intersection. However, a dataset showing single primary and multiple events with identical wavelet shape is not completely realistic. In practice, variations in wavelet shapes associated with uncorrelated time arrivals of primaries and multiple greatly contribute to the success of adaptation in a given time frame. Thus, better behavior of the proposed algorithm on actual data is expected. Firstly, it enhances primaries and multiples discrepancies in the wavelet domain. Secondly, it enables crosscorrelation reduction using larger window frames containing several signals, as demonstrated in the following field example. The examples shown in Figures 2a and3a are taken from a 3D real marine dataset, in common shot and channel gather, respectively. The main objective is to uncover primaries masked by strong multiple events using a multiple model. The 3D-SRME model [START_REF] Pica | 3D surface-related multiple modeling[END_REF] shown in Figures 2b and3b achieves good precision on the near and mid offset range. We compare the proposed method based on 1D complex unary filters in wavelet domain, with reference algorithms based on standard 1D and 2D adaptive filters. Figures 4 and5 represent subtraction results in common shot and channel gathers, respectively. Additional details are provided in Figures 6 (for improved visual amplitude assessment) and 7 (for accurate waveform comparison), to better illustrate the main features in data and differences between algorithms.

With a standard two-step method, the global step mainly contains waveform and time-shift corrections, combined with a global amplitude adaptation factor. In 2D, it is obtained with large 2D windows and long 2D filters. The global step is followed by a local adaptation, to account for local amplitude and phase discrepancies (in 2D, with small 2D windows and short 1D filters). Reference 2D algorithm parameters are taken from the best result obtained in an independent benchmark. Despite not being used in practice, the 1D adaptive filter is provided for completeness only.

With the proposed one-step method, a Morlet wavelet frame (with parameters ω 0 = 6.4 rad/s, j ∈ [START_REF] Abma | Comparisons of adaptive subtraction methods for multiple attenuation[END_REF][START_REF] Berkhout | Focal transformation, an imaging concept for signal restoration and noise removal[END_REF], v ∈ [0, 3] and b 0 = 1) is used to decompose the data and reconstruct the filtered signal. In the wavelet transformed domain, the optimum unary filter estimation for each sample uses a window of 0.636 s. Additionally, to provide robustness against delays higher than half signal period (at each scale), we have chosen the maximum normalized cross-correlation criterion within a range of ±12 ms. This subtraction operator may compensate for additional distortions as well, e.g., wavelet squaring introduced in convolution-based multiple models, since the estimation of the inverse wavelet is implicit in the estimation of an LSE filter.

single-pass, wavelet-based unary Wiener filters achieve slightly better non-coherent noise level (e.g., in Figure 6). As pointed out in ellipses in Figure 5, standard 1D lacks in terms of primary preservation and multiple attenuation. While the standard 2D approach exhibits best performance to this respect, as shown in the zoom of Figure 5, the proposed method exhibits close performance, despite not using additional information conveyed by neighboring traces. Figure 7 compares the waveforms of data, model, adapted and filtered results, for traces between indices 1470 and 1520 from Figure 6. We appreciate how the set of unary filters successfully reduces the differences between multiple model and actual multiple events, to attain sufficient multiple event attenuation to uncover main primary events. The trace overlay indicate the waveforms are well preserved.

To complete the picture, we provide additional results on near offset stacks. Figures 8 to 10 compare raw data and the different adaptive multiple subtraction methods. Note that stacking rea large portion of the multiples. The performance is thus evaluated on finer details, especially in the difference between raw data and filtered sections, as shown in Figure 10. The proposed 1D wavelet-based complex unary filter appears less aggressive with respect to the primaries previously masked by the multiples than the global adaptation step of a standard 1D method (after the local step, over-adaptation was even higher), compare the pointed signals in Figures 10a and10b. However, as previously inferred, the proposed method leaves a little more remnant multiples than the standard 2D method, as it can be appreciated by comparing figures around shot number 1200 and above.

The continuous wavelet frame provides an inherent scale-adaptive windowing that, when combined with the 1D complex unary filter, offers obvious improvements over standard 1D LSE timevarying frequency filters, and compares promisingly, at least away from isolated event crossings, with the industry standard 2D method on the common-channel gather, where a perfect removal of the source effects can still be demanding. Additional constraints may be added in the scale and in the space dimensions to improve the fractional delay estimation and reinforce the lateral coherence of the method, either with adaptive [9, 44] or low-redundancy multiscale transforms [START_REF] Chaux | Noise covariance properties in dual-tree wavelet decompositions[END_REF], which are robust to correlated noises and possess finer bandwidth selectivity than dyadic scales.

Conclusion

We propose a model-based multiple subtraction which combines complex Morlet wavelet frame with complex unary Wiener filters. This subtraction method counter-balances the redundancy provided by this complex wavelet frame with the freedom of designing fast filters which allow for simple non-stationary model adaptation. Complex unary filters with a low varying phase along scale can adapt fractional delays up to half the signal period at a given frequency. To add robustness against higher delays, we integrate an integer delay term in the kernel of the unary filters, which can be estimated globally or locally, depending on the application. Despite not competing with more advanced multiple removal methods, either regularized or in higher dimension spaces, the results of 1D unary filters in time-scale compare promisingly with standard 2D LSE methods. The structure of the algorithm allows for reduced memory footprint, higher code parallelization, and is potentially suitable to sparse acquisition. Further developments foresee additional diversity borrowed from neighboring traces, to improve apparent velocity separation and residual noise attenuation. The presented multiple removal technique aims at updating the portfolio of hybrid demultiple algorithms.

Figure 1 :

 1 Figure 1: Subtraction results on a synthetic dataset with two events (a) with multiple model (b), using: (c) standard 1D adaptation, (d) 1D complex unary filters in time-scale, and (e) standard 2D adaptation.

Figure 2 :

 2 Figure 2: Portion of common shot gather 1551. (a) Recorded data. (b) Multiple model.

Figure 3 :Figure 4 :

 34 Figure 3: Portion of first channel gather. (a) Recorded data. (b) Multiple model.

Figure 5 :

 5 Figure 5: Subtraction results on the first channel gather. (a) Standard 1D adaptive filter. (b) 1D complex unary filters in the time-scale domain. (c) Standard 2D adaptive filter in the time-space domain. See Figure 6 for a zoom on the delimited zone.

Figure 6 :Figure 7 :

 67 Figure 6: Zoom on traces between shot number 1400 and 1700 of Figure 5. (a) Recorded data. (b) Original model. Filtered data with (c) standard 1D adaptive filters, (d) 1D complex unary filters in the time-scale domain, and (e) standard 2D adaptive filters. P marks discrepancies on primary events.

Figure 8 :

 8 Figure 8: Near offset stack of the raw data. P and M mark primary and multiple events, respectively.

Figure 9 :

 9 Figure 9: Near offset stacks of filtered data using (a) standard 1D adaptive filters, (b) 1D complex unary filters in time-scale, and (c) standard 2D adaptive filters.

Figure 10 :

 10 Figure 10: Near offset stacks of the difference between the raw data and the filtered data using (a) standard 1D adaptive filters, (b) 1D complex unary filters in time-scale, and (c) standard 2D adaptive filters. P marks overadaptation to primary events.

Acknowledgments

The authors thank Statoil for allowing them to show the Norwegian Sea results. They also acknowledge IFP Energies nouvelles and CGGVeritas for the authorization to present this work. They are indebted to Anatoly Baumstein, Ramesh Neelamani, Tamas Nemeth and the anonymous reviewers for their constructive comments that helped improve significant parts of the paper. They finally recognize Bruno Lety (IFPEN) and Nicolas Quash (CGGVeritas) for their careful proofreading.


Jérémy Frey
email: jeremy.frey@inria.fr

Comparison of an open-hardware electroencephalography amplifier with medical grade device in brain-computer interface applications

Keywords: BCI, EEG, AMPLIFIERS COMPARISON, P300 SPELLER, WORKLOAD CLASSIFICATION

Comparison of an open-hardware electroencephalography amplifier with medical grade device in brain-computer interface applications

Jérémy Frey

INTRODUCTION

Brain-computer interfaces (BCI) are communication devices between humans and machines that rely only on brain activity (i.e. no muscular input) to issue commands or to monitor states [START_REF] Wolpaw | Braincomputer interfaces for communication and control[END_REF]. BCI is an emerging research area in Human-Computer Interaction that offers new opportunities for interaction, beyond standard input devices [START_REF] Tan | Brain-Computer Interfaces: Applying Our Minds to Human-Computer Interaction[END_REF]. In order to account for brain activity, portable and non invasive neuroimaging techniques are most commonly used, such as electroencephalography -EEG, which measures electrical current onto the scalp. The "interface" term covers many different areas of applications, for people with or without disabilities. However, while an increasing number of systems are being developed, from BCI aimed at controlling a cursor [START_REF] Wolpaw | Braincomputer interfaces for communication and control[END_REF] to adaptive systems (Zander and Kothe, 2011), more often than not the use of the technology is limited.

The price of the hardware is one of the main reasons that prevents the dissemination of non invasive BCI. Recently, more affordable EEG amplifiers ap-peared on the market, that could solve this issue. Among them, the OpenBCI board1 claims to bring BCI to the many. Enthusiasts and laboratories have started to use this board, but the quality of the recordings and the reliability of the resulting systems have yet to be assessed. In this study, we compare side by side the OpenBCI board with the g.tec g.USBamp amplifier2 , a device commonly used in BCI research. The price tag of the g.tec solution is around 20 thousands euros, 25 times more expensive than the 800 euros of 16 channels version of the OpenBCI board. Both OpenBCI and g.USBamp amplifiers can record up to 16 electrodes. This number of channels is sufficient to setup various BCI. We compared OpenBCI with g.USBamp for, on the one hand, a P300 speller application and, on the other hand, a workload monitoring application. Doing so, we could study respectively temporal and spectral features.

Note that here the question is not to assess which amplifier is the best device per se. Instead, we investigate if in a context of popular interactions -a nar-row scope compared to the possibilities that offers the g.USBamp -it is conceivable for researchers from the field or (well equipped) enthusiasts to make the leap. To which extend should we employ devices coming from the "DIY" ("Do It Yourself") community for actual BCI applications?

To answers this question, we adopted an approach somewhat different to what exists in the literature. Many papers deal with the comparison of electrodes, e.g. wet vs dry, with or without a conductive solution.

To do so, authors try to optimize the placement of both sets of sensors in order to get measures that originate from the same spots. However, no matter their efforts they could not merge sensors, and even clever montages, with electrodes of one sort positioned between electrodes of the other sort [START_REF] Tautan | Framework for evaluating EEG signal quality of dry electrode recordings[END_REF], are not ideal. It will produce a slight offset, hence a slight inaccuracy. Another alternative is to make separate measures by repeating the recordings with each system [START_REF] Nijboer | Usability of Three Electroencephalogram Headsets for Brain-Computer Interfaces: A Within Subject Comparison[END_REF], but once again the conditions could not be exactly the same.

In the present study we do not attempt to assess the quality of electrodes, but the behavior of amplifiers that are attached to them. Not a whole system, only the amplifiers. Therefore, we would not mind using the same electrodes during simultaneous recordings. This setup would ensure that the signal coming in each amplifier's inputs is exactly the same, avoiding any bias regarding the source of the measures.

We made that possible by crafting a dedicated adapter, one that basically splits in two the electrodes' wires. Such parallel measurement works because the amplifiers have high impedance circuits, that is to say that they are designed to not draw any amount of current from their source. As such, when one amplifier is connected, the readings of the other stay the same. Of course an infinite impedance cannot be achieved, and no matter the precautions this setup may cause a very slight difference compared to separate recordings. This is why in a second time we added to our adapter a circuit that prevents any interference between the two amplifiers, using ideal diodes to block current flows in one direction.

Using two different BCI applications, we investigated two types of EEG features. A task assessing workload aimed at assessing spectral information, and an oddball task sought temporal information. For each amplifier we measured the performance of a classifier based on those recordings, and additionally we compared both by correlating the signals that they recorded. No matter the financial aspects, the qualities of the g.USBamp amplifier make it the perfect baseline to gauge new challengers. This is also true for the electrodes developed by its manufacturer; in this study we are using g.tec wet and active (preamplified) electrodes.

FIRST EXPERIMENT: DIRECT CONNECTIONS

Experimental setup

We acquired 16 EEG channels using the active g.Ladybird electrodes from g.tec. In this system, the electrodes are attached to a box that powers their electrical components and retrieves the signal; the g.GAMMAbox. After studying the wiring of the g.GAMMAbox, we designed a printed circuit board (PCB) to connect both amplifiers. Our adapter plugs on one end to the D-sub 26 connector of the g.GAMMAbox. Thanks to a pinout composed of 2.54mm connectors that gave access to all the channels (16 EEG + reference + ground), we attached the OpenBCI board to the adapter -ground set to "bias" pin. On the other end of the adapter there was a Dsub 26 female connector, onto which we could plug the g.USBamp amplifier as if it were the regular end of the g.GAMMAbox. The schematics of the 2 layers PCB is presented in Figure 1. The EEG channels were positioned according to the 10-20 system at AFz, Fz, FCz, C3, C1, Cz, C2, C4, CPz, P3, Pz, P4, POz, O1, Oz and O2 -ground at FPz, reference on the left earlobe. Since the measures between both amplifiers were identical, only one recording session occurred, with one participant -there were no factors to counterbalance with repeated measures. The signals of both amplifier were acquired using OpenViBE 1.0.1 3 ; at 512Hz sampling rate for the g.USBamp and 125Hz sampling rate for the OpenBCI board.

The spectral features were investigated with a workload monitoring application, i.e. a BCI that is able to discriminate between several levels of mental effort. The system was trained using the N-back task, a well-known task to induce workload by playing on memory load and time pressure (see, e.g. [START_REF] Mühl | EEG-based workload estimation across affective contexts[END_REF]). The protocol we used was similar to [START_REF] Mühl | EEG-based workload estimation across affective contexts[END_REF], there were 360 trials presented during 6 blocks of alternate difficulty levels (0-back vs 2-back conditions). The recording session lasted approximately 12 minutes.

The temporal features were investigated using an oddball task directly implemented within OpenViBE with a visual P300 speller. P300 spellers are wellestablished BCI applications during which letters that randomly flash on the screen can be used to spell words with the sole brain activity. Indeed, when the letter that the user wants to spell flashes, a particular event-related potential (ERP) arises within the EEG, which possess a positive "peak" around t=300 ms after the stimulus onset. This is commonly referred to as the "oddball paradigm" since the occurrence of rare stimuli is used to elicit brain responses. During the recordings a matrix of 6 by 6 letters and digits was displayed in full screen on a 24-inch display. Only the calibration session occurred, during which one need to focus one's attention on a predefined sequence of letters. 32 letters composing a pangram were mentally "spelled" this way. The sentence was, without spaces, "pack my box with five dozen liquor jugs". Letters were flashing for 0.2s. There were 24 flashes per letter (12 times the row, 12 times the column), hence due to the matrix disposition there were in total 4608 trials, among which 768 were targets -"odd" trials, i.e. the letters of the target sentence were flashing. The recording session lasted approximately 30 minutes.

The acquisition of both amplifiers' signals and the P300 application occurred within the same Open-ViBE scenario (script). The recordings of each amplifier were synchronized with the appropriate events and exported in separate GDF files for later analyses. There was also only one scenario involved in the synchronization of all signals and events in the case of the N-back task; stimulation from the python script supporting this latter task were retrieved using the LSL protocol, a network protocol dedicated to physiological recordings which ensures accurate timings4 .

Signal processing

Two kinds of analyses were performed, using standard BCI signal processing pipelines. One aimed at assessing if and how the amplifiers differ in practice, when used for classification. The second then looked at the correlation between the acquired signals.

Classification

The signal processing of the data acquired during the N-back task is analogous to [START_REF] Mühl | EEG-based workload estimation across affective contexts[END_REF], i.e. 2s time windows, 5 frequency bands -delta (1-3 Hz), theta (4-6 Hz), alpha (7-13 Hz), beta (14-25 Hz) and gamma (26-40 Hz) -and spatial filters. We used common spatial patterns spatial filters to reduce the 16 channels to 6 "virtual" channels more discriminant between the workload conditions -see [START_REF] Ramoser | Optimal spatial filtering of single trial EEG during imagined hand movement[END_REF]. Additionally, we also tested a 3 frequency bands version of our pipeline, that consider only the lower frequencies, less prone to muscular artifacts -delta, theta and alpha.

Concerning the oddball task, we band-passed the signal between 0.5Hz and 40Hz, downsampled it by a factor 32 using the "decimate" Matlab function -by a factor 8 for OpenBCI because of the reduced sampling rate -, and applied regularized Eigen Fisher spatial filters -a spatial filter specifically designed for ERPs classification [START_REF] Hoffmann | Spatial filters for the classification of event-related potentials[END_REF] -to reduce channels' dimension from 16 to 5. We used 1s time windows after stimuli onsets -letters' flashesto epoch ("slice") our signal. However, in order to prevent data to overlap between consecutive stimuli due to the rapid pace of the flashes, after a first pass of epoching we discarded overlapping time windows from further analyses. This ensured that no part of the signal could be seen twice by the classifier between the training phase and the testing phase and bias the accuracy. The procedure was automatic, the first nonoverlapping epoch in order of appearance being kept. As a result, in the end we obtained 48 target trials and 240 distractor trials for classification, identical between the g.USBamp and the OpenBCI recordings.

Both for the workload and the P300 speller tasks, we used shrinkage LDA (linear discriminant analysis) for classification [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF]. To assess the classifiers' performance on the calibration data, we used 4-fold cross-validation. I.e. we split the collected data into 4 parts of equal size, used 3 parts to calibrate the classifiers and tested the resulting classifiers on the unseen data from the remaining part. This process occurred 3 more times so that in the end each part was used once as test data. Finally, we averaged the obtained classification accuracies. The accuracy was measured using the area under the receiver-operating characteristic curve (AUROCC). The AU-ROCC is a metric that is robust against unbalanced classes, as it is the case with oddball tasks. A score of "1" means a perfect classification, a score of "0.5" is chance. In order to make statistical comparisons between both amplifiers for each type of features that we studied, we ran 10 times the analyses -the trials were selected randomly for cross-validation.

Correlations

We compared, on the one hand, the frequency spectra associated to the different workload conditions and, on the other hand, the time course of the ERP that were caused by the flashing target letters. To do so, we used Pearson correlations, on par with the literature for similar analyses -e.g. (Zander et al., 2011). In order to ensure a 1-to-1 correspondence between our sets of data, the recordings from the g.USBamp were downsampled to 125Hz -same sampling rate as for the OpenBCI -using the "resample" function from Matlab R2014a signal processing toolbox.

Concerning the workload task, we first aggregated the 2s time-windows corresponding to each condition (0-back and 2-back). Then we used the "spectopo" function of the EEGLAB toolbox (version 13.4.4b) to compute the grand average power spectral between 1Hz and 40Hz, for each channel. The output of the function was then passed on to R (version 3.0.2) to compute correlations through the "rcorr" function from the "Hmisc" package.

For the oddball task, we first band-passed the signals between 1Hz and 8Hz -the approximate frequency band used for classification. Then we extracted time epochs starting 0.5s prior to the flashing of the target letters and ending 1s after stimuli onset. Contrary to what occurred for classification, we did not prune overlapping epochs in the oddball task when we compute the averaged ERP -there was no bias that could have been induced here. Finally, we averaged the ERP per channel before exporting the time points to the R environment.

Results

Classification

The results regarding classification accuracy are presented in Table 1, with the AUROCC scores for each one of the 10 repetitions, for both amplifiers and both tasks -including the 3 and 5 frequency bands pipeline for workload.

We tested for significance using Wilcoxon signedrank tests. There was a significant difference between amplifiers for the P300 tasks (p < 0.01). The AUROCC mean score for the g.USBamp was 0.961 vs 0.918 for the OpenBCI. There were however no significance but tendencies concerning the workload task, with mean AUROCC scores between 0.85 and 0.86 for the 3 bands pipeline (p = 0.095) and between 0.89 and 0.90 for the 5 bands task (p = 0.079) -see Table 1 for details.

Correlations

When we first analyzed our data to seek correlations regarding the oddball tasks, we realized that a shift occurred during the recordings, as denoted in Figure 2 by the grand average of the ERP for target trials across channels. This may have been caused by a software issue (see Discussion). In order to correct the shift and conduct proper comparisons between both amplifiers' measures, we used a cross-correlation to estimate the time shift, using the "ccf" function from the R "stats" package. We found a delay of 88ms between the two signals -11 data points at 125Hz, see Figure 2.

In Figure 3a, the averaged ERP were shifted by as much for each channel. Corresponding Pearson correlation R scores, that were computed using the "rcorr" function, are presented in Table 2. The mean R score is 0.9965 and is statistically significant (p < 0.001).

There was also a significant correlation (p < 0.001) for the spectral features, with a mean R score of 0.9983 for the 0-back condition and 0.9979 the 2-back condition (see Table 2 for details). Among the brain signals patterns that could be expected during the completion of a difficult task, the decrease nearby the alpha frequency band during the 2-task condition can be observed within per-channel spectra presented in Figures 3b and3c. Note that we did not correct time shifts prior to workload analyses due to the nature of the features -i.e. spectral and not temporal.

Discussion

The correlation between both temporal and spectral features tends so show that the signals acquired by the g.USBamp and the OpenBCI are, if not identical, very closely related. For every condition and channel tested, the Pearson R score was greater than 0.99.

There were however more dissimilarities in the classification accuracy obtained during the corresponding tasks. While there were hardly a difference between the AUROCC scores computed from both amplifiers with the N-back tasks, the g.USBamp performed significantly better than the OpenBCI during the P300 speller task. The time shift observed afterwards between the two amplifiers may partially explain this difference. Indeed, the detection of ERP is ) particularly sensitive to signals' latency, and a shift between events' timestamp and signal's acquisition could result in such degradation of performance when temporal features are involved.

AFz Fz FCz C3 C1 Cz C2 C4 CPz P3 Pz P4 POz O1 Oz O2 -2 0 2 4 -2 -1 0 1 2 -3 -2 -1 0 1 2 -2 -1 0 1 2 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 -2 -1 0 1 -2 -1 0 1 -2 -1 0 1 2 -3 -2 -1 0 1 -1 0 1 -2 0 2 -2 -1 0 1 2 -2 0
The radio transmission between the wireless OpenBCI board and the dongle plugged to the computer may be one of the cause of the situation. The problem could also originate from the software. As a matter of fact, the OpenViBE acquisition driver of the OpenBCI board was released no so long before our experiment, and was still labelled as "unstable" as for version 1.0.1 of the software. One "oddity" that may further highlight the youth of OpenBCI software integration: we realized during our analysis that the recorded signals were completely inverted on the Y axis. The voltage reported by the board were the opposite of what g.USBamp was claiming. Since on numerous occasions we acknowledged the accuracy of g.tec devices readings, it is the OpenBCI's signals that we inverted back to "normal" prior to correlation analyses.

Beside time shifts issues, as mentioned during the introduction we needed to strengthen those first insights by discarding the eventuality that both EEG signals may have influence each other due to the direct wiring with the electrodes.

SECOND EXPERIMENT: ISOLATED CONNECTIONS

The second set of recordings is very similar to the what was described during the first study. The second experiment only differs by the nature of the adapter that was employed. As such we will only discuss the changes that were made to the hardware and quickly dive into the results.

Ideal adapter

We modified the adapter that connects the amplifiers to the g.GAMMAbox -and by extent to the EEG electrodes. Instead of a direct connection between each amplifier's inputs and the EEG channels, we interposed "ideal" (or "super") diodes on the branches of the "Y" wiring. Diodes are electrical components that let the current flow in only one direction, the "forward" direction. Hence, this type of montage ensures that no current could travel directly from one amplifier to the other, contaminating the recordings. However, regular diodes cause a voltage drop. The voltage drop varies depending on the materials used for their construction, but it is at least 0.3V. Meaning that if the current coming in the forward direction is lesser than 0.3V, no signal will pass through. 0.3V is an order of Figure 4: Schematics of the adapter with the ideal diodes montage. Note that there is a set of ideal diodes on the ground channel, but they were shorted with a jumper during our experiment. magnitude superior to the range of EEG signal -approximately a thousand time, therefore regular diode could not be used.

To circumvent this problem, we utilized a particular montage that involved operational amplifiers (opamp). Op-amps are components widely used in electrical circuits, acting as sorts of "building blocks". Notably, in combination with a regular diode, one could use a precision rectifier configuration to obtain an "ideal" diode. This particular montage is also known as a "super" diode, since there will always be a slight voltage drop, but in this case, thanks to the gain of the op-amp, it becomes negligible.

We mounted 36 of such ideal diodes on the adapter. One on each end of the "Y" section associated to the 16 EEG channels, plus 2 for the reference. Due to the nature of the electrical recordings, only the ground was left without such circuit. We utilized Texas Instrument op-amps, model TLC2272ACPE4. The TLC227xA series are more indicated for precision application, and with 2 op-amps per chip we could limit the overall size of the adapter. The operational amplifiers were powered by an external circuit with regulated -2.5 / +2.5 voltage. The schematics of the adapter -also a 2 layers 2 layers PCBis presented in Figure 4. The ideal diodes montage, placed before amplifiers' inputs, prevented any current to flow in reverse direction from either amplifier to the adapter; it ensured that one set of recordings would not bias the other. One recording session occurred for each application and each condition.

Results

The signal processing and the analyses were strictly identical to the first experiment detailed above, refer to the previous section for related information.

Classification

As with the first study, the results regarding classification accuracy are presented in Table 3, with the AU-ROCC scores for each one of the 10 repetitions, for both amplifiers and both tasks -including the 3 and 5 frequency bands pipeline for workload. We tested for significance using Wilcoxon signed-rank tests. No matter the task there was no significant difference, although the 5% threshold was nearly reached for spectral features -p-value was 0.157 for the P300 task, 0.051 for the 3 bands version of the workload pipeline and 0.286 for the 5 bands version.

Correlations

Concerning the P300 oddball task, there was a offset of 88ms as well between the recordings of both amplifier with the isolated connection -see Figure 5 for the grand ERP average and the cross correlation. The perchannel averaged ERP are plotted in figure 6a. Corresponding Pearson correlation R scores are presented in Table 4. The mean R score is 0.8847 and is statistically significant (p < 0.001).

There was also a significant correlation (p < 0.001) for the spectral features, with a mean R score of 0.9976 for the 0-back condition and 0.9987 the 2-back condition (see Table 4 for details). The per-channel spectra are presented in Figures 6b and6c. As with the direct connection, the band frequency changes between the 0-back and the 2-back conditions can be observed in the spectra.

Discussion

The results with the isolated connections are not that different from what was obtained during the first experiment. This would suggest that directly connecting two high impedance amplifiers to the same EEG electrodes could be a viable montage for a side-by-side comparison.

Since with both types of connector there was only one set of recordings, we could not draw any conclusion about the lower classification accuracy obtained with the isolated montage. The vigilance level of the participant alone could explain these performances.

Thanks to signals' correlations, however, we may infer that noise was added to the system due to the presence of additional electrical components in the adapter. Indeed, while the spectra were once again strongly correlated, the averaged ERP achieved "only" a mean R score of 0.88. Here external factors such as the metal state of the participant or the quality of electrodes contacts could not have influenced one amplifier rather than the other. Since temporal features are more sensitive than spectral features to signal quality -e.g. one "peak" in the signal vs oscillatory patterns over several seconds -, it is instead more plausible that the difference with the first experiment comes from the adapter.

Nonetheless, even though the ideal diode montage did not produce ideal signals, those results still advocate for a close proximity between the g.USBamp and the OpenBCI. No device behaved "better" than the other, because no matter the lower correlation between averaged ERP, the classification accuracy is in practice comparable between both amplifiers. Each one probably endured different fluctuations since each had a dedicated set of ideal diodes.

CONCLUSION

During this preliminary study, we compared the OpenBCI board to the g.tec g.USBamp amplifier. We employed an original montage, based on the simultaneous recording of the same set of electrodes. While as a first approach we used a simple adapter with a direct connection between the amplifiers and the electrodes, in a second experiment we attempted to dis-card any possible interference that one amplifier could cause to the other.

To do so, we built an adapter that embedded "ideal" diodes, components that prevented electrical currents to flow "backward". This ensured that we could test both devices in isolation. We did not try to compare both adapters as the purpose was simply to gather more insights about the possibility of simultaneous recordings -this was a precaution to detect a possible bias. For all applications and conditions AUROCC scores were far beyond chancel level; the OpenBCI amplifier came close to the g.USBamp in terms EEG features and effective performance.

That is not to say that the OpenBCI could replace an equipment such as the g.USBamp, though. For instance, this open-hardware initiative does not aim at medical applications, hence it should be employed in sensitive contexts. It does not possess any certification; one reason why so many cheap EEG device are wireless is not for practicality, but to avoid any hazard due to power supply. Connecting somehow a body to the power grid requires extra precautions and a certified isolation, moreover when the impedance between the electrodes and the brain is intentionally lowered.

There are also few issues with the current state of the OpenBCI project. One concerns the sampling rate of the board. While 125Hz may be enough for our use-cases -no frequencies beyond 40Hz were used during the work presented here -it may not suffice others. The limitation of the sampling rate is caused by the wireless protocol used for data transmission. OpenBCI can deliver up to 250Hz signals to the computer, but only on 8 channels instead of 16. Note that this may be optimized in the future by updating the firmware or using alternate communications -as far as the board itself is concerned, the documentation of the ADS1299 chip ensuring analog-to-digital conversion claims a sampling rate up to 16,000Hz.

With those limitations in mind, overall the results suggest that the OpenBCI board -or a similar solution based also on the Texas Instrument ADS1299 chipcould be an alternative to traditional EEG amplifiers. Even though medical grade equipment possesses certification and still outperforms the OpenBCI board in terms of classification, the latter gives very close EEG readings. In practice, the obtained classification accuracy may be suitable for reliable BCI in popular settings, widening the realm of applications and increasing the number of potential users.

Figure 1 :

 1 Figure 1: Schematics of the direct adapter between electrodes and amplifiers.

Figure 2 :

 2 Figure 2: Left: Averaged ERP across channels of the target trials during the oddball task, before time shift correction. Right: Cross-correlation between the amplifiers. The computed lag of 11 data points corresponds to 88ms. (Direct connection.)

Figure 3 :

 3 Figure 3: Direct connection: averaged ERP for the target trials of the oddball task (a), averaged spectra for the 0-back (b) and 2-back (c) trials of the N-bak task.

Figure 5 :Figure 6 :

 56 Figure 5: Left: Averaged ERP across channels of the target trials during the oddball task, before time shift correction. Right: Cross-correlation between the amplifiers. The computed lag of 11 data points corresponds to 88ms. (Isolated connection.)

Table 1 :

 1 Classification accuracy (AUROCC scores) for the P300 and workload tasks studied during the first experiment -direct connection between the electrodes and the amplifiers. The 4-fold cross validations were repeated 10 times. Two pipelines are presented for the workload: 3 frequency bands ("WL 3b", δ + θ + α) as well as 5 frequency bands pipeline ("WL 5b", δ + θ + α + β + γ). Significance was tested using Wilcoxon signed-rank tests.

			1	2	3	4	5	6	7	8	9	10	Mean SD
	P300	g.USBamp 0.96 0.96 0.96 0.96 0.96 0.96 0.96 0.96 0.96 0.97 0.961 0.003
		OpenBCI	0.92 0.92 0.91 0.91 0.92 0.92 0.93 0.92 0.92 0.91 0.918 0.006
	WL 3b	g.USBamp 0.85 0.85 0.85 0.86 0.85 0.87 0.87 0.87 0.87 0.86 0.860 0.009
		OpenBCI	0.86 0.86 0.85 0.86 0.84 0.86 0.85 0.85 0.85 0.85 0.853 0.007
	WL 5b	g.USBamp 0.90 0.89 0.90 0.89 0.91 0.90 0.89 0.90 0.89 0.90 0.897 0.007
		OpenBCI	0.91 0.89 0.87 0.88 0.89 0.88 0.89 0.89 0.89 0.90 0.889 0.011
			AFz	Fz	FCz	C3	C1		Cz	C2	C4	CPz
	P300 target	0.998 0.997 0.997 0.997 0.997 0.994 0.996 0.992	0.994
	Workload 0-back 0.999 0.999 0.998 0.998 0.998 0.998 0.998 0.999	0.998
	Workload 2-back 0.999 0.999 0.998 0.998 0.998 0.998 0.998 0.999	0.997
			P3	Pz	P4	POz	O1		Oz	O2	Mean	SD
	P300 target	0.995 0.994 0.996 0.996 0.996 0.994 0.995 0.9965 0.0015
	Workload 0-back 0.998 0.998 0.998 0.999 0.998 0.998 0.998 0.9983 0.0003
	Workload 2-back 0.998 0.998 0.998 0.998 0.997 0.998 0.998 0.9979 0.0005

Table 2 :

 2 Pearson correlation R scores between g.USBamp and OpenBCI recordings at the 16 different electrode locations with a direct connection. The "P300 target" condition corresponds to temporal features (ERP averaged across trials) and the workload conditions to spectral features.

Table 3 :

 3 Classification accuracy (AUROCC scores) for the P300 and workload tasks studied during the second experiment -isolated connection between the electrodes and the amplifiers. The 4-fold cross validations were repeated 10 times. Two pipelines are presented for the workload: 3 frequency bands ("WL 3b", δ + θ + α) as well as 5 frequency bands pipeline ("WL 5b", δ + θ + α + β + γ). Significance was tested using Wilcoxon signed-rank tests.

		AFz	Fz	FCz	C3	C1	Cz	C2	C4	CPz
	P300 target	0.976 0.934 0.892 0.846 0.872 0.881 0.838 0.811	0.912
	Workload 0-back 0.999 0.998 0.998 0.997 0.998 0.998 0.998 0.996	0.998
	Workload 2-back 0.999 0.999 0.999 0.998 0.998 0.999 0.999 0.999	0.999
		P3	Pz	P4	POz	O1	Oz	O2	Mean	SD
	P300 target	0.849 0.823 0.896 0.910 0.878 0.982 0.853 0.8847 0.0483
	Workload 0-back 0.998 0.998 0.998 0.997 0.997 0.997 0.997 0.9976 0.0007
	Workload 2-back 0.999 0.999 0.999 0.999 0.999 0.998 0.998 0.9987 0.0004

Table 4 :

 4 Pearson correlation R scores between g.USBamp and OpenBCI recordings at the 16 different electrode locations with an isolated connection. The "P300 target" condition corresponds to temporal features (ERP averaged across trials) and the workload conditions to spectral features.

http://www.openbci.com/

http://www.gtec.at/

http://openvibe.inria.fr/

https://github.com/sccn/labstreaminglayer

ACKNOWLEDGMENTS

I thank Thibault Laine for his technical assistance during this project.