Iulia Dumitrescu
email: dumitrescuiulia@ymail.com

Smail Bachir
email: smail.bachir@univ-poitiers.fr

David Cordeau
email: david.cordeau@univ-poitiers.fr

Jean-Marie Paillot
email: jean.marie.paillot@univ-poitiers.fr

Mihai Iordache
email: iordache@elth.pub.ro

Modeling and Characterization of Oscillator Circuits by Van der Pol Model using Parameter Estimation

Keywords: Oscillator, Van der Pol Model, parameter estimation, Nonlinear system, modeling, continuous time domain

In this paper, we present a new method for the modeling and characterization of oscillator circuit with a Van-Der Pol (VDP) model using parameter identification. We also discussed and investigated the problem of estimation in nonlinear system based on time domain data. The approach is based on an appropriate state space representation of Van der Pol oscillator that allows an optimal parameter estimation. Using sampled output voltage signal, model parameters are obtained by an iterative identification algorithm based on Output Error method. Normalization issues are fixed by an appropriate transformation allowing a quickly global minimum search. Finally, the proposed estimation method is tested and validated using simulation data from a 1GHz oscillator circuit in GaAs technology.

Introduction

The interest in oscillators is motivated, among others, by their recent use in antenna array systems. Coupled oscillators are a simple and efficient method for phase control in microwave antenna arrays due to their synchronization properties [START_REF] Health | Beam steering of nonlinear oscillator arrays through manipulation of coupling phases[END_REF][START_REF] Liao | Beam Scanning With Coupled VCOs[END_REF][START_REF] Liao | A new phase-shifterless beam-scanning technique using arrays of coupled oscillators[END_REF][START_REF] Tohme | Analysis of the Frequency Locking Region of Coupled Oscillators Applied to 1-D Antenna Arrays[END_REF] . Nevertheless, for the design of an oscillator circuit or an array of coupled oscillators, one must be able to test and implement control laws on simple systems that reproduce well the dynamics of the real system. Classical Van der Pol (VDP) oscillators are usually used to model microwave oscillators. Indeed, in [START_REF] York | Nonlinear analysis of phase relationships in quasioptical oscillator arrays[END_REF][START_REF] York | Oscillator Array Dynamics with Broadband N-Port Coupling Networks[END_REF][START_REF] Lynch | York Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF] , R. York made use of simple Van der Pol oscillator to model an array of oscillators coupled through many types of circuits. Hence, this theory provides a full analytical formulation allowing to predict the performances of microwave oscillators arrays. Nevertheless, the identification of the parameters of the VDP oscillator is not always made easily.

The main difficulty in modeling an oscillating system is caused by its highly and complex non-linear behavior. The parameter computation for strongly non-linear devices is difficult and time consuming.

Parameter identification can be defined as mathematical modeling of a process in order to understand, predict and enhance its dynamic behavior using input and output set of data. Furthermore, such an identification can be carried out using an iterative identification algorithm based on Output Error (OE) method. Output Error techniques are based on iterative minimization of a quadratic criterion, also called cost function, by a Non-Linear Programming (NLP) technique [START_REF] Young | Parameter estimation for continuous-time models -a survey[END_REF][START_REF] Landau | Adaptive Control : The model reference aproach[END_REF][START_REF] Unbehauen | Identification of continuous systems, Systems and control series[END_REF] . This method requires much more computation and do not converge to unique optimum [START_REF] Unbehauen | A review of identification in continuous-time systems[END_REF][START_REF] Mensler | Experimental comparison of continuoustime model identification methods on thermal process[END_REF][START_REF] Tulleken | Grey-box modeling and identification using physical knowledge and bayesian techniques[END_REF] . Nevertheless, OE methods present very attractive features because the simulation of the output model is based only on the knowledge of the input so that the parameter estimates are unbiased [START_REF] Walter | Identification of parametric models from experimental data[END_REF][START_REF] Ljung | System identification: Theory for the user[END_REF][START_REF] Richalet | Identification des processus par la mèthode du modèle[END_REF] . Moreover, OE methods can be used to identify non linear systems. For these reasons, the OE methods are more appropriate in microwave systems characterization [START_REF] Bachir | New identification procedure for continuous-time radio frequency Power Amplifier model[END_REF] .

Due to these considerations, the aim of this paper is to present a tool allowing to extract the electrical and nonlinear parameters of the Van der Pol model, reproducing the behavior of a 1.025GHz GaAs FET oscillator circuit, using an identification technique based on Output Error (OE) method. This paper will be organized as follows. The principles of oscillator circuits and their representation with VDP model in state space will be described in section 2. In section 3, we will introduce and present the proposed method allowing to estimate the involved parameters using Output-Error technique. The application of this technique in the special case of an oscillator circuit will also be investigated. Two examples will then be studied in section 4: a theoretical one to analyze the performances of the method in term of robustness under noisy data and a practical one in order to see the efficiency of the study in a real case.

Oscillator design and modeling

The basic form of a harmonic oscillator consists in an electronic amplifier connected in a positive feedback loop to a frequency selective network like an RLC resonator. Such an architecture is presented in (Fig. 1) where the amplifier, also called the active part, compensates for the tank losses to enable a constant amplitude oscillation.

VDP model

The models used for real circuits need to take into consideration the energy dissipation and compensation. Thus, the circuit must contain a dynamic element with a negative resistance so that oscillation is maintained. Negative resistance circuit elements are purely theoretical and does not exist as discrete components. However, some elements present, in some parts of their operating characteristic, a negative resistance. The Van der Pol oscillator is considered as a widely used example in G β Fig. 1. Block diagram of a feedback oscillator the literature [START_REF] Van Der Pol | The nonlinear theory of electrical oscillations[END_REF] . It is highly nonlinear and it can exhibit both stable and unstable limit cycles. The main advantages of this model is the simplicity in simulation and high level in understanding and studying nonlinear dynamical phenomena, but also because it contains the necessary and sufficient elements to describe a real oscillator.

R u C (t) i N L (t) C L i L (t) G N L (t) = -a + b • u 2 C (t)

Fig. 2. The Van der Pol Model

The VDP model used contains a nonlinear conductance, noted G N L (t), and a RLC resonator circuit as shown in Fig. (2). The nonlinear conductance, G N L (t), models the active part of the oscillator and exhibits the necessary negative resistance region on its operating characteristic. The current i N L (t) to voltage u C (t) relationship describing the behavior of the nonlinear conductance is a cubic equation:

i N L (t) = -a • u C (t) + b • u 3 C (t) = G N L (t) • u C (t) (1)
Since, such an oscillator topology can be modeled by a quasi-linear representation allowing a simple analytical calculation [START_REF] Odynec | RF and Microwave Oscillator Design[END_REF] , the negative conductance presented by the active part in order to compensate for the tank losses R is equal to

-G N L (A) = -a + 3 4 • b • A 2
where A is the amplitude of the oscillations. Mathematically, a general form of the second order differential equation that describes the Van der Pol oscillator is:

üC (t) -µ (1 -u 2 C (t)) uC (t) + u C (t) = 0 (2)
where the parameter µ is a damping indicator, uC (t) =

VDP state space representation

The state equation representation of a physical system consists in describing the system as a set of inputs, outputs and state variables related by first-order differential equations. This modeling method, using matrix representation, can be applied for systems with multiple-inputs and multiple-outputs (MIMO) and the model includes the internal state variables and the output variables. Nevertheless, the most important advantage of this modeling form is that the representation consists in simple first order differential equations and provides directly a time-domain solution.

The state vector, noted x, composed of state variables, is a minimum set of variables that are fully describing the system and its response to any given set of inputs. These variables are related to the energy storage elements in a circuit. In this case, the state vector includes two state variables, the current through the inductor i L (t) and the voltage across the capacitor u C (t). Using elementary electrical rules (Kirchhoff's and Ohm's laws) and the nonlinear current expression given by (1), we can write the two following differential equations:

     C duC (t) dt + 1 R u C (t) + i L (t) = a • u C (t) -b • u 3 C (t) diL(t) dt = 1 L u C (t) (3)
For the Van der Pol oscillator, the state space representation of the system, written in matrix form, is the following set of first order differential equations:

ẋ(t) = A • x(t) + B(t) y(t) = C • x(t) (4)
where x(t) is the state vector such as:

x(t) =   u C (t) i L (t)   (5)
and the matrices are:

A =   a -1 R • 1 C -1 C 1 L 0   , B(t) =   -b C • u 3 C (t) 0   and C = 1 0
Let us note that the output signal y(t) is defined as the state variable of interest which is the output voltage u C (t) in our case.

where [.] T denotes the transposition operation.

The choice of this parameter vector is not random. Indeed, the parameters a and b of the non-linear conductance are not easily accessible and thus, need to be identified. Concerning the RLC resonator, designers have prior information on the resonator losses, i.e. the R parameter. Furthermore, the frequency of oscillation is determined by the LC tank. Hence, knowing one of these two parameters, the frequency of oscillation can be adjusted by varying the other one. In the practical case of a VCO(Voltage Controlled Oscillator), a variable capacitor is used to tune the oscillator. This has led us to the estimation of the parameter C. Let us note that this method is not restricted to this vector and can be generalized to the inductance L. ference between measured data and estimation. Thus, the output prediction error is defined as follows:

ε k = y * k -ŷk ( θ) (7)
where the predicted output ŷ is obtained by numerical simulations of the VDP model and θ is an estimation of the true parameter vector θ.

Parameter estimation with OE technique is based on minimization of a quadratic criterion defined as :

J = K k=1 ε 2 k = K k=1 y * k -ŷk ( θ) 2 (8)
Since ŷ(t) is nonlinear in parameters θ, the estimation is made iteratively using nonlinear programming techniques. Marquardt's algorithm [START_REF] Marquardt | An Algorithm for least-squares estimation of non-linear parameters[END_REF][START_REF] Trigeassou | Methods for Electrical Machines: Parameter estimation for knowledge and diagnosis of electrical machines[END_REF] was used to ensure an efficient and rapid convergence even with poor initialization. The parameters to be estimated are updated as follows:

θi+1 = θi -{[J ′′ θθ + λ • I ] -1 .J ′ θ } θ= θi (9)
J ′ θ and J ′′ θθ are respectively the gradient and the hessian such as:

J ′ θ = -2 K k=1 ε T k • σ y k,θ J ′′ θθ ≈ 2 K k=1 σ y k,θ • σ T y k,θ
with σ y k,θ = dy dθ the sensitivity function, λ the monitoring parameter and I the identity matrix.

Sensitivity computation

Sensitivity computation is an important point in the identification procedure. Indeed, the sensitivity functions express the effect of a parameter variation on the system output. We can define two types of sensitivity functions:

• σ y,θ = ∂y ∂θ : the output sensitivity function used in nonlinear programming algorithm;

• σ x,θ = ∂x ∂θ : the state sensitivity function.

The functions can be obtained directly using numerical differentiation methods, but they require high computation time and there is a risk of introducing errors caused by approximations. For this reason, it is recommended to solve, by simulation, the differential system that describes the dynamic of these variables. Thus, for each parameter θ i of the vector θ, the corresponding sensitivity functions are determined by partial differentiation of equation ( 4) so that:

     σx,θi = A σ x,θi + ∂A ∂θi x + ∂B(t) ∂θi σ y,θi = C σ x,θi (10)
For example, the sensitivity functions for the parameter C can be obtained by the simulation of the following differential model:

           σx,C =   a -1 R 1 C -1 C 1 L 0   σ x,C +   1 R -a • 1 C 2 1 C 2 0 0   x +   b C 2 u 3 C (t) 0   σ y,C = 1 0 • σ x,C (11)
It is important to note that all state space model like VDP model (Eq. 4) and sensitivity function models are solved using the 4 th order Runge-Kutta method [START_REF] Shampine | Numerical Solution of Ordinary Differential Equations[END_REF] .

Parameter normalization

If a hypothetical situation is studied, the values chosen for the different parameters are close. Nevertheless, in a real situation, the numerical values of the physical system parameters can be highly different, especially for RF systems. This can lead to convergence difficulties that are slowing or even make impossible the process to find the global optimum [START_REF] Bachir | New identification procedure for continuous-time radio frequency Power Amplifier model[END_REF][START_REF] Trigeassou | Methods for Electrical Machines: Parameter estimation for knowledge and diagnosis of electrical machines[END_REF] . The solution consist in normalizing the parameters values which means the normalization of the sensitivity functions.

The method is first explained for a single parameter θ n , then it will be extended to the entire vector of parameters θ to be estimated. Let us consider the parameter θ n with the initial value θ n0 so that θ n = θ n0 + ∆θ n . In this case, estimating the value for θ n is equivalent to estimate the ∆θ n parameter which represents the difference between the nominal and the initial value.

The ∆θ n expression depends on the parameter θ n0 in the following manner:

∆θ n = µ n θ n0 (12)
Thus, for the parameter θ n , we can write the following expression:

θ n = (1 + µ n ) • θ n0 (13)
Since the initial value θ n0 is known, the estimation of θ n can be made through a change of variables, i.e. by estimating µ n . The sensitivity function according to µ n is written as follows:

σ µn = ∂ ŷ ∂µ n = θ n0 • ∂ ŷ ∂θ n = θ n0 • σ θn (14)
where the sensitivity functions ∂ ŷ ∂µn are with the same magnitude order. In order to apply the normalization method for several parameters, all the sensitivity functions are normalized, leading to a well conditioned gradient and Hessian written as:

           J ′ µn = ∂J ∂µn = -2 • K k=1 ε k • σ k,µn J ′′ µnµm = ∂ 2 J ∂µn∂µm = 2 K k=1 σ k,µn • σ k,µm (15)
For the iterative estimation of all parameters of the vector µ, we apply the Marquardt's algorithm so that:

μi+1 = μi -{[J ′′ µµ + λ I] -1 • J ′ µ } μ=μ i (16)
Thus, (Eq. 16) allows to obtain the vector of parameters θi so that:

θi = θ 0 • (1 + μi .I) (17)
This method is very simple to implement because it does not change the structure of the algorithm, only the computation of the sensitivity functions changes.

Simulation results

The validation method consists in two main parts:

• a study of the robustness technique in presence of stochastic disturbances, • the characterization possibility of an oscillator using a CFY30 GaAs FET.

Validation of the technique and investigations

In order to validate the proposed technique, let us consider the identification of a 1.3GHz VDP oscillator of (Fig. ( 4)) simulated with Agilent's ADS software. C , the measured signal u C (t) is sampled with a period of T s = 0.01ns and the number of samples by realization is K = 10000. The proposed identification algorithm and data treatment is implemented on Matlab Mathworks software.

u C (t) + Measurement noise R = 200Ω L = 1.6nH C = 9.523pF G N L
In order to allow statistical analysis, a zero-mean white noise is added to the output data according to 30dB, 20dB and 10dB of SNR and 100 attempts were performed for each SNR value. Let us note that using noise free measurement, the estimated parameters are identical to the true values. It can be seen that the estimates are, on average, close to the exact value, whatever of the noise level. Regarding the accuracy of the estimation, let us note that once the value of SNR is decreasing, the uncertainty range increases for each parameter. Overall, the parameters identification, even in the presence of high noise level, gives satisfactory results.

In order to study the correspondence between the estimated output ŷ and the measured one y * and the parameters θ convergence to true values θ, we computed the Normalized Mean Square Errors (NMSE) based on the output data and on the

(ŷ-y * ) T •(ŷ-y * ) y * T •y * NMSE θ,dB = 10 log ( θ-θ) T •( θ-θ) θ T •θ (18)
The NMSE y,dB criterion is widely used by the radio frequency community, especially for describing the deviation level for time domain measurement. On the other hand, we have deliberately introduced the second criterion which reflects the deviation of the entire parameter space instead of the standard deviation which gives this information but only for a single parameter. Table (1) gives the identification results for a 10dB SNR. We can clearly see that the mean of the parameter θ imean is close to the true value with reduced standard deviation σ i . The normalized error NMSE y,dB on data and on parameters NMSE θ,dB , presented at the end of the table indicates a low error rate.

Simulation results for a 1.025GHz GaAs FET oscillator circuit

Oscillator Design

The schematic of the oscillator used in simulation is shown in Fig. (6). The active part of the oscillator provides the necessary negative resistance to compensate for the resonator's losses. The transistor used is a CFY30 GaAs FET. The total power consumption of the oscillator is 115 mW, under V D = 5 V. The passive part is represented by a RLC parallel resonator that determines the frequency of oscillation.

For the tank chosen in this case, the frequency of oscillation has a value of about 1.025GHz and the resulting quality factor of the resonator is 61.35.

A small signal study of the active part was first made. The objective of this part is to define the starting conditions of the VCO and the various constraints that can ensure it will work. The model used for the active part of the oscillator is shown in Fig. (7) where C gs and C gd are respectively the gate-source and the gate-drain capacitances. In steady-state, in order to establish and maintain a stable oscillation, the admittance presented by the active part, expressed as follows, for C gd = 0, should have a negative real part: Therefore, from this curve it can be deduced that, to ensure a negative real part of Y e , the value of the inductance L g must be chosen so that the resonant frequency of the tank L g C gs is less than the desired oscillation frequency. Thus, the ratio ω ωr should be greater than unity in order to produce and maintain a constant oscillation. Hence, to ensure proper start-up of the oscillator, the value fo L g was chosen to be equal to 7nH.

Y e = g m 1 -ω 2 L g C gs + jωC gs 1 -ω 2 L g C gs = i e V e (19)

Estimation results

The behavior of the oscillator presented in the previous section is simulated using Agilent's ADS software in Transient with a sampling period of 0.01ns. For a sake of clarity, the obtained simulated waveforms will be called the measured data in the following. In general case, the resulting output voltage u C (t) for VDP model simulation corresponds to a cosine-type function. To ensure a proper match between the measured data and the estimation during the identification procedure, we made a truncation in output data so that the measured voltage starts with the maximum value, corresponding to cosine type function.

In these conditions, the aim is to find the parameter vector θ so that the Van der Pol model describes the data collected from the real circuit as accurately as possible. Let us note that no prior knowledge is available for a and b parameters. After few iterations, the obtained values of parameters are: 9) shows the evolution of parameter estimates during identification procedure. We can see that the convergence to optimal values is obtained after 40 iterations. The comparison between the measured and the estimated responses is shown in figure (10). We can observe good agreement with a maximum error around 0.06V and a computed NMSE y,dB of -32dB.

   a = 6.2 • 10 -3 b = 1.3 • 10 -3 C = 9.5 pF

Conclusion

In this paper, the study and the implementation of an output error approach used for the identification of oscillating circuits is presented. The physical parameters of the corresponding Van der Pol model of an oscillator have been identified, using sampled time data. A brief description of the real system was made and the mathematical model of the considered system is also described. The mathematical equations used for the implementation of the method were developed and the improvement of the algorithm, for RF systems which have a highly different magnitude order of their parameters, was presented.

The simulation results, when the simulated circuit has the same structure as the model in the presence of noise, prove the algorithm convergence and show its robustness under high level of disturbances. A good agreement was found between the response of the circuit and the estimated model, showing the efficiency of such an identification technique. This approach constitutes a simple way to find the Van der Pol model for the large class of oscillating systems.

The investigated optimization technique behave very well and thus it is suitable for characterization of practical RF oscillator circuits with classical VDP model. The next objective is to study how to estimate the parameters of an array of coupled oscillators using the proposed method. Our first investigation in this field seem to be very promising.

3 .

 3 Parameter identification of the VDP model Parameter identification is based on the definition of a model. Once the previous mathematical model has been chosen, we can define the parameter vector to be estimated as follows: θ = [ a b C] T

Fig. 3 .6

 3 Fig. 3. Principle of Output Error Methods

8

  Oscillator Modeling and Characterization using Parameter Estimation

Fig. 4 .

 4 Fig. 4. Simulated circuit with added noise

Fig. 5 .Figure ( 5 )

 55 Fig. 5. Projections of estimated values for 100 realizations with differents SNR

10

  Oscillator Modeling and Characterization using Parameter Estimation parameters such as: NMSE y,dB = 10 log

Identification results θimean σ i a = 8 . 5 •

 85 10 -3 8.625 • 10 -3 2.77 • 10 -3 b = 71 • 10 -4 73.4 • 10 -4 5.35 • 10 -4 C = 9.523 • 10 -12 F 9.517 • 10 -12 1.874 • 10 -14 NMSE y,dB -10.06dB NMSE θ,dB -15.02dB

Fig. 6 .Fig. 7 . 1 √ 1 √= 1 √Fig. 8 .

 671118 Fig. 6. Oscillator schematic using a CFY30 GaAS FET

  Fig.(9) shows the evolution of parameter estimates during identification procedure. We can see that the convergence to optimal values is obtained after 40 iterations.

Fig. 9 .

 9 Fig. 9. Evolution of estimations during identification procedure

Fig. 10 .

 10 Fig. 10. Comparison between the measured data and the estimations.

Table 1 .

 1 Average results for 100 attempts for 10dB of SNR


M Ionita

D Cordeau

J M Paillot

S Bachir

M Iordache

A CAD TOOL FOR AN ARRAY OF DIFFERENTIAL OSCILLATORS COUPLED THROUGH A BROADBAND NETWORK

Keywords: broadband network, van der Pol differential oscillators, synchronization, coupled oscillators, design automation

50 words): A new expression of the equations describing the locked states of two oscillators coupled through a resistor is presented in this paper. This theory has led to the elaboration of a CAD tool which provides, in a short simulation time, the frequency locking region of two coupled differential oscillators.

I. Introduction

Arrays of coupled oscillators are the subject of increasing research activity due to their use in new applications such as power-combining techniques or beam steering of antenna arrays, by synthesizing a constant phase shift distribution between two adjacent oscillators [START_REF] Kuhn | Power combining by means of harmonic injection locking[END_REF][START_REF] Nogi | Mode analysis and stabilization of a spatial power combining array with strongly coupled oscillators[END_REF][START_REF] Health | Beam steering of nonlinear oscillator arrays through manipulation of coupling phases[END_REF][START_REF] Liao | Beam Scanning With Coupled VCOs[END_REF][START_REF] Liao | A new phase-shifterless beam-scanning technique using arrays of coupled oscillators[END_REF][START_REF] Toon | Evaluation of beam steering in circular planar array of coupled microwave oscillators[END_REF]. For a onedimensional oscillators array, the required inter-stage phase shift can be obtained by detuning the free-running frequencies of the outermost oscillators in the array [START_REF] Liao | A new phase-shifterless beam-scanning technique using arrays of coupled oscillators[END_REF]. Furthermore, it is shown in [START_REF] Liao | Beam Scanning With Coupled VCOs[END_REF][START_REF] York | Nonlinear analysis of phase relationships in quasioptical oscillator arrays[END_REF][START_REF] York | Injection and phase-locking techniques for beam control[END_REF] that the resulting phase shift is independent of the number of oscillators in the array.

In R. York's previous works, coupled microwave oscillators have been modeled as simple single-ended van der Pol oscillators coupled through either a resistive network or a broadband network [START_REF] Liao | Beam Scanning With Coupled VCOs[END_REF][START_REF] York | Nonlinear analysis of phase relationships in quasioptical oscillator arrays[END_REF][START_REF] York | Oscillator Array Dynamics with Broadband N-Port Coupling Networks[END_REF]. Unfortunately, these works are limited to the cases when the coupling network bandwidth is much wider than the oscillators' bandwidth. In these conditions, a generalization of Kurokawa's method [START_REF] Kurokawa | Injection locking of solid state microwave oscillators[END_REF] was used to extend the study to a narrow-band circuit allowing the equations for the amplitude and phase dynamics of two oscillators coupled through many types of circuits to be derived [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF].

Since these works, only a small number of articles have presented new techniques for the analysis of coupled-oscillator arrays [START_REF] Collado | Harmonic-Balance Analysis and Synthesis of Coupled-Oscillator Arrays[END_REF][START_REF] Suarez | Harmonic-Balance techniques for the design of coupledoscillator systems in both unforced and injection-locked operation[END_REF][START_REF] Georgiadis | New Techniques for the Analysis and Design of Coupled-Oscillator Systems[END_REF][START_REF] Suarez | Stability and Noise Analysis of Coupled-Oscillator Systems[END_REF]. In [START_REF] Georgiadis | New Techniques for the Analysis and Design of Coupled-Oscillator Systems[END_REF], a semi-analytical formulation is proposed for the design of coupled-oscillator systems, avoiding the computational expensiveness of a full harmonic balance synthesis exposed in [START_REF] Collado | Harmonic-Balance Analysis and Synthesis of Coupled-Oscillator Arrays[END_REF][START_REF] Suarez | Harmonic-Balance techniques for the design of coupledoscillator systems in both unforced and injection-locked operation[END_REF]. In [START_REF] Suarez | Stability and Noise Analysis of Coupled-Oscillator Systems[END_REF], a simplified closed-form of the semianalytical formulation proposed in [START_REF] Georgiadis | New Techniques for the Analysis and Design of Coupled-Oscillator Systems[END_REF] for the optimized design of coupled-oscillator systems is presented. Nevertheless, even if these new semi-analytical formulations permit a good prediction of the coupled-oscillator solution, it is only valid for the weak coupling case. On the other hand, the York's theory [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF] provides a full analytical formulation allowing to predict the performances of microwave oscillator arrays for both weak and strong coupling.

In the authors' previous work [START_REF] Tohmé | Analysis of the Frequency Locking Region of Coupled Oscillators Applied to 1-D Antenna Arrays[END_REF], mathematical manipulations were applied to the dynamic equations describing the locked states of the coupled oscillators proposed in [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF]. A reduced system of equations was obtained, leading to the elaboration of a CAD tool that provides, in a considerably short simulation time, the frequency locking region of the coupled oscillators, in terms of the amplitudes of their output signals and the phase shift between them. Nevertheless, although this CAD tool allows to obtain, in a short simulation time, the cartography of the locked states of the coupled oscillators, it is limited to the case of two single-ended van der Pol oscillators coupled through an RLC circuit. Furthermore, it does not allow an accurate prediction of the oscillators' amplitudes since the van der Pol model used by J. Lynch and R. York in [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF] is too simple. Moreover, since the theoretical limit of the phase shift obtained for an array of single-ended coupled oscillators is within the range of 90° [START_REF] York | Nonlinear analysis of phase relationships in quasioptical oscillator arrays[END_REF], it seems interesting to analyze the behavior of an array of differential oscillators. Indeed, in this case, the theoretical limit of the phase shift is within 360° due to the differential nature of the array, leading to a more efficient beam-scanning architecture for instance. Furthermore, differential oscillators are widely used in high-frequency circuit design due to their relatively good phase noise performances and ease of integration. Moreover, the use of a broadband coupling network, i.e. a resistor, instead of a resonant one, can lead to a substantial save in chip area.

Due to these considerations, the aim of this paper is to present a new CAD tool allowing to acquire, in a considerably short simulation time, the frequency locking region of two differential oscillators coupled through a resistive network. This paper will be organized as follows. A brief review of R. York's theory giving the dynamics for two van der Pol oscillators coupled through a resonant network will be presented in subsection II-A. In subsection II-B, a new expression of the equations describing the locked states of two van der Pol oscillators coupled through a broadband network using an accurate model allowing a good prediction of the oscillators' amplitude will be presented. Then, in subsection II-C, we will introduce a reduced system of equations with no trigonometric aspects allowing an easier numerical solving method. In section III, the developed CAD tool as well as its reliability test using Agilent's ADS software will be presented. The details concerning the variables estimation technique used to solve the system of equations describing the locked states of the two van der Pol coupled oscillators as well as the stability analysis of the locked states will be described in subsections III-A and III-B, respectively. Finally, in subsection III-C, the simulation results showing the reliability and the accuracy of the proposed tool will be presented.

II. Coupled van der Pol oscillators theory

A. Overview of R. York's theory

R. York's work considers the analysis of simple single-ended van der Pol oscillators coupled through either a resistive network or a broadband network that produce a constant amplitude and a phase shift between the oscillators [START_REF] York | Nonlinear analysis of phase relationships in quasioptical oscillator arrays[END_REF][START_REF] York | Oscillator Array Dynamics with Broadband N-Port Coupling Networks[END_REF]. This analysis is useful in showing the effects of the coupling network parameters but, unfortunately, it is limited to the broadband case. In these conditions, R. York used a generalization of Kurokawa's method [START_REF] Kurokawa | Injection locking of solid state microwave oscillators[END_REF] to extend the oscillators dynamics to the case of narrow-band coupling. Using a system made of two parallel resonant circuits containing nonlinear negative conductance devices and coupled through a series resonant circuit, as shown in Fig. 1, and starting from the admittance transfer functions (Y 1 , Y 2 , Y c ) binding the coupling current (I c ) to the oscillators' voltages (V 1 and V 2 ), J. Lynch and R. York described the oscillators' dynamic equations, as well as those for the amplitude and phase of the coupling current. Then, by setting the derivatives to zero, the algebraic equations describing the oscillators' frequency locked states were obtained as follows: ).

) ( sin 1 ) ( cos 1 ) ( sin 1 ) ( cos 1 2 1 0 2 0 02 1 0 2 2 2 2 0 1 2 0 2 0 01 2 0 1 2 1 2 0                                                                                            A A A A A A A A A A a c

B. New expression of the equations allowing an accurate prediction of the oscillators' amplitudes

The purpose of this subsection is to present a new formulation of York's equations describing the locked states of two coupled van der Pol oscillators using an accurate model allowing a good prediction of the oscillators' amplitudes. Indeed, in [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF], the negative conductance presented by the active part has the following expression:

            2 0 1 A G A G (2)
leading to a free-running amplitude of one van der Pol oscillator always equal to unity irrespective of the G 0 value. Thus, this result clearly shows the limitation of the van der Pol model used by Lynch and York in [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF]. In this case, in order to provide a more accurate model, let us consider the schematic of a van der Pol oscillator as shown in Fig. 2, for which the current in the nonlinear conductance

G NL is equal to       t bv t av t i 3   
with -a the negative conductance necessary to start the oscillations and b a parameter used to model the saturation phenomenon. Now, assuming a perfectly sinusoidal oscillation so that    

t A t v 0 cos .   , the expression of the current   t i
is the following:

      t A b t A bA a t i 0 3 0 2 3 cos 4 cos . . 4 3             (3)
Furthermore, let us recall that such an oscillator topology can be modeled by a quasi-linear representation allowing a very simple analytical calculation [START_REF] Odyniec | RF and microwave oscillator design[END_REF]. In this case, the expression in the first bracket in (3) represents the negative conductance presented by the active part. Now, in order to obtain the new system of four equations describing the locked states of two coupled van der Pol oscillators, let us consider the circuit of Fig. 3 made of two van der Pol oscillators coupled through a series RLC circuit. In this case, the frequency-domain equations, using the admittance transfer function binding the coupling current to the oscillators' voltages, can be written as follows:

                2 1 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 . . . , . . , V V Y I V Y A G V A Y I V Y A G V A Y I c c c L NL c L NL c                               (4)
where A 1 ,  1 , represents the amplitude and pulsation of oscillator 1, and so on and the subscripts "NL" and "L" denote respectively the nonlinear and linear portions of the circuit. The admittance transfer function for oscillator 1, assuming that the frequency of the oscillator remains close to its "free running" or uncoupled value 01  , has the following expression:

      1 01 1 1 1 1 2 1 ,                j C G A G G A Y L NL L (5)
where

C L 1 01 1  
is the tank resonant frequency of oscillator 1 and

  2 1 1 4 3 bA a A G NL    
is the negative conductance presented by the active part as expressed in (3) leading to sinusoidal oscillation.

Now, let us note

C G L a 2  
which represents the resonator "bandwidth" and using the substitution

    1 1 1 A S G A G L NL  
in order to simplify the notation, the expressions for the admittances of the oscillator 1 and 2, are respectively:

                  a L j A S G A Y     1 01 1 1 1 1 , (6)                   a L j A S G A Y     2 02 2 2 2 2 , ( 7
)
Thus, following the same steps as in [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF] but using ( 6) and ( 7) for the admittance functions of oscillators 1 and 2 and replacing 0 G by L G , the new system of equations describing the locked states of two van der Pol oscillators coupled through a resonant network is the following:

                                                    sin 1 cos sin 1 cos 2 1 0 2 0 02 1 0 2 2 0 2 1 2 0 2 0 01 2 0 1 2 0 1                                                  A A A A A S A A A A A S a c ac a a c ac a (8) where c L R G 1 0   , C G L a 2   and     L L i L i NL i G G bA a G A G A S 4 3 1 2      .
The case of a resistive coupling circuit can be deduced from ( 8) by letting the frequency bandwidth of the unloaded coupling circuit approach infinity, so that

  ac  leading to 1   and 0   . In these conditions, since       01 01    c and          02 02 c
, the equations describing the locked states can be written as follows:

                                sin cos sin cos 2 1 0 02 1 0 2 0 2 1 2 0 01 2 0 1 0 1 A A A A A S A A A A A S a a                       (9)

C. New expression of the equations describing the locked states of two van der Pol oscillators coupled through a resistor

Once more, a solution to [START_REF] York | Oscillator Array Dynamics with Broadband N-Port Coupling Networks[END_REF] indicates the existence of a frequency-locked state for the two van der Pol coupled oscillators. However, in this case, one can obtain the oscillators' amplitudes A ).

Nevertheless, due to the trigonometric and nonlinear aspect of equations ( 9), the solutions of this system are very difficult or impossible to obtain. In these conditions, mathematical manipulations are applied to the equations in order to obtain a new system of three equations with three unknowns A 1 , A 2 and ,  easier to solve. Indeed, from ( 9) we can write:

            2 1 0 02 1 0 2 0 2 1 2 0 01 2 0 1 0 1 sin cos sin cos A A A A A S A A A A A S a a                                            (10) Since 1 ) ( sin ) ( cos 2 2      
, we have:

  1 2 1 2 0 01 2 2 0 1 0 1                                       A A A A A S a       (11) and   1 2 2 1 0 02 2 1 0 2 0 2                                       A A A A A S a       (12)
Furthermore, according to [START_REF] Kurokawa | Injection locking of solid state microwave oscillators[END_REF], we obtain:

           0 1 01 ) cos( ) sin(         A S a (13) and             0 2 02 ) cos( ) sin(         A S a (14)
Now using the substitution

             0 1 01 1 1     A S , ω A X a and              0 2 02 2 2     A S
, ω A X a , then [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF] and ( 12) can be written as:

    2 2 2 0 1 2 1 2 0 1 2 1 1 A , ω A X A S A                  (15)
    2 1 2 0 2 2 2 2 0 2 2 2 1 A , ω A X A S A                  (16)
In the same way, ( 13) and ( 14) can be expressed as:

  , ω A X 1 1 ) tan(    (17)   , ω A X 2 2 ) tan(     (18)
And dividing [START_REF] Odyniec | RF and microwave oscillator design[END_REF] and ( 18) yields:

    0 1 ) ( ) ( 2 2 1 1 2 2 1 1      , ω A X , ω A X , ω A X , ω A X (19)
Hence, the new system of three equations with three unknowns A 1 , A 2 and ,  is presented as shown below:

        0 1 ) ( 1 ) ( 2 2 1 1 2 1 2 0 2 2 2 2 0 2 2 2 2 2 2 0 1 2 1 2 0 1 2 1                                 , ω A X , ω A X A , ω A X A S A A , ω A X A S A     (20)
where:

        a a A S , ω A X A S , ω A X                           0 2 02 2 2 0 1 01 1 1
Once this system of equations is solved, the solutions will then allow to determine the phase shift   between the two adjacent elements according to the pulsations 01  and 02  using the following expressions:

  ) ( arctan 1 1 , ω A X    (21) or   ) ( arctan 2 2 , ω A X     (22)

III. CAD Tool

In this section, the new CAD tool which allows to obtain, in a considerably short simulation time, the cartography of the locked states of two differential oscillators coupled through a resistive network is presented. Hence, the system of equations ( 20) was implemented on Matlab and solved using nonlinear programming technique.

A. Variables estimation technique

In this part, the problem of identifying the variables (or the roots) of the system of equations [START_REF] Marquardt | An Algorithm for least-squares estimation of non-linear parameters[END_REF] is discussed. The basic formulation of the general problem is the following: given a n th dimension system of nonlinear coupled functions

) ( k f
, we seek the value of the variables' vector  for which (

n 1, k , 0 ) ( f k    ).
In the case of the two van der Pol coupled oscillators model given previously, a new root finding procedure based on parameter estimation technique is developed. Let us consider the previous mathematical model given by equations [START_REF] Marquardt | An Algorithm for least-squares estimation of non-linear parameters[END_REF]. The vector of variables to be estimated is:

T a A A            
3 ) ( L a a L L G A G G bA a A f                                               (24) 2 2 1 2 0 2 02 2 2 0 2 2 2 2 2 1 4 3 ) ( L a a L L G A G G bA a A f                                               (25)                                                       0 2 1 02 0 2 2 01 3 1 4 3 1 4 3 ) (            L a a L a a G bA a G bA a f (26)
As a general rule, estimation with Output Error technique [START_REF] Ljung | System identification: Theory for the user[END_REF][START_REF] Bachir | New Identification Procedure for Continuous-Time Radio Frequency Power Amplifier Model[END_REF] is based on minimization of a quadratic multivariable criterion defined as :

   3 1 k 2 k ) ( f J  (27)
Unlike linear equations, most of the nonlinear equations cannot be solved in a finite number of steps. Iterative methods are being used to solve nonlinear equations. For criterion minimization, the variables values  can be driven iteratively to the optimum by Non Linear Programming techniques. Practically, we are using Marquardt's algorithm [START_REF] Marquardt | An Algorithm for least-squares estimation of non-linear parameters[END_REF] for off-line estimation:

    i J I J i i                 1 1 . . ˆ (28)
where

 ˆ is an estimation of the system solutions  ,         3 1 k , k k . f . 2 J J     is the gradient of criterion,         3 1 k T , k , k 2 2 . . 2 J J      
is the hessian,  is the monitoring parameter, and

      k k f , is the output sensitivity function or derivative of ) ( k f according to  .
It is important to note that the obtained  ˆ is an approximated solution if J ≈ 0, i.e.

1,3 , 0 ) (   k f k 
. This final value of the criterion corresponds to "small residual".

B. Stability analysis

For the stability analysis of the synchronized states, it is important to consider an initial differential system describing the amplitudes . Thus, dynamic equations developed in [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF] have been adapted to the studied case and give the following differential relations:

  02 01 . , , ,    a x x f  (29)
where

T cy cx A A A A x        2 2 1 1   and                                                                                                                                                                                                                                                           1 1 2 2 0 0 1 1 2 2 0 0 2 2 2
                                   A A A A A A A A A A A A A G G bA a A A A A A A G G bA a A x f

These differential equations are nonlinear in states

x , and unfortunately, the stability theory developed for the linear problem does not apply directly to this system. In practice, we typically linearize this system around a synchronized solution noted 0 x and consider the eigenvalues of the Jacobian matrix Hence, the obtained solution is asymptotically stable if all eigenvalues i  of the Jacobian matrix

      A so that:   x A x x f x x x    . . 0               (30)
x x x f A                is                                                                                                                                                     ac c
A A A A A A A A A A A A A A bA G a G A A A A A A A A A A bA G a G A                                                                   0 2 2 0 2 0 1 1 0 1 0 0 2 2 0 2 0 1 1 0 1 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1
      A satisfy        i eal  < 0 whereas the solution point is unstable if at least one eigenvalue i  satisfies 0         i eal  .
Let us note that in the developed CAD tool, the are obtained using the MATLAB software function eig().

C. Simulation results

In order to show the accuracy and the reliability of the proposed CAD tool, simulations using Agilent's ADS software will be performed. In these conditions, let us consider the design of two identical differential oscillators coupled through a resistor R c /2 of 200 Ω, as shown in Fig. 4. The oscillators' structure is based on a cross-coupled NMOS differential topology using a 0.35 μm BiCMOS SiGe process. The cross connected NMOS differential pair provides the negative resistance to compensate for the tank losses, and the tail current source is a simple NMOS current mirror which draws 28 mA. The frequency of oscillation is chosen to be close to 6 GHz and is determined by the LC tank at the drains. In these conditions, the inductance value, L 1,2 , is close to 0.8 nH and the capacitor value, C is close to 0.88 pF. The resistor value, R, is equal to 100 Ω so that the quality factor of the tank is equal to 3.3. A tail capacitor C T is used to attenuate both the high-frequency noise component of the tail current and the voltage variations on the tail node [START_REF] Hajimiri | Design Issues in CMOS Differential LC Oscillators[END_REF]. To ensure proper start-up of the oscillator, the transconductance of the NMOS transistor should be greater than R 1 . In these conditions, the sizes of NMOS transistors T 1 to T 4 are identical and chosen to be

m m L W g   35 . 0 70            .
Since the presented theory implemented in our CAD tool uses van der Pol oscillators to model microwave coupled oscillators, we performed the modeling of this structure as two differential van der Pol coupled oscillators as presented in [START_REF] Ionita | Analysis and Design of an Array of Two Differential Oscillators Coupled Through a Resistive Network[END_REF], using ADS simulation results for one differential NMOS oscillator at the required synchronization frequency. As a consequence, the two coupled oscillators of Fig. 4 can be reduced into two differential van der Pol coupled oscillators as shown in Fig. 5. In this case, let us note that the value of the coupling resistor on each path is equal to Rc/2 to match well with the theory based on the use of two single-ended van der Pol oscillators. For the modelling of the active part, the I = f(Vd 1 -Vd 2 ) characteristic of one differential NMOS oscillator of Fig. 4 at the required synchronization frequency was plotted leading to the typical cubic nonlinearity of a van der Pol oscillator. Hence, from this characteristic, the values of parameters a and b of the negative conductance presented by the active part of each oscillator were found to be respectively equal to 7.55 ּ10 -3 and 4 ּ10 -4 .

Then, knowing the parameters  0 ,  a , a and b, the proposed CAD tool provides the cartography of the locked states of the two differential coupled oscillators. For instance, for a synchronization frequency of 5.97 GHz, with a  = 5.68 ּ10 9 rad/s and a coupling constant 0  = 0.5, the cartography of the oscillators' locked states provided by the CAD tool is presented in Fig. 6.

Thus, this figure presents the variations of the phase shift,   , the oscillators' amplitudes A 1 and A 2 , and the synchronization frequency f s in function of f 01 and f 02 .

In order to validate the results provided by our CAD tool, we compared them to the simulation results of the two coupled differential oscillators of Fig. 4, obtained with Agilent's ADS software.

Let us note that with ADS, only a transient analysis of one point at a time of synchronization region allows to verify the synchronization results obtained with the CAD tool. The amplitudes of the oscillators' output voltages as well as the phase shift between them and the synchronization frequency are provided for each combination of f 01 and f 02 . For instance, let us now consider Fig. 6 where the point marked with an arrow in the four subplots represents a freerunning frequency f 01 = 5.73 GHz for oscillator 1 and f 02 = 6.21 GHz for oscillator 2. Let us note that these two free-running frequencies correspond to the equally spaced tunings line as presented in [START_REF] Lynch | Synchronization of Oscillators Coupled Through Narrow-Band Networks[END_REF]. The marked points lead to a phase shift of 32.39°, a synchronization frequency of the coupled oscillators of 5.97 GHz and an amplitude of 2.68 V at the output of each of the coupled oscillators.

In the same conditions, the two differential NMOS coupled oscillators of Fig. 4 simulated with ADS have lead to two sinusoidal waves at a synchronization frequency of 5.98 GHz, a phase shift of 31.23° and an amplitude of 2.72 V at the output of each oscillator, as presented in Fig. 7.

Thus, the errors regarding the phase shift and the amplitude predictions are only 3.71 % and 1.49 %, respectively.

Furthermore, Table I shows a comparison between the synchronization frequency, the phase shift and the amplitudes obtained in simulations with ADS for the two coupled NMOS differential oscillators, and with our CAD tool. Let us note that the free-running frequencies chosen for this table correspond to the equally spaced tunings line for which the synchronization frequency is equal to 5.97 GHz. Table II shows the same comparison as Table I, except that the free-running frequencies are not referred to the equally spaced tunings line. Hence, a good agreement was found between the results presented in both tables showing the reliability and the accuracy of the presented CAD tool. Nevertheless, as we approach the locking-region boundary, one can observe that the difference between the theoretical and simulated results is increasing. This is mainly due to the fact that the modeling of each NMOS differential oscillator as a differential van der Pol oscillator has been performed only at one frequency, i.e. the desired synchronization frequency.

Anyway, these results show that it is possible to adjust, with a relatively high accuracy, the freerunning frequencies of the two differential NMOS oscillators required to achieve the desired phase shift, amplitudes and synchronization frequency. The main advantage of this CAD tool is that, in an extremely short simulation time, one can obtain all the phase shifts, synchronization frequencies and amplitudes of the differential coupled system, in function of the free-running frequencies f 01 and f 02 of the two differential oscillators. Hence, since the inter-stage phase shift is independent of the number of oscillators in the array [START_REF] Liao | Beam Scanning With Coupled VCOs[END_REF][START_REF] York | Nonlinear analysis of phase relationships in quasioptical oscillator arrays[END_REF][START_REF] York | Oscillator Array Dynamics with Broadband N-Port Coupling Networks[END_REF], the proposed tool can also help the designer to rapidly find the free-running frequencies of the two outermost oscillators of the array required to achieve the desired phase shift. Indeed, without such a tool, a transient analysis of the oscillator array for different couples of oscillators' free-running frequencies must be performed in order to find the phase shift required for the targeted application. The proposed CAD tool can of course be used for other differential oscillators architectures and in strong and weak coupling cases.

IV. Conclusion

This paper presents a new CAD tool which provides, in a considerably shorter analysis time a cartography giving the phase shifts, synchronization frequencies and amplitudes of the differential output voltages of two differential oscillators coupled through a broadband network.

Starting from a new expression of the equations describing the locked states of two van der Pol oscillators coupled through a resistor, a simpler system of three equations with three unknowns was obtained. Solving this system on Matlab using nonlinear programming technique has led to the elaboration of the cartography of the locked states of the two differential oscillators coupled through a resistive network.

The reliability and the accuracy of the developed CAD tool was verified using Agilent's ADS software, where two differential NMOS oscillators coupled through a resistor were simulated.

The obtained results were in accordance with those generated by the presented CAD tool,

showing the usefulness of such a tool for the design of an array of differential oscillators coupled through a broadband network.

+ - R c /2 L 1 R C R L 1 C L 2 R C L 2 C R R c /2
G L G L -G NL (A 2 ) V 1 V 2 I c Y 1 Y c (V 2 =0) Y 2 R c /2

: 02 

 02 constant, where G 0 is the nonlinear conductance at zero voltage; the free-running frequencies or tunings of oscillators 1 and 2, respectively; us note that in the equations presented above, we refer the oscillators' freerunning frequencies or tunings, 01  and , and the synchronization frequency of the system,  , to the coupling circuit resonant frequency, solution to (1) indicates the existence of a frequency-locked state and allows to obtain the amplitudes A 1 and A 2 of the two oscillators as well as the inter-stage phase shift  

1 and A 2

 12 for any values of parameters a and b of the van der Pol non linearity under a sinusoidal assumption as well as the inter-stage phase shift   and the synchronization pulsation  for a

  from[START_REF] Marquardt | An Algorithm for least-squares estimation of non-linear parameters[END_REF] as:



  is the phase of the coupling current.

  asymptotically stable if and only if all the eigenvalues of the Jacobian matrix parts. Thus, after each estimation of variables values  , the resulting Jacobian matrix is evaluated and these six eigenvalues, noted i  , are computed.

Smaïl Bachir ,

 Bachir PhD received both his degree in Signal Theory from Polytechnic school of Algeria in 1997. He joined the scientific department of Leroy Somer Society and University of Poitiers in france where he received his Ph.D. degree in Automatic and Electrical Engineering in 2002. He is presently an Associate Professor at the university of Poitiers and a researcher in the department of XLIM/SIC of University of Poitiers. His research interest include signal processing, parameter identification, telecommunication and radiocommunication circuits. Prof. Eng. Mihai IORDACHE, PhD received the M.S. and Ph.D. degrees in electrical engineering from the Politehnica University of Bucharest, Romania, in 1967 and 1977, respectively. He is a Full Professor in the Electrical Department, Politehnica University of Bucharest, where he is working in the areas of circuit analysis and simulation, and in the Electrical Engineering Fundamentals. He is Doctoral Advisor at the Politehnica University of Bucharest, and the author or coauthor of more than 190 journal papers and 30 books. He is also a reviewer of different Scientific Conferences in the Analysis and Simulation Circuits. He was the recipient of the 2000 Romanian Academy Award.

Fig. 1 -Fig. 2 -Fig. 3 -Fig. 4 -Fig. 5 -Fig. 6 -Fig. 7 -

 1234567 Fig. 1 -Two single-ended van der Pol oscillators coupled through a resonant circuit

Fig. 2 G

 2 Fig.2

Fig. 4

 4 Fig.4

Fig. 5

 5 Fig.5

  Fig. 6

Fig. 7

 7 Fig. 7


Keywords: Design, Information Visualization, Incremental Visualization, Data Stream, Social TV

We present Visual Sedimentation (VS), a new design metaphor for visualizing streaming data inspired by the geological process of sedimentation. Our work started by early experiments visualizing political Twitter streams during the French 2012 presidential elections, and social interactions during a TV show. In both cases, the positive feedback we received expressed an unexpectedly high level of engagement from users, guiding our generalization of the metaphor.

This article explores VS and describes a new generative design space for Information Visualization. Geological sedimentation is our inspiration as it smoothly aggregates falling objects by compacting them into strata. We use this idea to visualize changing information in a new way by providing continuity between the representation of new and older data. The metaphor preserves an overall visual encoding while making it suitable for monitoring streaming data generated at unpredictable rates.

Towards Visual Sedimentation

INTRODUCTION

Visualizing streaming data, such as Tweets, is challenging. Among many reasons, the constant and unpredictable data updates make it difficult to allocate an optimal portion of the screen. Reflecting changes at different temporal scales, while keeping the visual mapping continuous for the viewer, is required to support such tasks as monitoring. Animation is probably the most common technique to convey updates, but it requires user attention, thus may not support real-time changes or emphasize enough recent and past updates.

CASE STUDIES

In this section, we present two prototypes of real-time micro-post visualizations based on early experiments. After assessing these two prototypes, we discuss a new design space that manages updating visualizations progressively.

Tweet Monitoring for Presidential Elections

Tweets can provide an informative pulse of what is being said on various public subjects. Even though tweets are short messages, their quantity, rhythms, and conversational structures are a way to understand a general opinion. One challenging aspect for monitoring near real-time updates is that we need to visually integrate the incoming data in the visual representation of data that has already arrived. This challenge arises for example during a meeting or a TV show, but also during a major event like presidential elections.

To visualize Twitter streams during one such events, the 2012 French presidential election, we developed a web-based visualization called Bubble-t [START_REF] Vuillemot | Bubble T : Et si twitter votait[END_REF]. It addresses the challenges of real-time visual update by using a particle system to fill in a bar chart (Figure 1 (a)). Each particle is a new piece of information that we refer to as a "token" (in this context, a Tweet). Each bar corresponds to a presidential candidate. Therefore, once a tweet is sent about a candidate, it is transformed into a token that is thrown into the candidate's bar. By accumulation, the columns are progressively filled in by the tokens. When a column is full, the n tweets that first arrived are flushed out of their column. We mapped the Twitter user's avatar to the token for additional informations and to further engage viewers. A static bar chart at the bottom captures all tweets during two different time periods (7days or 24hours).

The application received an award at the Google Dataviz 2012 challenge. Up to 81 273 unique visitors browsed the website in 6 months time. The average visit duration was 4min 21s, which is long. We received positive feedback from Twitter users with more than 1 800 tweets embedding the URL. This work also received attention from some national and authoritative media such as newspapers and radio. We also informally observed unexpected reaction from users. When the visualization was publicly presented, the audience in the room sent tweets, not only to test the system but also to show or see their own avatar inside the chart.

Despite these successes we also identified several issues:(p1) The transition between the bar holding the tokens and the static chart at the bottom is not continuous (Figure 1 (a) left). (p2) The time windows in which the token remain inside the bar is short due to the limited amount of screen space. (p3) The exploration of detail for aggregated tweets (in bars) is currently not possible.

Live Visualization for Social TV

With Bubble-t's popularity, a national TV broadcaster asked us to adapt it to monitor another class of live data: live Twitter poll results during a TV Show (Figure 1 (a) right). This mixed use of TV and social media is refers to a trend called "Social TV" where social media becomes the backchannel of an audience's sentiment. We developed a new tool called Bubble-TV to address the shortcomings of Bubble-t and to take the requirement of social TV into account.

At the TV show Bubble-TV was used as followed :

(1) At the beginning of the TV program, the presenters asked a question and the audience could reply by using Twitter with a specific hashtag. (2) During the show, every tweet sent by the audience was visually encoded into tokens which accumulated on top of the bar chart columns (Figure 1). The bar chart also encodes the aggregated values of the previous tweets. (3) Integration : The visualization was tightly integrated in the show: it was a part of the style of TV show's graphical elements. The TV show presenters frequently referred to the visualization to present trends, rather than reporting accurate numeric information. (4) The TV show's community manager also used it to monitor the Twitter activity and to present some tweets to the audience during the show.

This second iteration Bubble-TV differs from Bubble-t since a bar chart serves as a podium to structure the layout of tokens, trying to tackle (p1). It provides an overview of all tweets and details of last tweets. The community manager reported, after the show, that it was a useful to relate real-time Twitter activity to the content of the show. Bubble-TV has also significantly increased the participation on Twitter for the TV show: the number of tweets sent (5 300) was more than twice as during the previous show (2 000).

The positive part is that, despite our iteration, we did not break the attractiveness. Some problems remained: (p4) The relation between the bar chart and the tokens was better in this version but they were still disconnected. We observed the need to design a better transition between tokens and their aggregated state. (p5) The lack of visual cues for showing the density of tweet aggregated.

DESIGN GENERALIZATION

Trying to address the problems (p4) and (p5) from the previous section, inspired us to design a new framework for updating charts.

From Geological to Visual Sedimentation

Our assumption is that VS is highly familiar due to its physical equivalent (Figure 1 (b) left) that is referred as "sedimentation" and is defined in Wikipedia as: "the tendency for particles in suspension to settle out of the fluid in which they are entrained, and come to rest against a barrier.' This physical sedimentation surrounds us: mountains, hills or rivers are visible and familiar accumulations. The closest approach to ours is the use of the geological metaphor explored by Viegas et al. [START_REF] Viégas | Artifacts of the presence era: visualizing presence for posterity[END_REF], to compress and display data as strata over time. As far as we know, the full exploration of the sedimentation metaphor from falling objects, due to gravitational forces, to their aggregation into compact layers over time (Figure 1 (b) right) is novel. Our case studies also show it has a low learning curve.

Visual Sedimentation Design space

We call the visual extension of sedimentation Visual Sedimentation (VS). In its scope we defined three main states: suspension, deposition and aggregation; and two main transitions, decay and flocculation (Figure 2). To create a visualization using VS, continuous data streams are separated into discrete token elements. Barriers and walls have to be set up to provide containers for the tokens and define the general layout (bar chart, pie chart, . . . ). Tokens encode incoming data. They fall during suspension and land on other tokens and remain there during the deposition phase. Over time, the decay process changes the size of individual tokens. Finally, the flocculation process fluidly transforms token deposits into the aggregated chart area where the tokens adopt the aggregation state.

Generative Design

Given the previous guidelines, removing walls, changing the decay function or the attraction forces allows to recreate classic data charts and to produce new ones (Figure 3):

PERSPECTIVES

Throughout our case studies and our generalization, we presented a promising design space for aggregating streaming data using a smooth and continuous transition. Despite the positive user feedback we received, the technique still needs to be refined. Future work will consist in evaluating VS's applications to real data set and testing it on different space filling visualizations during an updating process.

Figure 1 :

 1 Figure 1: From left to right: (a) Our two case studies (Bubble-t and Bubble-TV), (b) illustration of sedimentary rocks [1] compared to our visual sedimentation process, and (c) three types of visual sedimentation charts: bar chart, pie chart and bubble chart.

Figure 2 :

 2 Figure 2: Details of the sedimentation process.

Figure 3 :

 3 Figure 3: Example of generated visualizations (with their walls and barriers as black lines).