Laurent Duraffourg

Eric Colinet

Sebastien Hentz

Eric Ollier

Philippe Andreucci

Bruno Reig

Philippe Robert

Compact and explicit physical model for lateral metal-oxide-semiconductor field-effect transistor with nanoelectromechanical system based resonant gate

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

Field effect transduction is certainly among the older techniques used in micro electromechanical systems (MEMS) [1]. However, it had rapidly been supplanted by capacitive techniques that use simple MEMS technologies. It is commonly admitted that capacitive detection with a low noise amplifier exhibits ultra low noise. Today, these noise levels are sufficiently low to reach the resolution needed for typical MEMS sensors such as inertial sensors, or pressure sensors [START_REF] Yazdi | Proceedings of the IEEE[END_REF]. Now, NEMS (Nano-Electro-Mechanical-Systems) devices are actively investigated because of their physical properties resulting from ultra miniature size elements [START_REF] Roukes | [END_REF]. Furthermore, they offer the opportunity to integrate mechanical structures and CMOS devices on the same die. NEMS advantages include ultra low power consumption, potential high resonant frequency, and high sensitivity to applied force, external damping or additional mass [4] [5].

Here, we study a detection technique consisting in a nano structure oscillating along the channel of the MOSFET. The lateral resonant gate MOSFET (LRG-MOSFET) and the CMOS circuit can be fabricated on the same die through thin SOI technology and SON (Silicon on Nothing) technology [6]. The same thin layer of single-crystal silicon is used for both mechanical structures and advanced ICs. In the resonant gate MOSFET, the gate moved along the channel width [7]. The drain current was then easily calculated by considering N basic transistors placed in parallel along its width. This paper is devoted to the modelling of LRG-MOSFET. As the suspended gate moves along the channel length, the developed approach in [7] is no longer usable and another electromechanical model has to be developed.

A self-consistent model computing both mechanical gate deflection and surface potential onto the channel is therefore presented. The electrical model is based on an explicit formulation of the surface potential to determine all charge densities in the MOSFET. The fringe effect and the mode shape of the beam are included in the computation of the electrostatic forces.

Theoretical results are finally compared with raw data measured on a device. The cubic term α 3 is the Duffing effect, Q is the quality factor, ω 1 is the free resonance frequency of the first mode. R F and L F correspond to the electrostatic forces acting on each

( ) L R F F z Q z z z - = + + + 1 1 3 1 3 1 2 1 1 & & & ω α ω (1)
side of the beam. Labels L and R are for left and right sides (see

Figure ). { } ∞ = , 1 ) ( n n x χ follows the normalisation condition: nm L m n dx x x S G δ χ χ ρ = ∫ 0 ) ( ) ( (2)
The electrostatic force densities are projected on

) ( 1 x χ : dx x z C V F E L L L L ∫ ∂ ∂ = 0 1 1 2 ) ( 2 χ with E g L V V V - = (3)
∫ ∂ ∂ = MOS L R R R dx x z C V F 0 1 1 2 ) ( 2 1 χ with )) ( ( x z V V V oxide g R - = (4)
V g , V E , V oxide are respectively the gate voltage, the electrode voltage and the voltage between the surface of the oxide and the ground. L C and R C are the air gap capacitances per length unit:

)

) ( ( ) ( 0 , x z g t C x C Si n R L ± = ε (5) 0 ε is the vacuum permittivity. n
C is the fringe effect factor [START_REF] Nishiyama | [END_REF]. R F is a bit more complex than L F because the potential difference depends explicitly on z(x) through V oxide , which depends on the capacitive bridge

R MOS g C C V + 1
. The classical MOSFET equations for each z-displacement are used to define MOS C .

First, the surface potential ) (x s ϕ is computed with an explicit model, which is obtained by asymptotic simplifications of the implicit equation for each operating regime [10]:

        + -       -         +       -         - + =         - - T k V q T k q q T k T k q q T k V V B Fi B s B B s B s s FB g ) 2 ( exp 1 exp 1 exp 2 φ ϕ ϕ ϕ γ ϕ (6) 2 0 2 i a Si C N q ε ε γ = , where 1 1 -         + = R ox ox i C t C
ε C i is the dielectric capacitance due to the oxide layer and the variable air gap.

( )

Fi s i ss i i ms FB C Q C Q V φ ϕ φ - - - =
is the flat band voltage. Q i and Q ss are respectively the charge density in the gate oxide and the charge density at the interface oxide / silicon. V(x) is the potential in the channel depending on V ds . q, N a , ε Si , k B , T and Fi φ are respectively the electron charge, the doping level of the channel, the relative dielectric permittivity, the Boltzmann constant, the temperature and the Fermi potential.

Once s ϕ known, the charge densities in the MOSFET and the MOS capacitance can be evaluated to determine oxide voltage and drain current:

( ) )) ( ( ) ( 1 ) ( exp ) ( ) ( x V V C x Q V V V V T k x q q T k x C x Q FB g FB g B s B S i dep - - - = - -       -         - + - = ϕ ϕ ϕ γ (7) 1 1 ) ( ; ) ( ) ( -         + = ∂ ∂ = ox ox sc MOS s sc sc t C x C x Q x C ε ϕ (8) Where, ) (x Q dep , ) (x Q inv and ) (x Q sc
are respectively the depletion charge density, the inversion charge density and the charge density in the gate. The drain current is the sum of the inversion charge density along the channel.

∫ - = d s V V inv eff dV V Q L W Id ) ( µ ( 9
)
µ eff is the effective mobility of electron (N-MOSFET). For each operating point (V g , V ds ), µ eff is approximated by the Mathiessen rule considering the phonon interaction, the Coulomb interaction and the surface scattering [11].

Finally, the model depends on the geometry and a set of physical parameters (µ eff , Q i , Q ss and N a ). The equations ( 1) to [START_REF] Landau | Theory of Elasticity[END_REF] have to be solved self consistently through a robust iterative algorithm based on the Runge-Kutta formulation.

To validate our theoretical model, we compare the simulation results with measurements reported in [12]. Devices were fabricated with a technology based on 200mm tools with a SON (Silicon on Nothing) approach [7]. This technology takes advantage of the single-crystal silicon for the NEMS structure without using SOI wafers. This technology is compatible with a front-end CMOS process. The Figure 2 shows a SEM micrograph of a typical device.

Let us consider a device (Figure 2) on which static measurement were carried out.

The measured characteristics are L G =16.1 µm, w=490 nm, g=107 nm, L MOS = 9 µm, W MOS =400 nm, t ox =2 nm, and N a = 5.10 15 at/cm 3 [12]. Figure 3 a) and b) give respectively the static characteristics I d (V g ) and g m (V g ) with V E =V S =0. The theoretical results are superimposed in each case. Notice that the bulk voltage (body potential) remains floating. To fit the model with the experimental results, the doping level N a , and the charge density in the gate oxide Q i are respectively fixed at 4.5.10 15 at/cm 3 and 5.10 10 cm -2 since the charge density Q os at the interface oxide / silicon is tuned. This parameter depends on the technology (etching step and oxide deposition) as well as the voltages. In the Figure 3 a), theory and measurement are well in agreement. Threshold voltages are quite similar around 2 V. Below the threshold, in the accumulation regime, the experimental device exhibits a quite strong leakage current from 0.1 µm to 1.3µA. This current variation depends linearly on V ds at V g =0V. It may be attributed to a photolithography misalignment when protecting channel from phosphorous implants generating a 2.1 MΩ short-circuit resistor. This process error can not be reproduced by the model and the resistance was added by hand in parallel of the transistor. For V ds =1.55V and 2.75V, the charge density Q ss is respectively fixed at 10 12 cm -2 and 10 13 cm -2 . The charge traps at the interface oxide / etched silicon may be gradually filled according to the surface potential, which depends both on V gs and V ds . This additional charge induces a modification of the flat band voltage and the vertical electrical field leading to both threshold voltage variation and slop variation. In the Figure 3Figure b), we observe a quite large dispersion of raw data which makes the comparison difficult. Anyway, the computed values are in the good order of magnitude. The maximum of transconductance, is around 1 µS for (V g =6.5V, V ds =2.75V) leading to an I d current of about 3 µA. For this operating point, µ eff was evaluated at 300 cm 2 /V.s.

Let us consider a device on which dynamic resonance were observed. The measured characteristics are L G =10 µm, w=165 nm, g=120 nm, L MOS = 6.8 µm, W MOS =400 nm, t ox =2 nm, and N a = 5.10 15 at/cm 3 [12]. RF characterizations were performed using a measurement bench with a vectorial network analyser Agilent 8753E (VNA) to measure the S 12 parameter (ratio of the transmitted power over the incident power in the logarithmic scale). The values of the bias voltages applied on the component were extracted from the static characteristics I d (V g ) and I d (V ds ). Electrode voltage V dc , gate voltage V g and drain voltage V ds were respectively V dc =10V, V g =4.6V, V ds =3V. The input power on the electrode was fixed at -41 dBm (ac-voltage around 4.5 mV RMS over 50 Ω). Pressure was fixed at 10 -6 Torr -The quality factor Q was evaluated to ~700 that seems to be a typical value for such a vibrating beam. A bulk pad is added to avoid any body effect. The gap reduction implies to take into account the Casimir force [13] that may have an impact at tiny gap. This force is already included in the model. Compared with other models, our model is more comprehensive since it computes the charges densities in the channel according to the mechanical oscillations. It may finally be suitable for any electromechanical components based on a field effect transduction.

Authors acknowledge financial support from the French ANR -NANORES contract.

1 Figure .

 1. Figure . Second Newton equation is applied in the Galilean referential frame of the

  Our AC-simulation is based on the model presented in precedent section considering the experimental operating point. The total noise floor of the measurement chain loaded with the component is calculated to know the theoretical background level of the S 12 parameter. This noise results from three main noise sources (uncorrelated noise assumption): the white mechanical noise of the beam, the thermal noise current of the MOSFET and the VNA noise. The mechanical noise is computed using theoretical characteristics of the beam. The electrical noise is computed from the transconductance g m , and the source gate capacitance C gs with our electrical model. g m , C gs and µ eff are respectively evaluated at 200 nS, 0.2 fF and 100 cm 2 /V.s. The VNA noise is known through an open loop calibration without device (~380 pA/√Hz). Theoretical results and measurements are shown in the Figure 4. The model is well in agreement with the experiment despite the only use of physical parameters. The theoretical resonance frequency is quite close to the experimental frequency. The difference between the theoretical level and the measured level is only of 1 dB higher. The background shown in the Figure 4 is only due to the VNA. A low distortion of the experimental resonance peak, which is not anticipated by the model, is also visible. Its origin is not yet explained.We detailed a unified electromechanical model for a LRG-MOSFET only based on physical assumptions. No any extra mathematical parameters to fit the behavior of the device are thus required. For static study, theory seems to be quite well in agreement with the measurement. The visible dispersion might be a poor control at the bulk potential because the body of the MOSFET was kept floating. Other reasons could be a strong surface density of charge traps and a strong roughness of the interface silicon/oxide inducing a large instability of the measurements. Next generation of component is in progress to improve the field effect and the responsivity to the beam movement. Thus, the gap and the channel length are scaled down.

Figure 1 :

 1 Figure 1: Schematic of the LRG-MOSFET -F L and F R are the electrostatic forces on each

Figure 2 :

 2 Figure 2: SEM microphotograph of the LRG-MOSFET [Ref 12]

Figure 3 :

 3 Figure 3: Static characteristics of the LRG-MOSFET -a) Id(Vg) -b) Transconduction

Figure 4 :

 4 Figure 4: Measured S 12 parameter and S 12 computed with our electromechanical model


Eric Colinet

Laurent Duraffourg

Sébastien Labarthe

Philippe Andreucci

Sébastien Hentz

Philippe Robert

Self-oscillation conditions of a resonant-nanoelectromechanical mass sensor

This article presents a comprehensive study and design methodology of co-integrated oscillators for nano mass sensing application based on resonant Nano-Electro-Mechanical-System (NEMS). In particular, it reports the capacitive with the piezoresistive transduction schemes in terms of the overall sensor performance. The developed model is clearly in accordance with the general experimental observations obtained for NEMS-based mass detection. The piezoresistive devices are much sensitive (up to 10 zg/√Hz) than capacitive ones (close to 100 zg/√Hz) since they can work at higher frequency. Moreover, the high doped silicon piezoresistive gauge, which is of a great interest for very large scale integration displays similar theoretical resolution than the metallic gauge already used experimentally.

Introduction

In the last several years, nano-electromechanical systems (NEMS) resonators have shown their strong potential as highly resolved nano-sensor for gas sensing and ultra low mass detection applications [1], [START_REF] Verd | Monolithic CMOS MEMS Oscillator circuit for sensing in the attogram range[END_REF], [3]. In these resonant applications, the adsorbed mass on top of the NEMS induces a shift of its resonant frequency that is continuously measured. To follow this frequency variation, the NEMS is embedded in a closed loop that can be either a phase-locked loop (PLL) or a self-oscillating loop [1]. In such architectures, the mechanical amplitude of the vibrating structure is usually kept constant to ensure the stability of the system. Our work is focused on the self-excited scheme because it is a compact co-integrated mass sensor system. Self-oscillating state is established if the electronic gain supplied by the loop is large enough to compensate the losses of the measurement chain. The electronics must also compensate the phase shift induced by the NEMS.

The mass-resolution of such a device is limited by different noise processes such as thermomechanical, Nyquist-Johnson, or adsorption-desorption noises. It is also limited by the electronics sensitivity. Other noise sources as thermal fluctuations and defect motion-induced noise can also be relevant when the targeted ultimate mass is close to some hundred Daltons.

The stability of the self-excited architecture must then be carefully studied to know the expected performance of the resonant mass sensor. Basically, the needed electronic gain to reach the stable state of the device depends strongly on the transduction technique. In this article, a theoretical study of a NEMS embedded in a self-excited loop with capacitive [START_REF] Verd | Monolithic CMOS MEMS Oscillator circuit for sensing in the attogram range[END_REF] [4], or piezoresistive (PZR) [5] detection principles is presented. Field effect transduction [6] that may also be used can be treated in the same manner as the capacitive detection.

The same NEMS geometry is used to estimate the performance of the two detection principles. There are many factors that determine the mass sensitivity for a given detection. To make a comprehensive comparison, the NEMS is interfaced at its optimal operating point, e.g. maximum bias voltage taking into account to the physical limitations (pull-in for the capacitive approach, self-heating of the gauges for the piezoresistive way).

The NEMS is based on a vibrating nano cantilever moving in the wafer plane [START_REF] Ph | French Patent EN07 55992[END_REF].

The cantilever is actuated with lateral electrodes. To perform the PZR detection, gauges are placed on each sides of the cantilever close to the anchor. The gauge are then strained or extended when the beam vibrates. For the capacitive detection is achieved with a second electrode.

The article is organized as follows: The first part the oscillator general topology is presented. The practical implementation of the mechanical structure is then detailed. This structure can be used in capacitive detection or in piezoresistive detection as well. The second part deals with the so-called "capacitive oscillator" where the detection and the actuation of the resonant NEMS are purely electrostatic. In the third part, the "piezoresistive oscillator" based on electrostatic actuation and piezoresistive detection is presented. Finally, ultimate performance (mass resolution) of each oscillator is evaluated.

General topology of an oscillator

Whatever the transduction principle, the oscillator can be described by the general topology shown in Fig 1 . The mechanical resonator displacement is converted into an electrical signal by an impedance variation Z through either capacitive or PZR detection. This signal is amplified and phase shifted to fullfil the Barkhausen condition [START_REF] Hajimiri | The Design of Low Noise Oscillators[END_REF] and then feedbacked to the NEMS through the electrostatic actuation.

Oscillation conditions

To realize a NEMS-based oscillator, it is critical to adjust the open-loop gain and phase changes to satisfy the Barkhausen criterion on resonant 1  ,

                 0 ) ( arg 1 ) ( 1 ) ( ) ( ) ( ) ( 1 1 1 1 1 1 j H j H j G j H j Vin Vout j H OL OL OL (1)
where OL H is the overall transfer function in open loop configuration.

    ) ( ) ( j J j H 
is the transfer function of the NEMS. J , ) (   j and  are the detection gain, the NEMS dynamic and the actuation gain respectively. G(j) is the electronic transfer function to be adjusted.

Whatever the mechanical resonator considered (beam, nanowire, disk…), its dynamical displacement can be reduced to an equivalent mass/spring/damper scheme as follows [START_REF] Lin | [END_REF],

Q j M B j K M j eff eff eff eff / / 1 1 ) ( 1 2 1 2 2 0                 (2)
where M eff , K eff and B eff are the effective mass, the effective stiffness and the effective damping coefficient respectively. The values of M eff , K eff , which depend on both NEMS geometrical parameters and anchorage conditions, define the eigen free frequency

eff eff M K   1
. B eff is related to the quality factor Q=M eff 1  /B eff . For NEMS Q is typically close to 1000 in vacuum and 100 in air.

Using (2), ( 1) is reduced the following expression,

                        90 ) ( arg ) ( 0 ) ( arg 1 ) ( 1 2 1 1 1 1 j G JQ M j G j H j H eff OL OL  (3)
where the product J  is considered positive.

Phase noise

The oscillator stability is characterized by its frequency noise   i.e. the frequency fluctuation around its eigen frequency 1  .   is related to the frequency noise spectral density  S via the Parseval-Plancherel theorem,

         ) ( 1 2 S ( 4
)
where  is the measurement bandwidth. The maximal bandwidth is set by the NEMS bandwidth Q 1  . To increase the performance, the signal must be integrated over several periods to have for instance 1 Hz-bandwidth resulting in a relative low response time. is the sum of all noise sources (white thermomechanical noise, Johnson noise, electronic noise) referred to the resonator mechanical displacement. y is the RMS value of the beam deflection, which can be controlled through electrical non linearity (e.g. saturation of the electronic gain or by controlling the loop gain).

Considering that the noises are not correlated the overall noise density of the open loop is the sum all noise densities, The mathematical expressions of each noise density have to be defined according to the transduction principle and will be described further in the next part. The equation (5) shows that the overall noise density is minimized for large displacement y. In other words, the dissipated power in the vibrating beam must be larger as possible to limit the frequency noise.

) ( ) ( ) ( ) ( 1 1 1 1         n n n n
The mass resolution (minimum mass detectable) M  can be deduced from the phase noise,

1 2      eff M M (7)
The effective mass must be low and the resonant frequency must be as high as possible. These characteristics are fulfilled by NEMS structures.

Practical implementation of the NEMS oscillator

Architecture

The considered structure is a cross beam as shown in Fig. 2. Table I gives the typical feature sizes of the device. This beam used as a lever arm is actuated by the electrostatic force. The detection of its motion is performed either by two perpendicular thin piezoresistive gauges or by the capacitive variation between the beam and a second electrode. This architecture is versatile and can be used in both configurations. The geometric parameters and the material properties remain the same in both configurations.

According to equation (3)  and J have to be defined for each detection case. The electronics transfer function will be designed according to  and J. The electronics is based on the same Pierce scheme as shown in Fig. 3.

Transduction modeling

Second Newton equation is applied in the Galilean referential frame of the substrate to the cross-beam. The cross-beam is assumed to follow Euler-Bernoulli equation [START_REF] Landau | Theory of Elasticity[END_REF]. The equation is reduced to its normalized lumped expansion on the first mode through the Galerkin method [START_REF] Ern | Theory and practice of finite elements[END_REF]. Here, we only consider the first eigen modal vector

) ( 1 x



, which observes the boundary conditions

) ( 1 1 l    = 0 ) ( 1 1  l  and ) ( ' ' 1 1 l l    = 0 ) ( ' ' 1 1   l l  . ) ( 1 x
 must also follow the mechanical continuity at the gauge position

) 0 ( 1    = ) 0 ( 1    and ) 0 ( ' 1    = ) 0 ( ' 1   
. The equivalent electromechanical response of the cross-beam in the "spring/mass/damper" model is expressed as,

) ( 1 ) ( 1 2 2 1 Q j M eff          (8) where 1 1 1 ) (    Sl l l a M eff  
is the effective mass of the resonant beam for the first eigen mode

2 1 1 k S EI    . b S Q 1   
is the quality factor of the structure (b is the damping coefficient).

     1 1 ) ( 1 1 2 1 l l a l l dx x l  
is a constant depending on the modal base and on the force repartition along the beam. S, , E and I are the section of the cross-beam, the density of the silicon (2330 kg/m 3 ), the Young modulus (169 GPa) and the quadratic moment respectively.

) (

1 1 l l  
is the value of the modal function at the beam end. k 1 is the wave vector that depends on the distance l 1 between the anchor of the lever beam and the gauge. This value has been set to enhance the stress inside the gauge sensors due to a mechanical lever effect and was fixed at l 1 =0.15l. In that case k 1 =2.12/l. The actuation is made through an electrostatic force, which is uniform along the electrode. Within the limit of small deflection compared to the gap g and for AC-voltages V ac smaller than the DC-voltage V dc , the force per unit length at =

1 is, ) cos( ) ( 1 2 0 t g V eaV t f ac dc    (9)
where  0 is the vacuum permittivity (8.85x10 -12 F.m -1 ).

From [START_REF] Lin | [END_REF], the actuation gain can be written as follows,

Capacitive oscillator

Detection gain

The displacement y()=()F() implies a capacitance variation between the lever arm and the electrode. The detection gain corresponding to the conversion of the displacement y to the capacitance variation C is then calculated as,

) ( ) ( ) ( 1 2 1 2 0 l l g el y C           (11)

Circuit modeling

The figure 4 shows the equivalent small signal circuit of the ensemble variable ). The fringe effect can be included using the results presented in Ref [START_REF] Nishiyama | [END_REF] and will increase the absolute value of some units. The electronic gain using Millman's theorem are defined at nodes 1 and 2,

) ( 0 2 1 C C C j jC V C j V V in fb fb dc         (12) L ds L fb m fb R g C C j V g jC V / 1 ) ( ) ( 1 2        (13)
At this stage some assumptions on the relative values of R L , C L , C fb , and C in have to be made. In particular, the input current must flow through the capacitance C in and the output current must flow through the capacitance C L ,

1 0 1 ) ( and ) 1 ( 1          C C C C C g g R R C in fb L fb m ds L L L L (14)
where L  is the cut-off frequency of the electronics.

From ( 13) and ( 14), one can deduce,

L m out C j V g V V  1 2    (15)
Using ( 12) and (15), G() is,

     ) ( ) ( ) ( ) ( 0 2 C C C jC V g C V G in fb L dc m      (16)
The minimal transconductance g m is defined from the oscillation conditions (3), and the expression ( 16),

1 3 2 2 0 2 1 0 4 3 1 3 1 0 ) ( ) ( ) (     l e QV l l C C C g C M Q J V M C C C C g dc in fb L eff dc eff in fb L m          (17)
The conditions ( 14) can be simplified using (17) and expressing the conductance g ds as a linear function of g m . Considering a MOSFET in saturation regime it is known that,

m E m th G ds g L V g V V g     2 (18)
V G , V th V E and L are the gate voltage, the threshold voltage, the Early voltage, and the channel length of the equivalent MOSFET respectively. For the 0.35 µm-CMOS technology, their typical values are 3.3V, 0.485V, 21V/µm and 5 µm respectively. Finally, the conditions ( 14) are,

1 0 2 1 1 , 1 ) ( 1         Q J V C C C M C R dc in fb eff L L   (19)

Overall noise

To evaluate the mass resolution, the noise sources must be considered (see equation ( 5)).

Here, the white thermomechanical noise and the electronics noise are estimated in a very straightforward manner.

Hence in (m²/Hz) From the equations ( 12) and ( 14), the noise density of the electronic (in m²/Hz) is as follows,

                     m B L m fb in fb dc y s Electronic g T k C g C C C C V S n 3 8 ) ( 1 ) ( 2 2 1 2 2 2 0 2 2   (22)         m B g T k 3 8
is the voltage noise density at the input of the electronics.

Regardless of the NEMS geometry, the critical parameters to reach self-oscillations (16) and optimal mass resolution (5), ( 22) are the DC-voltage, the gap, the capacitances and the deflection amplitude. Discussion on the system performance is done in the last part of the article.

Piezoresistive oscillator

Detection gain

The well known piezoresitive equation for the axial stress

 L is [15], ) 2 1 (        L d R dR (23) L L L L E d         (24)
L  is the longitudinal piezoresistive coefficient along the gauge. Substituting (24) into (23), we get the usual expression,

L L L L L E R dR           ) 2 1 ( (25)  L is the gauge factor.
The first part of (25) corresponds to the gauge elongation. The second one takes into account the intrinsic piezoresistive coefficient that is predominant for semiconductor gauge. Its large value permits to attain strong sensitivity. The drawback is its large resistance due to the small section of the gauge that is more or less doped (over a range of 10 15 cm -3 to 10 19 cm -3 ). This transduction might not be suitable for very tiny mass detection because of the strong Johnson noise floor. A metallic layers can also be used as the gauge layers (for instance coated thin layer on the silicon beam). In this case the elongation effect is only used. Its advantage is clearly the low resistance with a very low Johnson noise [16]. However, the gauge factor is only a function of the Poisson ratio and its value is reduced to some units.

According to equation (24), the relative resistance variation is proportional to the axial stress σ L exerted by the resonant lever beam. It is given by the difference of shear forces ) 0 ( and ) 0 (

        x V V x V V
. After few mathematical developments, σ L can be expressed as,

y s EIk ew V V j L 3 1        (26)
s is the section of the gauges.  is a constant that depends on the boundary and normalization conditions.

The displacement y induces the axial stress and implies a relative resistance variation leading to the PZR detection gain,

) ( / 1 1 3 1 l l s Ik y R R L         (27)
 depends on the mode shape of the lever beam, the geometry, and the gauge factor.

The gauge factor L  depends on the doping level and is between 10 and 100 [15].

Circuit modeling

Figure 5 shows the equivalent small signal circuit of the ensemble variable resistance and its electronics. The expression of the transfer function V out /R is straightforward using the Millman's theorem at nodes 1 and 2,  

) ( 2
) ( 2 2 2 1 in fb fb b C C j R C j V R R V V        (28) and, L L fb fb m out R gds C C j V C j V g V V 1 ) ( 1 1 2          (29)
R (Z in Fig. 1) is the resistance variation around its nominal value R. V b is the read voltage applied on the gauges.

At this stage some assumptions on the relative values of R L , C L , C fb , and C in have to be made. In particular, the input current must flow through the resistance R. In practice, typical resistance values of gauge obtained in this work are close to 3 k. C in is close to 10 fF in a co-integration approach to 1 pF when a stand-alone circuit is used. It means that the assumption is true up to 150 MHz for the worst case. Moreover the output current must go through the output capacitance C L . These conditions are summed up as,

)) 1 ( ( 1 1 L ds L L R g R C    and fb m C g 1   (30)
From ( 29) and (30), one can deduce:

L m out C j V g V V  1 2    (31)
Using ( 31) and ( 28

), G() is, ) 2 ( 2 ) ( ) ( m fb L b m out Rg C C j V g R R V G         (32)
The minimal transconductance g m is defined from the oscillation conditions (3), and the expression (32),

     3 1 0 1 1 2 3 1 3 1 ) ( Ik ea V QV l l s g C M Q J V C M g L b dc L eff dc L eff m      (33)
As in the capacitive case, the conditions (30) can be simplified using (33) and expressing the conductance g ds as a linear function of g m :

1 2 1 1 , 1 1       Q J V M C R dc eff L L   (34) and 1 2 3 1   R M Q J V C eff dc fb  (35)

Noise modelling

The thermomechanical noise remains unchanged (see equation [START_REF] Kacem | Twelfth Conf. on Nonlinear Vibrations, Dynamics, and Multibody Systems[END_REF]). The noise from the gauge is the Johnson noise. Its expression in (m²/Hz) is deduced from the equations ( 27) and (28),

2 2 2 2 1 4 ) (           L m fb dc B y Johnson C Rg C V TR k S n   (33)
where R is the sum of the gauge resistance and the lead resistance. The electronics noise in (m²/Hz) is for a MOSFET in strong inversion,

2 2 2 2 1 3 8 ) (           L m fb dc m B y s Electronic C Rg C V g T k S n   (34)
where g m is the transconductance according to the equation (32).

The model does not take into account the 1/f-noise due to low charge carrier density in the semiconductor material. However, we considered high doped silicon gauges (at least 10 19 cm -3 ) built with a typical CMOS process, which reduces the 1/f-noise. Discussion on the system performance is done in the last part of the article.

Notice that all characteristics driving the actuation gains are the same for the two oscillators. A priori, the only design parameters are R L and C L that should be adjusted for each mechanical resonant frequency to optimize the oscillation conditions. As shown in Tab. III, R L and C L have the same impact on the gain and phase conditions for each technique. These parameters are set at reasonable values so that

1 1 1   L L C R
up to 100 MHz. Clearly, the upper working frequency limit of the two oscillators could be optimized by fine tuning of these parameters. But this study is out of the framework of this paper.

For the capacitive transduction, the oscillation conditions (see table III) are shown in Fig. 6. The phase conditions are only observed if the oscillation frequency is below ~5 MHZ. Notice that the equation ( 19) follows a 1/f law at low frequency and a quadratic variation at high frequency. Gains J and  are independent of the frequency. The requested transconductance g m increases with the frequency. Stiffer mechanical structures demand larger energy level to reach self-oscillation. This simple law is well expressed by the equation (17). Over the usable frequency range, g m is always smaller than 2.10 -5 S. This value is rather easily reached with a co-integrated circuit (with a single aforementioned 0.35 µm-MOSFET). Similar results are presented in Fig. 7 for the semiconductor gauges. This transduction principle is more effective to get self-oscillation compared to the capacitive one. The useful transconductance g m for a semiconductor gauge varies over the 7-47 MHz frequency range from some 10 -4 S to 5.10 -4 S. Notice that the equation (34) follows a 1/f law at low frequency and becomes constant at high frequency. Gains J and  change according to a f 2 law.

Finally, the figure 8 gives the theoretical mass resolution for the two detection schemes. Table IV sums up all noise sources used in the model. Over 5 MHz-10 MHz frequency range, the resolutions are quite similar close to 10 zg for a 1 Hz-bandwidth.

For highest frequency up to 50 MHz, the expected resolution may reach 1 zg with the piezoresistive detection. These results are rather comparable with experimental mass performance already published [START_REF] Ilic | [END_REF] [21]. The difference between the two detection schemes (in particular the inflection point for the capacitive detection) is mainly due to the electronic noise. In particular, this noise is inversely proportional to the detection gain , which is constant with the capacitive detection and changes according to a f 2 -law with the piezoresistive detection.

Conclusion

We have done a theoretical comparison between two in-plane oscillating architectures based on piezoresistive and capacitive detections respectively. The sensors are embedded in a closed loop system with an electronic readout circuitry to define the mass-sensor. The overall measurement chains are studied through two specific models.

This work provides a comprehensive theoretical method that is valid for any capacitive or piezoresistive principles. It is a useful methodology to design compact oscillator circuit around NEMS application.

Considering that the capacitive and PZR oscillators are set in their optimal working PZR detection can work at higher frequency and offers therefore better mass resolution. This trend can be extended to any NEMS.

Moreover silicon gauge can reach similar mass resolution than metallic gauge commonly used elsewhere [21]. Although the Johnson noise is increased with the resistance value, the signal improvement is much higher than the noise worsening. This conclusion completes the analysis most widely held that considers the metallic gauge as the best candidate for ultimate mass resolution. The silicon piezoresistive detection is an alternative choice for the mass sensing application if the ambient temperature is not too low (freeze-out effect of carrier density).

To conclude, this work paves the way toward the very large scale integration of the piezoresistive NEMS with its co-integrated electronics.

Acknowledgement

Authors acknowledge financial support from the French Carnot -NEMS contract.

We also thank M.L. Roukes for fruitful discussions.

TAB. III. Oscillation conditions

Gain

Phase capacitive

Q J V M C C C C g dc eff in fb L m  3 1 0 ) (     1 0 2 1 1 , 1 ) ( 1         Q J V C C C M C R dc in fb eff L L   1 2 1   eff read fb M Q J V C  PZR Q J V C M g dc L eff m  3 1   1 2 1 1 , 1 1       Q J V M C R dc eff L L   1 2 3 1   R M Q J V C eff

From

  Robins's expression, the frequency noise spectral density is function of the open loop noise spectral density

  can now be studied for each detection principle.

  is the displacement of the cross-beam end in the Fourier space.

  capacitance and its electronics (seeFig 3.b). The electronics is reduced to an equivalent MOSFET transistor. The input capacitance is the gate/source capacitance in parallel with a certain coupling capacitance. The feedback capacitance is the gate/drain capacitance in parallel with the capacitance due to the actuation capacitance. C (Z in Fig.1) is the capacitance variation around its nominal value C 0

  and

  parasitic capacitances included in C in have to be minimized to optimize detection gain. The feedback capacitance is also a strong source of perturbation. C fb reduces the frequency range of the operation.

FIG. 1 .FIG. 5 .

 15 FIG.1. Closed loop representation of the self-exited system

Performance of piezoresistive and capacitive oscillators

We simulate the two kinds of oscillators according to the resonant frequency of the NEMS. The frequency is changed from 1 MHz to 50 MHz by varying the beam width w.

All other geometrical parameters are kept constant (Table I). To study in a proper way the capacitive and piezoresistive oscillators, both devices are settled in their optimal working point regime i.e based on their highest detection gain for each. For both oscillators, a monolithic integration is considered so that the parasitic capacitances are minimized [17].

The simulations are performed considering 0.35 µm-CMOS technology [START_REF] Verd | Monolithic CMOS MEMS Oscillator circuit for sensing in the attogram range[END_REF]. C in , and C fb are parasitic components. C fb will reasonably not exceed 5 fF (around 0.13 fF/µm for the 0.35 µm -CMOS technology) and C in will be close to 10 fF.

The actuation voltages V dc is the same for the two oscillators. We assume a deflection y =10 nm (~gap/10), well below than the non linear displacement limit and the dynamic pull-in [START_REF] Kacem | Twelfth Conf. on Nonlinear Vibrations, Dynamics, and Multibody Systems[END_REF]. For the piezoresistive detection, the readout voltage V j should be chosen as high as possible. In the same time, V j must be low enough to avoid the gauge breaking due to a large current flowing through [START_REF] Ch | Introduction to Solid State Physics[END_REF]. For the silicon gauge of ~3 k with a gauge factor of 100 for 10 19 cm -3 doping level [15], the voltage V j is set to 600 mV to keep the gauge temperature growth around 100 K. For the capacitive detection, the readout voltage is directly V dc and set to the maximum acceptable value for CMOS compatibility (V dc = 3.3 V). Table II gives the electrical parameters used in. the model.


Fabrice Muhlenbach
email: fabrice.muhlenbach@univ-st-etienne.fr

Stéphane Lallich
email: stephane.lallich@univ-lyon2.fr

A New Clustering Algorithm Based on Regions of Influence with Self-Detection of the Best Number of Clusters

Keywords:

Clustering methods usually require to know the best number of clusters, or another parameter, e.g. a threshold, which is not ever easy to provide. This paper proposes a new graph-based clustering method called "GBC" which detects automatically the best number of clusters, without requiring any other parameter. In this method based on regions of influence, a graph is constructed and the edges of the graph having the higher values are cut according to a hierarchical divisive procedure. An index is calculated from the size average of the cut edges which self-detects the more appropriate number of clusters. The results of GBC for 3 quality indices (Dunn, Silhouette and Davies-Bouldin) are compared with those of K-Means, Ward's hierarchical clustering method and DBSCAN on 8 benchmarks. The experiments show the good performance of GBC in the case of well separated clusters, even if the data are unbalanced, non-convex or with presence of outliers, whatever the shape of the clusters.

I. INTRODUCTION

Clustering is an unsupervised machine learning task used successfully in many fields, including image analysis, bioinformatics and marketing. Given a set of n data objects described by p attributes, clustering procedures aim at grouping the objects into clusters such that similar objects are in the same cluster and dissimilar objects are in separate clusters. In this paper, we are interested with crisp clustering methods which consider that each object belongs to exactly one cluster at the contrary of fuzzy clustering.

Partitioning clustering aims at directly decompose the set of objects into k disjoint subsets, where k has to be predetermined. K-Means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] is a popular iterative reallocation algorithm due to its simplicity and low complexity. Among the variants of K-Means, PAM [START_REF] Kaufman | Finding Groups in Data: an Introduction to Cluster Analysis[END_REF] is more robust against the noise and CLARANS [START_REF] Ng | Efficient and effective clustering methods for spatial data mining[END_REF] is more efficient and scalable.

Hierarchical algorithms produce a hierarchy of partitions. According to the direction of the procedure, at each step these algorithms either merge the most similar clusters (agglomerative procedure, e.g., using Ward's criterion [START_REF] Ward | Hierachical grouping to optimize an objective function[END_REF]) or split the most heterogeneous cluster (divisive procedure). Among others one can cite BIRCH [START_REF] Zhang | BIRCH: an efficient data clustering method for very large databases[END_REF] which concentrates on densely occupied portions of the space, CURE [START_REF] Guha | CURE: an efficient clustering algorithm for large databases[END_REF] which uses multiple representatives for each clusters to "capture" the shape of each one, and ROCK [START_REF]ROCK: A robust clustering algorithm for categorical attributes[END_REF], or its improved version Q-ROCK [START_REF] Dutta | QROCK: A quick version of the ROCK algorithm for clustering of categorical data[END_REF], that does not employ distances but links when it merges the clusters.

Density-based clustering algorithms are designed to discover arbitrary-shaped clusters. Clusters are considered as dense regions of objects in the data space that are separated by regions of low density. DBSCAN [START_REF] Ester | A densitybased algorithm for discovering clusters in large spatial databases with noise[END_REF] is a typical densitybased algorithm. The core point in DBSCAN is that for a fixed radius, the neighborhood of each object in a cluster has to contain at least a minimum number of other objects.

Grid-based algorithms partition the data space into a finite number of cells to form a grid structure, and determine the cells whose density exceed a prefixed threshold. Then they perform all clustering operations on the obtained grid structure (e.g., STING [START_REF] Wang | Sting: A statistical information grid approach to spatial data mining[END_REF] or CLIQUE [START_REF] Agrawal | Automatic subspace clustering of high dimensional data for data mining applications[END_REF]).

Graph-theoretic approaches are able to discover non usual data structures. The simplest algorithm is founded upon the minimum spanning tree (MST) [START_REF] Zahn | Graph-theoretical methods for detecting and describing Gestalt clusters[END_REF]. First the minimum spanning tree is built, where the nodes are the objects and the edges sizes are the Euclidean distances in the objects space. The edges having the greatest size are removed to construct disjoint connected components.

The purpose of our work is to extend Zahn [START_REF] Zahn | Graph-theoretical methods for detecting and describing Gestalt clusters[END_REF] and Urquhart [START_REF] Urquhart | Graph theoretical clustering based on limited neighbourhood sets[END_REF] graph-theoretic approach by showing an easy way of selecting the "ideal" number of clusters. In section II, we discuss the different state-of-the-art clustering methods and investigate how they process to choose the best number of clusters. We will show that no one can give an appropriate answer because if they do not ask explicitly the user to give them this information, they asked him on other parameters that he can not give a response without having an intuition on the domain. This problem is fixed in our new method, presented in section III, which is a divisive data set processing that tries to find automatically the good number of clusters. In the section IV we perform an extensive experimental evaluation with different clustering quality indices obtained on different clustering methods.

II. HOW TO CHOOSE THE NUMBER OF CLUSTERS?

On Table I, we list the characteristics of main clustering methods found in the data mining literature.

Algorithm Category

Method

Parameters / Properties Hierarchical clustering Ward's algorithm [START_REF] Ward | Hierachical grouping to optimize an objective function[END_REF] agglomerative algorithm MST Divisive based on graph theory Hierarchical Clustering [START_REF] Zahn | Graph-theoretical methods for detecting and describing Gestalt clusters[END_REF] Clustering Using each cluster is represented REpresentatives (CURE) [START_REF] Guha | CURE: an efficient clustering algorithm for large databases[END_REF] by a set of representatives RObust Clustering k : number of clusters using linKs (ROCK) [START_REF]ROCK: A robust clustering algorithm for categorical attributes[END_REF] Q-ROCK [START_REF] Dutta | QROCK: A quick version of the ROCK algorithm for clustering of categorical data[END_REF] θ: similarity threshold Hard Clustering k-Means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], [START_REF] Lloyd | Least squares quantization in PCM[END_REF] k : number of clusters Density-based clustering DBSCAN [START_REF] Ester | A densitybased algorithm for discovering clusters in large spatial databases with noise[END_REF] : distance to consider that 2 points are neighbors or not Sequential clustering Basic Sequential Algorithm Θ: dissimilarity threshold Scheme (BSAS) [START_REF] Theodoridis | Pattern Recognition[END_REF] k : max. number of clusters

Table I MAIN CLUSTERING METHODS

For the hierarchical clustering algorithms, the choice of the best number of clusters is not a simple task. Of course, a dendogram obtained with Ward's agglomerative method or the MST divisive algorithm can help the user, but in many cases it is only the knowledge of the domain that can be used to identify what is the best clustering. When no information is given to set k , the number of clusters, the user have to provide a specific parameter, like a similarity (or dissimilarity) threshold Θ, or to rely on some heuristic arguments, like a k = f (n) function that will provide a number of clusters depending on n, the size of the data set.

The density-based algorithms, like DBSCAN [START_REF] Ester | A densitybased algorithm for discovering clusters in large spatial databases with noise[END_REF], DB-CLASD or DENCLUE, consider that clusters are dense regions of points. These methods can handle arbitrary shaped clusters, outliers, and their time complexity is lower than other clustering methods. Unfortunately, the choice of the parameters required for these methods (e.g., for DBSCAN) may lead to different results and it is not easy to discover which one is the best without doing several experiments.

The clustering algorithms based on regions of influence are an extension version of the MST divisive algorithm [START_REF] Zahn | Graph-theoretical methods for detecting and describing Gestalt clusters[END_REF] aforesaid. These methods, based on the graph theory, are capable of detecting clusters for various shapes. These algorithms can provide good results when the clusters are well separated in the representation space. The key idea of these algorithms is to construct a graph by connecting with an edge each pair of points that are alone in a given region of influence, e.g. the relative neighborhood graph (RNG) [START_REF] Toussaint | The relative neighbourhood graph of a finite planar set[END_REF] or the Gabriel graph (GG) [START_REF] Urquhart | Graph theoretical clustering based on limited neighbourhood sets[END_REF], then removing the edges that are inconsistent compared with their neighboring edges. An edge is considered as inconsistent if it is greater of q standard deviations (typically q = 2) of the mean of a given number of its neighbors, which makes subjective this procedure. Nevertheless, an advantage of these algorithms is that the results do not depend on the order in which the data are considered and no initial conditions are required.

Finally an ideal method does not exist. In our opinion, a good clustering method has to be easy to implement and to use (no specific parameter values have to be asked to the user), and has to be suitable for various shapes of clusters.

III. GBC METHOD A. Neighborhood Graphs

The neighborhood graphs, which are special tools of the computational geometry, can be used in the clustering algorithms based on regions of influence, but also in many other data mining tasks, and especially in the supervised machine learning [START_REF] Muhlenbach | Multivariate supervised discretization, a neighborhood graph approach[END_REF], [START_REF] Muhlenbach | Identifying and handling mislabelled instances[END_REF], [START_REF] Zighed | A statistical approach to class separability[END_REF], [START_REF] Toussaint | Geometric proximity graphs for improving nearest neighbor methods in instance-based learning and data mining[END_REF]. Such neighborhood structure can be for example the k-nearest neighbor, the Delaunay triangulation, the MST, the relative neighborhood graph (RNG), the Gabriel graph (GG).

For each graph, a specific condition is required, depending on a region of influence, to link two points with an edge. For the MST, the condition is to connect all vertices together with the minimal size of edges; for the RNG, the region of influence is a lune, the intersection of two hyperspheres centred on each pair of points; and for the GG, the region of influence is an hypersphere with each pair as a diameter.

B. GBC, a New Clustering Algorithm

GBC, the new clustering method proposed in this paper, is conducted in 2 phases. The first phase (in 10 steps) consists in doing a list of µ values which will be used in the formula 1 to detect the appropriate number of clusters and is conducted as written on Table II.

1 construction of a neighborhood graph N G 2 descent sorting (by size) of the edge set E of N G 3 initialization step: k ← 2, Σ k ← 0, and ne ← 1 4 cutting the edge emax ∈ E of the (sub-)graph with the higher value 5 adding the size of emax to Σ k , the sum of the cut edges at the level k 6 testing if the (sub-)graph with the edges E -{emax} is still connected 7 if the graph is still connected, increasing the number of edges:

ne ← ne + 1 8 if the graph is not connected, modifying:

k ← k + 1 (new level for having k clusters) ne ← 1 (re-initialization of the number of cut edges) µ k ← Σ k /ne (µ k is the average of the sizes of the cut edges) Σ k ← 0 (re-initialization of the sum of cut edge sizes) 9

E ← E -{emax} (remove the bigger edge from the set E) 10 back to step 4 by using the next maximal edge emax ∈ E while k < n

Table II GBC ALGORITHM

After this first phase, we can calculate the δ values for each level k as follows:

δ k = µ k -µ k+1 µ k + µ k+1 , ∀k = 2, . . . , n -1 . (1)
The maximal value of δ k is used to select k δ max , the ideal number of clusters in the data set. The second phase is similar to the first one, but the loop runs until k = k δ max in the step 10 (instead of k < n).

Notice that we can equally use the MST, RNG or GG for the neighborhood graph on the step 1 of the algorithm.

In the experiments, we did not find significant differences in the results, but following [START_REF] Urquhart | Graph theoretical clustering based on limited neighbourhood sets[END_REF], we recommend the RNG of Toussaint [START_REF] Toussaint | The relative neighbourhood graph of a finite planar set[END_REF] which is a structure that overcome some problems encountered with the MST.

C. Main Characteristics of δ

The first phase considers all the data as a unique set and tries to separate this set in subsets as it is done for the descendant hierarchical clustering method, but in a much lower time complexity thanks to the graph structure. At each level k, the size average of the cut edges µ k is calculated.

Since the set of edges E is sorted by descending size, µ k is a monotone descending function of k and can not be used as it is for finding the more relevant number of clusters.

The value µ k -µ k+1 represents the decreasing of µ while getting from a partition of k clusters to k + 1 clusters. The parameter δ k allows to normalize this difference by taking into account µ k and µ k+1 . The splitting procedure continues while this normalized difference δ k is increasing.

D. Properties of the Method

First, GBC is very sensitive to the outliers [START_REF] Hawkins | The Identification of Outliers[END_REF], because an outlier will be far from the other data in the representation space, and it will be detected as an independent cluster.

Secondly, GBC is limited to hard clustering: it will fail to detect different clusters if there is a recovery between them. If the clusters are not well-separated in the representation space, there will not be an edge with a big size between the clusters in the neighborhood graph.

IV. EXPERIMENTATIONS

A. Methodology

For this study, we have compared GBC with 3 other clustering methods: K-Means [1], Ward's hierarchical clustering method [START_REF] Ward | Hierachical grouping to optimize an objective function[END_REF] and DBSCAN [START_REF] Ester | A densitybased algorithm for discovering clusters in large spatial databases with noise[END_REF]. For DBSCAN, we have used the value of given by the mean value of the 4-nearest neighbors (4-NN) of the edges of the data set [START_REF] Ester | A densitybased algorithm for discovering clusters in large spatial databases with noise[END_REF], and we have increased this value (indicated or ) for having better results (and less clusters).

Three quality indices have been calculated on Table III: Dunn's [START_REF] Dunn | Well separated clusters and optimal fuzzy partitions[END_REF], Silhouette [START_REF] Rousseeuw | Silhouettes: a graphical aid to the interpretation and validation of cluster analysis[END_REF], and the Davies-Bouldin's [START_REF] Davies | A cluster separation measure[END_REF] indices. In addition to these results, Table III indicates for each data set the numbers of individuals (Indiv.), of attributes (Attr.), of clusters discovered with GBC (k δ max ), of real classes or clusters for the domain when it exists (k r ). The number of clusters used for the test (k u ), and the respective sizes of the clusters obtained with these 3 different methods are also shown on this table.

B. Quality Indices in Clustering 1) Dunn's Index: This index, which tries to obtain wellseparated sets with examples very closely located to each others, is calculated with the dissimilarity measure between two clusters and the diameter measure of a cluster. Large values of Dunn's index correspond to good cluster partitions. The Dunn's index has two major drawbacks: its computation requires a considerable amount of time, and this index is very sensitive to the noise and to the outliers.

2) Silhouette Index: The global silhouette partition, which is the average of all example silhouettes, takes its value between -1 and 1, and the partition with the maximum global silhouette is considered as the optimal partition.

3) Davies-Bouldin's Index (D.-B.): D.-B. reaches to minimize the average similarity between the different clusters, so D.-B. is a function of the ratio of the sum of withincluster scatter to between-cluster separation. D.-B. exhibits no trends with respect to the number of clusters. The smaller D.-B. is, the better the partition is considered.

C. Data Sets

For these experimentations, we used 5 real benchmarks from the UCI Repository [START_REF] Asuncion | UCI Machine Learning Repository[END_REF], 3 artificial data sets (test-2c-2o, test-random, yin-yang), and ruspini [START_REF] Ruspini | Numerical methods for fuzzy clustering[END_REF].

The data sets were chosen to represent a variety within a specific class of data sets characterized by: (1) small number of true classes, which may or may not correspond to coherent clusters; (2) moderate number of observations; (3) moderate number of features; (4) numerical attributes (continuous or multivalued values) for calculating the distances. The observations with missing values have been removed from the data sets, all values have been transformed with the Milligan and Cooper method [START_REF] Milligan | A study of standardization of variables in cluster analysis[END_REF] and the Euclidian distance has been used for all the clustering methods.

D. Global Results

On Table III, the best values of the clustering validity indices are emphasized in a bold font. The experimental setup allows to compare the performances of GBC, K-Means and Ward's hierarchical clustering 11 times, because we used the 8 data sets described in the section IV-C, and 3 of them were used with two different numbers of classes (auto-mpg, e-coli and iris). On Table III, the performance of each algorithm on each data set is evaluated according to the 3 quality indices described in the section IV-B.

In the case of the Dunn index, GBC outperforms K-Means and Ward's hierarchical clustering. Compared to K-Means, GBC improves significantly the value of the Dunn index 8 times out of 11 and there are 3 equalities (p-value of the sign test = 0.008). Compared to Ward's hierarchical clustering, GBC improves the value of the Dunn index 7 times out of 11, with 3 equalities and 1 defeat (p-value = 0.070).

According the Davies-Bouldin index, GBC has better results than K-Means (7 wins, 2 defeats, 2 equalities), but this superiority is not significant (p-value = 0.180). The results of GBC and WHC are quite the same (5 wins, 4 defeats, 2 equalities).

The performances of the three algorithms according to the Silhouette index are very close (GBC wins 3 times, K-Means 4 times, WHC 2 times, and 2 equalities). It is not easy to compare the performances of the three algorithms with those of DBSCAN because DBSCAN determines automatically the number of classes in function of the values chosen for the two parameters of the algorithm. Only 4 comparisons can be achieved (iris, ruspini, test-2c-2o, yin-yang). For each of them, the results of DBSCAN are identical to those of GBC.

E. Detailed Results

Auto MPG is generally used for a regression task: the continuous class variable is to be predicted in terms of 3 multivalued discrete and 5 continuous attributes. Consequently, there is no "true real" number of classes to discover. GBC is not able to find an ideal number of clusters (k δ max = 388 392 = n) but a real number of cluster k r does not exist for this data set. GBC obtains however the best result with the Dunn index for k = 5 clusters.

Breast Cancer Wisconsin (Original) is a data set that contains 8 attributes and is used in a Boolean classification task. Here, GBC is not able to find an ideal number of clusters (336 seems to be irrelevant) but nevertheless, when it has to detect 2 clusters, it obtains good results with Dunn and D.-B. indices, just by isolating one outlier.

Ecoli Data Set contains 7 predictive attributes used to predict the cellular localization sites of proteins. They are 8 possible sites. With this data set, GBC does not detect the 8 classes but find 3 clusters; the 3 quality indices obtained with this method pass all the others.

Iris Data Set contains 4 predictive attributes with 3 classes of 50 instances each, where each class refers to a type of iris plant. One class is linearly separable from the other two; the latter are not linearly separable from each other. All the methods are able to discover the 2 clusters, but the GBC method gives an interesting result when it is forced to produce 3 clusters (one outlier detected, and best value with the D.-B. index).

Ruspini [START_REF] Ruspini | Numerical methods for fuzzy clustering[END_REF] is a data set consisting of 75 points in four groups with 2 variables giving the x and y coordinates. On this case, the results of the different methods are the same.

Test 2 Clusters with 2 Outliers (or Test-2c-2o) consists of 402 instances on 2 variables (the x and y coordinates) with 2 groups of 200 observations each and 2 outliers between the groups. Only GBC and DBSCAN are able to detect the 2 clusters and the 2 outliers.

Test with Random values (or Test-Random) is an artificial data set consisting of 1000 instances on 2 variables (the x and y coordinates) where the values have been randomly generated. The value of k δ max = 999 obtained is close to n = 1000, indicating that the method is not able to find a structure in the data set (which is an appropriate behavior with random data).

Yin and Yang is a non-convex artificial data set consisting of 1000 instances randomly extracted from the black points of a Yin and Yang picture (the Chinese symbol). The 2 variables correspond to the x and y coordinates. GBC and DBSCAN are able to find the appropriate clusters. The Silhouette index is not the best in this case because it penalizes a partition with unbalanced clusters.

V. CONCLUSION AND FUTURE WORK

The new clustering method that we proposed in this paper is in the line with the graph-theoretic approaches. GBC has been tested with 3 quality indices and compared to other clustering methods on 8 benchmarks. On artificial data sets where the clusters are well represented, GBC outperforms the other methods. On real data sets extracted from the UCI Repository, when clusters can intrinsically be found in the data, GBC detects automatically the correct cluster number and identify the relevant clusters. On the other hand, when the class clusters are not well separable, the method tends to "peel" the global clusters in isolated points.

Even if GBC is more time consuming than DBSCAN or K-Means, its computational complexity is equivalent to a hierarchical agglomerative clustering due to the neighborhood graph construction. We have also to remind that it is limited to the hard clustering and will fail if there is a recovery between the different clusters.

Nevertheless, this method has many advantages. First, we have shown that GBC obtains good results when the data are well-structured: it detects easily the well-formed clusters and the outliers, whether the cluster shapes are convex or not, whether the cluster sizes are homogeneous or not. Second, the main advantage of GBC is that the method does not need any parameter to perform on a data set. It is a considerable improvement compared with other clustering techniques developed in the data mining literature. Even for DBSCAN, which can discover the best number of clusters, the quality of the results is associated to one or more specific parameters (here, the size of the radius and the minimal number of neighborhood points for considering a core point).

Third, when they are some outliers in a data set, GBC can automatically find them, without any parameters, as well as DBSCAN when the parameter is well chosen. Fourth, GBC provides a dendogram -like the other hierarchical cluster algorithms-but can find automatically the ideal number of clusters. Fifth, it can be associated to a visualization method for neatly navigating into the data. This method is useful to prepare the data and to find the border points between clusters. Sixth, when there is no a-priori structure in the data of the set, GBC will indicate it by providing a value k δ max close to the number of examples in the data set. And seventh, while processing on a divisive way in a moderate computer complexity, GBC is better to represent the way of a human being processes (e.g., Zahn's approach has been driven by psychological considerations [START_REF] Zahn | Graph-theoretical methods for detecting and describing Gestalt clusters[END_REF]).

For our future work, we plan to add another quality index adapted to graph-based clustering, in a similar way of the indices proposed by [START_REF] Pal | Cluster validation using graph theoretic concepts[END_REF] which is an adaption of the Davies-Bouldin's index to neighborhood graphs. And finally we will make a more specific study of the outliers that GBC detects in a more selective way than DBSCAN.

Table III EXPERIMENTAL

 III RESULTS OBTAINED ON THE 8 DATA SETS.


Aurélien Hazan
email: aurelien.hazan@univ-paris1.fr

Michel Verleysen
email: verleysen@dice.ucl.ac.be

Marie Cottrell
email: cottrell@univ-paris1.fr

Jérôme Lacaille
email: jerome.lacaille@snecma.fr

Random model of vibrations for Foreign Object Damage detection in a civil aircraft engine

In this article we analyze several vibration time series measured on a real fan test rig before and after it is hit by a flying object. We show first evidence that a windowed autoregressive model may be used to detect the shock after it occurred. We compare these results with a second time series that measures the rotation period of the fan. Lastly we repeat this experiment using adaptive autoregressive estimation instead of the classical AR estimation algorithm. First results show that vibration data allow better prediction. Furthermore, the windowed AR estimation algorithm seems to perform better than the adaptive one.

Introduction

Foreign Object Damage (FOD) is a serious threat to civil aicraft engines, all the more given that it occurs and vanishes very quickly. One can detect FODs with the help of vibration sensors because they elicit a short and intense rise in vibration amplitude, as shown by Fig. 2 (top). However, for several reasons this burst may be hard to detect online, or may trigger false alarms when buried in real engine noise.

In this article we consider detecting FODs even after the shock, by monitoring vibrations sensors and the duration of each shaft revolution, and modeling them as random signals. This approach is commonplace in structural health monitoring [START_REF] Poulimenos | Parametric Time-Domain Methods for Non-Stationary Random Vibration Modelling and Analysis: A Critical Survey and Comparison[END_REF] but to the best of our knowledge was not applied so far to FOD detection. 1

The test system is not a complete aircraft engine, but a partial fan with a small number of blades used by Snecma for testing, as depicted by Fig. 1. The shaft angular position can be computed thanks to a keyphasor that emits a pulse at the end of each cycle, and whose raw output is noted tacho(n). Furthermore, the amplitude of vibration is measured with a accelerometer mounted on the test rig, and that outputs vib(n). At a specific time instant, a cannon fires a projectile at a controlled speed, aimed at a precise blade.

Figure 1: Test rig composed of a single fan with a small number of symmetric blades, and equipped with a tachometer and a vibration sensor.

The available data are provided by the Health Monitoring Department of Snecma. It includes vibrations and raw keyphasor for three different experiments, sampled at 25000 Hz, during several hundred seconds. Raw vibration data is kept almost unchanged, since it shows no trend, variance fluctuation nor slow periodic variation. Because our aim is to detect the shock after its occurrence, we remove the short time interval that corresponds to this event, represented by Fig. 2 (top). The beginning of the shock is found by thresholding the amplitude with a fixed arbitrary threshold, the length of the event is set to 2s. Meanwhile, raw keyphasor data is preprocessed with a local polynomial interpolation and simple root-finding so as to get precise zero-crossing time, and hence the duration of successive periods, given in seconds. The result is shown by Fig. 2 (bottom), where the shock is easily seen. In a similar way, we remove the time interval where the shock occurs. Moreover, we detrend and deseasonalize the data, to get rid of a periodic modulation of the signal. The obtained time series are modeled as autoregressive (AR) stationary processes, whose coefficients are estimated by two algorithms. First we use the classical Burg algorithm [START_REF] Brockwell | Introduction to Time Series and Forecasting[END_REF] with nonover-

Results

We first present results relative to windowed AR estimation and projection, in the case of vibrations. Each vibration time series is divided into 5 blocks of equal length. The first four blocks correspond to pre-shock data, the last block to post-shock data. The order of AR models is set to 10, after a first model selection step that is not discussed here. The dimension of the reduced space obtained by PCA is set to 2. For each of the three separate FOD experiments, Fig. 3 shows the resulting scatter plot. Pre-shock (resp. post-shock) dots are in blue (resp. red). In all three cases, post-shock data are represented by an isolated point. However, in cases 1 and 3 there is another distant point, whose presence is unexplained so far. Now we compute adaptive autoregressive coefficients [START_REF] Ljung | System Identification: Theory for the User[END_REF], and project them in 3-dimensional space. We discard the beginning of the obtained time series to get rid of convergence transient, and plot the three trajectories in Fig. 4. The green part of the curve stands for the post-shock time interval. In cases 1 and 3 we notice a sudden change in the direction of the tangent to the curve. This lack of frank discontinuity may result from the very nature of the adaptive algorithm, which uses past data to overcome the underdetermination of the estimation problem.

Lastly we have made the same experiment with the time series per(t), but must conclude that so far no interesting result was found. Our first intuition was that the aftershock could last longer with the period duration, but it turns out that the order of magnitude is roughly the same, as evidenced by Fig. 2.

Conclusions and future work

In this article we tackled the issue of detecting an FOD after the shock has occured. We have presented preliminary results obtained with a test rig that is different from a real aircraft engine. However some evidence is shown that elementary filtering and dimension reduction could be sufficient for such detection.

As more in-depth investigations are needed, we first plan to fine-tune the model selection and estimation of the autoregressive model coefficient. The process type is so far restrained to AR model, but may be extended to several classical processes (ARMA, ARIMA, . . . ) as well as others more likely to be found in vibration and health monitoring litterature (TARMA, . . . ). The order choice must be discussed in detail, as well as goodness-of-fit tests results. Similarly, the tradeoff between the window length and the number of available points will be debated on statistical grounds in order to balance the lack of data. Lastly, the physical justification for modelling normal vibrations with such AR process will be traced back in the relevant litterature 1 ; then the mechanical effect of a FOD on this model will be derived, so as to explain the effect on the estimated AR coefficients.

In the much more complex case of a real engine, all the previous steps will be taken successively so as to decide whether or not FOD can be detected after their occurrence. We remarked already that vibrations have a much more harmonic content in that case than with the test rig, which could impose further constraints on the AR model.

Figure 2 :

 2 Figure 2: Time series of (top) vibration sensors vib(t) and (bottom) period duration.

Figure 3 :

 3 Figure 3: Scatter plot of AR coefficients projected in a bidimensionnal reduced space, for the three successive vibration data blocks. The last point -in red-corresponds to the post-shock part of the time-series.

Figure 4 :

 4 Figure 4: Trajectory plots of adaptive AR coefficients projected in a 3-dimensional reduced space, for the three successive vibration data blocks. The last part -in green-corresponds to the post-shock part of the time-series.

it is easy to see why an approximate solution of the deterministic differential equation x ′′ = sin(ωt) can be found in an AR form, by simply replacing x ′′ by its finite difference counterpart x(t + 1) -2x(t) + x(t -1) and considering the sine as a random process.


Layane Abboud #
email: layane.abboud@supelec.fr

Andrea Cozza #
email: andrea.cozza@supelec.fr

Lionel Pichon
email: lionel.pichon@supelec.fr

Utilization of Matched Pulses to Improve Fault Detection in Wire Networks

Keywords: wire networks, fault detection, time reversal

A new concept to fault detection in wire networks, based on the properties of time reversal, is presented. The method, called the Matched Pulse approach (MP), propose to adapt the testing signal to the analyzed network, instead of using a predefined signal, as opposed to existing reflectometry methods. Through mathematical study and numerical simulations, we show the benefits of this technique. A physical interpretation is also presented to better understand the proposed approach.

I. INTRODUCTION

Location of faults on wires has become an area of international concern [START_REF] Furse | Down to the wire[END_REF]. Wired networks can be found everywhere: in transportation systems, industrial machinery, buildings, nuclear facilities, power distribution systems, etc., so that faulty electrical wiring has been considered among the most significant potential cause of failure and maintenance cost in these structures [START_REF]Review of federal programs for wire system safety[END_REF].

Away from the classical methods of wire testing such as visual inspection and impedance testing [START_REF] Furse | Down to the wire[END_REF], reflectometry methods are widely used today to help detecting and locating wire faults. The basic concept of these methods is presented in the next section. Generally, hard faults (open and short circuits) are detectable through reflectometry, while soft faults (damaged insulation, etc.) are more critical to detect, especially in complex wire networks.

In this paper, we propose an approach to fault detection in wire networks, called the Matched Pulse approach (MP). This proposed method relies on the idea of adapting the testing signal to the wire network to be analyzed, instead of injecting a predefined testing signal, as opposed to reflectometry methods. This method will be shown to be always beneficial when compared to standard Time Domain Reflectometry (TDR) [START_REF] Furse | A critical comparison of reflectometry methods for location of wiring faults[END_REF], as it permits to increase the detection probability of the fault, especially when it comes to identifying soft faults.

The paper is organized as follows: in the next section we give a brief reminder of the principle of reflectometry methods, then in the third section introduce the MP approach, based on the matched filter principle. The fourth and fifth sections present the investigation tools we used for comparing the MP and TDR approaches : by first defining an echo-energy gain, followed by a parametric study which shows the benefits of the new approach. This mathematical study is followed by a physical interpretation of the obtained results, based on the properties of time reversal [START_REF] Fink | Time reversed acoustics[END_REF].

II. REFLECTOMETRY METHODS

Reflectometry methods are widely used today for fault detection and location. These methods send testing signals down the wire network to be analyzed. They include TDR, which uses a fast rise time step or pulsed signal as the testing signal, Frequency Domain Reflectometry [START_REF] Furse | Frequency-domain reflectometery for on-board testing of aging aircraft wiring[END_REF] which uses multiple sinusoidal signals, Sequence TDR [START_REF] Smith | Analysis of spread spectrum time domain reflectometry for wire fault location[END_REF] which uses pseudo noise, and Spread Spectrum TDR which uses pseudo noise modulated onto a sinusoidal carrier signal. Several other methods also exist, such as Noise Domain Reflectometry [START_REF] Lo | Noise-domain reflectometry for locating wiring faults[END_REF], Mixed Signal Reflectometry [START_REF] Tsai | Mixed-signal reflectometer for location of faults on aging wiring[END_REF], etc. However, the principle of these methods remains the same: we first inject a predefined signal into the system to be examined, and then we analyze the reflected signal. The first step in the analyzing process is to decide whether or not a fault exists; this is the detection phase, where the reflected signal is analyzed to try to determine the echo from an eventual fault. Once detected, the second step is to determine the position and nature of this fault.

As we already mentioned, hard faults are easily detected using reflectometry methods, while soft faults are generally not. Additional complexity arises when using wire networks, because faults can be masked by the presence of several reflection sources (such as junctions, branches, mismatched loads, etc.); this produces multiple echoes, making it difficult to identify the one corresponding to the actual fault. Furthermore, those same reflections disperse the energy of the interesting echo and make the detection process even more difficult.

The problem of fault detection in wire networks can be compared to the problem of target detection encountered in radar, when the target to be detected is embedded in a cluttered medium [START_REF] Jiang | Antenna array detection in highly cluttered environment using time reversal method[END_REF]. In order to enhance the probability of detection of the target and reduce the false-alarm rate, the idea is to maximize the energy of the echo from the target to be detected, thus increasing its signal to noise ratio. Such approach has been studied in [START_REF] Fink | Time reversed acoustics[END_REF] [START_REF] Fink | Time reversal of ultrasonic fields -part 1: basic principles[END_REF], by exploiting the properties of timereversal.

H(f) H(f) I (f) TDR E (f) TDR I (f) TDR E (f) TDR E (f) MP I (f) MP * (a) (b) H(f)
In this paper, we introduce the MP approach, based on Time Reversal techniques, and compare it with the classical TDR. The analyzing method is the following: we will consider that the testing signal (i.e., the injected signal) is the input to our system, and the reflected signal is the output; the transfer function of this system in the absence of the fault is denoted as H 0 (f ) and in the presence of the fault is H F (f ). Generally, we can either analyze the reflected signal directly [START_REF] Furse | A critical comparison of reflectometry methods for location of wiring faults[END_REF], or take the difference of the two reflected signals (with and without the fault), especially when considering soft faults [START_REF] Griffiths | The invisible fray: a critical analysis of the use of reflectometry for fray location[END_REF], so that the echo from the fault would be more easily detected. this latter method is the one we will be using, assuming that our system is linear time invariant. Analyzing the difference of the two reflected signals is consequently equivalent to analyzing the output of an equivalent system whose transfer function H(f ) is defined as follows

H(f ) = H F (f ) -H 0 (f ) (1)
This equivalent system will be referred to as the difference system.

III. MATCHED PULSE

In this section we introduce the basic concept of the MP approach. Instead of injecting a predefined signal into the network to be analyzed, as for existing reflectometry methods, the idea is to adapt the injected signal to the system, so that the energy of the echo from the fault is maximized, thus increasing its detection probability. Such a tailor-made signal, or MP, can be synthesized considering the matched filter approach used in signal processing [START_REF] Papoulis | Signal Analysis[END_REF]. According to this theory, a matched filter is obtained by correlating a known signal, or template, with an unknown signal to detect the presence of the template in the unknown signal. This is equivalent to convolving the unknown signal with a time-reversed version of the template. The matched filter is the optimal linear filter for maximizing the signal to noise ratio (SNR) in the presence of additive white noise.

By examining Figure 1, the same idea is encountered. Let us remind that the testing signals are I T DR (f ) and I M P (f ), and H(f ) is the transfer function. In the TDR scheme, we simply have the predefined signal injected into the system, thus resulting in the echo E T DR (f ). As for the MP, if we consider for example that we are injecting a Dirac Pulse (in the time domain), then the TDR echo would be the transfer function of the system. By sending the time reversed version of this echo into the same system, we are effectively auto correlating it. This is the same as doing the matched filtering of the TDR echo to obtain the MP echo E M P (f ). Therefore we have a higher efficiency in the MP case, because the signal is adapted to the studied system. This fact will also be shown in section V.

IV. SIMULATION SETUP

In order to better illustrate the effects of introducing the MP approach, we chose to compare it with the TDR. The simulation method and analyzed configuration are presented in this section.

A. Analyzed Configuration

The configuration we chose to study is represented in Figure 2. It consists of a lossless transmission line with a 75 Ω characteristic impedance, as for coaxial cables. The line is terminated by an open circuit at its far end. Two parallel impedances of, respectively, 140 Ω and 200 Ω, stand for equivalent loads, such as for two further parallel branches. The fault to be detected is located between those two impedances. We chose to study such a configuration as it is difficult in practice to analyze it with the existing reflectometry methods, due to the number of echoes resulting from the network discontinuities. The fault to be detected is modeled as a parallel impedance. In this paper we will only examine such types of faults.

B. Simulation Method

We used the transmission line theory as presented in [START_REF] Paul | Analysis of Multiconductor Transmission Lines[END_REF]. Briefly, let us consider a uniform scalar transmission line, as illustrated in Figure 3; V S is the injected voltage, and Z S is the impedance of the source, considered equal to the characteristic impedance of the line Z c . The impedance of the load at the end of the line is denoted as Z L . The voltage and current at any point of the line which position is z can be calculated according to the following formulas:

V (z) = V + 0 exp(-jγz) + V - 0 exp(jγz) (2)
Z s V s Z L z = 0 V - 0 V + 0 Fig. 3.
A model illustrating the transmission line theory used in the paper: V + 0 and V - 0 are the modal variables to be determined in order to characterize the voltage and current propagation along the line.

I(z) = I + 0 exp(-jγz) -I - 0 exp(jγz) (3)
V + 0 and V - 0 are the modal voltages in the frequency domain, and I + 0 and I - 0 are the modal currents. γ is the propagation constant. We can write:

Z c = V ± 0 /I ± 0 .
In order to determine the voltage and current at any point of the line, we need to determine the values of two modal variables. The advantage of this method is therefore obvious: we have a full description of the evolution of the voltage and current along the transmission line just by obtaining the values of two modal variables, which can be determined by imposing the boundary conditions at the source and the load. The fact that this analytic model describes physical phenomenons made it interesting to use in our software implementation.

V. GAIN OF THE NORMALIZED ECHO ENERGIES OF THE MP AND TDR APPROACHES

In order to test the ability of the proposed method in delivering higher echo energies, we now introduce the tools we used to compare the MP and TDR approaches. This is followed by a parametric study and numerical simulations that prove the advantages of our new approach.

A. Definition

We will first define a gain G(Z F ) as the ratio of the normalized energy of the MP echo to the normalized energy of the TDR echo. We note that the normalization is done according to the energy of the testing signal in each case, that is the normalized energy of the TDR echo is the ratio of the energy of this echo to the energy of the gaussian pulse, and the normalized energy of the MP echo is the ratio of the energy of this echo to the energy of the TDR echo (which is the time reversed version of the testing signal in this case). This normalization makes possible the comparison of the two methods.

The normalized energy of the TDR echo is and the normalized energy of the MP echo is

E TDR (f ) = |E T DR | 2 df |I T DR | 2 df = |H(f )| 2 |I T DR (f )| 2 df |I T DR | 2 df (4)
E MP (f ) = |E M P | 2 df |I M P | 2 df = |H(f )| 4 |I T DR | 2 df |H(f )| 2 |I T DR | 2 df ( 5
)
The gain is the ratio of these normalized energies, therefore we can write

G(Z F ) = E MP (f ) E TDR (f ) (6)
Using the Cauchy-Schwartz inequality [START_REF] Papoulis | Signal Analysis[END_REF], one can easily prove that G(Z F ) is always greater than 1, and consequently the MP approach is always beneficial compared to TDR.

B. Parametric Study

We analyze the previously defined gain in terms of the nature and position of the fault. But first, it is interesting to be able to predict the variation of the gain according to the nature of the fault. Therefore, we chose to express this gain in terms of the reflection coefficient from the fault Γ F , and the scattering matrix S of a two ports system. This representation would allow us, for a single numerical characterization (corresponding to a given network configuration), to directly predict the system's response for any fault impedance Z F .

To better illustrate this idea let us examine Figure 4. We considered the setup of Figure 2. The system is described by the scattering matrix S, linking its inputs to its outputs: the inputs are the incident voltage waves from the two ports, and the outputs are the reflected waves. The impact of introducing a fault to the system is described by Γ F . The transfer function of the system may now be written as follows

H(f ) = S 12 S 21 (1 -Γ F ) (1 -S 22 )(1 -Γ F S 22 ) (7)
By introducing this expression into Equation ( 6), we can predict how the gain changes according to the impedance value of the fault, in the previously illustrated configuration (see Figure 2).

C. Numerical Results

We consider the numerical results obtained by changing the value and position of the fault in the system illustrated in Figure 2, and computing the corresponding value of G(Z F ). In our experiment, we consider that the impedance value of the fault is expressed as

Z F = R + 1 jωC ( 8
)
where R is the real part and -1/ωC the imaginary one. The values of this impedance are chosen in order to simulate as much as possible the real faults [START_REF] Griffiths | The invisible fray: a critical analysis of the use of reflectometry for fray location[END_REF]. The values of G(Z F ) in terms of R and C are shown in Figures 5a and5b. By examining the gain's variation in the two figures, we can see that for very small values of the capacitance, the value of the gain is about 3.5. Note also that in this case, the value of R does not have a big influence on G(Z F ), and consequently the softness or hardness of the fault does not affect the gain.

As the capacitance increases, the gain decreases around 1.2 for a big value of C. Furthermore, when closely examining those results, we can notice that it is the constant τ = RC which influences the most G(Z F ). When the value of 1/τ (i.e., the cutoff frequency of Z F ) starts to reduce the bandwidth of the injected signal, the gain starts to decrease. This fact is not similar to the one we will be examining in Figure 7, because in this case the injected energy is still the same, but the bandwidth of the echo is reduced due to the interactions with the loads. Here it is also important to note the effect of the dispersion introduced by the presence of the capacitance C: when its value increases, the gain decreases; In fact the state equation of reactive elements does not meet the requirement of having a time reversible system, unlike the equation of the wave propagation. Therefore the dispersion introduced by such an element would result in decreasing the effectiveness of the time reversal operation, which corresponds to decreasing the value of G(Z F ).

The next step is to observe the variation of G(Z F ) according to the variation of the fault's position. In our case we considered the previously studied configuration, and limited the variation of the fault's position between the two impedances ( which values are 140 Ω and 200 Ω ). The origin of the position is the 140 Ω impedance, and the fault is being moved toward the other impedance. We chose an arbitrary value of the fault, where R = 50 Ω and C = 1 nF. The numerical results are shown in Figure 6, where it is clear that the position of the fault influences the gain G(Z F ), keeping the advantage of the MP over the TDR. Note that the values of the gain change within 1.329 and 1.15. We notice that G(Z F ) decreases to reach its lowest values when the fault is close to the 140 Ω and 200 Ω impedances. In these cases, we lose the multiple reflections, which will eventually decrease the efficiency of the time reversal process.

Another interesting aspect is to observe the influence of the bandwidth of the injected signal on G(Z F ). Here, we remind that the choice of the sampling step directly affects the resolution capacity of the TDR.

According to Equation (6), it is obvious that the bandwidth of the testing signal limits the integration interval, and also the spectrum of the injected signal plays the role of a weighting function for the transfer function of the difference system. Given that the theoretical prediction of the influence of this variation is difficult, we computed the values of G(Z F ) in terms of the bandwidth of the injected signal, for the same

0.5 1 1.5 2 2.5 3 3.5 x 10 7 1.1 1.2 1.3 1.4 1.5 1.6
Bandwidth of the injected signal (Hz) Gain Fig. 7. The Gain G(Z F ) in terms of the bandwidth of the injected gaussian pulse. The system is the same as illustrated in Figure 2.

given value of Z F (R = 50 Ω and C = 1 nF). The results are shown in Figure 7. We can notice that while increasing the bandwidth of the injected signal, we are averaging the integrand, so after a certain bandwidth (about 35 MHz in our example of Figure 7) the gain's value gets steadier. Note that we are not able to predict the optimal value of the bandwidth because it is directly related to the transfer function H(f ).

VI. STUDY OF A PARTICULAR CASE

In this section we focus on a particular case where the fault is soft, in order to observe and compare the injected and received signals in the TDR and MP cases, keeping the same configuration as before (see Figure 2). If we observed the equivalent scheme of Figure 4, we will notice that a soft fault can be modeled by a small variation of the reflection coefficient compared to the case where the fault doesn't exist. In other words, a value of the fault close to the case where we have an open circuit. The value we chose is 600 Ω, and consequently |Γ F | = 0.059. Thus, we have a very small reflection on the fault.

A. Testing Signals and Echoes

Let us examine Figure 8a. The testing signal in the TDR method is a gaussian pulse. In the MP case, the testing signal is the time reversed version of the TDR echo. The two testing signals are normalized to the same energy. We can notice the shape of the testing signal in the MP case, containing the information from the analyzed system; the different echoes with positive and negative peaks finally contribute to maximize the energy of the echo from the fault to be detected. This will be observed in the next section. We also notice that the highest peak has its amplitude lower than the amplitude of the gaussian pulse, a fact that might be useful when examining electromagnetic interference problems.

As for the echoes obtained in both cases, they are shown in Figure 8b. It is clear that the highest peak of the MP echo is greater than the TDR one, thus increasing the detection probability of the fault. It is also observable that the ratio of the amplitude of the highest peak to the amplitude of the secondary peak is greater in the MP case (about 1.532 in the MP case compared to 1.379 in the TDR case); this is important because it reduces the ambiguity of detection in the case of multiple echoes. Further discussion is presented in the next section, when examining the propagation of the two echoes.

B. Physical Interpretation

To better understand physically the reason why the MP approach results in a higher efficiency, we first illustrate the propagation of the echoes in the TDR and MP cases, according to the diagrams of Figures 9a and9b. Those diagrams represent the propagation of the voltage along the transmission line, both in space (horizontal axis) and time (vertical axis). The testing signals are normalized to the same energy. Here we considered the difference system, so that is equivalent to subtracting the reflected signals obtained with the presence and absence of the fault. The obtained signals contain the information from the fault to be detected. The three vertical black lines indicate respectively the positions of the three impedances: the 140 Ω impedance, the fault, and the 200 Ω impedance.

In the TDR case, we observe that the echo from the fault to be detected results mainly from the reflection of the testing signal on the fault. We should note that the testing signal is not seen on the diagram at the injection time, because it is the same with and without the fault; the subtraction of the reflected signals thus eliminate it. In the MP case, there is an increasing number of echoes compared to TDR, due to the shape of the testing signal; those echoes contribute to maximize the energy from the fault to be detected, as we can observe in Figure 9b at the fault's position: we can notice that in the TDR case, in the first reflection on the fault, the energy is almost equally distributed to the left and to the right of the fault, while in the MP case the energy is much more concentrated to the left of the fault (i.e., to the source side where we receive the echo). The echoes obtained in Figure 8b can be verified by looking at the position z = 0 in the two diagrams, along the time axis. We can observe the origins of the highest peaks in both cases. Also note that the 140 Ω impedance reduces the energy of the echo, as observed in the two diagrams, while the 200 Ω impedance is high enough not to significantly affect the echo propagation.

The results we observed in the case of the MP are due to the properties of time reversal. The Time Reversal is a self adaptive technique that compensates for the distortions due to the propagation through inhomogeneous media, in other words it takes into account the information recorded from the medium thus resulting in a higher efficiency. When we examined the previously studied diagrams, we saw how we used the TDR echoes (normally a nuisance in fault detection) to obtain an echo with a higher energy in the MP case. The Time Reversal allowed to synchronize these echoes, and thus made it possible to maximize the energy of the MP echo. This was not observed in the TDR case, where the reflected signal only takes into account the time delay to the fault.

VII. CONCLUSION

This paper presented an approach to fault detection in wire networks, based on the properties of time reversal: the Matched Pulse approach. This method was compared to the classical Time Domain Reflectometry, through mathematical analysis where we defined and studied a gain in order to show the benefits of this approach. The MP was found to be always beneficial compared to TDR, as it permits to increase the detection probability of the fault to be detected. A physical interpretation based on the properties of time reversal allowed us next to better understand the proposed method.

Fig. 1 .

 1 Fig. 1. Comparison of the TDR and MP methods: (a) is the TDR block diagram, and (b) is the MP one.

Fig. 2 .

 2 Fig. 2. The equivalent wire network analyzed in this paper.

Fig. 4 .

 4 Fig. 4. The equivalent illustration of the analyzed network, used in the simulations.

Fig. 5 .

 5 Fig. 5. The values of G(Z F ) in terms of the impedance of the fault, for several values of R (a) and C (b). The studied system is illustrated in Figure 4.

Fig. 6 .

 6 Fig. 6. The values of G(Z F ) in terms of the fault's position. The studied system is represented in Figure 2. The fault is being moved from the 140 Ω impedance towards the 200 Ω impedance.

Fig. 8 .

 8 Fig. 8. The injected signals in the TDR and MP cases, normalized to the same energy (a); The echoes obtained in both cases (b). The system is the same as illustrated in Figure 2.

Fig. 9 .

 9 Fig. 9. The Space-Time diagram of the TDR echo showing its propagation in the system illustrated in Figure 2 (a); The Space -Time diagram of the MP echo (b).