Hideo Aochi

Julien Rey

John Douglas

NUMERICAL SIMULATION OF WAVE PROPOGATION IN THE GRENOBLE BASIN: BENCHMARK ESG2006

. As the resolution of the basin model is not enough good and the inelastic attenuation is not fully taken into account, it is still difficult to quantitatively compare signal durations. However it is visible that the basin resonates a long time after the passage of body waves.

Introduction

This article presents the ground motions computed for the three earthquake scenarios of the ESG2006 benchmark on ground motion simulation. First simulations from two small real earthquakes (M 2.8-2.9, named W1 and W2) are compared with the velocity timehistories observed at three different stations. The other two scenarios consist of hypothetical pure strike-slip earthquakes of M w 6.0, occurring close to the Grenoble sediment-filled valley, named S1 and S2, whose hypocenters are the same as W1 and W2 respectively. Ground velocities at 40 (9 on rock and 31 within the valley) stations are discussed (Figure 1).

The first section of this paper briefly introduces the simulation method used, both in terms of the kinematic rupture model adopted and the computation method used for wave propagation. The following section presents the simulated ground velocities for the four scenarios. In particular, comparisons are made between the simulated PGV and those predicted by recent empirical ground motion estimation equations.

Figure 1. 2D map (S-wave velocity on the ground surface) simulated in this study. Map scale is in meters.

Simulation method

For the simulations here the following simulation scheme was used, which is a similar scheme to that used by [START_REF] Douglas | The importance of crustal structure in explaining the observed uncertainties in ground motion estimation[END_REF] to model ground motions in three French regions. A standard equally spaced staggered finite difference method (FDM), which is fourth order in space and second order in time [START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF][START_REF] Olsen | Simulation of three-dimensional wave propagation in the Salt Lake Basin[END_REF], is used.

Our code was previously used for testing dynamic rupture propagation and the consequent wave propagation along finite segmented faults for the 1999 Izmit earthquake [START_REF] Aochi | The 1999 Izmit, Turkey, earthquake: Nonplanar fault structure, dynamic rupture process, and strong ground motion[END_REF]. Other aspects of this type of FDM are detailed in [START_REF] Graves | Simulating seismic wave propagation in 3D elastic media using staggered grid finite differences[END_REF]. Furthermore, our procedure is characterised by the implementation of a finite fault model [START_REF] Olsen | Hybrid Modeling of Curved Fault Radiation in a 3D Heterogeneous Medium[END_REF] and the PML absorbing condition [START_REF] Collino | Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media[END_REF].

The lower depth of the basin shape is given within a file at every 250m. In our simulation, we linearly interpolate it where necessary. We truncate the northern end of the given model by a few km to improve calculating efficiency. Our modelling volume is then 30km x 30 km x 10 km. The medium heterogeneity is averaged in the manner of [START_REF] Graves | Simulating seismic wave propagation in 3D elastic media using staggered grid finite differences[END_REF].

The qualify factor of the attenuation Q is taken as infinity in the bedrock and 50 within the basin. The station location is approximated at one of the closest FD grid points after converting the given WGS coordinates into geographical Lambert coordinates. As in Figure 1, our model globally conforms to the true situation, but in detail, for example, station 1 is located on the bedrock in the field but is on sediment in the simulation. Though station 8 (541.5 m depth) is located correctly in the bedrock medium, the lower depth of the basin is embedded at around 480-500 m depth in our simulations. The output at this station is strongly affected by this medium boundary. The border area of the strong contrast in the medium is not suitable for quantitative comparison because of the numerical treatment in the program as well as the inherent numerical instability of the FDM.

Due to the strong contrast of the basin sediment and its non-negligible dimension, it is necessary to calculate with a fine grid spacing. Preliminary simulations are carried out using a grid of 100 m (time step 0.005 sec), with which our regional scale simulations are often executed. However, our target earthquakes generally have a magnitude of more than 4 and the assumed slowest S-wave velocity is around 1000 -2000 m/s, that is, the usual target is a few Hz at rock sites. Finer simulations with a grid of 50 m (time step 0.003) are tested here with help of our Computational System and Technology division (STI). This size is required by a simple analogy of the need of 5 grid points for a target wavelength (1 Hz wave for Vs=300m/s and the thickness of sediment layer of a few hundred meters). This may allow observation up to 1 Hz in the basin (on the sediment) and 10 Hz at rock sites at best. Through a brief study, the waveforms calculated with different resolutions are similar at rock sites. The early phase of waveforms at sediment sites seems similar, while the significant differences appear in later phases (resonance of the basin). However, even with these finer simulations, it is difficult to physically discuss those later phases, because the given structure model and the grids may not be sufficient and because the method of wave attenuation is too simple.

Concerning the seismic sources used as input, we use a locally supported B-spline function of order 4 (degree 3) for describing a smoothed piece-like function instead of the differential error function. It is an advantage that the B-spline function has a local support (finite duration) so that the origin time is always zero in our simulations. When this function is used as a source time function of a point source, it is known that the frequency context is not the one expected (i.e. flat at low frequencies in acceleration) [START_REF] Douglas | The importance of crustal structure in explaining the observed uncertainties in ground motion estimation[END_REF]. However this is sufficient to discuss the wave form in the time domain such as PGV, if we choose a reasonable duration according to the magnitude. When this is used as a slip function in a finite fault model, the form of the function becomes less important than the other factors such as rupture velocity, directivity, dimension etc. One has to introduce a complex asymmetric function for the slip history for more realistic simulations.

In any case, in our simulation, we do not introduce any further factors generating complexity in the waveforms. Some further remarks should be added. Though a rise time of 0.03 s is required for the small earthquakes, but it is in fact too short to realise because it may lead to unnecessary numerical oscillations. We limit the rise time to 0.2 s (5 Hz) for the tests. The simulation results are comparable not with the amplitudes and frequency contents of the observations but with the arrival times and with relative comparisons among different stations. For the large earthquakes, the proposed values of a rise time of 1.116 s, a uniform final slip of 1.116 m and a rupture velocity of 2.8 km/s as well as the given fault location and dimension are followed except for the mechanism of S2 (which was originally stated to be right lateral faulting).

Results

This section presents the results obtained for the 4 scenarios (W1, W2, S1 and S2), which are compared with the observed ground motions at the given stations, within the simulations with different numerical procedures and with empirical equations.

Comparison with the observations (W1: M w 2.9)

Figure 2 compares the simulated ground velocities for the seismic source W1 at the stations 1, 6 and 8, which correspond to OGMU, OGFH and OGFB, respectively. This simulation is carried out with a grid of 100 m, which is not fine enough for wave propagation within the basin. In the simulation, station 8 is just located on bedrock (the output is sampled at a depth of 600m) like the corresponding station OGFB so both should be comparable. Characteristic waveforms are commonly observed, but the arrivals are systematically shifted (they are earlier in the simulation). This may be because the event time or location and/or the medium properties are not well constrained, unless the problem is due to numerical problems. Station 6 shows the surface ground motion at the same location. The waveforms are not similar because of the lack of resolution. We recall again that station 1 is located on sediments in the simulation unlike the corresponding station OGMU. The

Duration of event (W2: M w 2.8)

Figure 3 compares the ground motion at three different stations on rock between two simulations for the seismic source W2. One of the simulation uses a rise time of 0.2 s as in all the other simulations, while the other adopted 0.003 s as proposed. Briefly it may be estimated that the resolution is up to 5 Hz at most with respect to wave propagation in the bedrock (Vs of about 3000 m/s, grid size of 100m). As clearly observed, there are very strong spikes after the arrivals of body waves in the raw output (left column). They must be from inevitable numerical disturbances within the FDM because they do not disappear until the application of a filter with a cut off up to 5 Hz. In this case, the information on the high frequency generation (source duration shorter than this period) is lost and the obtained synthetic waveforms become identical. It leads to the conclusion that short period behaviors of the seismic source become unimportant when compare with the numerical resolution of the FDM.

Figure 3. Comparison of ground motions at stations on bedrock for the seismic source W1 using different event durations: 0.2 s in red, and 0.03 in green.

Large event (S1: M w 6.0)

Figure 4 shows the simulated ground velocity for six different stations for the seismic source S1 from two simulations using different grid sizes. As expected, both waveforms are similar for the rock sites (numerical resolution is sufficient). From the waveforms at these rock sites, it can be seen that the rupture process is very simple without any complexity regardless of its finite dimension. On the other hand, within the basin on the sediments, the early phases seem to be similar in both simulations, but the later phases differ especially for stations in the near field and in the horizontal direction due to a lack of resolution. It confirms that high resolution should be required for discussing the ground motions in the basin. In general, the simulations well reproduce the amplification and long shaking duration in the basin. The same tendency is observed in the simulations obtained for the source S2.

Comparison to empirical ground-motion estimates

Similarly to the approach adopted by [START_REF] Douglas | The importance of crustal structure in explaining the observed uncertainties in ground motion estimation[END_REF], the predicted ground motions from the M6.0 earthquakes (S1 and S2) are compared to empirical groundmotion models for the estimation of PGV. These five recent equations for the estimation of PGV from shallow crustal earthquakes have been selected: Campbell (1997), [START_REF] Sadigh | Updated relationships for horizontal peak ground velocity and peak ground displacements for shallow crustal earthquakes[END_REF], [START_REF] Chapman | On the use of elastic input energy for seismic hazard analysis[END_REF], [START_REF] Tromans | The attenuation of strong-motion peaks in Europe[END_REF] and [START_REF] Pankow | The SEA99 ground-motion predictive relations and new peak ground velocity relation[END_REF]. For the equation of [START_REF] Tromans | The attenuation of strong-motion peaks in Europe[END_REF], which uses surface-wave magnitude (M s ) rather than moment magnitude (M w ), the conversion formulae of Ambraseys & Free (1997) has been used, which in fact gives a M s of 6.0 for a M w of 6.0.

The different distance metrics needed by the different models have been calculated for the source and site geometries of the benchmark. The empirical models have been evaluated for the site classes corresponding to the site conditions at the receivers. For the models of Sadigh & Egan (1998), [START_REF] Chapman | On the use of elastic input energy for seismic hazard analysis[END_REF], [START_REF] Tromans | The attenuation of strong-motion peaks in Europe[END_REF] and [START_REF] Pankow | The SEA99 ground-motion predictive relations and new peak ground velocity relation[END_REF] this means either rock or soil or rock or soft soil. However, for the model of [START_REF] Campbell | Empirical near-source attenuation relationships for horizontal and vertical components of peak ground acceleration, peak ground velocity, and pseudo-absolute acceleration response spectra[END_REF], which seeks to model the effect of sediment depth, account was also taken for the depth of the sediments at receivers within the basin.

Residuals defined as the difference between the common logarithm of the simulated PGV minus the common logarithm of the PGV predicted by the empirical equation are calculated. Figure 5 shows the simulated peak ground velocity (PGV) for the two scenarios against rupture distance. It shows that, as expected, the ground motions decay with distance for rock sites. However, the amplification of ground motions due to the basin and the generation of surface waves mean that PGV values at sites within the basin do not show much decay with distance. It is remarkable that the PGV at station 4 (at about 11km) on rock for scenario 1 is large (42 cm/s). This station is located on the extension of the earthquake fault, so this large value may be attributable to strong directivity and a maximum within the radiation pattern for shear waves.

Figure 6 shows the computed residuals for both scenarios using the model of [START_REF] Campbell | Empirical near-source attenuation relationships for horizontal and vertical components of peak ground acceleration, peak ground velocity, and pseudo-absolute acceleration response spectra[END_REF] with respect to rupture distance. The pattern of residuals is similar for the other empirical models although the absolute value of the residuals changes. Table I lists the mean and standard deviation of the residuals for the five empirical models for both rock and soil stations. The figure and the table show that, in general, the simulated PGV values are reasonably well predicted, however, ground motions in the basin are, in general, underestimated by an average factor of two to three times. As a continuation of the work present here, it would be interesting to compare the simulated long-period motions within the basin with those predicted by the recent models of basin effects developed for use in California (e.g. [START_REF] Joyner | Strong motion from surface waves in deep sedimentary basins[END_REF][START_REF] Choi | Empirical model for basin effects accounts for basin depth and source location[END_REF]. These models use the depth to the bedrock and, for some models, the distance from the edge of the basin to the site.

Conclusions

Numerical simulations of wave propagation in the Grenoble basin are carried out using a FDM for a few earthquake scenarios. The obtained ground motions are compared with the observed ground motions at the given stations, within the simulations for different numerical procedures and with empirical equations. The equations derived by [START_REF] Campbell | Empirical near-source attenuation relationships for horizontal and vertical components of peak ground acceleration, peak ground velocity, and pseudo-absolute acceleration response spectra[END_REF] better match the simulations for the rock sites outside the basin than for the soil sites within the basin. It is remarkable that PGVs at some basin sites show changes of several times at the same distance from the earthquake fault. We also observe strong directivity effects from rupture propagation along the fault. The obtained scatter is compatible with other simulations which we have performed recently [START_REF] Douglas | The importance of crustal structure in explaining the observed uncertainties in ground motion estimation[END_REF][START_REF] Douglas | The importance of crustal structure in explaining the observed uncertainties in ground motion estimation[END_REF] so that this is considered to be inherent due to the complexity in the rupture process and the crust structure models.

Figure 2 .

 2 Figure 2. Simulated (red) and observed (green) ground velocities for W1.

Figure 4 .

 4 Figure 4. Simulated ground velocities with different grid sizes: 100m (red) and 50m (green) at six different stations. The left two columns are for stations on bedrock and the right two columns are for stations in the basin (on sediment).

Figure 5 .

 5 Figure 5. Simulated peak ground velocities for scenario 1 (blue) and scenario 2 (red) against rupture distance.

Figure 6 .

 6 Figure6. Residuals using model of[START_REF] Campbell | Empirical near-source attenuation relationships for horizontal and vertical components of peak ground acceleration, peak ground velocity, and pseudo-absolute acceleration response spectra[END_REF] for scenario 1 (blue) and scenario 2 (red) against rupture distance.

Table I .

 I Mean residuals for the five empirical models for rock and soil sites

	Model	Mean and standard	Mean and standard
		deviation of	deviation of
		residuals (rock)	residuals (soil)
	Campbell (1997)	0.36 (0.18)	0.27 (0.25)
	Sadigh & Egan (1998)	0.10 (0.19)	0.40 (0.25)
	Chapman (1999)	0.21 (0.19)	0.35 (0.25)
	Tromans & Bommer (2002)	0.23 (0.18)	0.43 (0.23)
	Pankow & Pechmann (2004)	0.21 (0.19)	0.45 (0.24)

Acknowledgements

We first appreciate the effort of the organizers at LGIT Grenoble for preparing this benchmark test. The basic simulations of this study have been carried out as a part of the European project Interreg SISMOVALP (Seismic Movements of Alpine Valleys, 2003-2006) and the French national project SISMODT (2004SISMODT ( -2006) ) funded by MEDD. A part of the comparison between the simulated ground motions with empirical equations is a contribution to the national project QSHA (Quantitative Seismic Hazard Analysis, 2006-2008) funded by ANR. The high performance simulations with fine grids were undertaken by the Computational System and Technology division (STI) of BRGM in the framework of the national project NUMASIS (2006NUMASIS ( -2008) ) funded by ANR. This work is also a general contribution to the seismic risk research project of BRGM (SEISMORISK).


L Campagnolo

S Roman

J Perchoux

S Lorthois

A new optical feedback interferometer for measuring red blood cell velocity distributions in individual capillaries: a feasibility study in microchannels

Keywords:

Introduction

The dynamics of blood flow in microvascular networks is of great importance in the exchange of nutrient and waste substances between blood and living tissues. In these vessels of diameter less than 100µm, the development of quantitative methods for measuring the velocity of red blood cells (RBCs) is still challenging. In this context, the Dual-Slit (DS) technique, a temporal correlation technique first introduced by Wayland and Johnson [START_REF] Wayland | Erythrocyte velocity measurement in microvessels by a two-slit photometric method[END_REF], is commonly used. Our group has recently shown that, provided that several operational conditions are fulfilled, this technique can provide a precise measurement of the transverse velocity profile of RBCs, i.e along the x axis in Fig. 1. The measured velocities are maximal velocities in the depth of the channel, i.e. y direction [START_REF] Lorthois | An optimized technique for red blood cells velocity measurement in microvessels[END_REF]. However, this technique requires expensive equipment (high speed camera) and the related data-treatment is time consuming. Moreover, the duration of the acquisition is long (typically 2 to 40s). For these reasons, on line measurements, especially in cases of transient regimes, are not possible. The aim of the present work is to determine whether Optical Feedback Interferometry, a new optical technique in the microfluidic domain which is based on the optical feedback effect in laser diodes, can be used to perform quantitative measurements of RBCs velocity in channels of size less than 100µm. Optical Feedback Interferometers (OFI) are indeed compact, low cost and simple sensors. They are known for providing much shorter response times than DS, currently lower than 10ms, while keeping the precision of traditional Laser Doppler Interferometers. This technique has been previously applied and validated to measure velocity profiles in 300µm diameter channels [START_REF] Nikolic | Flow profile measurement in micro-channels using changes in laser junction voltage due to Self-mixing effect[END_REF], but its feasibility in smaller channels is still to be demonstrated. For that purpose, the velocity profiles obtained by OFI and DS are compared, using PDMS microchannels and spheric monodisperse particle suspensions in set-up configurations where the DS has been previously validated [START_REF] Lorthois | An optimized technique for red blood cells velocity measurement in microvessels[END_REF].

Methods

Particle suspensions (0.1%(v/v)) are prepared from a concentrated stock solution of 4µm diameter latex microspheres (In Vitrogen). The flow is controlled in rectangular PDMS microchannels of size 20µmx100µm using a syringe pump (PHD 22/2000 Harvard Apparatus) and Hamilton Gastight glass syringes. First, the DS technique is applied. In this technique, the channel under study is transilluminated and two photo-sensors (slits) are positioned, separated by a known distance, L s , along the vessel axis (z direction in Fig. 1). In our case, the slits are regions of interest (1x3pixels²) on images of particle flows recorded using a high speed video camera (PCO Dimax). The time modulation of light, produced by the passage of the particles flowing through the channel, is deduced by performing the sum of grey levels in both slits at each time step. For various distances L s between the two slits, a cross-correlation velocity, V ds =L s /T d , is obtained, where T d is the time delay for which the cross-correlation between the two signals is maximum. For an optimal L s [START_REF] Lorthois | An optimized technique for red blood cells velocity measurement in microvessels[END_REF], the final velocity at a given point on the x axis represents the maximal velocity in the y direction. Finally, the slits are successively positioned in the x direction to obtain a profile of RBCs maximal velocities. *Corresponding author. Email: Lucie.Campagnolo@alumni.enseeiht.fr Then, the OFI technique is applied to the same microfluidic set-up. The optical feedback effect occurs when a portion of the light emitted from a laser is reflected from an external target and reenters the laser cavity [START_REF] Petermann | External optical feedback phenomena in semiconductor lasers[END_REF]. For flow measurement, light scattered from a single moving particle that is suspended in a fluid is shifted in frequency by the Doppler Effect. The light suffering from optical feedback mixes with the initial light, causing optical interferences inside the laser cavity. By amplifying and analyzing the laser junction voltage in the frequency domain, the velocity of the particle can be calculated [START_REF] Riva | Laser Doppler Measurements of Blood-Flow in Capillary Tubes and Retinal Arteries[END_REF]. The OFI dedicated experimental set-up was as follows. The laser was a commercial Firecomms Vertical-cavity surfaceemitting laser (VCSEL) lasing at 667 nm. A duallens system has been used to minimize the sensing volume, and therefore improve spatial resolution. The optical head was mounted on a computer controlled 3-axis motorized micrometer translation stage, and the focus point was stepped through the channel along the x axis to reconstruct the flow profile. The VCSEL terminal voltage was ACcoupled to a low noise amplifier module and then fed into an analog-to-digital converter for digital signal processing on computer. For practical purpose, optical system collects light backscattered from a number of different particles. Therefore, the frequency spectrum contains a distribution of frequencies which corresponds to a distribution of particle velocities. For this configuration, the maximum velocity in the sensing volume can be associated with the 3dB cutoff frequency of the Doppler spectrum [START_REF] Riva | Laser Doppler Measurements of Blood-Flow in Capillary Tubes and Retinal Arteries[END_REF]. Both experimental set-ups are represented on Fig. 1 As shown by Roman et al. [START_REF] Lorthois | An optimized technique for red blood cells velocity measurement in microvessels[END_REF], in the conditions of this experiment for the DS technique, the measured velocity at each point in the x direction corresponds to the maximal velocity in the depth of the channel (y direction). Similarly, due to the optical configuration, the maximal velocity in the depth of the channel is also measured by the OFI technique.

Results and Discussion

In order to validate the OFI technique, the laser spot was first focused in the middle of the channel, and the pumping rate was set from 0.1 to 1 µL/min. The cutoff frequency position varies linearly with the flow rate, showing good agreement with theory (data not shown). Both reconstructed profiles, as well as the newtonian theoretical profile [START_REF] Patzek | Shape factor and hydraulic conductance in noncircular capillaries I. One-phase creeping flow[END_REF], normalized by the maximal measured velocity, are plotted on Fig. 2. The profile shapes are similar for the two techniques showing that OFI is suitable for measuring velocity distributions in such channels. These preliminary results demonstrate the feasibility of OFI measurement in microchannels of size below 100µm. The next step will be to use both sensing techniques simultaneously in the same experiment in order to perform quantitative comparisons. The technique will then be applied on RBCs flow. This promising low-cost technique brings the possibility of real time applications, such as in-vivo blood flow measurements. Moreover, OFI technique may have the potential for simultaneous measurement of hematocrit (RBCs volume fraction).

Figure 1 :

 1 Figure 1: OFI (left) and DS (right) set-ups.

Figure 2 :

 2 Figure 2: Normalized profiles of maximal velocity measured in a rectangular PDMS channel (100µmx20µm), using OFI (blue curve) and DS (red curve).

Acknowledgments

This work has been partly supported by BQR 2011 INPT VLAM and GDR CNRS 2760. We also thank David Bourrier from "Plate forme technologique du LAAS-CNRS".


J Juyon

I Burciu

T Borr

S Thuries

E Tournier

A low spurious level fractional-N frequency divider based on a DDS-like phase accumulation operation

Keywords: Phase accumulator, fractional-N frequency divider, phase error, frequency division ratio, Direct digital synthesis

This paper presents a novel architecture of a low spurious level fractional-N frequency divider. It has already been shown in previous work that the use of a direct digital synthesizer (DDS) represents a promising solution to the mitigation of the fractional spurs that appear at the output of fractional frequency dividers, but with the drawback of causing a non-linear control of the division ratio. The DDS-based architecture proposed in this paper recovers the benefit of having a linear tuning of the division ratio, while still having similar performance in terms of fractional spurs.

I. INTRODUCTION

The Phase Lock Loop (PLL) structure implements a feedback loop technique that forces the phase of a generated signal to follow the phase of a reference signal. In the radiocommunication domain, the major advantage of implementing this technique relies on the improvement of the phase noise (PN) performance of a voltage controlled oscillator (VCO). Studies show that the PN of a VCO based frequency synthesizer is significantly improved when using PLL architecture in order to lock the phase of the VCO to that of the output of a low frequency reference [START_REF] Skworks | Basic of Dual Fractional-N Synthesis[END_REF].

Moreover, in order to be able to take profit of the low frequency reference PN at a high frequency, the PLL architecture acts as a frequency multiplier by including a frequency divider in the feedback loop (Figure 1). The phase of the signal present at the divider output is then compared with the one of the reference signal by the Phase Frequency Detector (PFD).

The polarity of the phase difference measured by the PFD is used to charge and discharge the filter with the Charge Pump (CP). This amount of charge is proportional to the phase difference measured by the PFD and will turn into voltage that controls the VCO. When the loop reaches the lock conditions, the phase of the VCO output signal will be the same than that of the reference signal. In order to obtain a frequency tuning of the VCO output, the commonly used solution relies on the use of a frequency divider with a programmable division ratio.

When speaking of Integer-N Frequency synthesizer, the programming of the frequency divider is done only by using integer values for the division ratio. As a consequence, the synthesizable frequencies will be multiples of the frequency of the reference signal. The main constraint of the integer-N synthesizers is the fact that the minimum frequency step size is equal to the frequency of the reference signal. In order to improve frequency resolution of PLL based frequency synthesizers, we have two types of solutions:

• lowering the frequency of the reference signal while using an integer-N frequency divider in the feedback loop.

• keeping constant the frequency of the reference signal while using a fractional-N frequency synthesizer in the feedback loop [START_REF] Swaminathan | A Wide-Bandwidth 2.4 GHz ISM Band Fractional-N PLL With Adaptive Phase Noise Cancellation[END_REF].

Each of these two solutions has its drawbacks. The decrease of the reference signal frequency implies an increase of the PLL lock time, while the use of fractional-N frequency synthesizers generates parasitic spurs (also called fractional spurs) in the output signal of the CP that will automatically degrade the PN of the global frequency synthesizer [START_REF] Razavi | RF Microelectronics[END_REF], [START_REF] Banerjee | PLL Performance: Simulation and design[END_REF].

In this paper we are assessing the conception and the design of frequency fractional dividers as well as the impact that these can have on the global performance of a PLL. In the second section, a presentation of different types of fractional frequency dividers will be made. The third section will be dedicated to a novel architecture we propose in order to realize a fractional frequency divider based on a DDS-like phase accumulator operation. In the fourth section, we present several simulation results that were done in order to validate the theoretical study concerning the novel architecture we propose. Finally, conclusions about this study are drawn and the follow-up to this work is provided.

II. FRACTIONAL-N FREQUENCY DIVIDERS

A. Digital Fractional-N Frequency Synthesis

One the most common used methods for achieving a digital fractional N-frequency synthesis is based on division ratio averaging over time. In order to obtain a fractional frequency division ratio, an output signal is generated from an input signal by alternating two or more integer division ratio over time. The global division ratio is therefore a noninteger which value depends on the ratios used and on their appearance duration as shown in:

N.F = N 1 P 1 + N 2 P 2 + ... + N k P k (1)
Where N.F denotes the average division ratio, and Ni and Pi are each integer division ratio and their appearance duration.

The upper part of the schematic presented in Figure 2 shows the architecture of a digital fractional-N frequency divider based PLL. In the lower part we present different architectures that are commonly used in order to implement the digital fractional-N frequencydivider.

For the first two implementation scenario, the fractional frequency divider architecture is based on a technique that involves the use of a multiple modulus integer frequency divider. The division ratio is dynamically switched over time through the use of a simple modulus controller. The implementation of this block varies depending on the number of integer division ratios that can be used. In fact, a simple digital phase accumulator (first order sigma delta modulator) can be used for a scenario where the division ratio can vary between only two integer values [START_REF] Skworks | Basic of Dual Fractional-N Synthesis[END_REF]. If the scenario implies the use of several integer division ratios, the most suited modulus controller architecture is considered to be that of a high order (superior to one) Sigma Delta Modulator (SDM) [START_REF] Kozak | Rigorous analysis of delta-sigma modulators for fractional-N PLL frequency synthesis[END_REF].

The third type of architecture that allows implementing the fractional-N frequency divider relies only on the use of an accumulator. The accumulator is clocked by the input signal while the output is generated by the register Most Significant Bit (MSB) [START_REF] Wu | A Quantization Error Minimization Method Using DDS-DAC for Wideband Fractional-N Frequency Synthesizer[END_REF]. The frequency control word (FCW) imposes the increment step. As shown in [START_REF] Wu | A Quantization Error Minimization Method Using DDS-DAC for Wideband Fractional-N Frequency Synthesizer[END_REF], if the value chosen for the FCW is not a power of two, the time between two successive overflows of the register is not constant over time. As a consequence, the output signal will have a "period" that varies between two values. Therefore, thisaccumulator based architecture is naturally switching over time the frequency division ratio between two integer values. The averaged frequency of the output signal fout is defined by:

f out = f in • 2 n FCW (2)
Where f in is the frequency of the input signal, n is the number of bits used for the register design and FCW is the control word.

B. Fractional Spurs

If the division ratio is constant over time for the ideal fractional frequency divider scenario, the averaging technique used to implement the fractional-N frequency divider imposes a variable frequency division ratio. The period of the output is therefore oscillating between two consecutive multiples of the input signal's period. As a consequence, a phase error appears when comparing the divider's output and the ideal divided signal. The schematics shown in Figure 3 show the evolution of this error over time: each of the arrows represent the end of a period of the divider's output. The phase error signal is represented by the evolution over time of the ε quantity. For an ideal fractional frequency divider this phase error is constantly equal to zero while it evolves for an average technique based divider.

From the PLL point of view, the phase error signal in the output signal of the fractional frequency divider is detected by the PFD. As a consequence, this phase error gives rise to a parasitic current signal at the output of the CP. This parasitic current is transformed into a voltage parasitic signal by the loop filter. If a periodicity characterizes the phase error waveform, the parasite voltage waveform will share this periodicity. Therefore, the spectrum of the VCO input signal will contain parasitic components that depend on the phase error periodicity.

Consequently, spurs appear in the PLL output as the VCO is modulated periodically. The magnitude of the periodic waveform is large compared to random jitter in the divider based on the average technique. As presented in [START_REF] Razavi | RF Microelectronics[END_REF], the fractional spurs generated by the phase error are yielding at typical values of 20 or 30 dBc. Various methods of mitigating the influence of this phase error have been proposed [START_REF] Banerjee | PLL Performance: Simulation and design[END_REF]. Three types of correction can be mentioned: current injection based fractional compensation, sigma delta noise shaping and use of Interpolation Zero Crossing method applied to a DDS output.

The current compensation of the phase error is realized by injecting a current pulse train to the integrating capacitor in the loop filter [START_REF] Huh | A CMOS dual-band fractional-n synthesizer with reference doubler and compensated charge pump[END_REF]. This current pulse train has the same width but opposite sign than the current signal generated by the phase error. This method is accurate but sensitive to the temperature variation.

The sigma delta noise shaping correction is based on a spread of the energy in the low frequency of the phase error signal towards higher frequency. This noise shaping is realized by destroying the periodicity of the phase error when using a high order sigma delta modulator in order to switch the division ratio between several integer values. As a consequence, the loop filter will be able to filter easily the frequency components present in the noise shaped phase error waveform as they are rejected at high frequency [START_REF] Kozak | Rigorous analysis of delta-sigma modulators for fractional-N PLL frequency synthesis[END_REF].

Finally, the third method of realizing a spur compensated fractional frequency division is based on the integration of a phase accumulator like in direct digital synthesizer architecture. This method is developed in [START_REF] Wu | A Quantization Error Minimization Method Using DDS-DAC for Wideband Fractional-N Frequency Synthesizer[END_REF]: the phase information generated by the accumulator is converted into an waveform by using a predefined table and a digital to analog converter (DAC). The analog signal is then filtered by a low pass filter. As a consequence, the amplitude of the phase error signal at the output of this component is greatly reduced.

III. DIGITAL N-FRACTIONAL FREQUENCY DIVIDER BASED ON A DIRECT SIGNAL GENERATION

As mentioned in the previous section, one promising solution of realizing a high frequency resolution PLL relies in the use of a DDS based fractional frequency divider. Nevertheless, if this solution offers a better performancecomplexity tradeoff compared with PLL using sigma delta compensated dual modulus divider, its major drawback relies in a non-linear evolution of the PLL output frequency depending on the evolution of the FCW injected at the input of the DDS. As mentioned in [START_REF] Wu | A Quantization Error Minimization Method Using DDS-DAC for Wideband Fractional-N Frequency Synthesizer[END_REF], when implementing a phase accumulator fractional frequency divider, the division ratio (DR) depends on the number N of the accumulator's bits and on the FCW as follows:

DR = f clk f out = 2 n FCW (3)
Where f out is the frequency of the output signal and f clk is the frequency of the signal that is driving the accumulator. The control of the division ratio is realized by changing the

f PLLout = f ref • DR = f ref • 2 n FCW (4)
Knowing that the reference frequency and the number of bits of the accumulator are constant, it becomes obvious that the evolution of the output frequency of the PLL is not linear depending on the evolution of the FCW.

In order to tackle this linearity evolution issue, this paper aims at proposing and at assessing a novel phase accumulator based fractional frequency divider architecture. The authors consider that this architecture is meant to replace the single phase accumulator of the DDS used to realize a frequency divider similar to that presented in [START_REF] Wu | A Quantization Error Minimization Method Using DDS-DAC for Wideband Fractional-N Frequency Synthesizer[END_REF].

Figure 4 shows the phase accumulator based architecture that is proposed in order to obtain a linear evolution of the DR.

The main innovation of this proposed architecture is the fact that the FCW that controls the phase accumulator is no more constant over time but is switching between two FCW. The control of this switching operation is realized by implementing a loop through a counter, a comparator and a switch.

As shown in the schematic model, the clock input of the counter is driven by the MSB of the phase accumulator. This corresponds to a counting of the phase accumulator overflow.

The output of the counter is tied to the input of the dual comparator. Depending on the result of the comparisons between the value forwarded by the counter and the A and B values, the two outputs of the comparator will be enabled. As long as the counter output exceeds the A or B value, the comparator will enable the output corresponding to A or B, respectively. While the output corresponding to the comparison with B is used to drive the reset input of the counter, the output that corresponds to the A value comparison is used to drive the switch between two FCW.

Based on the number of the phase accumulator's overflows, the FCW that controls the step increment of this block is changed over time. Moreover, the global system output is generated using the comparison with the B value through the use of a flip-flop.

When analyzing the output signal of the proposed architecture, the evaluation of its period's mean value has to take into account the fact that the FCW injected into the phase accumulator is switched twice during this period. As a consequence, when evaluating the time between two consecutive overflows of the phase accumulator, two scenarios appear depending on the value of the FCW that is injected.

As long as the comparison between the number of overflows doesn't exceed the A value, the mean value of the time between two consecutive overflows of the phase accumulator can be evaluated as DR 1 periods of the input signal. Based on (3), the DR 1 is defined by:

DR 1 = 2 n FCW 1 (5)
In a similar way, as long as the counter outputs a number between the A and B value, the mean value of the time between two consecutive overflows of the phase accumulator can be evaluated as DR 2 periods of the input signal, where:

DR 2 = 2 n FCW 2 (6)
Therefore, based on ( 5) and ( 6), the mean global frequency division ratio can be defined by:

DR GLOBAL = 2 A • 2 n FCW 1 + (B -A) 2 n FCW 2 (7)
For the proposed fractional frequency divider the tuning of the division ratio is controlled through the A input signal; the B, FCW 1 and FCW 2 input signals are constant over time. It can be observed that the DR evolution becomes linear, at the opposite of the simple phase accumulator architecture.

Therefore, the resolution of the division ratio is defined by the relation shown in (8). This is realized by incrementing A and by developing the difference between the two DR thus obtained:

∆DR GLOBAL = 2 n+1 • FCW 2 -FCW 1 FCW 1 • FCW 2 (8)
IV. RESULTS

In order to validate the results of these theoretical studies, simulations were done using realistic models of a "classical" phase accumulator, and the proposed phase accumulation dedicated architecture. The software used in order to realize these simulations is the Cadence system simulation and analysis [8]. Moreover, the implementation of the two compared structures is modeled by using the IBM BiCMOS 0.13 µm technology.

The first series of simulations aims at validating the linear division ratio evolution regarding the control value. The different parameters characterizing the design of the electronic component and the test set-up are shown in Table

I.

The Figure 5 shows that the evolution of the output frequency is linear and that the frequency resolution is equal to 2 kHz for the considered scenario. Concerning the proposed phase accumulation dedicated structure, the parameters used during the modeling stage are presented in Table II. A second fractional frequency division structure with a single 20 bits accumulator was chosen. This accumulator size of 20 bits is chosen in such a manner that the division ratio performance is similar to those obtained when using the proposed architecture. The value of the FCW injected in this simple phase accumulator was chosen to be equal to 3753. For both compared structures, the frequency of the input signal is equal to 10 GHz. The output signal's frequency is of 35.7891 MHz and 35.7914 MHz for the implementation scenario of the proposed architecture and the single 20 bits accumulator respectively.

The phase error signals generated when using each of the two phase accumulation dedicated architectures are cosimulated by Cadence/Matlab. The time difference is evaluated between the rising edges of the output signals of each of the two fractional frequency dividers and reference signals having constant frequencies of 35.7891 MHz and 35.7914 MHz respectively. The phase error is analyzed by using a probabilistic approach and by a comparative study of the frequency spectrums. The probabilistic approach is done by evaluating the Cumulative Distribution Function (CDF) realized over a large population of phase error values. The CDF function can be therefore defined by: F (x) = CDF (err) = P robability(err < x)

where err represents the phase error variable.

The graphics of Figure 6 show the CDF of phase error populations issued of implementation scenarios of the two analyzed fractional frequency dividers.

It can be seen that the time evaluated values of the phase error are linearly distributed between the same levels for the two implementation scenarios.

Moreover, the frequency spectrums of the signals corresponding to the phase error generated by the use of each of the two frequency dividers are presented by the charts of Figure 7. It can be seen that the power levels of the spurs generated by the phase error are similar for the two scenarios.

V. CONCLUSION

This paper presents a novel DDS-based fractional frequency divider with a phase accumulation dedicated structure. Unlike a single accumulator architecture that is commonly implemented in the DDS structure in order to realize the phase accumulation operation, the integration of this block allows a linear control of the division ratio when using the DDS as a fractional frequency divider. It has been shown that the performance of the proposed architecture is similar to that of a single accumulator structure but with a significantly reduced complexity. It can be concluded that the proposed architecture offers a much better performancecomplexity-power consumption trade-off when compared to a single accumulator structure.

Figure 1 .

 1 Figure 1. Phase-locked loop

Figure 2 .

 2 Figure 2.Fractional-N frequency Synthetiser and the different existing fractional frequency divider architectures.

Figure 3 .

 3 Figure 3. Temporal evolution of the phase error for three different scenarios depending on the implemented fractional frequency divider architecture.

Figure 4 .

 4 Figure 4. Phase accumulator based architecture

Figure 6 .

 6 Figure 6. CDF of the phase error corresponding to a use of the simple 20 bits phase accumulator (up) and to the proposed phase accumulation dedicated structure (down).

Figure 7 .

 7 Figure 7. Frequency spectrums of the phase error corresponding to a use of the simple 20 bits phase accumulator (up) and to the proposed phase accumulation dedicated structure (down).

  Evolution of the frequency at the output of the proposed strucutre depending on the evolution of the control value A. The freqeuncy of the input signal is equal to 10 GHz.

	Table I
	DIVISION RATIO EVOLUTION TEST PARAMETERS USED TO MODEL THE
	PROPOSED PHASE ACCUMULATION DEDICATED STRUCTURE
	Number of bits of the accumulator	9
	FCW 1	50
	FCW 2	51
	A	64
	B	1 to 64
	Frequency of the input signal	10 GHz
	Figure 5.


S Mercier
email: sabine.mercier@univ-rouen.fr

J J Daudin
email: daudin@inapg.inra.fr

Exact distribution for the local score of one i.i.d. random sequence

Keywords: P-value, sequence analysis, local score, Markov chain

Let X 1 :::X n be a sequence of i.i.d. positive or negative integer valued random variables and H n = max 0 i j n (X i + ::: + X j ) the local score of the sequence. The exact distribution of H n is obtained using a simple Markovchain. This result is applied to the scoring of DNA and protein sequences in molecular biology.

INTRODUCTION

The assessment of the statistical signi cance of scores of DNA and protein sequence is an important stage in the work of molecular biologists. Let A 1 :::A n be a nucleic or protein sequence. In order to identify interesting patterns, appropriate scoring values can be assigned to each residue. Scoring assignments for nucleotides or amino acids may arise from a variety of considerations likebiochemical categorization, physical properties, or association with secondary structures, (see [START_REF] Kyte | A simple method for displaying the hydrophatic character of a protein[END_REF][START_REF] Karlin | Methods for assessing the statistical signi cance of molecular sequence features by using general scoring schemes[END_REF] for examples). The local score of the sequence A 1 :::A n according to a scoring scheme is de ned as follow.

H n = max 1 i j n j X k=i (A k ) ! :
The local score is a very useful tool for biological sequence analysis in order to identify unusual sequence pattern or similaritythatmay re ect biological significance. It is desirable to know whether interesting patterns can arise by chance.

We are therefore interested in the distribution of H n under the null hypothesis of only random variation, so that wemay judge the statistical signi cance of a local score of a real biological sequence. For surveys of this subject, see [START_REF] Waterman | Introduction to Computational Biology[END_REF][START_REF] Durbin | Biological Sequence Analysis[END_REF]. Let X = (X i ) (i=1 ::: n) be a sequence of independent and identically distributed random variables of positive or negative integers. Let S k = X 1 + ::: + X k be the partial sums and S 0 =0. The local score of X is de ned as follows H n = max 0 i j n (S j ; S i )= max 0 i j n (X i + ::: + X j ) :

(1)

The local score has been already studied a lot, and asymptotic approximations have been obtained but its exact distribution is unknown untill now. Following [START_REF] Iglehart | Extremes values in the GI/G/1 queues[END_REF][START_REF] Karlin | Methods for assessing the statistical signi cance of molecular sequence features by using general scoring schemes[END_REF][START_REF] Karlin | Limit distributions of maximal segmental score among Markov-dependent partial sums[END_REF] proved the next limit when E X] < 0andfor X non lattice

lim n!+1 P H n log n + x] = exp ; ;K e ; x (2)
where and K depend only on the probability distribution of the X i . For X lattice, there is no limit for P H n log n + x] and ( 2) is replaced by exp (;K ; e ; x ) lim inf n!+1 P H n log n + x]

(3)

lim sup n!+1 P H n log n + x] exp (;K + e x )
where K ; and K + depend on the distribution of the X i too.

This result applied to sequence alignment, (see Dembo et al. 1994(b)), is implemented in BLAST, (see [START_REF] Altschul | Basic Local Alignment Search Tool[END_REF]), for gapless alignmentbetween two sequences.

In this paper we give the exact distribution of the local score H n for X lattice independently of the sign of E X]. This article is organized as follows. The result and its proof is given in Section 2. Section 3 is devoted to a comparison between our exact result and the approximation used by [START_REF] Karlin | Methods for assessing the statistical signi cance of molecular sequence features by using general scoring schemes[END_REF] on the numerical examples given by these authors. Section 4 deals with the ungapped alignment problem. A conclusion is given in section 5.

2. EXACT DISTRIBUTION OF H n Let T k be the time of the k th successive minimum of the process fS k g T 0 = 0 and T k+1 =inffi>T k : S i ; S T k < 0g and m(j) = supfk : T k jg, the number of successive minimum over the time from 1 to j. The T k may be in nite for k 1 if E X] > 0. These stopping time have already been studied by [START_REF] Karlin | Limit distributions of maximal segmental score among Markov-dependent partial sums[END_REF] to obtain their approximation of the distribution of the local score. See also Karlin and Taylor 1981 (Chapter 17). Let U be the process de ned by : U j = S j ; S T m(j) for j 1 and U 0 =0 Lemma 1 U possesses the following properties:

1. 8k =1 ::: m(n) U T k =0, 2. 8j > 0 U j 0, 3. 8j > 0 U j =(U j;1 + X j ) + . Proof 1
. By de nition of U j and m(j): 2. If U j < 0 we should have S j < S T m(j) which contradicts the de nition of T m(j) .

3. U j = S j ; S T m(j) = S j;1 + X j ; S T m(j) . There are two possible cases:

T m(j) = j or T m(j) <j.

T m(j) = j: On the one hand U j = 0 (see the rst property), and on the other hand U j;1 + X j = S j;1 ; S T m(j;1) + X j = S j ; S T m(j;1) = S T m(j) ; S T m(j;1) :

The last term is negativeby the de nition of T i . Thus (U j;1 +X j ) + =0 and we get U j =(U j;1 + X j ) + . T m(j) <j. Then T m(j;1) = T m(j) and S T m(j;1) = S T m(j) and we have U j = S j ; S T m(j) = S j;1 + X j ; S T m(j;1) = U j;1 + X j :

As U j 0 (property 2), U j;1 + X j 0. Therefore (U j;1 + X j ) + = U j;1 + X j and U j =(U j;1 + X j ) + .

The process U and the local score are linked by the following lemma Lemma 2

H n = max 0 k n U k : (4) Proof
With j held xed and by de nition of T m(j) wehave S j ; S i S j ; S T m(j) for all i j , with equality for i = T m(j) . Then H n = max 0 i j n (S j ; S i )= max 0 j n max i j (S j ; S i )] = max 0 j n (S j ; S T m(j) )= max 0 j n U j :

(5)

Let a be a positiveinteger and a the dual variable of U j : a = inffj 1 U j ag : Using (4), we get f a >n g = fmax k n U k <a g = fH n <a g : So the distribution of the local sore H n is given by P H n <a ]=P a >n ] : (6) Let U j = U j for j < a and U j = a for j a . The Lemma 1 implies that the process U is a Markov chain whose states are f0 1 : : : a g. Let =( h`) 0 h ` a be the probability matrix transition of U : h `= P U n = `=U n;1 = h] for 0 h ` a : By construction, a is an absorbing state, thus we have a a = 1 and a `=0for`=0 ::: a ; 1 :

(7) Moreover, h a = P U n = a=U n;1 = h] = P U n a=U n;1 = h] = P (U n;1 + X n ) + a=U n;1 = h] = P U n;1 + X n a=U n;1 = h] = P X n a ; h=U n;1 = h] :
As the X k are iid and U n;1 only depends on the (X k k n ; 1), we obtain h a = P X n a ; h]=P X 1 a ; h] for h =0 ::: a ; 1 : (8) For the state 0, we have h 0 = P X 1 ;h] pour h =0 ::: a ; 1 (9) and for the general transition from h, h = 1 : : : a ; 1 to `, `= 1 ::: a ; 1, we obtain h `= P X n = `; h]=P X 1 = `; h] : (10) Let f (k)=P X i k] and p(k)=P X i = k] for 0 h ` a. According to ( 7), ( 8), ( 9) and ( 10 The distribution of U n , denoted by P n = (P U n = 0] : : : P U n = a]) is given by an elementary result on Markov chains, (see [START_REF] Freedman | Markov chains[END_REF].

P n = P 0 n with P 0 =(1 0 0 ::: 0) : (11) Then P a >n ]=1; P a n]=1; P U n = a] : (12) Taking into account (11), ( 12) and ( 6), we obtain the following theorem Theorem 1 The c.d.f. of the local score is given by P H n <a ]=1; P n (a)=1; P 0 n P 0 a with P a =(0 0 0 ::: 1) : (13) ii) and c) iii)) where a is the highest among all examples. Generally the relative error is positive, a fact which is in agreement with KA] remark about the use of the constant K + which leads to conservative P-values. However the error is negative in two cases which proves that the KA] approximation is not always conservative. This formula is good for the tail of the distribution. One can use two di erent ways in order to judge of its accuracy.

First, when the P-value obtained is small enough (< 10 ;3 ), the asymptotic bounds (see ( 3)) can be use with con dence.

The second way stands on x, where x = a ; log(n)= . The ratio H n = log(n)

converges almost surely to 1= where depends on the distribution of X (see Dembo et al. 1994(a), [START_REF] Arratia | Critical phenomena in sequence matching[END_REF], and x corresponds to the deviation between the calculated local score a and the almost surely \limit". Larger is x, better are the conditions. For a safe use of (2), empirical tests lead to consider that x>3i sasu cient condition if n is enough high and E X ] << 0.

The two negative relative errors in the Table 1, c)(i) and e), are linked to the too large P-value : 3 56:10 ;2 for c)(i) and 9:5 ;1 for e). We are not in the tail of the distribution and the upper bound given byKarlinet al. can be less than the exact P-value. The deviation x equals to 2,39 for c)(i) and -14,88 for e). However note that the case d)(ii) is an example of high P-value but (3) is correct.

These numerical examples are only illustrations about the quality of the KA] approximation. A more complete numerical study must be made in order to give more general results and to produce advices about its safe use specially with short sequences.

Table 1 4. UNGAPPED ALIGNMENT PROBLEM

In order to compare two biological sequences A and B , scores can be assigned to couple of amino acids or nucleotides re ecting the similaritybetween the two components : (A i B j ). We can de ne the local score of two sequences as

H n = max 1 i j n j X k=i (A i B j ) ! :
The statistical signi cance of this problem is already solve by the case of one sequence (see previous sections). But this approach of evolution between two sequences implies that they have the same length, and does not consider shifted couple of segments. It does not consider as well the fact that components can be inserted or deleted, what we call indel, or gap. For the ungapped alignment problem, but considering shifts, see R Mott and R. Tribe 1999, we can use the exact probabilities in the following way. Let us consider two sequences of respectivesizem and n with m n. We can align the rst sequence with m+n;1 sequences obtained by slipping the second sequence (see Figure 1). The local score is then de ned as Figure 1 H n m = max 0 ` min(n m);1 1 i n;1 j m;`X k=0

(A i+k B j+k )

As most authors we do not take into account the dependence between them. Their weak dependence is ignored. Therefore the probability that the local score does not exceed a is approximately the product of the corresponding probabilities for each of the m + n ; 1 couple of sequences.

The length of each couple is equal to the common length, between 1 and n. There are 2 couples of length i for i =1 ::: n; 1 (see the rst case and the third cases of Figure 1) and m;n+1 couples of length n (see the middle case of the gure). Therefore the probability of exceeding a is approximatively equal to P(H n m a)=1; Y i=1 ::: n;1 P(H i <a) 2 P(H n <a) m;n+1

Note that we are able to take into account the size of each couple of sequences, which is not possible when an asymptotic method is used.

CONCLUSIONS

Our exact method is very easy to implement. A MATLAB computer program (available upon request) contains only thirty lines of code. It is not computationally intensive if one possesses a good subroutine for computing n th -powers of a matrix. For example the computation of the P-value a 100,000-long sequence with a = 100 is obtained in less than one second using a Pentium-II PC. Note that the KA] approximation also needs a computer program in order to obtain the constant K . The imprecision due to computation errors can be controlled using the fact that the sum of the probabilityofthe states of the Markov chain must be equal to one. The rst wrong nger in the above example is the 12th and we think that the computation errors are small for P-values greater than 10 ;10 . In our mind it would be useful to implement the computation of the exact P-values in any scoring sequence software. It would give exact P-values even if n is small and if E X] is near from zero, which is not the case for the asymptotic expression (3). Some work remains to deal with the alignment scoring problem with gap and with the case of Markovian dependent sites.

Figure 1 :

 1 Figure1: There is a total of m + n ; 1 di erent positions between the two sequences that can be regrouped in three main cases.

  3. COMPARISON BETWEEN THE EXACT AND THE ASYMPTOTIC RESULTS ON EXAMPLESKarlin and Altschul 1990 ( KA]) have presented representative examples of maximal scoring segments of protein sequences with scores based on charge, hydropathy, cysteine clusters and amino acid similarity. We use the same examples and give the exact result in comparison with the approximation given by KA]. The following table (Table 1) summarize the results. The indexing of the examples is the one used in KA]. Even if the relative error may be large in some case, one can see that the KA] approximation is correct and does not lead to dramatic errors of decision. However the relative error is high in two cases (Examples a)


A dramatic expansion of research in the area of electrochemical energy storage (EES) during the past decade has been driven by the demand for EES in handheld electronic devices, transportation, and storage of renewable energy for the power grid (1)(2)(3). However, the outstanding properties reported for new electrode materials may not necessarily be applicable to performance of electrochemical capacitors (ECs). These devices, also called supercapacitors or ultracapacitors ( 4), store charge with ions from solution at charged porous electrodes. Unlike batteries, which store large amounts of energy but deliver it slowly, ECs can deliver energy faster (develop high power), but only for a short time. However, recent work has claimed energy densities for ECs approaching ( 5) or even exceeding that of batteries. We show that even when some metrics seem to support these claims, actual device performance may be rather mediocre. We will focus here

True Performance Metrics in Electrochemical Energy Storage

Y. Gogotsi 1 and P. Simon 2 Exceptional performance claims for electrodes used in batteries and electrochemical capacitors often fail to hold up when all device components are included.

on ECs, but these considerations also apply to lithium (Li)-ion batteries.

Typically, the performance of both batteries and ECs is presented by using Ragone plots (see the fi gure) that show the relation between energy density (how far an electric car can go on a single charge) and power density (how fast the car can go). A commercial EC can harvest or release more energy than a typical Li-ion battery can deliver on time frames of less than 10 s, and it can be used for an almost unlimited number of charge and discharge cycles ( 4). A near-term application will be storing energy for car starter motors to allow engine shut-offs when stopped ( 6) and harvesting braking energy.

Increasing the energy density of ECs usually comes at the cost of losses in cyclability ( 5) or power, which are the most important properties of ECs and without which they become mediocre batteries. A major effort has been directed toward increasing the energy density of ECs by either increasing the capacitance of the material, C, or the operation voltage window, V, or both, since the energy stored is proportional to CV 2 . Some recent publications on graphene and nanotube-based materials have used Ragone plots to argue that supercapacitors can achieve the energy density of batteries. Those claims are summarized in the gray area in the upper right corner of panel A in the fi gure.

Reporting the energy and power densities per weight of active material alone on a Ragone plot (panel A) may not give a realistic picture of the performance that the assembled device could reach because the weight of the other device components also needs to be taken into account. ECs are similar to Li-ion batteries in that they contain current collectors, electrolyte, separator, binder, connectors, and packaging, in addition to carbonbased electrodes. Because the carbon weight accounts for about 30% of the total mass of the packaged commercial EC, a factor of 3 to 4 is frequently used to extrapolate the energy or power of the device from the performance of the material. Thus, the energy density of 20 Wh/kg of carbon will translate to about 5 Wh/ kg of packaged cell.

However, this extrapolation is only valid for electrodes with thicknesses and densities similar to those of commercial electrodes (100 to 200 µm or about 10 mg/cm 2 of carbon fi lm). An electrode of the same carbon material that is 10 times thinner or lighter will reduce energy density by three-to fourfold (from 5 down to 1.5 Wh/kg based on cell weight, see panel A), with only a slight increase in power density. Our ability to predict the performance of a 200-µm-thick electrode by testing a 1-µm fi lm ( 7) or a small amount of material in a cavity microelectrode ( 8) is still very poor. Experimental data show that there may be an additional drop in the capacitance by a factor of 2 to 3 when the thickness of the nanoporous carbon electrode increases ( 7).

Much of this uncertainty stems from reporting gravimetric, rather than volumetric, energy and power densities of materials and devices. Many nanomaterials, such as nanotubes or graphene, have a low packing density (<0.5 g/cm 3 ), which leads to empty space in the electrode that will be fl ooded by electrolyte, thereby increasing the weight of the device without adding capacitance. An extreme case would be the use of a carbon aerogel with 90% porosity. The volumetric energy of such an electrode will be 20% that of a carbon electrode with just 50% porosity.

If we consider a low-density graphene electrode (0.3 g/cm 3 ) with an extremely high gravimetric energy density of 85 Wh/kg (gray area in panel A of the fi gure), its volumetric density will be 25.5 Wh/liter for the electrode and ~5 Wh/liter for the device (panel B), which is a typical value for commercial ECs with activated carbon. If a 2-µm fi lm of the same material is used in the device, a much greater drop occurs, which is why "paper batteries" or thinfi lm ECs cannot be used for storing large amounts of energy.

The gravimetric energy density is almost irrelevant compared to areal or volumetric energy for microdevices and thin-fi lm ECs, because the weight of the active material used in a micrometer-thin fi lm on a chip or a nanotube coating on a smart fabric is negligible. These systems may show a very high gravimetric power density and discharge rates, but those characteristics will not scale up linearly with the thickness of the electrode ( 7), i.e., the devices cannot be scaled up to power an electric car. Ragone plots are only one measure of a device; they do not show other important properties, such as the device's cycle lifetime, energy effi ciency, self-discharge, temperature range of operation, or cost. They may also provide misleading information for fl ow and semisolid batteries ( 3,9), where energy and power densities are decoupled.

By presenting energy and power densities in a consistent manner, we can facilitate introduction of new materials and fi nd solutions for EES challenges the world faces. National and international testing facilities should be created for benchmarking electrodes and devices similar to the facilities that exist for benchmarking photovoltaics. Clear rules for reporting the performance of new materials for EES devices would help scientists who are not experts in the fi eld, as well as engineers, investors, and the general public, who rely on the data published by the scientists, to assess competing claims.

A B

Device with 1 mg of carbon per cm 2 Device with a 2-µm film (0.07 mg of carbon per cm 2 ) Device with a 120-µm film (4 mg of carbon per cm 2 ) Device with 10 mg of carbon per cm 2 Energy density (Wh/kg) Energy density (Wh/liter) Power density (Wh/liter) Power density (Wh/kg) A tale of two plots. One way to compare electrical energy storage devices is to use Ragone plots [START_REF]Details for the fi gure plots are[END_REF], which show both power density (speed of charge and discharge) and energy density (storage capacity). These plots for the same electrochemical capacitors are on a gravimetric (per weight) basis in (A) and on a volumetric basis in (B). The plots show that excellent properties of carbon materials will not translate to medium-and large-scale devices if thin-fi lm and/or low-density electrodes are used [START_REF]Details for the fi gure plots are[END_REF].


Anna Jezierska

Caroline Chaux

Hugues Talbot

Jean-Christophe Pesquet

Image quantization under spatial smoothness constraints

Keywords: Convex optimization, combinatorial optimization, proximal methods, graph cuts, image coding, segmentation, denoising, entropy

INTRODUCTION

Quantization is a fundamental task in digital image processing. It plays a prominent role in early processing stages such as image digitization, and it is essential in lossy coding. It bears close resemblance to high level tasks such as denoising, segmentation, and data classification. In particular, quantizing a grey scale image in Q levels can be viewed as a classification or segmentation of the image in Q areas following an intensity homogeneity criterion. Each segmented area then corresponds to a decision class of the quantizer.

A classical solution for designing an optimal quantizer of an image is provided by the celebrated Lloyd-Max (LM) algorithm [START_REF] Max | Quantizing for minimum distortion[END_REF][START_REF] Lloyd | Least squares quantization in PCM[END_REF]. An extension to the general vector case is the LBG algorithm [START_REF] Linde | An algorithm for vector quantizer design[END_REF]. The LM algorithm proceeds iteratively by alternatively optimizing quantization levels and decision levels so as to minimize a flexible quantization error measure. It is known to present good convergence properties in practice. However, one drawback is the lack of spatial regularity of the quantized image. Spatially smooth properties may THIS WORK WAS SUPPORTED BY THE AGENCE NATIONALE DE LA RECHERCHE UNDER GRANT ANR-09-EMER-004-03.

be useful in low-rate compression when using advanced coding algorithms (e.g based on run length, differential or multiresolution techniques). It may also be of interest for quantizing images featuring noise. In the latter case, quantization can be viewed as a means for denoising discrete-valued images that are piecewise constant.

Since the LM algorithm is closely related to K-means, which are widely used in data classification, enforcing spatial smoothness of the quantized image could be achieved by resorting to fuzzy C-means clustering techniques and their extensions [START_REF] Chuang | Fuzzy c-means clustering with spatial information for image segmentation[END_REF]. These algorithms are however based on local measures of smoothness.

In this paper, we propose a quantization method that enforces some global spatial smoothness. This is achieved by introducing an adjustable regularization term in the minimization criterion, in addition to a quantization error measure. Similarly to the LM algorithm, the optimal design of the quantizer is performed iteratively by alternating the minimization of a label field i D and of a quantization level vector r. The latter minimization reduces to a convex optimization problem whereas the former is carried out by efficient combinatorial optimization techniques.

The problem is formulated in Section 2 and the notation used throughout the paper is introduced. Section 3 describes the proposed quantizer design algorithm. Section 4 provides more details on the combinatorial optimization step. Finally, some simulation results are provided in Section 5 to show the effectiveness of the proposed quantization method before a conclusion is drawn in Section 6.

PROBLEM

In order to define a quantizer, we introduce the following variables: Q is a positive integer, D = (D k ) 1≤k≤Q is a partition of {1, . . . , N } × {1, . . . , M } and r = (r 1 , . . . , r Q ) is vector in C Q , where C Q is the closed convex cone:

C Q = (s 1 , . . . , s Q ) ∈ R Q s 1 ≤ • • • ≤ s Q . (2.1)
The partition D can be characterized by the label image i D (n, m) 1≤n≤N,1≤m≤M , defined as: for every (n, m) ∈ {1, . . . , N } × {1, . . . , M } and k ∈ {1, . . . , Q},

i D (n, m) = k ⇔ (n, m) ∈ D k .
(2.2)

A scalar quantized image over Q quantization levels r 1 , . . . , r Q associated with the partition D is then given by

q i D ,r = (r i D (n,m) ) 1≤n≤N,1≤m≤M ∈ {r 1 , . . . , r Q } N ×M Let f = f (n, m) 1≤n≤N,1≤m≤M ∈ R N ×M denote the original image. An "optimally" Q-level quantized im- age q i D ,r of f is usually obtained by looking for (i D , r) ∈ {1, . . . , Q} N ×M × C Q solution to the following problem: minimize i D ,r ϕ(q i D ,r , f ) (2.3)
where ϕ is some measure of the quantization error. Standard choices for ϕ are:

• the weighted ℓ p norm measure (p ∈ [1, +∞[) ∀g = g(n, m)) 1≤n≤N,1≤m≤M ∈ R N ×M ϕ(g, f ) = N n=1 M m=1 ω n,m |g(n, m) -f (n, m)| p where (ω n,m ) 1≤n≤N,1≤m≤M ∈ [0, +∞[ N ×M .
As special cases, the mean square error criterion is found when p = 2, and the mean absolute error criterion when p = 1.

• the sup norm measure

∀g = g(n, m)) 1≤n≤N,1≤m≤M ∈ R N ×M ϕ(g, f ) = max (n,m)∈{1,...,N }×{1,...,M } |g(n, m) -f (n, m)|.
One of the drawbacks of this approach is that it does not guarantee any spatial homogeneity of the resulting quantized image. To alleviate this shortcoming, we propose to solve the following problem:

minimize i D ,r ϕ(q i D ,r , f ) + ρ(i D ) (2.4)
where ρ is some regularization function which is used to promote the spatial regularity of the label image. Typical choices for ρ that can be made are the following

• isotropic variation functions ρ(i D ) = µ N -1 n=1 M -1 m=1 ψ( ∇i D (n, m) ), µ ≥ 0 (2.5) where ∇i D (n, m) = i D (n + 1, m) -i D (n, m), i D (n, m + 1) -i D (n, m) is the discrete gradient of i D at location (n, m).
• anisotropic variation functions

ρ(i D ) = µ N -1 n=1 M m=1 ψ(|i D (n + 1, m) -i D (n, m)|) + N n=1 M -1 m=1 ψ(|i D (n, m+1)-i D (n, m)|) , µ ≥ 0. (2.6)
In the above two examples, ψ is a function from [0, +∞[ onto [0, +∞[. When ψ is the identity function, the classical isotropic or anistropic total variations are obtained. If ψ = (•) 2 , then a Tikhonov-like regularization is performed. Another interesting choice of ψ is the binary cost function

(∀x ∈ [0, +∞[) ψ(x) = 0 if x = 0 1 otherwise.
(2.7)

PROPOSED OPTIMIZATION METHOD

Even if ϕ(•, f ) and ρ are convex functions, Problem (2.4) is a nonconvex optimization problem due to the fact that i D belongs to a nonconvex set of discrete values. In order to solve numerically this problem, we propose to use the following alternating optimization algorithm:

Fix Q ∈ N * and r (0) ∈ C Q . For ℓ = 0, 1, . . . i (ℓ) D ∈ Argmin i D ∈{1,...,Q} N ×M ϕ(q i D ,r (ℓ) , f ) + ρ(i D ) r (ℓ+1) ∈ Argmin r∈C Q ϕ(q i (ℓ) D ,r , f )
It is worth noticing that this algorithm constitutes an extension of the LM algorithm which would correspond to the case when ρ is the null function. At each iteration ℓ, the determination of i (ℓ) D given r (ℓ) is a combinatorial optimization problem for which there exist efficient solutions for some choices of ϕ and ρ, as shown in the next section.

In turn, if ϕ(•, f ) is a convex function, the determination of r (ℓ+1) given i (ℓ) D is a conic constrained convex optimization problem the solution of which can be determined numerically. For example, in the case of the weighted ℓ p -norm criterion, this problem is equivalent to find

r (ℓ+1) ∈ Argmin r∈C Q Q k=1 (n,m)∈D (ℓ) k ω n,m |r k -f (n, m)| p (3.1)
It is thus possible to solve (3.1) through existing convex optimization approaches [START_REF] Boyd | Convex Optimization[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF].

Note that, in the proposed iterative algorithm, it may happen that a decision class D (ℓ) k0 , for some k 0 ∈ {1, . . . , Q} becomes empty at iteration ℓ. In this case, a random value of r (ℓ+1) k0 is drawn, so that r (ℓ+1) remains in C Q . It can then be shown that the resulting algorithm has convergence properties similar to the LM algorithm.

OPTIMAL COMBINATORIAL PARTITIONING

We now consider combinatorial optimization methods for finding

i D ∈ Argmin i D ∈{1,...,Q} N ×M ϕ(q i D ,r , f ) + ρ(i D ) (4.1)
for a given value of r ∈ C Q . Here we seek to use standard methods in combinatorial optimization which have proved useful in applications to denoising. In this context, a common form for regularization problems is the following:

minimize i D ∈{1,...,Q} N ×M φ(i D , f ) + ρ(i D ), (4.2)
where φ is a data fidelity function, ρ a regularization function, f the initial image and i D the resulting discrete one. To formulate our problem in this framework, we need to introduce the auxiliary function

χ r : {1, . . . , Q} N ×M → {r 1 , . . . , r Q } N ×M i D → q i D ,r .
Then, our problem becomes

minimize i D ∈{1,...,Q} N ×M ϕ(χ r (i D ), f ) + ρ(i D ). (4.3)
Note that χ r is monotonic but nonlinear. Note further that the set {r 1 , . . . , r Q } changes at each iteration of the complete algorithm. However, during the regularization step this set is fixed, so we are free to choose any optimization method that can be framed as in (4.2).

In this article we use graph cut-based algorithms, which have been shown to be effective in the context of smoothing, denoising and segmentation [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]. More specifically, a solution can be found at each regularization step using the graph proposed by Ishikawa [START_REF] Ishikawa | Mapping image restoration to a graph problem[END_REF], with some small modifications. Namely, i) labels ι take values from 1 to Q, ii) data edges for each vertex u ι,j have capacity ϕ(r ι , f (n j , m j )), and iii) penalty edges have capacity µ. In case of few quantized levels Q, the Ishikawa framework is very efficient. This graph method is known to yield globally optimal solutions for convex functions ϕ and ρ as in (2.6) with ψ the identity function. Alternatively, the convenient though suboptimal alphaexpansion algorithm [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] can be used for all combinations of ϕ and ψ, where ρ is the anisotropic TV of (2.6), as described in Section 2. In order to solve Problem (4.3) with the αexpansion algorithm, we propose to define the capacities of edges as described in Table 1. Note that α takes values from 1 to Q.

Other optimization methods can also be used. When minimizing isotropic TV as in (2.5), one might want to use for instance Chambolle's algorithm [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. Similarly to the Ishikawa framework, we can obtain the global optimum in this case. In the following, due to lack of space, we illustrate our proposed method only using graph-cut based regularization methods.

SIMULATION EXAMPLES

In this section, we present two experiments in order to demonstrate the practical performance of our method in case of: i) low resolution quantization, and ii) quantization in the presence of noise. edge capacity condition

[α, u] +∞ ∀u, i u = α [α, u] ϕ(r iu , f u ) ∀u, i u = α [α, u] ϕ(r α , f u ) ∀u [u, t {u,v} ] ψ(|i u -α|) [t {u,v} , v] ψ(|α -i u |) ∀{u, v}, i u = i v [α, t {u,v} ] ψ(|i u -i v |) [u, v] ψ(|i u -α|) ∀{u, v}, i u = i v
We compare our approach with the LM method [START_REF] Lloyd | Least squares quantization in PCM[END_REF]. The algorithm performance is measured by the SNR between the original and quantized images and also by the the Shannon entropy of order (2, 2) (that is the entropy over image blocks of size 3 × 3). In the following experiments, we used 8 bit microscopy images of size 512 × 512 from public domain (source: http://remf.dartmouth.edu).

First, we consider grey-scale image quantization over 8 levels. The original image is shown in Fig. 1(a). Here we minimize the energy function defined by (2.4), with function ϕ defined as the ℓ 2 norm and ρ defined by (2.6) where ψ is the identity. The regularization parameter µ was hand-optimized to 400. The global optimum solution to the problem can be found with the modified Ishikawa graph shown in Section 4. Both methods, LM and ours, were initialized with uniform decision levels. In order to solve (3.1), the FISTA algorithm [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problem[END_REF] was used. As expected, Fig. 1 shows that the images have a better spatial smoothness with our method than with LM. This is also verified by inspecting the entropy value, which in our case is equal to 0.56 bpp and in case of LM to 0.84 bpp. In this example, we have shown that, in case of quantization with high reduction of levels, our method provides smaller entropy rate while maintaining the desired fidelity.

Next, we consider the problem of grey-scale image quantization over 32 levels in the presence of noise. To generate a noisy image (Fig. 2(b)) from the original one (Fig. 2(a)), we added a zero-mean Laplacian noise with variance 202 (this corresponds to an initial SNR between original and noisy image equal to 18 dB). We applied our algorithm minimizing energy (2.4), where ϕ is the ℓ 1 norm, and ρ is defined by (2.6) where ψ is the binary cost-function (2.7). The associated regularization parameter µ was experimentally chosen to 15. Both methods, LM and ours, were initialized with decision levels computed on a cumulative histogram. To minimize the energy, we used a modified alpha-expansion graph as described in Section 4. One can observe that there is almost no noise in our result (Fig. 2(d)), while the LM (Fig. 2(c)) pre- serves noise in the image. This is also verified by inspecting SNR values, which are equal to 27.6 dB for our method and 18.3 dB for LM. The difference is even greater in terms of entropy: we obtained 0.49 bpp for our method and 1.64 bpp for LM. In this example, we have shown, that in case of quantization in the presence of noise, our method reconstructs the original image, while quantization is being performed.

CONCLUSION

In this paper, we have proposed a new quantization approach while enforcing spatial smoothness in the resulting image. We have shown that this approach is robust to the presence of noise. These features may be interesting in the context of image compression, in particular to avoid sparing bitrate in coding noise corrupting the image in some applications. We have shown that our method may outperform standard quantization design methods, both in terms of entropy and visual quality.

In future work, we will extend our work to multichannel images and vector quantization.

Fig. 1 :

 1 Fig.1: Figures (a,b,c) illustrate fragment of original image, LM and our result, respectively.

Fig. 2 :

 2 Fig. 2: Figures (a,b) illustrate original and noisy images, respectively. Images after quantization are presented in Figure (c) for LM algorithm and (d) for our method.

Table 1 :

 1 Capacity for the α-expansion graph ([x1, x2] denotes an edge between nodes x1 and x2, iD(nu, mu) is denoted by iu, f (nu, mu) is denoted by fu, {u, v} denotes a pair of neighbouring nodes u and v).