C Siviero

Member, IEEE R Trinchero

Fellow, IEEE, I. S. Stievano S Grivet-Talocia

Senior Member, IEEE G Signorini

Member, IEEE M Telescu

Constructive Signal Approximations for Fast Transient Simulation of Coupled Channels

Keywords: Modeling, digital integrated circuits, LVDS buffers, differential drivers, signal integrity, high-speed signaling, equalization, pre-emphasis, de-emphasis, feed-forward equalization

This paper addresses the problem of fast transient simulation of high-speed coupled channels driven by transceivers that include Transmitter Feed-Forward Equalization (TX-FFE). Practical circuit hardware implementations of such drivers may not be correctly represented by industry standard algorithmic modeling interfaces, which are usually based on ideal digital finite impulse response filters. In particular, slow transients may arise when aggressive FFE is activated and particular switching sequences significantly stress local voltage regulators integrated in the output buffers. In this paper, we model this effect using a constructive hierarchical approximation of the analog waveforms that form a linearized vector Thevenin or Norton representation of the driver. This representation enables the accurate simulation of coupled channels including both differential and common-mode signals. Analog waveforms at the receiver input are readily obtained by standard time-frequency transformations, while including channel characteristics and any additional linear equalization block, if present. This proposed framework is compatible with industry-standard Algorithmic Modeling Interfaces (AMI), extending their scope to coupled channels with non-ideal TX-FFE circuit blocks. Eye diagrams and transient waveforms corresponding to millions of bits are computed in few seconds using a non-optimized software implementation. Random jitter and crosstalk are seemlessly included as in standard AMI frameworks. The results of various numerical experiments on commercial transceivers are reported in this paper, confirming both the accuracy and the efficiency of the proposed framework.

I. INTRODUCTION

The design of high performance electronic equipment, such as smartphones, tablets and PCs, requires a careful characterization of the non-ideal behavior of the entire end-to-end electrical interconnects, in order to evaluate any detrimental effect on the overall system reliability. This is even more important in the most recent and modern devices, whose high-speed wired communication interfaces (e.g., between the application processors and other specialized integrated circuits or pheripherals) have been constantly improved over the years in order to support higher and higher data rates. Faster signaling inevitably leads to more likely and stronger C. Siviero, R. Trinchero, S. Grivet-Talocia and I. S. Stievano are with the EMC Group, Department of Electronics and Telecommunications, Politecnico di Torino, 10129 Torino, Italy (e-mail: claudio.siviero@polito.it, riccardo.trinchero@polito.it, stefamo.grivet@polito.it, igor.stievano@polito.it).

G. Signorini is with Intel Corporation, Munich, Germany. M. Telescu is with Université de Bretagne Occidentale; CNRS, UMR 6285 Lab-STICC; Brest, France. distortion effects due to interconnect nonidealities, like losses, attenuation, dispersion, impedance mismatches and return path discontinuities, to name just a few. Such effects can significantly affect the reliability of hardware implementations, whose first-pass quality has become a must in order to comply with very short time-to-market requirements.

In this context, a key role is played by full-channel analog simulations, aimed at assessing the link signal and power integrity [START_REF] Swaminathan | Designing and Modeling for Power Integrity[END_REF], [START_REF] Fan | Signal integrity Design for High-Speed Digital Circuits: Progress and Directions[END_REF]. For data rates in Gb/s range (e.g., 8, 10, 16 Gb/s or more), attenuation and dispertion introduced by interconnets can be partly compensated using equalizers in high-speed transceivers. For these devices, macromodel-based simulations have been proven to be the only viable solution, since the computational burden required by a full-detail transistor-level description cannot be handled in realistic simulation times. Macromodels avoid this problem by replacing the detailed models with much more compact and efficient surrogates, without significant loss in accuracy.

In the past twenty years many efforts have been spent by both industry and academia to provide efficient macromodeling solutions with tunable accuracy, ranging from the standard approach based on the Input/Output Buffer Information Specification (IBIS) [START_REF] Sharawi | Circuit Modeling in High-Speed Designs[END_REF], [START_REF]I/O Buffer Information Specification version 5[END_REF] to more elaborate techniques based on either circuit-inspired equivalents of blind blackbox representations [START_REF] Stievano | Reliable eyediagram analysis of data links via device macromodels[END_REF]- [START_REF] Dghais | Reduced-Order Parametric Behavioral Model for Digital Buffers/Drivers With Physical Support[END_REF]. For the latter set, more recent enhancements and applications are available in [START_REF] Zhu | Accurate and Scalable IO Buffer Macromodel Based on Surrogate Modeling[END_REF]- [START_REF] Dghais | Mixed-Signal Overclocked I/O Buffers Model Abstraction for Signal Integrity Assessment[END_REF].

All the above approaches, however, share the same limitation of dealing with nonlinear model representations and even if the proposed surrogates are more the one order of magnitude times faster than transistor-level models, the efficiency improvement is not enough to process millions of bits in a reasonable amount of time. The only available solution widely recognized by the industry and which has been successfully implemented in most of the commercial state-of-the-art tools is represented by the Algorithmic Modeling Interface (AMI) extension of IBIS. IBIS-AMI offers a linear modeling and simulation tool where Transmitter Feed-Forward Equalization (TX-FFE) is accounted for by a suitable post-processing of the differential signal flowing through the communication channel. Even if the simulation speed allowed by this technique is appropriate for modern CAD flows, the limited accuracy in reproducing the rich dynamical effects of recent transceivers, including the common-mode signal in coupled channels, demand for improvements.

This paper extends the preliminary results presented in [START_REF] Siviero | A Fast Channel Simulation Framework Based on Hierarchical Waveform Representations[END_REF] and discusses a novel high-speed serial link simulation method based on a different linear modeling framework. The proposed approach aims at evaluating the transient waveforms at the receiver side of the channel, for arbitrarily long bit sequences, and in just a few seconds. Both differential-and commonmode waveforms are predicted with a tunable accuracy depending on the number of hierarchical basis functions, which are introduced in the model to account for possible broadband and even long-term memory effects caused by switching mechanisms within the transceiver circuits. This phenomenon may be hard to model within the IBIS-AMI framework, which relies on time invariance since based on weighted precursors and postcursors that are just rescaled versions of the same elementary pulse. The proposed signal representation overcomes this limitation through successive hierarchical refinements that are shown to converge to the reference responses.

The proposed signal construction process is nonlinear and adaptive, depending on the particular binary sequence being transmitted. However, linear superposition fits in the analog part of the signal processing flow, thus enabling to compute the transient responses via either frequency-domain (e.g., Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform (IFFT)) or time-domain approaches (e.g., recursive convolution). We adopt the former frequency-domain approach, which naturally plugs into the existing IBIS-AMI framework by suitably extending its scope and applicability.

The paper is then organized as follows. Section II derives the proposed model representation for differential transceivers from the well known nonlinear model structures used so far in the literature. Section III presents the step-by-step procedure for model generation, via a set of simple simulations of the reference transistor-level models of devices connected to predefined loads. Model extension involving the above mentioned hierarchical bases for TX-FFE representation is illustrated in Sec. IV. Section V discusses the fast channel simulation via a mixed time and frequency-domain approach and Sec. VI collects the numerical results based on the application of the proposed technique to two commercial devices, demonstrating its feasibility and strength. Final remarks and conclusions are given in Section VII. Our proposed modeling framework is based on a generic multiport (differential) driver structure, as depicted in Fig. 1, where the individual single-ended output currents i 1,2 and voltages v 1,2 are defined. These are readily converted, if desired, into differential and common-mode quantities, as

II. MODELING MULTIPORT TRANSCEIVERS

v d = v 1 -v 2 , v c = v 1 + v 2 2 . (1)
Similarly, conversion to single-ended and/or differential scattering representations in terms of incident and reflected transient wave signals can be considered too, leading to a pure scattering formulation as in [START_REF] Olivadese | Macromodel-Based Iterative Solvers for Simulation of High-Speed Links With Nonlinear Terminations[END_REF]. In the following, we base our derivations to single-ended voltages and currents for simplicity of illustration, but any other representation that can be linearly related to such variables can be adopted as a starting point as well. Our numerical implementation is based on single-ended scattering variables for improved numerical robustness.

With reference to Fig. 1, we propose the following model structure

i(t) = -Gv(t) + j(t) + g(v(t)) (2)
where

vectors i = [i 1 , i 2 ] T and v = [v 1 , v 2 ]
T collect the port variables. In (2) we recognize various terms as follows:

• the static behavior of the output port currents is captured by a constant conductance matrix G; • the dynamic behavior of the driver when excited at its output by the voltages v is represented by a dynamical Linear Time-Invariant (LTI) submodel g; • the time-varying source term j(t) = [j 1 (t), j 2 (t)] T represents the switching activity of the driver as a time-varying analog waveform. This switching activity results from the particular bit pattern imposed by the inner stages of the driver, which are represented in Fig. 1 by the voltage v in (t). The model structure ( 2) is indeed a particular case of two well-established formats for behavioral transceiver representation, namel the IBIS standard and the MpiLog framework [START_REF]I/O Buffer Information Specification version 5[END_REF], [START_REF] Stievano | Mπlog, macromodeling via parametric identification of logic gates[END_REF], [START_REF] Signorini | Present and Future of I/O-Buffer Behavioral Macromodels[END_REF], [START_REF] Signorini | Macromodeling of I/O Buffers via Compressed Tensor Representations and Rational Approximations[END_REF], see also [START_REF] Diouf | Simplified topology for IC buffer behavioural models[END_REF]. We recall that basic Mpilog models are written as time-varying combinations of two-piece dynamical submodels as

i(t) = w H (t)F H (v(t)) + w L (t)F L (v(t)) (3)
where F ν with ν ∈ {H, L} are nonlinear dynamical multivariate relations accounting for the static and the dynamic behavior of the driver in fixed logic states, and w ν are time varying switching functions. A very similar structure holds for IBIS [START_REF]I/O Buffer Information Specification version 5[END_REF], with some simplifications. The proposed model structure (2) can be readily derived from (3) under the following assumptions:

• the two switching functions describing the analog transitions L → H and H → L are constrained to be symmetric through w H (t) = 1-w L (t), which effectively represents switching behavior through a single timedomain signal; • the F ν submodels, which can be highly nonlinear and dynamic in some applications [START_REF] Stievano | Mπlog, macromodeling via parametric identification of logic gates[END_REF], [START_REF] Mutnury | Macromodeling of nonlinear digital I/O drivers[END_REF], [START_REF] Cao | A New Training Approach for Robust Recurrent Neural-Network Modeling of Nonlinear Circuits[END_REF], are here approximated by LTI compact models through

F ν (v) = -G ν v + I ν + g ν (v), ν ∈ {H, L} (4)
where the conductance matrices G ν and the static bias current vectors I ν define linear multivariate static characteristics (in Norton form) approximating the static behavior of the device in each logic state (see Fig. 2 (top panel), which considers the behavior of a commercial 90-nm low-power driver available as a transistor-level encrypted netlist [21]); we denote with g ν the two LTI submodels that reproduce the dynamical response of the device. Using the two above simplifications in (3) leads to

i(t) = w H (t) [-G H v + I H + g H (v)] + +(1 -w H (t)) [-G L v + I L + g L (v)] (5)
Introducing now the additional assumption that the driver is symmetric in both its static and dynamic parts, the representation in ( 5) can be recast as in [START_REF] Fan | Signal integrity Design for High-Speed Digital Circuits: Progress and Directions[END_REF], where:

G = G H = G L , g = g H = g L (6)
and

j(t) = [w H (t)I H + (1 -w H (t))I L ] . (7)
Assumption ( 6) may appear overly restrictive, but extensive numerical tests on various benchmark cases demonstrated that the accuracy loss is negligible when modeling high-speed differential drivers. In particular, this assumption is crucial for ensuring computational efficiency.

III. CONSTRUCTING TRANSCEIVER MODELS

Model parameters can be readily obtained via a simple estimation procedure based on a set of Transistor-Level (TL) simulations. Specifically, the 2 × 2 matrix G and the vectors I ν , defining the approximated linear static characteristics are computed via simple least-squares fits from the real static multivariate surfaces of the device response when operating in a fixed logic state. These responses are extracted by a set of double DC sweeps computed, e.g., in SPICE. The proposed linear approximation corresponds to two planes, depicted for a test driver in Fig. 2 (top panel) in gray solid color. Note that these two planes corresponding to H and L static characteristics are parallel, since we enforce G H = G L = G, while the DC current offsets I H and I L define the reference bias levels of each logic state. The figure shows that the approximation of the actual DC characteristics with parallel planes is excellent, as highlighted by the error surfaces shown in the bottom panel of Fig. 2. The error in magnitude turns out to be always less than 0.5 mA in the full considered region of the port voltages v 1 and v 2 . In particular, the above mentioned error is one order of magnitude less (i.e., < 50 µA) in the zone explored during normal operation of the device (see the black stars in the top panel), thus enabling the proposed modeling assumption.

The switching sources j(t) are instead computed from (2) from a set of current and voltage responses of the driver switching on a given load, according to the setup in the top panel of Fig. 3. The device is driven to produce a set of complete transition events, i.e., a "010" bit pattern, allowing to generate the basic switching sources associated to the corresponding up or down transitions. The load is chosen to cancel the effects of the dynamical contribution g. For the i-v representation (2), the above mentioned load corresponds to two ideal short circuits connected between the positive and negative output terminals of the device and ground, forcing v 1 = 0 and v 2 = 0. For the scattering representation which is indeed adopted in the present work, two single-ended 50 Ω resistors are considered instead. The switching term j(t) plays in our framework the same role of the weighting functions w ν (t) in standard IBIS/Mpilog models, providing in particular an effective behavioral representation of the TX-FFE blocks. Fig. 2. Top Panel: trajectory of the current response i 1 (t) of an example driver during a real simulation test with a lossy channel (black stars). The trajectory is superimposed to the real and the approximated static characteristics of the driver in the fixed logic H or L states, i.e., i sH and i sL , respectively (gray color corresponds to the planes of the linear approximation). Bottom panel: error computed between the actual and the approximated static characteristics at the H (blue surface) and L (red surface) fixed logic state.

For the dynamical submodel g, a multiport state-space structure is assumed, with parameters computed through a Time-Domain Vector Fitting (TDVF) process [START_REF] Grivet-Talocia | Package Macromodeling via Time-Domain Vector Fitting[END_REF] applied to the transient responses of the currents i 1 (t) and i 2 (t) and voltages v 1 (t) and v 2 (t), collected when the driver switches on a transmission line load (see the testbench shown in the bottom panel of Fig. 3, where the far-end load consists of two resistors connected to an ideal voltage source forcing the nominal common mode voltage V C,nom ). We remark that the TDVF algorithm is able to identify a LTI system from its transient responses excited by arbitrary inputs. This feature is particularly useful here, since both inputs and outputs used for the identification of submodel g are known only implicitly and are obtained by extracting port voltages and currents from the results of a dedicated transient analysis testbench. As an example, Fig. 4 shows the transient voltages and current i 1 (t) for the above mentioned test case. Model prediction is included as well, highlighting its good accuracy in reproducing the rich dynamical behavior of the driver.

Differently from IBIS and Mpilog [START_REF]I/O Buffer Information Specification version 5[END_REF], [START_REF] Signorini | Macromodeling of I/O Buffers via Compressed Tensor Representations and Rational Approximations[END_REF], the proposed 3). Bottom panel: dynamical part of the response of the current i 1 (t) for the same test case, compared with the response of the dynamical submodel g for the first terminal current. The above waveforms, toghether with the companion response of i 2 (t) are used for the identification of g.

model is completely linear, thus enabling frequency-domain methods for fast simulation. Differently from the widelyadopted IBIS-AMI standard (which is also based on linearity), the proposed model is better grounded in circuit theory and natively multivariate, thus allowing for an accurate representation of both the differential and the common mode signals.

IV. MODELING LONG-TERM MEMORY EFFECTS OF TX-FFE BLOCKS

We now consider embedding an appropriate representation of FFE blocks in the proposed driver model structure. Essentially, TX-FFE circuitry boost the driver strength when switching occurs, in order to induce a high-pass behavior that should compensate the inevitable low-pass characteristic of the channel due to metal and dielectric losses. Modern CAD tools for channel simulation and eye diagram construction embed effective and simple-to-use FFE blocks, which are commonly represented as Finite-Impulse-Responses (FIR) filters. Given a single-bit pulse p(t) that would be applied without FFE activated, the enhanced pulse with FFE is represented as

p(t) = r r=-r α r p(t -rT B ), (8)
where T B is the bit time, α r are suitable weights used to design the FIR frequency response, and r, r denote the number of precursors and postcursors, expressed as translated and rescaled identical replicas of the main pulse p(t). When the boosted elementary pulse p(t) is used for transmitting a given bit pattern {B k ∈ {0, 1}, k = 1, 2, . . . }, the resulting analog waveform is still expressed as a superposition of identical replicas of the initial pulse p(t),

k B k p(t -kT B ) = k r r=-r B k α r p(t -(r + k)T B ) (9)
as a consequence of the assumed translation invariance. Evidence collected from several numerical experiments challenges the above representation of the pre-emhpasis behavior. Let us consider Fig. 5, where the switching patterns of the example driver is depicted with and without its (1-tap) FFE activated. The analog waveform without FFE is welldescribed by a superposition of switching events with the same analog rising and falling fronts and can thus be represented as a weighted superposition of successive translations of a single and well-define pulse shape p(t). Conversely, when FFE is enabled, we observe that the expected boosting effect is characterized by some slow transients: the maximum peakto-peak amplitude of the waveforms is obtained after a few successive switchings, with a slow dynamic saturation effect. This effect makes each individual switching front different from each other, based on the number of earlier consecutive switching events. This effects has a well-defined extent in time, which can be estimated in Fig. 5 to last for about 4-5 bit times T B for this particular device. We conclude that a linear FIR approximation as in [START_REF] Signorini | Present and Future of I/O-Buffer Behavioral Macromodels[END_REF] is not adequate for this driver, since the shape of the switching fronts is not translationinvariant. What is instead required is a nonlinear description of the TX-FFE block, which is able to reproduce the switching front dependence on the number of immediately preceding consecutive switchings.

A. Hierarchical functions for FFE representation

The non-linear and adaptive representation of analog switching fronts that depend on the particular bit pattern is here obtained by a simple constructive approximation process. The main idea shares a common background with widely adopted adaptive approximation framework as those enabled by wavelet decompositions or image and video compression. An initial coarse-level approximation is derived, which provides a gross approximation of the response. Successive details are then added where a correction is required, in order to refine approximation. Corrections are added hierarchically through levels, until a sufficiently small approximation error is attained. The main enabling factor for proposed approach is the finite memory length in time of the TX-FFE analog behavior (see Fig. 5), which in turn induces sparsity, so that corrections need to be applied only locally and with a sparsity pattern that increases with increasing hierarchical levels.

Let us consider each of the two driver source components (n = 1, 2) in ( 7), which we detrend from their initial condition as ξ n (t) = j n (t) -j n (0), and which we represent through the following sparse hierarchical superposition

ξ n (t) = L =0 k∈Ω ( ) n,u ϕ ( ) n,u (t -kT B ) + k∈Ω ( ) n,d ϕ ( ) n,d (t -kT B ) .
(10) The signal ξ n (t) is expressed as a superposition of individual L → H and H → L transitions, denoted respectively with the subscripts u, d. In [START_REF] Cao | A New Training Approach for Robust Recurrent Neural-Network Modeling of Nonlinear Circuits[END_REF], the index ranges through hierarchical approximation levels, up to a maximum level L. The construction of these levels is discussed below.

Let us first assume that no FFE is applied, so that no successive approximations are required. In such case, we set L = 0 (index is removed). The signal ξ n (t) is thus constructed as a sequence of individual transitions L → H and H → L, centered respectively at time instants k ∈ Ω n,d (t), respectively. If translation-invariance holds, then these two basis functions are sufficient to reconstruct all switching events. For instance, a single-bit pulse p(t) can be represented with only two switching events as

p(t) = ϕ (0) n,u (t) + ϕ (0) n,d (t -T B ) (11)
Consider now the case where FFE is activated, so that the shape of a switching event at t = kT B changes depending on the presence of another switching event at the immediately preceding time t = (k -1)T B . We collect in sets k ∈ Ω n,d such switching instants. It is clear that since there are switching events that are preceded by at least two non-switching intervals. We define two level-1 basis functions ϕ

Ω (1) n,ν ⊂ Ω (0) n,ν , ν ∈ {u, d}, (12)
n,u (t) and ϕ

n,d (t) as the local corrections that need to be applied to the level-0 approximation, in order to account for the difference in switching waveforms due to the presence of an immediately preceding switching event. We use these basis functions, centered at kT B with k ∈ Ω n,d (t) are used to represent the correction that is required locally at those instants. As a result, the hierarchical approximation [START_REF] Cao | A New Training Approach for Robust Recurrent Neural-Network Modeling of Nonlinear Circuits[END_REF] is obtained, where less and less basis functions are needed when the level increases, since by construction

Ω ( ) n,ν ⊂ Ω ( -1) n,ν , ν ∈ {u, d}, = 1, . . . , L. (13)
In addition to sparsity, the signal representation (10) has another very convenient feature. As depicted in Fig. 6, the amplitude of such basis functions decreases with increasing levels . This makes the proposed hierarchical representation more and more accurate as the levels increase, in a predictable and controllable manner. It is expected that uniform convergence holds to the exact signals, provided of course that no other sources of errors is present (e.g., the driver behaves exactly as a linear structure, which is not true in general).

As already mentioned, the size of the index sets Ω ( ) n,u

and Ω ( ) n,d decreases when increasing , since the probability of occurrence of consecutive switching events decreases. These index sets are determined by a digital preprocessing of the logic bit sequence and they are exact, so that the computatioal complexity in the construction of the hierarchical approximation [START_REF] Cao | A New Training Approach for Robust Recurrent Neural-Network Modeling of Nonlinear Circuits[END_REF] is known in advance. Therefore, the signal construction is based on a nonlinear and adaptive approximation, which is however based on a sequence of iterative (through levels) linear superpositions of known basis functions.

B. Estimation of hierarchical basis functions

We now consider the construction of the basis functions ϕ ( ) n,u (t) and ϕ ( ) n,d (t). These are estimated for each level through a recursive peeling procedure, starting from a training response pair ξ( ) n,u and ξ( ) n,d computed using the reference TL driver model excited respectively by a pair of bit sequences S ( ) u , S ( ) d that include all combinations of required switching patterns. Estimation is performed via incremental subtraction, while truncating the duration of transient sequences (the basis functions) at a time T , determined when the expected steadystate is reached. In particular:

• at level = 0, only isolated switching events are used to define the basis functions as

S (0) u = {0,1,1,1,...} ϕ (0)
n,u (t) = ξ(0) n,u (t) for the ν = u transition and

S (0) d = {1,0,0,0,...} ϕ (0) n,d (t) = ξ(0) n,d (
t) for the ν = d transition.

• at level = 1, two consecutive switching events are used to define the basis functions as

S (1) u = {0,1,0,0,0,...} ϕ (1) n,u (t) = ξ(1) n,u (t) -ϕ (0) n,u (t) -ϕ (0)
n,d (t -T B ) for the ν = u transition and

S (1) d = {1,0,1,1,1,...} ϕ (1) n,d (t) = ξ(1) n,d (t) -ϕ (0)
n,d (t) -ϕ (0) n,u (t -T B ) for the ν = d transition.

• at level = 2, three consecutive switching events are used to define the basis functions as

S (2) u = {0,1,0,1,1,1,...} ϕ (2) n,u (t) = ξ(2) n,u (t) -ϕ (0) n,u (t) -ϕ (0) n,d (t -T B ) -ϕ (0) n,u (t -2T B ) -ϕ (1) n,u (t) -ϕ (1)
n,d (t -T B ) for the ν = u transition and

S (2) d = {1,0,1,0,0,0,...} ϕ (2) n,d (t) = ξ(2) n,d (t) -ϕ (0) n,d (t) -ϕ (0) n,u (t -T B ) -ϕ (0) n,d (t -2T B ) -ϕ (1)
n,d (t) -ϕ (1) n,u (t -T B ) for the ν = d transition.

• at higher levels > 2, the process is iterated by considering + 1 consecutive switching events and defining the basis functions by removing all contributions required by the approximation of the training signals based on the basis functions at levels up to -1. A graphical illustration of the estimated basis functions spanning four levels = 0, 1, 2, 3 for the example driver of Fig. 5 is provided in Fig. 6.

V. COMPLETE CHANNEL SIMULATION AND RESULTS

Once all basis functions ϕ ( ) n,ν (t) are estimated from suitable identification signals, and the index sets Ω ( ) n,ν in the hierarchical approximation [START_REF] Cao | A New Training Approach for Robust Recurrent Neural-Network Modeling of Nonlinear Circuits[END_REF] are available for a given bit pattern, the transient switching waveforms j(t) of the driver model ( 2) are completely determined. The driver model is then completed with the static conductance G and DC biasing terms I ν , as discussed in Sec. III.

We can now use this model to establish a fast transient simulation framework for a complete high-speed link, including a complete channel description together with the analog input impedance of the receiver (also assumed to behave as a LTI system) and its Continuous-Time Linear Equalization (CTLE) block, if present. As Fig. 7 shows, the entire system from driver up to the analog receiver part is a LTI system subject to a transient input j(t). Let us then consider the desired output (voltage) waveforms η m (t) = v m (t) -v m (0), where also in this case a DC shift by the initial operating point has been applied to setup a transient simulation with zero initial conditions. Due to linearity, the transfer function H mn (jω) between the n-th input component j n (t) and the m-th output component η m (t) can be determined by a standard frequency-domain analysis, by suitably cascading all frequency-dependent characterizations of the various blocks. This setup is the same as in standard AMI tools.

We now exploit the linear superposition in [START_REF] Cao | A New Training Approach for Robust Recurrent Neural-Network Modeling of Nonlinear Circuits[END_REF]. Although this model is obtained as an adaptive nonlinear approximation, the latter is predetermined and statically computed once the bit pattern to be simulated is known. Therefore, we can consider [START_REF] Cao | A New Training Approach for Robust Recurrent Neural-Network Modeling of Nonlinear Circuits[END_REF] as a superposition of known elementary transient waveforms, one for each of the 2(L + 1) basis functions. A standard Fast Fourier Transform (FFT) processing is then applied to evaluate the response at the output due to a single individual basis function, which is expressed as

ψ ( ) mn,ν (t) = F -1 H mn (jω) • F ϕ ( ) n,ν (t) (14)
where F is the discrete Fourier operator (see Fig. 8). Thanks to superposition, the received voltages are then expressed as

η mn (t) = L =0 k∈Ω ( ) n,u ψ ( ) mn,u (t-kT B )+ k∈Ω ( ) n,d ψ ( ) mn,d (t-kT B ) (15
) Note that this expression is identical to [START_REF] Cao | A New Training Approach for Robust Recurrent Neural-Network Modeling of Nonlinear Circuits[END_REF] but uses different (known and predetermined) basis functions. In particular, it inherits sparsity due to the hierarchical multilevel expansion. The index sets Ω ( ) n,ν are unchanged. The constructive approximation at the receiver is illustrated in Fig. 9. The top panel shows the received voltages at the end of a coupled lossy channel driven by the driver of Fig. 5 with FFE enabled. The different curves demonstrate the increasing approximation quality as the number of levels is augmented. When only level = 0 is included, only single and isolated switching events are correctly represented. Including also level = 1 provides an accurate representation of any pair of bits that are consecutively switching. Adding more levels leads to convergence for all possible switching sequences. We conclude this Section with a few remarks. First, the proposed constructive approximation is able to estimate also common-mode signals in coupled channels, an example will be provided in next section. Second, nothing prevents application to multiconductor channels with more than two coupled signals, as far as adequate coupled models are available for all system parts. Finally, since the proposed constructive approximation aims at releasing the translation-invariance of the algorithmic TX-FFE part typical of the IBIS-AMI framework, it is expected that with minimal (and possibly no) modifications, the same approach can approximate well the responses of drivers operating in overclocking conditions, where similar effects due to single bit pulse responses extending beyond bit time are observed. A detailed analysis of this scenario is left to a future report.

VI. RESULTS

All numerical results presented in this section were obtained using a Workstation equipped with an Intel Core i9-7900X CPU @ 3.30 GHz and 32 GB RAM. For the same example driver considered in the previous sections (henceforth testcase #1) transmitting at a data rate of 3.125 Gbps and a real channel, Figure 10 compares the received differential and commonmode voltages computed using the true TL driver model with the results predicted by the proposed approach and the standard IBIS-AMI model. The responses predicted by the proposed linear model in [START_REF] Fan | Signal integrity Design for High-Speed Digital Circuits: Progress and Directions[END_REF] are obtained by post-processing the received waveforms [START_REF] Cunha | Validation by Measurements of an IC Modeling Approach for SiP Applications[END_REF], which are conveniently converted into the corresponding transient responses. This comparison highlights the superior accuracy introduced by the proposed approach. The bottom panel of Fig. 10 lacks the prediction for the common-mode voltage from the IBIS-AMI simulation, since this waveform is not supported by the current IBIS-AMI model structure. The common-mode voltage prediction of proposed approach is instead fairly accurate, with an approximation error that is likely due to the residual nonlinearity of the driver, which is here approximated as a pure LTI system.

Figure 11 shows the eye diagram (including jitter and crosstalk) obtained for a PRBS-31 pattern of one million bits (this bit stream was selected just to avoid any periodic repetition in the bit pattern). Thanks to linearity, common approaches for inclusion of deterministic and random jitter, as well as crosstalk, can be used in a post-processing phase, as in standard IBIS-AMI flows. Our prototypal non-optimized MATLAB implementation provided this result in 24 seconds. This performance should be compared to a state-of-the-art IBIS-AMI solver, when applied to the same testcase. The IBIS-AMI was simulated on the same hardware used to run the proposed model. The cumulative runtime was 1.7 seconds for the statistical simulation (which is not comparable in scope with our full transient simulation). Conversely, the bit-bybit simulation million bits) with the IBIS-AMI solver required 49 seconds. These times include parsing and all preprocessing steps by the adopted commercial solver. Therefore, we see that the full transient simulation based on the proposed approach is even faster than a high-end commercial software, with the additional advantages of common mode estimation and accurate prediction of non-ideal TX-FFE effects.

In terms of accuracy, both the IBIS-AMI model and the proposed hierarchical solver are seen to underestimate Eye Width (EW) and Height (EH) with respect to the TL results. However, the relative errors of proposed method (10% and 17%, respectively on EW and EH) are significantly smaller than what is achieved by IBIS-AMI (30% and 42%). We remark that these errors are exceedingly large for allowing reliable use of such models in channel qualification, and should be regarded as a worst-case scenario due to the particularly strong non-ideal TX-FFE behavior of this particular TL device. The above results are thus a clear indication that any model should be used only within its limits of applicability.

As a second validation, Fig. 12 and Fig. 13 collect similar comparisons for a more recent commercial 40-nm low-power driver [21] transmitting at a data rate of 2.5 Gbps, denoted as testcase #2. Both figures report a comparison of the proposed approach with a reference TL simulation and with a prediction achieved by a nonlinear Mpilog macromodel [START_REF] Stievano | Mπlog, macromodeling via parametric identification of logic gates[END_REF]. This second test further confirms the excellent accuracy of the proposed approach in reproducing the rich dynamical behavior of the device. Concerning the simulation times, the TL simulation takes 1146 seconds to simulate 4500 bits, the Mpilog model takes 739 seconds, whereas the same simulation based on the proposed hierarchical linear model requires less than 1 second. For this driver, the eye pattern shown in Fig. 13 exhibits a wider opening and reduced dynamical effects with respect to the one obtained for the previous and older device (see Fig. 11). Similar to the previous testcase, the CPU-time required for the eye diagram construction with one million bits is about 27 seconds with proposed method, and about 50 seconds with the IBIS-AMI bit-by-bit solver. In terms of accuracy, the proposed solver allows estimating the EW and EH with a relative error of 6% and 8%, respectively, which is comparable to the accuracy achieved by the IBIS-AMI model (i.e., 5% and 9%) and it is exactly the same of the one provided by the Mpilog model (i.e., 6% and 8%). Thus, we see that all models are nearly equivalent in the estimation of eye diagram metrics, although Fig. 12 clearly illustrates that, for the displayed bit sequence, the IBIS-AMI waveforms appear less accurate than those obtained by Mpilog and by proposed hierarchical solver. Summarizing, the proposed hierarchical linear model structure performs in all cases with a similar accuracy as the fully nonlinear Mpilog model but requiring simulation times that are comparable to the IBIS-AMI framework, albeit with some small overhead due to the implementation of the constructive hierarchical superposition.

VII. CONCLUSIONS

This paper presented an improved modeling and simulation methodology for the fast simulation of high-speed transmission over differential channels. On one hand, the state-of-theart and de facto industrial standard IBIS-AMI is extended by proposing a compact linear macromodel of a differential driver which embeds the TX-FFE mechanism via an ad hoc hierarchical-based expansion of its source terms. The proposed model allows the accurate modeling of both the differential and common mode device responses and a tunable accuracy which depends on the number of basis functions accounted for. On the other hand, channel simulation is done in a mixed time and frequency domain, yielding reduced simulation times and allowing to compute eye patterns for millions bits in few seconds. This is accomplished net of any code optimization, using a simple prototypal MATLAB implementation.

Fig. 1 .

 1 Fig. 1. Typical structure of a differential driver with its relevant port variables.

Fig. 3 .Fig. 4 .

 34 Fig.3. Testbenches used for the collection of device responses required for model parameters estimation. Top panel: .TRAN analysis used for the computation of the switching terms j(t) = [j 1 (t), j 2 (t)] T ; bottom panel: .TRAN analysis used for the computation of the dynamical submodel g. See text for details.

Fig. 5 .

 5 Fig. 5. Switching pattern of the differential voltage v d (t) = v 1 (t) -v 2 (t) with and without FFE for the example driver.

  . The rising and falling fronts are represented through dedicated basis functions ϕ

Fig. 6 .

 6 Fig. 6. Basis functions ϕ ( ) n,u (t) and ϕ ( ) n,d (t), for n = 1 and = {0, 1, 2, 3}.

  and k ∈ Ω

  , respectively, to obtain a refined approximation.The process is then iterated through levels = 2, . . . , L, by identifying the sets k ∈ Ω ) n,u and k ∈ Ω

  that are preceded by at least consecutive switchings. Approprate basis functions ϕ ( ) n,u (t) and ϕ ( )

Fig. 7 .Fig. 8 .

 78 Fig. 7. Channel response computation.

Fig. 9 .

 9 Fig. 9. Received voltages computed through superposition of increasing hierarchical levels (top) and error (bottom).

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Test case #1. Received differential (top) and common-mode (bottom) voltages, computed by different models.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Test case #2. Received differential voltages (top), computed by a reference transistor-level model, the proposed hierarchical simulator, an Mpilog [6] model and by IBIS-AMI. Bottom panel reports the transient error waveforms of the two macromodel-based simulations with respect to the TL reference signals.

Assistant with the Electromagnetic Compatibility Group, Politecnico di Torino, where he works on algorithm development for optimization and order reduction of large-scale dynamical systems, and behavioral modeling of nonlinear circuit elements, with a specific application to the characterization of digital integrated circuits. His current research interests include macromodeling of electrical interconnects for electromagnetic compatibility and signal integrity problems.

Riccardo


Ashvaany Egambaram

Nasreen Badruddin

Vijanth S Asirvadam

Eric Fauvet

Christophe Stolz

Tahamina Begum

Unsupervised Eye Blink Artifact Identification in Electroencephalogram

Keywords: Electroencephalogram, EB Artifacts, Automated Threshold

I. INTRODUCTION

The Electroencephalogram (EEG) signal has been in use to interpret cognitive processes and physiological activity of the brain for medical purposes such as diagnosing epilepsy, sleep disorders, coma, encephalopathies and brain death. The EEG signal is also being extensively used for various research purposes, for example in neuroscience, cognitive science, cognitive psychology, neurolinguistics and psychophysiological research [START_REF] Urigüen | Eeg artifact removal Ůstate-ofthe-art and guidelines[END_REF]. In addition, EEG signal is being progressively investigated for real time applications such as continual epilepsy monitoring and brain computer interface (BCI). Hence, it is important to extract meaningful neurological information from EEG signal to facilitate its interpretation. However, EEG signals are often contaminated with undesired non-neurological information which causes deviations in the signal of interest. These undesired signals are termed artifacts. Superimposition of these artifacts with EEG signal could potentially mystify EEG's interpretation. This is particularly relevant in medical field where EEG signals are being used as a sole diagnostic source, thus failing to recognize artifacts may severely affect clinical decisions. Therefore, artifact identification in EEG signal processing is the first and most crucial step prior to artifact removal. The most common types of artifacts contaminating EEG signals are the cardiac artifact, the muscle artifact and the eye blink artifact. The muscle artifact is induced by muscle movement and contraction during EEG recording which takes place when the subject talks and swallows. The pattern of the artifact purely relies on the degree of muscle movement and contraction. Cardiac artifact on the other hand is an effect of electrical activity of the heart. The influence of cardiac artifact on the scalp is typically low, with low amplitude, repetitive and are of regular pattern. The eye blink (EB) artifact is the strongest and most prominent type of artifact recorded along with EEG, hence this research will focus on identifying it. EB artifacts appear as spikes with amplitudes of around 10 times greater than the actual brain signals, noticeable in the delta wave range and can last up to 200ms to 400ms [START_REF] Klem | The ten-twenty electrode system of the international federation[END_REF], [START_REF] Sanei | EEG signal processing[END_REF].

A lot of methods have been developed in the past for automatic identification or detection of EB artifacts in EEG signal. One of the easiest and preferred ways for EB artifact recognition is by simply using an amplitude threshold [START_REF] Nolan | Faster: fully automated statistical thresholding for eeg artifact rejection[END_REF]. This method determines if an EB artifact present in an EEG segment if the amplitude in the segment exceeds a predetermined amplitude threshold. However, the amplitude of eye blinks may vary depending on the blinking strength of an individual, hence this method may not identify EB artifacts from an individual who exhibits gentle blinks which are lower than the threshold value. The other common way of detecting EB artifacts is through feature based identification which determines the presence of EB artifact in an EEG segment after extracting certain features. Some of the common features used are the kurtosis, maximum absolute value, entropy-based features and second-order difference in [START_REF] Nolan | Faster: fully automated statistical thresholding for eeg artifact rejection[END_REF], [START_REF] Barbati | Optimization of an independent component analysis approach for artifact identification and removal in magnetoencephalographic signals[END_REF], [START_REF] Mognon | Adjust: An automatic eeg artifact detector based on the joint use of spatial and temporal features[END_REF], [START_REF] Klein | A reliable statistical method to detect eyeblink-artefacts from electroencephalogram data only[END_REF]. As a general rule, these artifact detection features require a certain threshold value to classify or make a binary decision whether or not an EEG segment is contaminated by EB artifact. As elaborated earlier for the amplitude threshold, applying a fix threshold value for the features discussed may lead to detection errors due to the individual variance in blinking pattern and blinking strength. As a result, the threshold values may need to be tailored for every individual which will be impractical in real time applications. This study proposes a new approach to automatically identify EB artifacts with varying threshold without any supervision on the EEG signal. The performance of the proposed approach is measured by validating if the proposed approach is accurate in identifying EB artifacts in comparison with the use of a fix threshold.

The paper is organized as follows: Section II illustrates the methodology, Section III elaborates the proposed algorithm, results and discussions are included in Section IV and V. Finally, Section VI provides the conclusion.

III. PROPOSED APPROACH A. Unsupervised EB Artifact Region Detection

In an EEG signal recording, the EB artifacts are primarily captured in the frontal electrodes, Fp1 and Fp2. This is because the frontal electrodes are positioned close to the eyes. Another logical point to note here is, both eyes of any individual blinks simultaneously, hence the Fp1 and Fp2 electrodes should theoretically unveil high correlation whenever eyes blink. To validate this theory, the correlation between Fp1 and Fp2 channels of an EEG signal is computed in segments of 500 samples ( 1.95 seconds) per segment. This window size is chosen so that at least one EB artifact can be captured in this window. Figure 1 shows an example of Fp1's and Fp2's recordings plotted out and their corresponding correlation coefficient values in each segment.

The test has revealed that segments of Fp1 and Fp2 without EB artifact produces correlation below than 0.7, whereas segments containing EB artifact results in higher correlation, usually more than 0.9. Thus, existence of EB artifact in a particular segment or window is supported by the high correlation coefficient value, which proves the theory. Following this, segments of the EEG signal contaminated with EB artifacts are identified. However, EB artifact components should be identified for subsequent analysis or artifact removal, which requires a threshold. Next section will discuss on the automated threshold level determination on every window that was identified containing an EB artifact.

B. Automatic Thresholding to Identify EB Artifacts

Displacement or deviation of amplitude from the mean is chosen as a threshold criteria to classify the EB artifact's onset point. The displacement of amplitude is chosen as EB artifacts are higher in amplitude in nature relative to that of the EEG or brain signal, hence the EB artifacts are expected to yield higher displacement from mean amplitude in comparison with uncontaminated EEG potentials. A threshold value is required in the case of recognizing the starting event of an EB artifact. First, the windows of Fp1 channel that were confirmed containing EB artifacts identified through steps illustrated in section III-A are subjected to attain displacement distribution from the mean amplitude using:

Displacement[t] = |[X[t] -µ]| (1)
where, X[t] is the EEG signal's amplitude at sample t = n until t = n+500, n is any sample point of the signal where the window starts, and µ is the mean of that window or segment.

An EB artifact event detection algorithm within a window or segment confirmed containing an EB artifact is then designed. The onset of EB artifact is assumed 100 samples (0.39 seconds) before the first sample of amplitude displacement more than the threshold, within the window. The reason setting the onset point in advanced of 100 samples before the threshold is to provide a buffer for any subsequent analysis. The end point of the EB artifact is found by aligning a 1 second (256 samples) window apart from the first sample of amplitude displacement more than the threshold. A 1 second window is chosen, because an eye blink can last up to 0.8 seconds (205 samples) in duration which can fit well into this frame (100+256=356 samples).

An experiment was conducted to define classification criteria to classify EEG potentials and EB artifacts. The threshold value for EB artifact for windows containing EB artifacts are checked and determined by this experiment. From the displacement distribution, a mean (µ) and a standard deviation (σ) are acquired. Figure 2 shows an example of displacement distribution by setting the threshold value to be any displacement values beyond 1σ from the mean, while figure 3 shows EB artifacts components plotted in red when threshold is set beyond 1σ.

Displacement Distribution for Threshold greater than 1

EB Locations (Threshold = + 1 )

EEG signal EB Artifact Region

EEG detected as EB artifact

Fig. 3. EB Artifacts for Threshold greater than 1σ

From Figure 3, it can be seen that some of the EEG potentials are identified as EB artifacts when the threshold is set beyond 1σ from the mean.

Figure 4 shows the displacement distribution of the same segment by setting the threshold value to be any displacement values beyond 2σ from the mean, while figure 5 shows EB artifacts components plotted in red when threshold is set beyond 2σ.

From Figure 5, it can be seen that all EB artifacts are identified as EB artifacts when the threshold is set beyond 2σ from the mean.

Figure 6 shows the same displacement distribution by setting the threshold value to be any displacement values beyond 3σ from the mean, while figure 7 shows EB artifacts components plotted when threshold is set beyond 3σ.

From Figure 7, it can be seen that the EB artifacts are not identified as EB artifacts when the threshold is set beyond 3σ from the mean.

Therefore the threshold for EB artifact dominating the (2)

Figure 8 summarizes the EB artifact detection algorithm in a flowchart.

IV. RESULTS

This section aims to measure the accuracy of the proposed approach in detecting EB artifacts compared to the use of a constant or common threshold. The constant threshold is evaluated by fixing the threshold values of more than 10uV, 20uV, 30uV, 40uV and 50uV, whichever amplitude displacement that exceeds these thresholds are considered to be EB artifacts.

The binary prediction is used to test the detection of EB artifacts, which produces below outcomes and the confusion matrix as in Table I:

Accuracy: ACC is the opposite of error rate, where it is the ratio of correct EB artifact and EEG detections by total number of detections. The best accuracy is 1.

ACC = T P + T N T P + T N + F N + F P (4)
The results obtained from the proposed method in comparison with constant thresholds in detecting EB artifacts correctly, applied on the Fp1 channel of 10 EEG signals are tabulated in Table II.

V. DISCUSSIONS

As elaborated in Section III-A, an automated EB artifact detection algorithm is designed. The proposed approach detects EB artifacts without any human supervision. First EEG segments containing EB artifacts are recognized with the concept of correlation between EEG electrodes, Fp1 and Fp2. Secondly, to address the issue of fixing a constant threshold, achieved through fixing thresholds between 10uV to 50uV of amplitude displacement. The individual accuracy value for each EEG signal of the proposed automated threshold reveals the accuracy has not fluctuated much, between 94% to 99% compared to the constant thresholds, between 89% to 99% for 10uV, between 91% to 99% for 20uV, between 75% to 100% for 30uV, between 27% to 99% for 40uV, and between 21% to 99% for 50uV. The constant threshold has achieved 100% accuracy for EEG set 6, which means the constant threshold of more than 30uV for EB artifact detection suits well for this data set. On the other hand, the 30uV threshold has only achieved 75.9% of accuracy for EEG set 4. This indicates the performance of fixing a constant threshold is purely dependent on the nature of the EEG signal and is not consistent across all EEG signals. While the automated threshold changes for every window, whereby the threshold will correspond to the varying nature of the EB artifacts. Hence, threshold is automatically determined for every window without setting any specific value.

VI. CONCLUSION

Conventional EB artifact detection algorithms depends on constant thresholds or constant features to make a binary decision to recognize if an EEG segment contains EB artifact or not. In this paper, an unsupervised EB artifact detection algorithm in EEG signal is proposed and evaluated. The algorithm relies on the concept of correlation between two EEG electrodes and the amplitude displacement range when there is an event of EB activity. Based on the results and discussions, it is apparent that the automated and unsupervised EB artifact detection algorithm proposed in this paper is accurate in identifying EB artifact events in an EEG signal. On the other hand, the algorithm is also consistent in detecting EB artifacts across different EEG signals compared to a conventional algorithm which is fixed with a constant threshold to detect EB artifacts. Thus, it can be concluded that the proposed algorithm is a reliable solution in detecting EB artifacts across all types of EEG signals which may have individual variance due to blinking pattern and strength.

  II. METHODOLOGY A. Evaluation on Real EEG Signal 10 sets of EEG signals collected at Universiti Teknologi PETRONAS, (UTP) are used to evaluate the proposed algorithm. EB artifacts can be clearly captured in the frontal channels, Fp1-Fp2 electrodes of the EEG recordings. Hence, the proposed algorithm is evaluated on the frontal channel, Fp1, of these EEG signals. All recorded signals are of different durations, with a sampling frequency of 256 Hz. The evaluation was performed using MATLAB R2016b in Windows 7 Professional (64 bit OS) with a 4GB RAM.

Fig. 2 .

 2 Fig. 2. Displacement Distribution for Threshold greater than 1σ

Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .

 4567 Fig. 4. Displacement Distribution for Threshold greater than 2σ

•

  Fig. 8. Flowchart of the Automated EB Artifact Identification Algorithm

ACKNOWLEDGMENT

The authors would like to thank the Ministry of Education, Malaysia for supporting this research through the Fundamental Research Grant Scheme, FRGS (FRGS/2/2014/TK03/UTP/02/1) and the Higher Institution Centre of Excellence (HICoE) Scheme.


Thierry Contaret
email: thierry.contaret@im2np.fr

Jean-Luc Seguin
email: jean-luc.seguin@im2np.fr

Philippe Ménini
email: philippe.menini@laas.fr

Khalifa Aguir
email: khalifa.aguir@im2np.fr

Dr Anna G T Mignani

J.-L Contaret

K Seguin

Aguir

Physical-Based Characterization of Noise Responses in Metal-Oxide Gas Sensors

Keywords: Adsorption-desorption noise, low-frequency noise, metal-oxide gas sensors, noise spectroscopy

The noise level in the gas microsensors is a tool for characterizing the electrical conduction under various gases and a means to improve selectivity. Metal-oxide gas microsensors with WO 3 sensitive thin film are characterized using a low-frequency noise technique. The spectral form of the noise responses measured using our specific systems is similar for tested gases (ozone and nitrogen dioxide). We observe a clear Lorentzian behavior according to adsorption-desorption (A-D) noise theory. To identify the detected gas, a physical-based characterization model of A-D noise source is proposed and compared with the empirical flicker noise model. We show that the excess noise is due to the A-D processes on the surface of the sensors sensitive film. The Lorentzian parameters depend on the nature of the gases and the noise-level dependence with gas concentration is clearly demonstrated. This confirms the interest on noise spectroscopy to improve the selectivity of gas sensors.

I. INTRODUCTION

N OISE spectroscopy is of great interest for gas detection since it is expected to allow enhanced selectivity of gas microsensors. Indeed, the noise level in the gas microsensors is an indicator of performance but also a tool for characterizing the electrical conduction in presence of gas. Some experimental results have been shown that, using noise spectroscopy, an initially non-selective sensor can be used to analyze the composition of a number of different gases [START_REF] Kish | Extracting information from noise spectra of chemical sensors: Single sensor electronic noses and tongues[END_REF]- [START_REF] Kish | Detecting harmful gases using fluctuation-enhanced sensing with Taguchi sensors[END_REF].

In 2005, we published the first attempt to develop a model of adsorption-desorption noise in metal oxide gas sensors [START_REF] Gomri | Modelling on oxygen chemisorption induced noise in metallic oxide gas sensors[END_REF], based on the simplest adsorption isotherm derived from Langmuir's theory. Next we improved this model using Wolkenstein's theory to describe the noise generated by free electron's density fluctuations in the sensing layer due to instantaneous fluctuations in the number of the adsorbed molecules [START_REF] Seguin | Adsorption-desorption noise in gas sensors: Modelling using Langmuir and Wolkenstein models for adsorption[END_REF] and adding the effect of mobility fluctuations [START_REF] Gomri | A mobility and free carriers density fluctuations based model of adsorption-desorption noise in gas sensor[END_REF]. The first noise measurements made with different gases in metal oxide gas microsensors with Sn0 2 sensitive film were used to validate the consistency of our models [START_REF] Contaret | Noise spectroscopy measurements in metallic oxide gas microsensors[END_REF]. We observed experimentally that the noise is superimposed on the intrinsic noise and could be described by a sum of Lorentzian component. In this paper, we show new results of low frequencies noise characterization of metal oxide gas microsensors with WO 3 sensitive thin film. It aims to show the interest of noise spectroscopy to characterize measured noise spectra in metal oxide gas sensors in order to improve their selectivity. We will also show that the developed model can explain the measured noise origin.

II. NOISE MODELING IN GAS MICROSENSORS

In metal-oxide gas sensors, noise sources have similar behavior to those present in the devices based on intrinsic semiconductors. On the other hand, the physical origin of noise sources is different and is related to the chemical processes to the surface of the sensitive layer of the sensor.

A. Noise Sources in Semiconductors 1) From a microscopic approach, the fluctuations of number of free carriers in a semiconductor generate three noise components: The shot noise due to fluctuations of the carriers that create the conduction current. At low frequencies, the spectral density S I conduction of conduction current I conduct ion is constant and is written:

S I conduction = 2 • q • I conduct ion . ( 1
)
2) Thermal noise due to carrier diffusion. It presents a white spectrum. The spectral density S I conduction of diffusion current flowing through the measuring electrodes of the semiconductor is proportional to the conductance G between these electrodes:

S I thermal = 4 • k • T • G (2)
with k the Boltzmann constant, T the thermodynamic temperature; 3) The generation-recombination (G-R) noise caused by the carriers' number fluctuations associated with the processes of generation and recombination. It's an excess noise in addition to fundamental noise sources (shot noise and thermal noise). The spectral density of current fluctuations associated with the capture and emission of carriers by a single trap in the bulk semiconductor is given by the relation:

S I = K L F 1 + 4π 2 f 2 τ 2 (3)
with K L F the current noise spectral density at low frequencies proportional to the G-R induced variation on the current, and τ the average time associated with a trap.

When the generation-recombination involves interface traps, the phenomenon is called Random Telegraph Signal (RTS). If the fluctuations come from a large number of generationrecombination process, the spectral density approximately varies as 1/f in the low frequencies band [START_REF] Wong | Low-frequency noise study in electron devices: Review and update[END_REF]. The 1/f noise or flicker noise is predominant over the thermal noise or shot noise. It involves various energy traps. It is located in the volume [START_REF] Hooge | Experimental studies on 1/f noise[END_REF] or at the surface of semiconductor or at the interface between semiconductor and electrical contact [START_REF] Mcwhorter | Semiconductor Surface Physics[END_REF]. It is known that the characterization of excess noise in microelectronics devices is a measure of performance and quality. It thus provides relevant information on the traps density, type and location [START_REF] Contaret | Low frequency noise characterization and modeling in ultrathin oxide MOSFETs[END_REF].

B. Noise Sources Related to the Chemical Environmental

The main noise sources in the metal-oxide gas sensors are related to the properties of the metal-oxide of the sensitive layer and the chemical environment. In metal oxides, the noise depends strongly on the oxygen stoichiometry and displacement of oxygen atoms. The electronic charge transport is associated with the presence of oxygen at the surface and in the crystal lattice sites of the metal-oxide. The adsorptiondesorption process of oxygen atoms and the presence of inhomogeneities, stresses and grain boundaries in the metal oxide cause fluctuations of the oxygen density and, thus fluctuations of the electrical conductance [START_REF] Raychaudhuri | Measurement of 1/f noise and its application in materials science[END_REF].

In a gaseous environment, conductance fluctuations of microsensors due to free carrier's number and mobility fluctuations are thus related to concentration and distribution fluctuations of chemical species. The measured noise spectra result from the superposition of three noise sources contributions.

1) The adsorption-desorption (A-D) noise due to the adsorption-desorption of gas molecules. It is similar to the G-R noise in intrinsic semiconductor.

2) The diffusion noise related to the diffusion of molecules adsorbed on the surface adsorption, 3) The shot noise due to the current through the potential barriers at grain boundaries in the sensitive layer.

C. Adsorption-Desorption Noise Model

The dominant noise source in microsensors is the adsorption-desorption noise [START_REF] Contaret | Noise spectroscopy measurements in metallic oxide gas microsensors[END_REF]. The modeling of the spectral behavior A-D noise source is based on Wolkenstein's theory [START_REF] Wolkenstein | Electronic Processes on Semiconductor Surfaces During Chemisorption[END_REF]. Electrical conductivity fluctuations due to free carrier number and mobility fluctuations are proportional to the density fluctuations of adsorbed molecules:

S δσ ( f ) = σ 2 0 S δ N ( f ) (4)
where S δσ ( f ) is the power spectral density (PSD) of electrical conductivity fluctuations. σ 0 depends of free carrier density and mobility at adsorption-desorption equilibrium. S δ N ( f ) is the PDS of adsorbed molecules density fluctuation which can be written [START_REF] Wolkenstein | Electronic Processes on Semiconductor Surfaces During Chemisorption[END_REF], [START_REF] Djuric | Adsorption-desorption noise in micromechanical resonant structures[END_REF]: where

S δ N ( f ) = H 0 1 1 + f f c 2 (5)
H 0 = 4δ N 2 τ f c = 1 2πτ (6)
δ N 2 is the mean square value of the adsorbed molecules density fluctuation δ N. The time constant τ is

τ = dd dN eq - da dN eq -1 (7)
where a and d are the number of adsorbed molecules and the number of molecules desorbed per unit area, respectively. Whether using Langmuir's theory or Wolkenstein's theory, calculating the time constant shows dependence with the gas parameters (mass of adsorbed molecule, adsorption energy, gas partial pressure …). Details of the calculation of τ and δ N 2 can be found in Reference [START_REF] Gomri | A mobility and free carriers density fluctuations based model of adsorption-desorption noise in gas sensor[END_REF]. So, H 0 and f c are function of gas species. For a sensor of resistance R crossed by a current I 0 , the PSD of the fluctuations of the terminal voltage is

S δV ( f ) = I 2 0 R 4 s L 2 σ 0 H 0 1 + f f c 2 (8)
where s and L are respectively, the cross section area and length of the sensing layer.

III. NOISE SPECTROSCOPY TECHNIQUE

A. Experimental Setup

The specific noise spectroscopy system which we developed to measure the sensor noise response in air, nitrogen dioxide and ozone is presented in Fig. 1.

The power spectrum of the amplified current across the sensor is measured in the frequency range of 0.1 Hz-100 kHz using low noise current preamplifier. The sensors are placed in a stainless steel measurement chamber where mixed gases can be admitted and evacuated via distribution valves. To obtain a linear flow and homogeneous mixture, the incoming gases concentration is controlled by precision electronic flow meters. The gases are then mixed by a diffuser. The gas chamber, the preamplifier, the sensors biasing and the power supply for sensor local heating are incorporated in a Faraday cage. All these devices are battery powered. The configuration and control of the spectrum analyzer with LabVIEW software provide a uniform spectrum in the range of measurement. The gas flow controllers are also supported by LabVIEW.

B. Noise Measurements Results

Experimental noise measurements have been performed on metal-oxide gas microsensors with WO 3 sensing layer. The studied gas-sensitive structures come from the development of a multi-sensors technology involving four cells on one chip.

Each sensor consists of a suspended membrane of Si0 2 /SiN x . The membrane includes two conductive platinum layers: one as a heater, and a second for the sensing electrodes separated from each other by SiO 2 . A four-microsensor array was developed from this design by duplicating the microhotplate and adding bonding pads.

The WO 3 sensing layer was deposited by IM2NP-CNRS laboratory in Marseille on the micro-hotplate fabricated by the LAAS-CNRS lab in Toulouse [START_REF] Chalabi | Development of integrated smart nose based on wo3 gas sensors[END_REF].

To precisely characterize the measured noise responses under gas, it is necessary to evaluate the minimum noise level of the measurements system. This minimum level, often called system noise, is equal to the noise of preamplifier. For a sensibility of 10 -8 A/V, the noise spectrum of SR570 low-noise current preamplifier is reported in Fig. 2 (curve 4).

The spectrum varies as 1/f from 0.1 Hz to 10 Hz and shows a bandwidth up to 200 Hz. Taking into account the sensibility, the current noise is equal to 0.1 pA/ √ Hz at 1 Hz and corresponds to a very low noise level (see section IV-A). In figure 2, we also present the voltage noise spectra under dry air with a constant flow of 100 ml/min (curve 3), nitrogen dioxide (NO 2 ) with a high concentration of 90ppm (curve 2) and ozone (0 3 ) with a concentration of 0.2 ppm (curve 1). All measured noise spectra are higher by one or more orders of magnitude compared to amplifier noise level and have a Lorentzian-like behavior.

All series of measurements indicated that the studied sensors were much more sensitive to ozone than to nitrogen dioxide. This explains the low concentration of ozone and those higher for nitrogen dioxide. Moreover, the listing of maximum noise level for each gas presented in Fig. 2 shows that noise decreases when ozone and nitrogen dioxide concentration increases. So, the noise responses clearly vary with the nature and concentration of detected gas.

IV. CHARACTERIZATION OF MEASURED NOISE RESPONSES

The physical-based characterizations of noise responses in gas microsensors consists of two parts, first the acquisition of low frequency spectra in the presence of gas using specific noise measurements system, and, second the identification of each noise source with a spectral decomposition method.

A. Frequency and Noise Analysis of Amplifier

To extract useful information from the measured noise across the sensor in the presence of gas, the noise voltage spectral density is reduced to a noise current spectral density by taking into account the characteristics of low noise amplifier using (1):

S I sensor = S V measured (A L N A ( f )) 2 (9)
where S I sensor is the intrinsic current spectral density of sensor, and S V measured is the noise voltage spectral density [START_REF] Gomri | A mobility and free carriers density fluctuations based model of adsorption-desorption noise in gas sensor[END_REF]. Indeed, the noise generated due to the adsorption of a gas on an adsorption site with a given adsorption energy has a Lorentzian spectrum. For both studied gases, we observe that two Lorentzians dominate the low frequency spectrum. The white noise over 1 KHz is mainly the contribution of thermal noise due to sensitive layer resistance and the thermal noise of the amplifier due to the feedback resistor.

Taking into account the contribution of observed noise sources, the sensor current noise spectral density can be decomposed using the electrical model given by ( 2): [START_REF] Mcwhorter | Semiconductor Surface Physics[END_REF] where K i and f c i are the parameters of each Lorentzian and are related to gas properties and adsorption-desorption process as seen in section II. In particular, K i depends on the number of adsorbed molecules per unit area at adsorptiondesorption equilibrium, and f c i depends on the time constant of adsorption-desorption process. Fig. 4 presents a significant example of spectral decomposition by using the modeling of all observed noise source. As Lorentzian parameters of Table I, equation [START_REF] Mcwhorter | Semiconductor Surface Physics[END_REF] enables to accurately describe the spectral behavior of noise response. In particular, the correct description of the both Lorentzians show that the adsorption-desorption noise is dominant and A study of adosrption sites on WO 3 surfaces would be useful to confirm our opinion, but no data was found in the literature.

S I sensor = 2 i=1 K i 1 1 + f f c i 2 + 4kT R sensor + S I ampli f ier
To show the consistency and validity of our method, a comparison between the adsorption-desorption noise model and the empirical flicker noise model is given in Fig. 5. The empirical model is described by using (3):

S I empirical = K f β + 4kT R sensor ( 11
)
where K is a constant proportional with Hooge's constant and β the coefficient of spectrum slope. The values of noise parameters for both gases measured are reported in Table II. Whatever the detected gas, the empirical model approximately describes the behavior of noise spectra. The empirical flicker noise model does not highlight the observed physical origin of the A-D noise due to some adsorption sites. Such a description would mean that the A-D noise comes from the A-D process on a large number of adsorption sites. So, by identifying the origin of different noise sources, only the modeling based on the spectral decomposition method can correctly describe the behavior of the noise spectral response of the sensor in the presence of gas.

B. Gas Dependence of Noise Responses

In order to analyze the influence of the detected gases on the noise responses of the studied sensors, we applied the spectral decomposition method to characterize the measured spectra under different concentrations of O 3 and NO 2 . Figure 6 shows the results of noise modeling for five significant concentrations of ozone and three concentrations of nitrogen oxide. For all concentrations of both tested gases, the modeling of measured spectra clearly identifies two distinct Lorentzians as mentioned in the previous paragraph.

The extracted Lorentzian parameters' are reported in Tables II and III for ozone and nitrogen oxide, respectively. The cut-off frequencies of Lorentzians not significantly depend on the gas concentration but have some variations depending on the gas species. This is consistent with our theoretical approach presented in part II. With regard to LF levels of Lorentzians, we observe dependence with the nature and concentration of gas in according to developed model. Whether under ozone or nitrogen dioxide, the variation of the noise current spectral density is consistant with the variation of the gas concentration. However, the variation range of the spectral density seems approximately the same for both tested gases.

To discriminate the detected gas, we must take an interest in the evolution of LF levels depending on the gas concentration. In Fig. 7, we report the normalized LF noise levels of the both Lorentzians of each modeled spectrum for all concentrations of 0 3 and NO 2 . The relative LF noise levels are calculated by using following expression:

K i relative = K i gas K i air ( 12
)
where K i gas is the LF noise level (for 0 3 and NO 2 ) of the Lorentzian with the cutoff frequency f i gas and K i air is the LF noise level for air of the Lorentzian with the cutoff frequency f i air .

Either with ozone or nitrogen dioxide, the variation of noise levels with the concentration is similar for the two Lorentzians of each spectrum. The dynamic range is greater for the Lorentzian with the lowest cutoff frequency. On the other side, the variation of the levels is very different for the two gases and can therefore be a criterion of selectivity of gas sensors.

In addition, we also observed that the maximum dynamic range of noise levels is roughly the same for ozone and nitrogen dioxide and seems to be rather a characteristic of the sensors used. The ozone concentrations differ from those of nitrogen dioxide due to different sensitivity of the sensor for the both gases.

To highlight the possibility of noise spectroscopy to improve the selectivity of gas sensors, we extracted the value of sensors resistance from the thermal noise level observed from 10 KHz. The values of sensors resistance are reported in Table IV for each ozone and nitrogen dioxide concentration and, are consistent with the resistance values directly measured in a previous study [START_REF] Chalabi | Development of integrated smart nose based on wo3 gas sensors[END_REF]. Ozone and nitrogen dioxide being oxidizing gases, we verify an increase in sensor resistance when the concentration increases.

From the value of sensor resistance under air, we plotted in Fig. 8 the relative resistance depending on the ozone and nitrogen dioxide concentration. If we compare the figures 7 and 8, we immediately see that variation range of relative sensor resistance is less than a decade while the noise level variation range is about three decades. The detection of gas by measuring the variations in the noise levels would therefore be more sensitive and selective than sensor resistor measurement.

V. CONCLUSION

In this paper, a physical-based method to characterize the noise responses of gas sensor has been presented. Experimental noise measurements have been performed on metal oxide gas microsensors with a WO 3 thin-film and a platinum microhotplate. Using specific measurements system, the sensor noise response under gas exposure is measured in low frequencies range.

To analyze the power spectrum of the amplified current fluctuations across the sensor, we developed a noise model based on Wolkenstein's theory. Using spectral decomposition, the modeling of noise current spectral density enables to explain the origin of excess noise and gives more information on the detected gas. Lorentzian behavior of the measured spectra is correctly taken into account by using the description of adsorption-desorption process of gases. We have shown that the dependence of noise parameters with the nature and concentration of gas is consistent with the theory of adsorption-desorption noise. The most relevant result of our study is the largest dynamic range of noise levels variation relative to sensor resistance variation.

So, these results show the growing interest in the development of noise spectroscopy technique to improve the selectivity of gas sensors.

Fig. 1 .

 1 Fig. 1.General diagram of measurements system used to characterize microsensors noise responses.

Fig. 2 .

 2 Fig. 2. Voltage noise spectra under air, nitrogen dioxide (NO 2 ), ozone (O 3 ), and system noise level.

Fig. 3 .

 3 Fig. 3. SR570 low-noise current amplifier for a sensibility of 10 -8 (a) Frequencies characteristics. (b) Current noise.

Fig. 4 .

 4 Fig. 4. Spectral decomposition of a noise response under ozone (0.2 ppm). Curve 1: measurements. Curve 2: extracted noise model of low-noise current preamplifier. Curve 3: white noise proportional to the resistance of the sensor. Curve 4: adsorption-desorption noise modeled by lorentzians. Curve 5: total noise.TABLE I EXTRACTED NOISE PARAMETERS

Fig. 5 .

 5 Fig. 5. (a) O 3 gas. NO 2 gas. Curve 1: comparison between measurement. Curve 2: modeling using the adsorption-desorption noise model. Curve 3: empirical flicker noise model.

Fig. 6 .

 6 Fig. 6. Modeling of noise current spectral density. (a) Several ozone concentrations. (b) Nitrogen dioxide concentration.

NO 2 K 1 (A 2 / 15 Fig. 7 .

 212157 Fig. 7. Relative LF noise level for the both Lorentzian of each concentration of NO 2 and O 3 .

TABLE II EXTRACTED

 II NOISE PARAMETERS FOR O 3

TABLE III EXTRACTED

 III NOISE PARAMETERS FOR NO 2

TABLE IV EXTRACTED

 IV RESISTANCES FOR NO 2 AND O 3

	NO 2 (ppm)	R NO 2 (M )	O 3 (ppm)	R O 3 (M )
	40	110	0.2	66
	70	184	0.4	83
	90	236	0.6	110
			0.8	166
			1.0	331

Fig. 8. Relative sensor resistance for NO 2 and O 3 concentrations.

ACKNOWLEDGMENT

The authors would like to thank A. Combes from IM2NP/CNRS UMR 6242/Aix-Marseille University, Marseille Cedex, France, for technical support. They would also like to thank LAAS-CNRS, Toulouse, France, for their collaboration in the development of the advanced microsensors used in this work.

Thierry Contaret was born in Saint Germain-en-Laye, France, in 1975. He received the Ph.D. degree from the University of Montpellier II, Montpellier, France, in 2003, with a thesis on noise modeling and simulation in submicron metal oxide semiconductor transistors.

He joined the MicroSensors Group, Institut Matériaux, Microélectronique, Nanoscience de Provence, Marseille, France, in 2006, after a Post-Doctoral Research with the Institut de Microélectronique Electromagnétisme et Photonique, Grenoble, France, on noise characterization of advanced metal oxide semiconductor field-effect transistor. His current research is on improving the selectivity of gas sensors by using noise spectroscopy. He is a Senior Researcher with Paul CEZANNE -Aix-Marseille III University, Marseille. He is a Senior Lecturer on electronics with the Institute of Technology of Marseille, Marseille. He specializes in thin films preparation and characterization for applications in microsystems. Since 1997, he has been interested in gas microsensors and he developed a selective ammonia sensor based on CuBr mixed ionic conductor. He is currently with the Institut Materiaux Microelectronique Nanosciences de Provence, Marseille, on WO 3 gas sensors and selectivity enhancement strategies, including noise spectroscopy. He is a Professor with Aix Marseille University, Marseille, France. He is currently the Head of Sensors Group, Institute of Materials Microelectronic Nanosciences of Provence, Aix-Marseille University, Marseille. His current research interests include metal -oxide (WO 3 , SrTiO 3 , CuO, CeO 2 ) and organic thin films for gas sensors, flexible gas sensors, microsystems, selectivity enhancement strategies including surface modification of sensors, signal treatment, adsorption -desorption noise spectroscopy, modelling of sensor responses analysis, and low noise amplifier design.

Jean-Luc

Philippe Ménini was born in