R Wilson

F Dalaudier

H Luce

Atmospheric Measurement Techniques

 to evaluate the climatology of atmospheric turbulence through the detection of overturns in the free atmosphere by applying a Thorpe analysis on relatively low vertical resolution (LR) profiles collected from standard radiosoundings. Since then, several studies based on this idea have been published. However, the impact of instrumental noise on the detection of turbulent layers was completely ignored in these works. The present study aims to evaluate the feasibility of overturns detection from radiosoundings. For this purpose, we analyzed data of two field campaigns during which high-resolution (HR) soundings (10-20 cm) were performed simultaneously with standard LR soundings. We used the raw data of standard meteorological radiosondes, the vertical resolution ranging from 5 to 9 m.

A Thorpe analysis was applied to both LR and HR potential temperature profiles. A denoising procedure was first applied in order to reduce the probability of occurrence of artificial overturns, i.e. potential temperature inversions due to instrumental noise only. We then compared the empirical probability density functions (pdf) of the sizes of the selected overturns from LR and HR profiles.

From HR profiles measured in the troposphere, the sizes of the detected overturns range from 4 to ∼1000 m. The shape of the size pdf of overturns is found to sharply decrease with increasing scales. From LR profiles, the smallest size of detected overturns is ∼32 m, a similar decrease in the shape of the pdf of sizes being observed. These results suggest Correspondence to: R. Wilson (

1 Introduction [START_REF] Thorpe | Turbulence and mixing in a Scottish Lock[END_REF] proposed an elegant and simple-to-use method allowing to identify and to characterize overturns produced by turbulent patches from in situ measurements in oceans and lakes. The method is based on the comparison between an observed vertical profile of potential density and a reference profile corresponding to a minimum state of available potential energy. The reference profile is built by sorting in increasing order the observed potential density profile. It corresponds to an adiabatic re-arrangement of the observed fluid parcels. This method can be applied to atmospheric data as well by considering potential temperature for non saturated air instead of potential density. Since a vertical profile of potential temperature in a stably stratified atmosphere is a monotonic function of altitude, overturns can be easily identified by comparing the observed and sorted profiles. Such overturns are the signature of convectively unstable regions which can be due either to convective or dynamical instabilities of the flow (prior to the onset of developed turbulence) or to inversions produced by turbulent eddies. The convective and dynamical instabilities giving rise to overturns are expected to break down and to produce small-scale Published by Copernicus Publications on behalf of the European Geosciences Union. R. Wilson et al.: Small-scale turbulence from standard radiosondes turbulence. The detected overturns are thus believed to be associated with the occurrence of turbulent events. However, we cannot distinguish from the only in situ measurements, overturns resulting from the rolling up of air due to instabilities before turbulent breakdown from those produced by turbulent eddies in a mixing layer.

Simultaneous observations of convectively unstable regions from in situ and additional instruments (such as radars, lidars or aircrafts, ...) remain very sparse in the free atmosphere. [START_REF] Luce | Strong mixing events observed near the tropopause with the MU radar and hight-resolution Balloon techniques[END_REF] reported observations of strong and isotropic radar echoes (suggesting the presence of turbulence) probably related to several unstable or nearly neutralized layers which were detected by high-resolution balloon measurements. These observations support the hypothesis of turbulent events related to overturns in the potential temperature profile. However, such a comparison cannot be fully conclusive since the air masses encountered by the balloon passed at distances of few kilometers from the radar, 1.5 to 16 km, with a time delay of about 10-16 min. Further experimental validation is still desirable.

In principle, the Thorpe method can be applied to standard sounding data, either conductivity-temperaturedensity (CTD) in the oceans or lakes, (e.g. [START_REF] Thorpe | Turbulence and mixing in a Scottish Lock[END_REF][START_REF] Galbraith | Identifying overturns in CDT profiles[END_REF][START_REF] Ferron | Mixing in the Romanche fracture zone[END_REF][START_REF] Alford | Observation of overturning in the thermocline: The context of ocean mixing[END_REF] or pressure-temperature (PT) in the atmosphere [START_REF] Luce | Strong mixing events observed near the tropopause with the MU radar and hight-resolution Balloon techniques[END_REF][START_REF] Gavrilov | Turbulence parameter estimations from high-resolution balloon temperature measurements of the MUTSI-2000 campaign[END_REF][START_REF] Clayson | On Turbulence and Mixing in the Free Atmosphere Inferred from High-Resolution Soundings[END_REF][START_REF] Wilson | Turbulent Patch Identification in Potential Density/Temperature Profiles[END_REF]. It has been recently proposed by [START_REF] Clayson | On Turbulence and Mixing in the Free Atmosphere Inferred from High-Resolution Soundings[END_REF] to apply Thorpe analyses to the huge data-base of meteorological radiosondes (hereafter noted RS), in order to infer the space-time variability of atmospheric turbulence in the free atmosphere. Since then, several such studies were published [START_REF] Alappattu | First observations of turbulence parameters in the troposphere over the Bay of Bengal and the Arabian Sea using radiosonde[END_REF][START_REF] Nath | Turbulence characteristics over tropical station Gadanki (13.5 • N, 79.2 • E) estimated using highresolution GPS radiosonde data[END_REF]. Also, [START_REF] Balsley | On the use of Slow Ascent Meter-Scale (SAMS) Radiosondes for observing overturning events in the free atmosphere[END_REF] proposed to use slow-ascent radiosondes for observing turbulence in the atmosphere.

However, the issue of instrumental noise was not addressed at all in these aforementioned works. From our point of view, the issue of noise is a key point of the Thorpe method, that must be considered with great care. Oceanographers long-time know that instrumental noise can generate artificial overturns which are very difficult to distinguish from real ones (e.g. [START_REF] Thorpe | Turbulence and mixing in a Scottish Lock[END_REF][START_REF] Galbraith | Identifying overturns in CDT profiles[END_REF][START_REF] Ferron | Mixing in the Romanche fracture zone[END_REF][START_REF] Gargett | Determining Thorpe Scales from Ship-Lowered CTD density profiles[END_REF]. Surprisingly, except a few works [START_REF] Gavrilov | Turbulence parameter estimations from high-resolution balloon temperature measurements of the MUTSI-2000 campaign[END_REF][START_REF] Wilson | Turbulent Patch Identification in Potential Density/Temperature Profiles[END_REF], the issue of noise is generally ignored in most studies dealing with atmospheric turbulence arguing that: "it is less of an issue in the free atmosphere, where the mean overall stratification is quite strong..." [START_REF] Clayson | On Turbulence and Mixing in the Free Atmosphere Inferred from High-Resolution Soundings[END_REF]. From our point of view, as shown by [START_REF] Wilson | Turbulent Patch Identification in Potential Density/Temperature Profiles[END_REF], this assertion does not hold, neither for standard meteorological measurements, nor for HR atmospheric measurements, because nearly neutral layers are very common, even at the vertical resolution of the standard radiosondes.

Hereafter, according to the terminology used by [START_REF] Johnson | Effect of noise on Thorpe Scales and run lengths[END_REF], an inversion is defined as a localized decrease of potential temperature versus height, whatever its origin may be, either induced by turbulent motions or by noise. (The present definition should not be confused with the standard definition of an inversion in meteorology which refers to a layer for which temperature increases whith height, for instance at the top of the planetary boundary layer.) The term overturn will specifically refer to an inversion resulting from atmospheric motions (developed turbulence or convective overturns produced by instabilities in the flow). The selection of artificial inversions as overturns, may result in a dramatic overestimation of their frequency of occurrence, thus leading to overestimates of both the turbulent fraction of the atmosphere and the turbulence energy. The space-time inhomogeneity of turbulence can give rise to diffusion coefficient estimates ranging over several orders of magnitude [START_REF] Wilson | Turbulent diffusivity in the free atmosphere inferred from MST radar measurements: a review[END_REF]. It is therefore crucial to apply a quantitative procedure allowing to discriminate overturns from noise-induced inversions.

A quantitative method for selecting overturns in a potential temperature profile has been recently proposed by [START_REF] Wilson | Turbulent Patch Identification in Potential Density/Temperature Profiles[END_REF]. The method is based on a hypothesis test applied to the data range within the detected inversions. The range of a data sample is defined as the difference between the maximum and the minimum values in that sample (the term "sample" refers here to a subset of a dataset). The detected inversions which cannot be distinguished from artificial inversions, i.e. inversions induced by instrumental noise, are not selected as overturns. The selection method will be described in more details later.

The present paper addresses the following two questions:

1. Can we detect small scale turbulence (or more exactly overturns) from meteorological radiosondes in the troposphere and lower stratosphere?

2. If yes, what are the sizes of the detectable overturns considering the vertical resolution and the noise level of radiosonde measurements?

In order to address these issues, we analyzed the datasets of two field campaigns during which high-resolution (HR) balloon soundings, 10 or 20 cm vertical resolution, and standard RS with relatively low-resolution (LR), 5 to 9 m vertical resolution, were simultaneously performed: the SFT campaign (1998) and the MUTSI campaign (2000). (SFT stands for Structure Fine de Temperature) and MUTSI for MU radar, Temperature Sheets and Interferometry). Seven (fourteen) RS were launched during the SFT and MUTSI campaigns respectively. One of the seven LR flights of the SFT campaign is not used in this study as the ascent speed was abnormally slow.

The paper is structured as follows. The datasets and data processing methods are described in Sect. 2. The analysis methods of the potential temperature profiles are presented The first difference of altitudes z inferred from raw measurements (cyan) and by using the approximated pressures (black).

in Sect. 3, a quantitative criterion for the selection of the overturns being described. In Sect. 4, we present the distribution in size of the selected overturns from both LR and HR profiles in the troposphere and in the lower stratosphere. Discussion and conclusions are given in the final section.

2 Dataset and data processing

The datasets

Meteorological RS provide temperature, pressure and relative humidity measurements. With the aim of detecting small scale motions, we used the raw RS data during the SFT and MUTSI campaigns. The RS data were obtained with Vaisala radiosondes RS80G (SFT) and RS90G (MUTSI) at an initial sampling frequency of 0.7 Hz. The altitude resolution is variable due to the variations of the vertical ascent velocity of the balloon.

For each of the campaigns (SFT and MUTSI), three of the balloons carried both HR and LR sensors. The HR profiles were acquired with a sampling frequency of 25 Hz for the SFT flights and 50 Hz for the MUTSI flights. The reader can find more details about the HR measurements in [START_REF] Luce | Strong mixing events observed near the tropopause with the MU radar and hight-resolution Balloon techniques[END_REF] and [START_REF] Gavrilov | Turbulence parameter estimations from high-resolution balloon temperature measurements of the MUTSI-2000 campaign[END_REF]. The HR and LR sensors were on the same gondola, and the data were acquired simultaneously.

Construction of potential temperature profiles from raw data of meteorological radiosondes

For the RS80 and RS90 Vaisala radiosondes used in this study, the measured temperature and pressure are quantized with a resolution of 0.1 K and 10 Pa respectively. As previously mentioned, the measurements are irregularly sampled with altitude.

In order to perform a Thorpe analysis, it is necessary to have data sampled at a regular vertical step, the reference profile being a vertical rearrangement of the observed profile. The measured profiles need to be resampled at a regular vertical interval δz. The first task is to estimate the altitude of the measurements. By assuming hydrostatic equilibrium, the i-th altitude difference, z M i = z M i+1 -z M i , inferred from the measured pressure and temperature at level i, P M i and T M i respectively, is given by

z M i = - R T M i g i P M i P M i ( 1
)
where R is the specific gas constant for air (J kg -1 K -1 ) and g i the acceleration of gravity at level i.

P M i = P M i+1 -P M i
is the pressure difference between successive measurements.

The superscript "M" stands for "measured". The i-th altitude level is simply

z M i = z 0 + i-1 k=0 z M k (2)
Due to the measurement noise, the first differences of measured pressure, i.e. the series of the differences between successive measurements, are not always negative (as it should be during the ascent of the balloon) especially for high altitude levels, i.e. above 20 km altitude. The top panel of Fig. 1 shows the first differences of pressure measurements (cyan curve). It is observed that quantization effects are increasingly important with increasing altitudes. As a consequence, the calculated raw altitude, z M , from the measured pressure is not monotonically increasing during the ascent of the balloon: the altitude differences are quantized, the vertical steps growing with increasing altitudes (Fig. 1, lower panel), some differences being negative. Consequently, the raw altitudes z M cannot be used for building the vertical profile.

A monotonic pressure profile is desirable for overcoming this difficulty. For that purpose, we estimate the pressure, P A , which, if quantized, would give the measured pressure, P M . The pressure P A is evaluated through a least square cubic spline approximation (superscript "A" stands for "approximation"). The black curve of the upper panel of Fig. 1 shows the first differences P A which are now negative for all levels. The altitude of measurements z A is estimated from Eq. (1) with P A instead of P M . In the lower panel of Fig. 1, the black curve shows the first differences of altitudes z A estimated from P A and T M . The vertical steps are found to range from 5 to 9 m, due to the variations in the balloon ascent velocity.

For this specific flight (MUTSI-LR10) an average height interval of 7 m is found. The measured temperature T M and approximated pressure linear interpolation or a spline interpolation, the differences between these two methods being found to be negligible. The resampled temperature and pressure, T R and P R , respectively are the data analyzed in the present work. The method by which the instrumental noise is estimated will be described in Sect. 2.3. The pressure, temperature and altitude are not independent quantities as they are related through the hydrostatic equilibrium relationship (Eq. 1). The knowledge of two of the three quantities fully determines the third one. For instance, the knowledge of T R regularly sampled at step δz allows to calculate the vertical profile of pressure for the same altitudes. The calculated pressure from T R i is almost the same as the resampled pressure P R i although T and P are independently measured: they differ by less than 0.02 Pa in the average, giving us strong confidence in the data processing method.

After re-sampling at a constant vertical step, the vertical resolution of LR temperature and pressure profiles is ∼ 6-7 m, depending on the flight. The vertical resolution of the resampled HR profiles is ∼ 20 cm (10 cm) for the SFT and MUTSI datasets, respectively.

The potential temperature θ is inferred from the resampled temperature and pressure, T R i and P R i , by using the relation for ideal diatomic gas:

θ i = T R i 1000 P R i 2/7 (3)
where T is the temperature (K), and P the pressure (hPa). The right panel of Fig.

(2) shows the potential temperature plus and minus two standard deviations of noise. Several thick nearly neutralized layers can be seen in the troposphere, around 4 km altitude and above 10 km altitude as well as in the stratosphere, below the altitude of 20 km altitude.

Estimation of instrumental noise

A key point of the selection method (presented in the next section) is to estimate the noise on potential temperature resulting from instrumental noise on T and P . The noise standard deviations for the LR and HR data are experimentally estimated from both the temperature and pressure resampled profiles. After comparing several methods based on spectral or on wavelet analyses for some of them, we finally retained a simple-to-use method based on the estimate of the variance of the first differences of the data (T or lnP ). One important outcome of these tests is that the instrumental noise on temperature and pressure depends on altitude for both the LR and HR profiles. Therefore, we must estimate a vertical profile of the instrumental noise level. The method used for estimating the noise variance was made as follows:

1. Splitting the entire profile in segments of about 200 m in length (i.e. 32 or 2000 points according to the vertical resolution).

2. Removing a linear trend for each segment.

3. Calculating the variance of the first differences. This variance is an estimate of twice the noise variance.

4. Smoothing the resulting noise profile.

Figure 3 shows the estimated level of instrumental noise for the temperature and pressure measurements for flight MUTSI-LR10. The stairs step curve shows the noise estimates for each data segment, the smoothed curve displaying Atmos. Meas. Tech., 4, 795-804, 2011 www.atmos-meas-tech.net/4/795/2011/ the noise level profile after smoothing (using a least square spline approximation). For the pressure measurements, the noise is estimated from the logarithm of the measurements in order to apply a linear detrending (step 2).

For the twenty LR flights of the SFT and MUTSI campaigns, the averaged noise standard deviation of the resampled temperature profile is found to be 27 mK. For the six HR flights, the noise standard deviation of temperature data is found to be 4 mK in the average. The noise level of the LR temperature measurements is observed to decrease slightly with altitude (Fig. 3). The noise standard deviation of the resampled pressure profile is found to be negligible: 3 × 10 -2 Pa in the average. Consequently, the noise on the potential temperature essentially depends on the instrumental noise of temperature measurements. From Eq. ( 3), the standard deviation of noise of potential temperature, σ θ is estimated by

σ θ = θ σ T T 2 + 2 7 σ P P 2 ≈ θ σ T T = 1000 P 2/7 σ T (4)
where σ T (σ P ) are the noise standard deviation of temperature and pressure respectively. The noise level of potential temperature grows with altitude, due to the pressure term of Eq. ( 4): it increases from ∼35 to ∼70 mK (2.6 to 3.5 mK) for LR (respectively HR) profiles for heights ranging from the ground level up to 27 km.

3 Analysis method of the potential temperature profiles

Selection method of overturns

In order to evaluate the filtering effects due to both the vertical resolution and the level of instrumental noise, we analyzed the profiles by using an identical selection criterion. For discriminating overturns from inversions induced by noise, we used a hypothesis test on the range of the data. The method is described in detail by [START_REF] Wilson | Turbulent Patch Identification in Potential Density/Temperature Profiles[END_REF].

Here, we shortly describe the three important steps of the method.

1. Determination of the optimal vertical resolution. We first estimate vertical resolution that should be used by taking account for the mean stratification and the level of the instrumental noise. Such a resolution can be inferred from the bulk trend-to-noise ratio (tnr), ζ of the profile. The bulk tnr, ζ , is a measure of the average of first differences of the data (potential temperature) τ , scaled by the noise σ θ :

ζ = τ σ θ (5)
with

τ = θ i+1 -θ i = 1 n -1 (θ n -θ 1 ) (6)
where θ i is the potential temperature measured at altitude z i (1 ≤ i ≤ n). The average of first differences τ depends on both the stability of the stratification and the vertical resolution since

τ = ∂θ ∂z δz (7)
where z is the vertical coordinate and δz the (constant) vertical step. The bulk tnr provides an indication on the minimum size for the overturns to be reliably detected [START_REF] Wilson | Turbulent Patch Identification in Potential Density/Temperature Profiles[END_REF]. The smaller ζ be, the larger the data sample must be in order to be distinguishable from a pure noise sample. If ζ is too small (smaller than 1 typically), a preliminary denoising procedure is required. The denoising procedure is based on both a filtering and an undersampling of the data. It aims at increasing the tnr while preserving the independence of the noise of each data bin. The tnr of the filtered and undersampled data increases by a factor m 3/2 where m is the undersampling factor (i.e. the number of bins of the running filter).

2. Detection of the inversions. The potential temperature profile (denoised if needed) is sorted. The profile of Thorpe displacements can then be estimated. The (artificial and real) inversions are found from the cumulative sum of the Thorpe displacements since the sum of Thorpe displacements is null within an inversion.

3. Selection of the overturns. An inversion is identified as an overturn if the range of the potential temperature data belonging to the inversion exceeds a large prescribed percentile of the range of a noise sample of same length, with standard deviation equal to that of the instrumental noise. For practical purpose, we tabulated the moments of the range, as well as various percentiles, for normally distributed random variables as a function of the sample size (ftp.aero.jussieu.fr/pub/os/WN.txt).

The vertical profile of tnr ζ LR for the MUTSI-LR10 flight is shown in Fig. 4a. The bulk tnr ζ LR is ∼0.9 (dashed line), the local values showing minima smaller than 10 -1 around 11 km altitude. For such a bulk tnr, one cannot expect to identify overturns smaller than n ∼ 60 bins, i.e. ∼420 m (see Fig. 5 of [START_REF] Wilson | Turbulent Patch Identification in Potential Density/Temperature Profiles[END_REF]. In order to improve the detection threshold, a denoising procedure is required. For the region with tnr of 0.1, an undersampling m = 2 allows to reach a local tnr of 2 √ 2 × 10 -1 ≈ 0.3, allowing to detect overturns of 20 bins, i.e. 20 × 14 = 280 m.

Overturns selection

The average tnr being less than 1 for both LR and HR profiles in the troposphere, a denoising procedure is applied on the two datasets. The undersampling factor is 2 for the LR data, and is 5 ( 11 respectively. The bulk tnr reaches 2-3 for both datasets. Consequently, the vertical resolution is ∼12 m for the LR profiles, and is 1 m (SFT) or 1.1 m (MUTSI) for the HR profiles.

The tnr values of both LR and HR data are substantially larger in the stratosphere than in the troposphere due to the higher stability of the stratosphere. Consequently, the undersampling factor required to get a tnr of ∼2-3 is smaller than in the troposphere, allowing to reach better vertical resolutions. The vertical resolution of the LR data is ∼6 m (no undersampling required) whereas the vertical resolution of the HR data is ∼0.5 m.

In order to evaluate the impact of the selection method (including the denoising procedure) we compare the results of Thorpe analyses without and with the selection procedure. Figure 5 shows histograms of the number of occurrences of inversions (blue-filled) and overturns (red-transparent) as a function of their sizes for all the LR profiles. The distribution of the detected inversions (blue-filled histogram) results from the original potential temperature profiles (vertical resolution ∼6 m), without applying any denoising and selection procedure. The distribution of the selected overturns for the same LR data (red-transparent histogram) results from the denoised (undersampled and filtered) profiles at a vertical resolution of ∼12 m. The selection procedure shows that only 11 % of the detected inversions can be regarded as overturns. Consequently, due to both the noise level of RS data and the average stability of the troposphere, most of the detected inversions (91 %) are only due to instrumental noise! It is therefore necessary to apply a selection procedure in order to identify real overturns, at least in the troposphere.

Size distribution of the detected overturns

We analyzed the LR and HR profiles by using the same selection method for detecting overturns. As mentioned in the previous section, the bulk tnr ζ of all profiles ranges from 2 to 3 after applying the denoising procedure. We present in this section the empirical probability density functions (pdf), i.e. the relative frequencies, of the sizes of the selected overturns from both LR and HR datasets in the troposphere and in the stratosphere.

Results for tropospheric data

Figure 6 shows the number of occurrences of Thorpe displacements (in absolute value) as a function of their sizes for the six HR profiles. The two histograms result from all displacements found within the selected overturns in the troposphere. The corresponding histograms from the HR and LR profiles are shown in green-filled and red-transparent shadings, respectively. Due to the large dynamic of the observed values, the histograms are shown in a log-log scale. The widths of the intervals (abscissa) being constant on a log scale, they increase as a power law. We simply used a dyadic scale, the maximum value of each interval being equal to twice the minimum value (i.e. 1-2 m, 2-4 m, and so on.). The number of occurrences is not directly comparable because the vertical step of LR profile is in the present case about 12 times larger than that of HR profiles (12 m vs. 1 m). The considered displacements range from 1 m to ∼ 500 m and from 12 to ∼1000 m for the HR and LR profiles, respectively. As expected, the smallest displacements correspond to the vertical resolution, ∼1 m (within the 1-2 m interval) for HR data and ∼10 m (within the 8-16 m interval) for LR data.

On the other hand, about 20 displacements observed from LR profiles are larger than the largest displacement obtained from HR profiles. Such larger displacements in the LR profile are likely the consequence of the poorer vertical resolution. Indeed, the LR profiles do not always allow to identify distinct overturns within weakly stratified regions. Such regions may appear as a single overturn in the LR profile, whereas the HR profile reveals several distinct overturns.

Several examples of such merged overturns can be seen in the paper by [START_REF] Wilson | Turbulent Patch Identification in Potential Density/Temperature Profiles[END_REF]. The largest LR displacements likely belong to these merged overturns. This effect is clearly a major source of bias when estimating Thorpe lengths and turbulence parameters from LR profiles in weakly stable regions.

The two histograms of Fig. ( 6) are proportional to the empirical pdf of Thorpe displacements. They reveal a tendency for the probability densities of displacements to decrease with increasing sizes. The pdf from LR measurements appears similar in shape but shifted by about 10 m toward larger scales.

Next, we consider the size of the selected overturns. The empirical pdf's of the overturn sizes for both LR and HR profiles are shown in the left panel of Fig. 7. The right panel shows the occurrence number of overturns according to their size. The empirical pdf's are of course proportional to the number of occurrences. However, the plot of the occurrence numbers (right panel) allows a direct comparison of the numbers of events independently detected from HR and LR data during the six HR/LR flights. The size of the shallowest overturns from the HR data lies between 4 and 8 m whereas the shallowest overturns from LR data are within the 32-64 m interval. It should be noted that the inner scale of turbulence, which ranges from ∼0.1 to 10 cm in the free troposphere and lower stratosphere [START_REF] Eaton | Preliminary estimates of the vertical profiles of inner and outer scales from White Sands Missile Range, New Mexico, VHF radar observation[END_REF], is much smaller than the smallest scale detectable from the HR profiles. Consequently, the decreasing number of occurrences for the shallowest overturns in the HR profile is a cutoff effect similar to the one observed at much larger size for the LR profiles. The relative frequencies and the cumulative distributions of the overturn sizes are detailed in Table 1. It is found that 93 % of the overturns identified from the HR profiles are smaller than 32 m, i.e. smaller than the shallowest overturns detected from the LR profiles. These results demonstrate that it is only possible to identify overturns which size belongs to the tail of the pdf from the LR profiles in the troposphere. The selected overturns from the LR profiles represent ∼7 % of the overturns found from HR profiles. On the other hand, the number of selected overturns from HR and LR profiles compares quite well, except for the smallest overturns in the LR profiles (the number of overturns in the HR profiles is about 3 times the number of overturns in the LR profiles for the 32-64 m interval). Although a large number of small overturns is missed, this result suggests that the largest overturns are correctly identified from the LR profiles.

Figure 8 shows the empirical pdf of the sizes of overturns for all the flights, i.e. the six HR and twenty LR profiles of the MUTSI and SFT campaigns, in the troposphere. Again, the smallest overturns of the LR profiles are within the 32-64 m interval. Also, the shape of the two probability densities compare quite well, as both densities decrease with a comparable rate. This observation suggests that the tail of the distribution in size of the overturns is correctly sampled from the LR profiles (except maybe for the extreme values, i.e. for the first and last intervals).

Results for stratospheric data

The same selection method was applied to the HR and LR profiles collected in the stratosphere. Figure 9 shows the distribution of the number of occurrences of Thorpe displacements for the selected overturns (HR and LR) in the stratosphere for the six HR profiles. Here, the Thorpe displacements range from 0.5 to ∼50 m for the HR profiles, and from 6 to ∼200 m for the LR profiles. Again, the few larger displacements from LR data (less than 20 occurrences) are likely due to the much lower resolution of LR profiles, leading to a poor definition of the overturns within weakly stable regions.

The left panel of Fig. 10 shows the empirical pdf of the sizes of overturns in the stratosphere for the six HR profiles. From the HR profiles, the size of overturns now ranges from 2 to ∼100 m. The relative frequencies of overturn size are found to sharply decrease with increasing size. From the LR profiles, the overturn size ranges from ∼25 m (16-32 m) to ∼200 m (128-256 m). Table 2 gives the cumulative frequencies of the overturn size. It is observed that 96 % of the selected overturns from the HR profiles have a size smaller than 16 m (Table 2). The right panel of Fig. 10 shows the number of occurrences of the overturn size for these six HR profiles. A striking feature is the small number of selected overturns from the six LR profiles: 85 events, i.e. about fourteen for each profile for all the sounded stratosphere (roughly from 10 to 25 km). As in the troposphere, the number of occurrences of selected overturns according to their sizes is very similar for the LR and HR profiles (except for the smallest class, i.e. 16-32 m, for the LR profiles). One can again concludes that, within the stratosphere, the overturns exceeding ∼25 m in size are correctly sampled from the LR profiles. But such overturns are rather rare in the stratosphere since they only represent about 4 % of all the overturns detected from the HR profiles.

Figure 11 shows the empirical pdf of the overturn size in the stratosphere for all the (six HR and twenty LR) profiles. The sizes of overturns from the LR profiles now range from ∼25 to 512-1024 m. Both pdf's are sharply decreasing with increasing scales. The shape of the two distributions is very similar, giving extra credence to the detection of the larger overturns from the LR profiles. Again, the scales for which LR overturns are detected and selected represent only 4 % of the HR overturns since the selection procedure does not allow to retain overturns smaller than 16 m from LR data.

Conclusions

We have analyzed a dataset of balloon soundings obtained during two field campaigns. The dataset includes six HR profiles (10-20 cm) and twenty LR profiles (6-7 m). The HR and LR profiles were obtained simultaneously, the sensors being on the same gondola. The LR profiles result from raw data collected from standard meteorological radiosondes vertically resampled at a constant step (6 or 7 m). After applying a denoising procedure, we performed a Thorpe analysis. We applied a selection method allowing to reject artificial inversions, i.e. inversions that cannot be distinguished from noise induced inversions. We can only retain 11.4 % (25 %) of the detected inversions as overturns from LR (respectively HR) profiles in the troposphere. In the stratosphere, the situation is worse as the noise on potential temperature increases whereas the sizes of the turbulent overturns diminishes. Only 7.9 % (20.8 %) of the detected inversions are selected as overturns from the LR (respectively HR) profiles. Consequently, for detecting atmospheric turbulence from a Thorpe analysis, the impact of instrumental noise is a key issue which must be carefully considered in the troposphere as well as in the stratosphere.

A noticeable result of the present work, the interpretation of which is beyond the scope of this paper, is that the atmosphere is unstable for more than 56 % in the troposphere, and more than 37 % in the lower-stratosphere for vertical scales larger than one meter (from the HR profiles). The former value is in good agreement with the one reported by [START_REF] Cho | Characterizations of tropospheric turbulence and stability layers from aircraft observations[END_REF] who found that 54 % of the troposphere is unstable when measured down to scale of ∼ 10 m.

An important conclusion of the present work is that, despite the above mentioned limitations and difficulties, it is possible to detect overturns resulting from atmospheric turbulence from the raw data of standard meteorological radiosondes sampled at a rate of 0.7 Hz. However, we showed that only the deepest overturns corresponding to the tail of the size pdf, can effectively be detected from standard meteorological radiosondes, in the troposphere as well as in the stratosphere. From our data set, the detectable overturns from RS in the troposphere correspond to the 7 % largest events that are observed from HR profiles. In the stratosphere, the detected overturns from RS correspond to about the 4 % largest events. The latter conclusion opens interesting perspectives about the means of observing turbulence in the free atmosphere. Indeed, provided that raw data of radiosondes are available, and by applying a rigorous treatment for the instrumental noise, the huge data base of standard meteorological radiosoundings can be used in order to carry out a climatology of atmospheric turbulent events through the detection of overturns, at least for the most energetic events.

  Fig. 1. (a) The first differences of measured pressure P (cyan) and approximated pressures (black) for thr MUTSI-LR10 flight. (b)The first difference of altitudes z inferred from raw measurements (cyan) and by using the approximated pressures (black).

Fig. 2 .

 2 Fig. 2. Temperature profile (MUTSI-LR10 flight) (left). An expended view is shown for the 4500-5000 m altitude domain. Potential temperature profile (right).

  Figure 2 (left panel) shows the temperature profile of flight LR10 of the MUTSI campaign. The curves show the temperature plus and minus two standard deviations of noise (only visible on the in the expended views in each panel).

Fig. 3 .

 3 Fig.3. Vertical profiles of the estimated level of noise for the MUTSI-LR10 flight of the temperature (left), the pressure (middle), the potential temperature (right). The noise level on potential temperature is inferred from the estimated noise on temperature and pressure. The blue curves show the estimates from data segments, the red curves showing smoothed profiles of these noise estimates obtained by spline approximations (see text).

Fig. 4 .

 4 Fig. 4. MUTSI-LR10 flight: trend-to-noise ratio (left); Thorpe signal (middle); Thorpe lengths for the detected inversions (right). The selected overturns are highlighted.

Fig. 5 .

 5 Fig. 5. Number of occurrences of the size of the detected inversions for all the LR profiles (blue-filled). Number of occurrences of the size of the selected overturns after applying a denoising and a selection procedure (red).

Fig. 6 .

 6 Fig. 6. Empirical distribution of Thorpe displacements within the selected overturns in the troposphere for the six HR flights of the MUTSI and SFT campaigns. The two histograms show the distributions obtained from HR profiles (green-filled) and LR profiles (red-transparent).

Fig. 7 .

 7 Fig. 7. Empirical pdf of the size of overturns for the six HR flights of the MUTSI and SFT campaigns (left). Number of occurrences of the size of overturns for the same six flights (right).

Fig. 8 .

 8 Fig. 8. Empirical pdf of of the size of the overturns for all the flights of the MUTSI and SFT campaigns (6 HR and 20 LR flights).

Fig. 9 .

 9 Fig.9. Same as Fig.6, but for the stratosphere.

Fig. 10 .

 10 Fig.10. Same as Fig.7, but for the stratosphere.

Fig. 11 .

 11 Fig. 11. Same as Fig. 8, but for the stratosphere.

Table 1 .

 1 Number of occurrences and cumulative frequencies (brackets) of the size of selected overturns in the troposphere, for the MUTSI and SFT campaigns (6 HR and 20 LR profiles).

	Sizes (m)	2-4	4-8	8-16	16-32 32-64 64-128 128-256 256-512 512-1024
	HR	245	1307	668	280	131	38	14	5	3
		(0.09) (0.58) (0.82) (0.93) (0.98) (0.992)	(0.997)	(0.999)	(1)
	LR	0	0	0	0	43	48	19	4	4
	(HR flights)	(0)	(0)	(0)	(0)	(0.36)	(0.76)	(0.92)	(0.95)	(0.98)
	LR	0	0	0	0	153	166	78	22	12
	(all flights)	(0)	(0)	(0)	(0)	(0.35)	(0.76)	(0.92)	(0.97)	(0.99)

Table 2 .

 2 Number of occurrences and cumulative frequencies (brackets) of the size of selected overturns in the stratosphere, for the MUTSI and SFT campaigns (6 HR and 20 LR profiles).

	Sizes (m)	2-4	4-8	8-16	16-32	32-64	64-128 128-256 256-512 512-1024
	HR	1960	1517	583	188	36	8	0	0	0
		(0.46) (0.81) (0.95) (0.990) (0.998)	(1)	(1)	(1)	(1)
	LR	0	0	0	33	39	9	4	0	0
	(HR flights)	(0)	(0)	(0)	(0.39)	(0.85)	(0.95)	(1)	(1)	(1)
	LR	0	0	0	419	115	50	16	11	1
	(all flights)	(0)	(0)	(0)	(0.68)	(0.87)	(0.95)	(0.98)	(0.998)	(1)

Atmos.Meas. Tech., 4, 795-804, 2011 www.atmos-meas-tech.net/4/795/2011/

www.atmos-meas-tech.net/4/795/2011/ Atmos. Meas. Tech., 4, 795-804, 2011

Acknowledgements. The authors wish to acknowledge CNES (Centre National d' Études Spatiales) for operating the balloons in Japan and the Research Institute for Sustainable Humanosphere (RISH) for their contribution to the MUTSI campaign.

Edited by: S. Buehler

The publication of this article is financed by CNRS-INSU.


Maryline Chetto
email: maryline.chetto@univ-nantes.fr

Damien Masson
email: d.masson@esiee.fr

Serge Midonnet
email: serge.midonnet@univ-paris-est.fr

Fixed Priority Scheduling Strategies for Ambient Energy-Harvesting Embedded systems

The new generation of embedded systems will have the capability to harvest energy from the environment. The electrical energy which is available to power these devices changes over time and is limited by the size of the energy storage unit such as battery or capacitor and the size of the harvester such as a solar panel. In order to cope with this limitation, the system has to dynamically decide when to be active and when to sleep in order to provide the best quality of service without wasting the harvested energy. In this paper, we study this problem for a uniprocessor architecture where periodic tasks have to execute with deadline constraints according to a preemptive fixed priority rule. We evaluate and compare several scheduling approaches by means of simulation.

I. INTRODUCTION

Up to now, Wireless Sensor Networks (WSN) have used batteries as portable energy source. However, systems continue to become smaller. As less energy is available on board, this leads to a short run-time for the devices. Unfortunately, there are emerging embedded applications where sensors are required to operate for very long durations after they are deployed. Extended life of these electronic devices is of particular importance when they have limited accessibility. Ambient energy harvesting is also known as energy scavenging. It is defined as the process where energy is drawn from the environment and then converted and stored for use in electronics applications. We generally use this terminology for low power and small autonomous devices such as wireless sensor networks, and portable electronic equipments. Consequently, ambient energy harvesting has generated a lot of interest for research. It will permit that autonomous embedded systems be powered perpetually. When compared with energy stored in classical storage units such as batteries, the environment represents an infinite source of available energy. Furthermore, an important advantage is eliminating replacement and periodic recharging of batteries that constitute a major part of service and maintenance. Many environmental sources can be exploited, including thermal, optical, mechanical, fluidic, etc. The selection of energy sources must be considered according to the application characteristics. Self powered sensors for medical implants for health monitoring and embedded sensors in structures such as bridges, buildings for remote condition monitoring are typical new generation embedded systems called Energy-Harvesting (EH) embedded systems. A classical EH system is composed of three parts. The harvester scavenges the energy from ambient surroundings and converts it to usable electrical power. The battery or/and capacitor stores the energy. And the computing system further uses this energy to run software. In EH systems, we have to make the best use of the available energy which is highly dependent on the environment. The energy consumption of the system should be adjusted to maximize the performance instead of minimizing the energy consumption as in classical battery powered systems. A new role of the operating system is to manage correctly the activity of the processing unit such that, at every time there is sufficient energy in the storage unit to satisfy all the constraints. Most environmental energy sources do not deliver a constant power over time. The energy generated may arrive in bursts and has to be stored so that the device can still be operated at a later moment. The main challenge for an EH system lies in the optimization of its performance while respecting the time-varying amount of energy without wasting or exhausting the energy stored in the battery or capacitor. We then say that the system operates in an energy neutral mode by consuming only as much energy as harvested [START_REF] Kansal | Harvesting aware power management for sensor networks[END_REF]. In most applications, the EH system is embedded and has real-time requirements expressed in terms of deadlines attached to the execution of programs, generally called tasks. Then, the performance of such system is measured by the deadline miss rate and heavily depends upon the stored energy and the energy harvested from the environment. Unfortunately, the scavenging power is time-varying and thus very unstable. Therefore, the accurate modeling for energy source plays a key role in designing a good policy to schedule the tasks and reduce the deadline miss rate. So, in that paper, we will describe how to dynamically manage the activity of a single processor system in charge of executing real-time tasks which are periodic and have to be scheduled according to a fixed priority driven policy. We must be able to answer questions like the following: When should the system use energy? When should it be idle and recharge the energy storage? Specifically, we explore different scheduling heuristics and compare their performance. This paper focuses on the same problem studied in [START_REF] Chetto | Real-time scheduling of periodic tasks in a monoprocessor system with rechargeable energy storage[END_REF] but considers variants of the Fixed Priority (FP) scheduler instead of Earliest Deadline First (EDF) scheduler. Simulations illuminate the merits of our heuristics.

The remainder of the paper is organized as follows: We first present the necessary background relative to real-time scheduling with energy harvesting constraints in the remainder of this Section. The model is described in Section II. We give fundamental definitions and introduce scheduling heuristics in Section III. A method for computing the slack time in a fixed priority environment is presented in Section IV. Section V is concerned with an experimental study for comparing the scheduling heuristics with different criteria. Finally, we conclude the paper in Section VI.

Related Work

Only in the past decade, researchers started to address power and scheduling issues with the objective of either minimizing power usage under timing constraints or maximizing the system performance under the energy constraints. Nevertheless, they did not consider the rechargeability of the batteries. For example, the well known Earliest Deadline First scheduler has been extended to variable-voltage processors. The idea is to save power by slowing down the processor just enough to meet the deadlines of the tasks. But solely applying these techniques has limitations in energy harvesting systems because they minimize CPU power, rather than they dynamically manage power according to the profiles of both available energy and processor workload. In the work by Allavena et al. in [START_REF] Allavena | Scheduling of frame-based embedded systems with rechargeable batteries[END_REF], power scavenged by the energy source is constant and all tasks consume energy at a constant rate. Later in [START_REF] Moser | Real-time scheduling for energy harvesting sensor nodes[END_REF], Moser at al. propose LSA (Lazy Scheduling Algorithm) to optimally schedule tasks with deadlines, periodic or not. In that work, the total energy consumption of every task is directly connected to its execution time through the constant power of the processing device. But in a real application, instantaneous power consumed by tasks may vary along time depending on circuitry and devices required by the tasks. Very recently, in [START_REF] Chetto | Real-time scheduling of periodic tasks in a monoprocessor system with rechargeable energy storage[END_REF], we relax the restrictive hypothesis that links together energy requirement and execution time of tasks. We present an on-line scheduling scheme, EDeg (Earliest Deadline with energy guarantee) that is a variant of EDF which relies on two fundamental concepts, namely slack time and slack energy. EDeg dynamically applies a slack-based method for making the processor inactive as long as possible without injuring deadlines feasibility. Performance evaluations have been performed and show the efficiency of this scheduler. However, EDeg is a clairvoyant scheduler that needs knowledge of future task arrival times to dynamically build an optimal schedule, which seriously limits its scope. To the best of our knowledge, all works for scheduling energy-harvesting systems so far consider dynamic priority scheduling based on the famous Earliest Deadline First rule. In this paper, our contribution will be to discuss the same question for fixed priority systems which are in the majority among current real-time embedded applications.

II. MODELS AND ASSUMPTIONS

This paper explores an embedded system which is equipped with an energy harvesting device. Our system model then includes the application model, the energy source model and the energy storage model. a) Application model: We consider here a set of independent periodic tasks that can be denoted as follows:

Φ = {τ i , i = 1, . . . , n}. A four-tuple (C i , E i , D i , T i ) is associated with each τ i .
In this characterization, task τ i makes its initial request at time 0 and its subsequent requests at times kT i , k = 1, 2, ... called release times. The least common multiple of T 1 , T 2 , . . . , T n (called the hyperperiod) is denoted by T LCM . Each request of τ i requires a Worst Case Execution Time (WCET) of C i time units and has a Worst Case Energy Consumption (WCEC) of E i . We assume that the WCEC of a task has no relation with its WCET. A deadline for τ i occurs D i units after each request by which task τ i must have completed its execution. We assume that 0

< C i ≤ D i ≤ T i for each 1 ≤ i ≤ n.
The tasks are scheduled on a single processor system based on a fixed priority driven scheduler.

b) Energy source model: We assume that ambient energy is harvested and converted into electrical power. We cannot control the energy source but we can predict the expected availability with a lower bound on the harvested source power output, namely P r (t)). For sensor nodes deployed in a certain environment, a predictor is used to predict the amount of harvested energy in the future [START_REF] Raghunathan | Design considerations for solar energy harvesting wireless embedded systems[END_REF]. P r (t) is then the instantaneous charging rate that incorporates all losses caused by power conversion and charging process.

c) Energy storage model: Our system uses an energy storage unit that has a nominal capacity, namely C, corresponding to a maximum energy (expressed in Joules or Wattshour). The stored energy, denoted as E s (t) at a given time t, has to remain between two boundaries E min and E max .

Consequently, E = E max -E min . E s (t) can be measured with reasonable accuracy, used at any time later with no leak over time. We assume that energy production times can overlap with the consumption times.

III. FIXED PRIORITY SCHEDULING

Schedulability analysis of periodic task sets can easily be performed both under fixed and dynamic priority assignments. In particular, a lot of work has been done for the Rate Monotonic (RM) and the Earliest Deadline First (EDF) algorithms [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real time environment[END_REF]. RM scheduling is the optimal fixed-priority (or staticpriority) scheduling policy for periodic tasks for the case where task deadlines are coincident with the end of a task's period. That means that if any other fixed-priority scheduling policy can meet deadlines, so can RM. In [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real time environment[END_REF], the authors also derived a sufficient condition for the schedulability of task sets that use a RM priority assignment. Leung and Whitehead in [START_REF] Leung | On the complexity of fixed-priority scheduling of periodic, real-time tasks[END_REF] later showed the optimality of a Deadline Monotonic (DM) priority assignment for the case where periodic tasks have deadlines that are at or before the end of their period. However, these well known analysis methods do not take into account energy constraints. In this example, τ 1 has a priority higher than τ 2 . We can see that there is no solution for producing a valid schedule τ 1 if τ 2 is not paused before it consumes all the energy. However the task set is feasible if τ 2 is paused before is has consumed all the remaining energy. Note that this result also applies to dynamic priority assignments rules such as EDF.

A. Feasible Systems

In this paper, we focus on preemptive Fixed Priority (FP) scheduling strategies such as RM or DM that are applied at run time in an EH system as one described in the previous section. We are interested with online algorithms i.e. schedulers that must make their decisions at run-time. In what follows, we will use the following definitions: A schedule for a set of periodic tasks Φ is valid if the deadlines of all tasks of Φ are met, starting with an energy storage unit fully charged. A periodic task set Φ is timely-feasible if there exists a valid schedule for Φ without considering its energy constraints. A periodic task set Φ is feasible if there exists a valid schedule for Φ. A scheduling algorithm is optimal if it finds a valid schedule whenever one exists.

B. Scheduling Heuristics

We propose here to define and evaluate four scheduling heuristics, all using the same FP rule. The first one is pretty straight forward and natural: it consists to make the system idle for a fixed arbitrary amount of time whenever the energy storage unit becomes empty i.e. at time t where E(t) reaches E min . Our objective is to demonstrate that such a situation may be avoided or at least be delayed. The second heuristic comes from the idea that the current energy level in the battery has to be considered dynamically. In one hand it is not necessary to let the system paused for a long time if a short interval is sufficient to refill a significant fraction of the battery capacity. In the other hand, if the energy level has not sufficiently increased, there is a high probability to quickly reach again the empty state (E(t) = E min ). This intuition is also enforced by the fact that realistic harvester/battery models do not permit a time linear refill strategy. The third heuristic considers that the longer the battery is recharging, the better it is. The problem then amounts to compute the length of the longest interval during which the system can be paused while guaranteeing no deadline miss. That necessarily implies for the scheduler to dynamically determine the slack time at run time. Explanations relative to that computation will be given in Section IV. However, in order to avoid wasting energy if the storage unit fully replenishes during the paused period, we put the processor in idle state only if there is slack time and the energy level is not maximum. We will show that none of these online scheduling heuristics is optimal. In other words, we can always find at least one feasible task set that is not feasibly scheduled by the heuristic. Indeed, we can construct a scenario where all these heuristics fail to produce a valid schedule while such a schedule exists through adequate idling periods. This scenario is represented in Figure 1. For that reason, we propose an additional heuristic using a threshold value for the energy level: whenever the threshold value is reached, we let the processor idle. However, with an arbitrary value for this threshold, system performances cannot be significantly optimized, even in the average case.

IV. SLACK TIME IN FIXED PRIORITY SYSTEMS

We define the slack-time of a task set Φ at time t as the maximum amount of time the system can be suspended from time t without missing any deadline. It is denoted S(t). The slack-time of a task τ i at time t, denoted S i (t), corresponds to the maximum amount of time τ i can be delayed from time t without missing its deadline. This value represents the unused CPU time units at priorities higher than or equal to i. We have S(t) = min ∀i S i (t). These notions were introduced in [START_REF] Lehoczky | An optimal algorithm for scheduling soft-aperiodic tasks fixed priority preemptive systems[END_REF] in order to address the issue of scheduling mixed traffic composed of hard real-time periodic tasks and soft real-time aperiodic ones. The slack-time computation then officiates as an on-line admission control test for aperiodic traffic. Algorithms to dynamically compute the exact value and upper bounds in fixed priority driven systems are given in [START_REF] Davis | Scheduling slack time in fixed priority pre-emptive systems[END_REF], [START_REF] Davis | On exploiting spare capacity in hard real-time systems[END_REF], [START_REF] Masson | Userland Approximate Slack Stealer with Low Time Complexity[END_REF].

A. Processing exact slack time S i (t)

To calculate S i (t), we consider the system between t and the next deadline of the i-level priority task (τ i ). During this interval the system can be viewed as a succession of ilevel busy periods (periods where the processor is servicing priorities higher or equal to i) and i-level idle periods (periods where the processor is idle or periods where processor is servicing priorities lower than i). The slack time value S i (t) is then the sum of the i-level idle period lengths. To process the exact slack value we use the two equations to compute the end of a busy period starting at instant t and to compute the length of an idle period starting at time t. The process steps can be summarized as follows: 1) initializing the slack value S i (t) at time 0, 2) processing w i (t), the length of the busy period starting at time t, 3) processing vi(t, w i (t)), the length of the idle period starting at time t + w i (t), 4) incrementing the slack value by the length of the idle period, 5) restarting the process at step 2 or stopping it if the deadline of the task τ i is reached. This algorithm could not be processed at any time instant because of its computational complexity. So, in order to simplify its dynamic computation, the slacktime value is obtained at time t from the slack-time value computed at time t. Two general cases have to be studied: none of the hard periodic tasks ends its execution in [t, t ) or one periodic hard real-time task (τ i ) ends its execution at time t ∈ [t, t ). In the first cas, then there are two possibilities for the processor activity between t and : the processor is is executing hard periodic tasks. In the first case, the slack is reduced by (t -t) for all priorities. However, if the processor is executing the hard real-time task τ i , then the system is idle for higher priorities (k < i), and the slack is reduced by (t -t) only for these priorities. In the second case, then all i-level idle times present in [t, d i (t)) will be present in [t, d i (t) + T i ) = [t, d i (t )), which is the new interval to consider for the i-level slack times computation. Therefore, the τ i termination can only increase S i (t) but never decreases it. Consequently, S i (t) has to be recomputed each time τ i ends a periodic activation. So, assuming that there is a time t where the S i (t) was up to date for all tasks, the algorithm to compute S(t ) is :

1) if none of the periodic hard real-time tasks ends in [t, t ) a) if the processor is idle or executing soft aperiodic requests

∀j : S j (t ) = S j (t) -(t -t) (1)
b) if the processor is executing hard periodic task τ i ∀j < i :

S j (t ) = S j (t) -(t -t) (2)
2) if hard real-time task τ i ends at time t ∈ [t, t ), S i (t ) has to be computed using the recursive analysis described at the beginning of this section.

V. EXPERIMENTAL RESULTS

We developed a discrete event simulator to evaluate the effectiveness of the scheduling heuristics previously introduced, all using a simple and easy to implement fixed priority rule. For all schedulers, we consider the system as failed as soon as one deadline is missed. The simulator is available under the GNU Public License at http://dajam.fr/Softwares/, with sufficient material for reproducing the experiments.

A. Simulation restrictions

In our simulations, we assume that P r (t) is constant along time. Moreover, we assume that all tasks linearly consume their energy budget over time, and their consumption rate, given by E i /C i , is always greater than P r (t) (every task execution leads to discharge the battery).

B. Formal definition of Scheduling strategies

Five Energy Harvesting in Fixed Priority (EHFP) scheduling policies have been evaluated:

1) EHF P 1 : All tasks are processed as soon as possible according to the Fixed Priority (FP) rule until the storage unit is empty or there are no more task ready to be executed. If the storage unit is empty and there is at least one ready task, the processor is put into sleep mode during x units of time where x is an input of the scheduler. During that period, the energy storage unit replenishes. Default value for x is 1.

2) EHF P 2 : All tasks execute as soon as possible according to FP. When the storage unit becomes empty, the processor is put into sleep mode until the energy level reaches a threshold value, namely E th , given as an input of the scheduler.

3) EHF P 3 : All tasks execute as soon as possible according to FP. When the storage unit becomes empty, the processor is put into sleep mode as long as the slack time remains positive.

4) EHF P 4 : All tasks execute as soon as possible according to FP. When the storage unit becomes empty, the processor is put into sleep mode as long as the slack time remains positive and the energy level remains less than E max .

5) EHF P 5 : There are two threshold parameters, namely E thmin and E thmax . All tasks execute as soon as possible according to FP. When the energy level reaches E thmin , the system is put into sleep mode as long as the slack time remains positive and the energy level remains less than E thmax .

C. Simulation methodology 1) Task sets generation: Groups of periodic task sets are generated with utilization factors respectively equal to 30, 50, 70 and 90%. Results presented in this section are averages over groups of ten task sets. Each task set contains six realtime periodic tasks. Periods are randomly generated with an exponential distribution in the range [40-2560] time units. Costs are randomly generated with an uniform distribution in the range [1-period]. In order to test task sets with deadlines less than periods, we also randomly generate deadlines with an exponential distribution in the range [cost-period]. Nevertheless, these results are not presented here since they do not bring any additional conclusions. Priorities are assigned to tasks according the Deadline Monotonic policy and solely timely-feasible task sets are evaluated.

2) Energy Parameters Generation: Energy consumption of tasks (WCEC) are randomly generated, but constrained by E max and P r values. Simulations have been performed with different values for E max and P r . Energy consumption per time unit of each task, E i /C i , is generated with an uniform distribution in the range (P r , E max /3]. E i is then obtained by multiplying this value by C i .

3) Simulations Description: Every simulation starts with a fully charged energy unit i.e. E(0) = E max , in the worst case timing scenario that corresponds to the synchronous activation of all periodic tasks at time t = 0. Every simulation terminates when either a deadline is missed or time t = 100000 is reached. The different measures are: 1) the number of feasible task sets (on the interval [0-100000]), 2) the average finishing date of the schedule corresponding to the first deadline violation whenever one exists, 3) the average time spent with maximum energy level in the storage i.e. E(t) = E max , 4) the average time spent with minimum energy level in the storage i.e. E(t) = E min , 5) the average duration of busy-periods (corresponding to discharging phases of the energy storage) and 6) the average duration of idle-periods (corresponding to recharging phases of the energy storage or energy storage fully charged) Finally, we compare these measures for different scenarios. We evaluate the impact of increasing/decreasing E max i.e. the capacity of the energy storage unit and P r . For policies EHF P 1 and EHF P 2 , we also evaluate the impact of varying the input parameter of theses heuristics. In both cases, we observe a correlation between the number of failures (in other terms the ratio of non-feasibly scheduled task sets) and the length of the simulation. Moreover, we observe that the performances are constant until x reaches a value (10 for 2(a) and up to 5 for 2(b)). There after this threshold, the pause time is too big and a higher number of deadlines are missed.

Conclusion: This set of simulations enables us to conclude that the length of the pause intervals has an impact on system performance. A value in a reasonable small range permits to increase the mean system life time, but there is no optimal value. 5) Effects of parameter E th on EHF P 2 : Simulations are performed first with E max = 100 and P r = 10, and second with E max = 200 and P r = 40. Results are respectively presented on Figure 3(a) and Figure 3(b). First, we observe that in the worst and best cases, EHF P 2 behaves similarly to EHF P 1 . As for the x parameter of EHF P 1 , the threshold parameter attached to EHF P 2 gives the same ratio of performances until a critical value (10% for 3(a) and 20% for 3(b)). Identically to EHF P 1 we cannot draw an optimal value for that parameter.

Conclusion:

We demonstrate through these simulations that basing the paused time duration on a criteria that depends on the battery size permits to enlarge the range of values that offer the best average performances. However, there is no optimal value for this threshold in the general case, and EHF P 2 does not provide a better performance than EHF P 1 , in terms of deadline missings.

D. Comparison with fixed E max and P r

Simulations are now performed by making vary E max and P r . We pay attention to compare the different heuristics when the energy constraints vary. We set x = 6 for EHF P 1 and E th = 10% for EHF P 2 and EHF P 5 . Simulation results are presented on Figure 4. Complete and detailed results can be consulted on-line at http://dajam.fr/greenCom2011. Due to its very bad performance, policy EHF P 3 is eliminated from this section. That demonstrates that a policy, for providing acceptable performance, has to take into account the maximal upper bound for the battery charge. In consequence EHF P 3 is not competitive.

1) Performance evaluation in terms of feasibility: One crucial criteria for measuring the performance of any realtime scheduling strategy is the ratio of feasible task sets. Experiments empirically demonstrate that, using this criteria, EHF P 4 outperforms the other policies in most cases. However, even if this ratio is never less than for other policies, task sets which fail may not be the same. Moreover the gain is generally quite small. It is highly probable that no online scheduler can be optimal. In other terms, only a clairvoyant algorithm (i.e. an algorithm which knows in advance all the future) can produce a valid schedule whenever one exists.

2) Performance evaluation in terms of busy/idle period duration: Every schedule can be characterized by the average number of idle time periods and so the average duration of idle time periods. This criteria is of particular importance in systems using a DPM (Dynamic Power Management) mechanism. Many architectures provide the equivalent of a halt instruction that reduces CPU power during idle periods. However, there are usually significant latencies and overheads for entering and exiting these states. Consequently, the number of state switchings will influence the resulting efficiency of the system. Higher will be the average duration of the idle time periods, lower will be the energy and time overheads incurred by implementation of a DPM mechanism. Moreover, this number will have an impact on the lifetime of the energy storage unit, in general a battery or capacitor. Charging a battery is not a linear process and longer it is paused, the more energy it harvests. As for the previous performance criteria, we cannot conclude on the existence of an optimal strategy regarding the average number of state switching. No policy appears as the best one amongst the heuristics under evaluation. Nevertheless, we observe that, in comparison with other heuristics and in most cases, EHF P 4 maximizes the average idle/busy period duration. Only in very few cases, EHF P 4 is outperformed by an other heuristic with a non significant gap.

VI. CONCLUSION

Ambient energy harvesting can provide an extended lifespan and support to conventional electronics systems including embedded systems with real-time constraints. This paper has explored how to optimize the performance of these systems which are provided with a preemptive fixed priority scheduler. In our model, there is no restrictive assumption on the profile of the energy source as well as tasks can consume energy with variable power. We show that the conventional greedy algorithm called FP do not work very well. We propose several variants of FP to derive more efficient solutions. The main conclusion that emerges from our experiments is the highly probability of the non existence for an optimal algorithm if the recharging rate is not taken into account both to decide when to stop the system and for how-long. However we proposed an heuristic, EHF P 4 , simply named slack-time based heuristic which offers pretty good average performance in a very large range of situations. Nevertheless, in the context of energy and time constrained systems, constructing a valid schedule for every feasible task set requires clairvoyance and cannot be performed with an online scheduling algorithm. Future work will then consist in proposing an optimal clairvoyant scheduling algorithm. Our objective is to adapt to fixed priority systems, the concept of slack energy initially proposed in [START_REF] Chetto | Real-time scheduling of periodic tasks in a monoprocessor system with rechargeable energy storage[END_REF] for a dynamic priority driven system based on EDF scheduling. How to dynamically compute the slack energy at run time will be the key point of the scheduler as well as proving its optimality.

Fig. 1 .

 1 Fig.1. Non optimality of strategies that do not anticipate energy shortage. In this example, τ 1 has a priority higher than τ 2 . We can see that there is no solution for producing a valid schedule τ 1 if τ 2 is not paused before it consumes all the energy. However the task set is feasible if τ 2 is paused before is has consumed all the remaining energy. Note that this result also applies to dynamic priority assignments rules such as EDF.

  (a) for Emax = 100 and Pr Emax = 200 and Pr = 40

Fig. 2 .Fig. 3 . 4 )

 234 Fig. 2. Effects of parameter x on EHF P 1

  (a) Emax = 200 and Pr = (b) Emax = 200 15 (c) Emax = 200 and Pr = 20 (d) Emax = 200 and Pr = 25 (e) Emax = 200 and Pr = 30 (f) Emax = 200 and Pr = 35 (g) Emax = 400 and Pr = 10 (h) Emax = 400 and Pr = 15 (i) Emax = 400 and Pr = 20 (j) Emax = 400 and Pr = 25 (k) Emax = 400 and Pr = 30 (l) Emax = 400 and Pr = 35

Fig. 4 .

 4 Fig. 4. Performance evaluations


Jacques Le ´vy

V E ´hel

Khalid Daoudi
email: daoudi@bora.inria.fr

Generalized IFS for Signal Processing

Introduction

Several methods have been proposed to estimate the Ho ¨lder exponents of signals [START_REF] Mallat | Characterization of Signals from MultiScale Edges[END_REF][START_REF] Gonzalve | Bilinear Time-Scale Analysis Applied to Local Scaling Exponents Estimation[END_REF]. In this paper, we propose a new approach, based on a generalization of iterated functions system (IFS), which is well adapted to irregular continuous 1D signals. We also use these generalized IFS to build parsimonious models of complex signals and to perform segmentation on them. This paper is organized as follows : in section 2 we recall the definition of generalized iterated functions systems (GIFS) and the definition of the Ho ¨lder exponent of a nowhere differentiable continuous function. In section 3, we recall a result obtained in [START_REF] Daoudi | Construction of continuous functions with prescribed local regularity[END_REF] concerning the Ho ¨lder exponents of the attractors of GIFS. We then propose a method to solve the inverse problem for GIFS and we use it to estimate the Ho ¨lder exponents of discrete signals. We present results obtained on generalized Weierstrass functions. In section 4, we develop a synthesis scheme for a given signal. The parameters of the GIFS, along with the Ho ¨lder exponents at different resolutions of the signal, allow to give a functional representation of discrete data. This in turn permits signal segmentation. An application on a residual speech signal is presented.

Background on generalized IFS

Consider a collection of sets IN , where each is a non-empty finite set of contractions on a 0 ; 1 , for 0 1, where 1 is an integer, which denotes the cardinal of . We denote by the contraction ratio of , for 0 1, and IN . Consider the operator , defined on every compact 0 ; 1 , by : 1 0 Then, under some conditions [START_REF] Andersson | Recursive Construction of Fractals[END_REF], there exists a unique compact such that : In the following we consider the special case of IFS where the 's are affine functions and the number of functions at step is , for some integer 2. Let be defined as the set of affine transformations 0 represented in matrix notation by: 1 0

We suppose 0 1 , 1 1 and that each is operating only on 1 ; 1 1 if is even, and on 1 1 ; 1 if is odd. Every maps to ; 1 . Suppose, also, that we want to interpolate the points , for 0 , and IR. Then, it was shown in [START_REF] Daoudi | Construction of continuous functions with prescribed local regularity[END_REF] that if is function from 0 ; 1 to 0 ; 1 , which is a lower limit of a sequence of continuous functions, there exists 's such that the attractor of the IFS defined above is the graph of a continuous function that verifies : of non negative integers strictly smaller than , are uniquely determined by : if

1 1 else if 0 1 else
This yields a non parametric method to estimate the Ho ¨lder function of a signal. Indeed, if for a signal, we are able to determine a GIFS whose attractor is the graph of the signal (inverse problem), then, using proposition 1, we get an expression of the Ho ¨lder exponent at every location of the signal. However, there still remain some limitations. The 's belong to 1 ; 1 , which implies that the Ho ¨lder exponents computed with this method belong to 0 ; 1 . This means that the attractor of a GIFS is necessarily the graph of a nowhere differentiable continuous function. The scheme we propose below to solve the inverse problem thus yields better results on irregular signals. Let 0 2 be a given discrete signal. For 1 , consider the set defined by : 2 0 2

Namely, is the signal sub-sampled with a step equal to 2 . For a given 1 1 , the set 0 2 1 is simply obtained as the set of contraction ratios of the 2 affine functions that permit to map the polygon defined by to the one defined by 1 . To clarify, we show in Fig. 5 an example with 3. The signal samples are represented by circles. For 1 2, the GIFS coefficients are obtained as follows : . This function is a generalization of the Weierstrass function and was shown [START_REF] Daoudi | Construction of continuous functions with prescribed local regularity[END_REF] to have Ho ¨lder function . Fig. 2 displays the corresponding estimated Ho ¨lder function using our method and a method based on timefrequency energy distributions [START_REF] Gonzalve | Bilinear Time-Scale Analysis Applied to Local Scaling Exponents Estimation[END_REF]. The distribution used here is a continuous wavelet (Morlet wavelet) transform. The theoretical Ho ¨lder function is the thick line, the GIFS estimation is the thin line and the continuous wavelet transform estimation is the dotted line. The results obtained using our method, in this case, are satisfactory. Other examples and comparison to other methods will be shown in the extended version.

Signal modeling using GIFS

In this section, we show how the 's obtained by solving the inverse problem corresponding to can be used to develop a parsimonious synthesis model of the data. The idea is to keep only some of the 's, which correspond in the time domain to some sub-sampling of the signal. We then determine the missing 's from the Ho ¨lder function of the signal at different scales. Remark that, for a given 1 1 , the set ; 1 0 2 1 determines the GIFS whose attractor at scale is 1 , when starting with 1 . Let 0 be a fixed integer, and suppose that we want to keep only ; 1 0 0 2 1 . This corresponds to keeping one sample out of 2 0 samples. Then, for every 0 1 0 we compute the Ho ¨lder function of the attractor at scale . We then perform a moving average on each in order to obtain a smooth function which can be defined with a few parameters. Then, we put :

2 2 0 2 1
This is justified by the fact that, if is continuous function from 0 ; 1 to 0 ; 1 , then a GIFS whose contraction ratios are of the form : 2 2 has an attractor which is the graph of a continuous function with Ho ¨lder exponent at each point equal to . Hence, we build a signal model which is all the more parsimonious than 0 is large. The choice of 0 depends on the desired synthesis quality. Fig. 3 shows a residual speech signal obtained by subtraction of the harmonic part to the original speech signal. The harmonic part was extracted using the HNS model [START_REF] Laroche | HNS: Speech Modification Based on a Harmonic + Noise Model[END_REF] developed by J. Laroche, Y. Stilyanou and E. Moulines. Fig. 4 shows the reconstructed residual signal using our synthesis scheme with 0 2. The method can be generalized to perform signal segmentation based on a split and merge scheme. The associated energy is a function of the 's. This will be detailed in the extended version of the paper.

  call the attractor of the GIFS IN and 1 its attractor at scale (which depends on ).

Proposition 1

 1 ¨lder exponent of at the point , and is defined as follows : lim inf 0 log log In the rest of the paper we will call the Ho ¨lder function of . The attractor of the IFS defined above is the graph of a continuous function such that : and where, for any positive integer , the -uples 1 and 1

  Fig. 1 displays the function defined on 0 ; 1 by :

Estimation of Ho ¨lder exponents using GIFSIn this section, we develop a technique to estimate the Ho ¨lder exponents of discrete signals. This technique requires an investigation on the inverse problem of GIFS. For this purpose, we need to recall a result proved in[START_REF] Daoudi | Construction of continuous functions with prescribed local regularity[END_REF] :


Multiple threshold voltage for glitch power reduction

Mariem SLIMANI,Philippe Matherat

Institut TELECOM, TELECOM-ParisTech, LTCI-CNRS 46, rue Barrault, 75634 -Paris Cedex 13, France Email: {mariem.slimani, philippe.matherat}@telecom-paristech.fr Abstract-We address the problem of circuit-level design for low power. We describe a new method for glitch power reduction based on threshold voltage adjustment. The proposed method achieves both dynamic and leakage power reductions. We develop an optimization algorithm that transforms the circuit netlist in an optimized one achieving glitch energy reductions without affecting the overall circuit delay requirement. Applying the algorithm to C17 benchmark circuit implemented in a 65nm industrial Low Power CMOS process, we have achieved 14% of total energy savings and 78% of leakage energy savings at the expense of just 5% of delay increase.

Index Terms-Glitch power reduction; Circuit-level Design; threshold voltage variation

I. INTRODUCTION

Glitches are spurious transitions generated by a logic gate when its input signals arrive at a different time due to paths of varying delay. These transitions consume an important amount of energy that can be as much as 40% of the total power in a standard synchronous circuit. [START_REF] Raja | Variable input delay CMOS logic for low power design[END_REF]. Many techniques are proposed to avoid glitches. All of them try to compensate the differential delay of input signals either by adding extra delay to early arriving signals or by increasing the inertial delay of the corresponding gates such that the glitch will be filtered (hazard filtering) [START_REF]Low-power design by hazard filtering[END_REF]. The first method is called path delay balancing. In [START_REF] Kim | New path balancing algorithm for glitch power reduction[END_REF], this technique inserts delay buffers in the path in which extra delay is desired. In [START_REF] Raja | Variable input delay CMOS logic for low power design[END_REF], authors add an extra resistance to introduce the extra delay on selected paths. For the hazard filtering method, many ideas emerged to increase gate's delay. In [START_REF] Wroblewski | Automated transistor sizing algorithm for minimizing spurious switching activities in CMOS circuits[END_REF], enlarging transistors was proposed as a solution to achieve equal path delays in the circuit. But, since large transistors means high values of parasitic capacitances, a multiobjective algorithm is developed to optimize both differential path delay and capacitive power consumption by a judicious choice of transistor's size. Like transistor sizing, gate sizing technique searches for the optimal size of the gate [START_REF] Raja | Variable input delay CMOS logic for low power design[END_REF], modeled as an equivalent inverter. All techniques above mentionned come with some drawbacks: In path delay balancing, since delays are added only to the fast paths, critical path is not affected. However, extra elements inserted consume themselves energy and hence reduce the energy saved. In hazard filtering technique, the overall delay of the circuit may increase. So, critical path should be always verified to still meets the delay constraints. Furthermore, techniques belonging to this method, like transistor or gate sizing, consume an extra capacitive power that reduces the achievable power savings. Adding delay buffers or enlarging transistors will add an extra area to the circuit. But, these methods are especially used for low power applications with low demand in area requirements. In this paper, we develop a new technique for hazard filtering based on threshold voltage variation approach. Recently, this approach has attracted a strong interest for low power design. To our knowledge, no previous work for glitch power reduction has adopted such approach. Techniques in this approach are used especially to achieve leakage power reduction [START_REF] Mudge | Combined dynamic voltage scaling and adaptive body biasing for lower power microprocessors under dynamic workloads[END_REF] [START_REF] Agarwal | Leakage in nanoscale technologies: mechanisms, impact and design considerations[END_REF]. As threshold voltage increases, sub-threshold currents decrease with an increase in the propagation delay of the gate. Thus, to reduce the power consumed while maintaining the performance of the circuit, high threshold voltage are generally used in transistors off-critical path or in standby mode [START_REF] Wei | Design and optimization of dual-threshold circuits for low-voltage low-power applications[END_REF] [8]. It's the supposed negative point of high threshold voltage that motivate this work. The increase of the propagation delay due to threshold voltage increase will be exploited to eliminate spurious transitions. Threshold voltage of gates with eventual glitch appearance will be adjusted to compensate different path delays of its input signals. Using the adequate threshold voltage in the adequate place, our technique reduces both glitch power consumption and leakage power dissipation while still maintaining the delay constraints. This paper is organized as follows. Section II briefly reviews glitch filtering technique. In Section III, the basic idea of multiple threshold voltage for glitch filtering is presented and techniques to achieve threshold voltage variation are analyzed. Finally, section IV presents the proposed algorithm and provides simulations results.

II. GLITCH FILTERING APPROACH

Glitch filtering approach was first proposed in [START_REF]Low-power design by hazard filtering[END_REF]. The basic idea is to adjust the inertial delay of the gate, where glitch can be generated, to be greater than the differential path delay (DPD) of its input signals. Figure 1 illustrates this technique. The principle is simple : the transition of input signal A to logic state 1, will start the charging of output signal C. When the delay of the gate is smaller than DPD, the output will found enough time to be completely charged before being dicharged due to the change of input signal B. Suppose now that the delay of the gate is increased by a factor DPD. Normally, the output will be at logic state 1 at time t1+D+DPD=t2+D. But, the change of input signal B to logic state 0, will start the discharge phase of output signal at time t2. Since the output is not being charged yet, the glitch will be filtered and the output remains at logic state 0 all the time. This technique may increase the overall delay of the circuit. For instance, consider the AND gate, with two rising transitions in its inputs. The output will be at the desired logic state after t+D, where t is the latest arrival time of input signals and D is the delay of the gate. Obviously, if the delay of the gate increases, the output will be delayed by a factor corresponding to the delay increase.

III. HAZARD FILTERING BY THRESHOLD VOLTAGE

VARIATION

In this paper, we develop a new technique for hazard filtering based on threshold voltage variation. The idea is by reason to the relation between the delay of a gate and threshold voltage that can be expressed as [START_REF] Luo | A novel gatelevel NBTI delay degradation model with stacking effect[END_REF]:

t pd = C L V dd I d = C L V dd L ef f µC ox W ef f (V gs -V th ) α (1)
where V gs and V th are respectively the gate-source voltage and the threshold voltage of the transistor, µ is the mobility, C ox is the oxide capacitance and α is the velocity saturation index. L ef f and W ef f are respectively the channel length and width of the transistor and C L is the equivalent capacitance of internal parasitic capacitances that consume the same power. The variation of the gate propagation delay resulting from a threshold voltage variation can be given in Equation 2 by diffrentiating Equation 1 with respect to V th .

∆t pd = α∆V th V gs -V th .t pd0 (2)
Where t pd0 is the original gate delay without V th variation. The idea of our technique exploits this relation. It consists to adjust the threshold voltage in order to increase the delay of the desired gate by an optimal factor that achieve glitch energy savings while maintaining a certain performance requirement on the overall circuit delay.

To achieve multiple threshold voltage, different techniques can be used. In [START_REF] Wei | Low voltage low power CMOS design techniques for deep submicron ICs[END_REF], authors discuss these techniques. Some of them are integrated in the fabrication process like ion implantation or oxide thicknesses (t ox ) variation. The others, like multiple channel length setting or adaptive body biasing technique can be used by the designer in the transistor level.

Transistor length variation increases gate's capacitances which induce more power consumption. Thus, in this paper, we will use adaptive body biasing technique to achieve high threshold voltage. It's a well known technique used for leakage energy reduction [START_REF] Agarwal | Leakage in nanoscale technologies: mechanisms, impact and design considerations[END_REF]. This technique consist to apply a voltage between the source and the bulk of the transistor that influence threshold voltage value. In fact, in BSIM4 model [START_REF] Xi | BSIM4. 3.0 MOSFET Model User; s Manual[END_REF], the threshold voltage of a MOSFET can be expressed as follows:

V th = V th0 + γ 2φ F + V sb -2φ f (3)
Where V th0 is the threshold voltage at zero substrate bias, φ F is the Fermi potential, γ is the body bias coefficient and V sb is the source-bulk voltage.

A variation in the source-bulk voltage results in a threshold voltage change by a factor of : γ √ 2φ F + V sb -2φ f . This is especially used to increase threshold voltage by reverse biasing source to bulk junction [START_REF] Layton | Low-Voltage Analog CMOS Architectures and Design Methods[END_REF]. Since this technique require a triple well process, we are aware that achieved energy savings will be at the expense of a certain area increase. Figure 2 shows the variation of threshold voltage for a standard and a high threshold voltage low power (SVTLP and HVTLP) NMOS transistor resulting from a V bs variation. As we can see, a reverse bias of about 0.5V results in an increase of threshold voltage by about 90 mV. Since transistor models don't take into account a V bs variation beyond V dd + 10%, a maximum reverse biasing of 1.32V is allowed. A marge of 0.27V raging from 0.58V to 0.85V is possible.

IV. OPTIMIZATION PROCEDURE FOR GLITCH POWER

REDUCTION

In this section a motivation example for glitch filtering applying adaptive body biasing technique is presented. The proposed algorithm is then described and experimental results are shown.

A. Motivation Example

Consider the cicruit with a glitch formation shown in figure 3. Glitches appear due to the difference of input arrival times of the AND gate. As mentioned in section II, To eliminate glitches, threshold voltage will be adjusted by reverse biasing source to bulk junction. Spice simulations are carried out in 65nm process. Figure 4 shows the output of the AND gate for different threshold voltages achieved using SVTLP and HVTLP transistors with and without substrate biasing. Here, the V bs represent the value of the bulk-source voltage of PMOS transistors and the source-bulk voltage of NMOS ones. One can observe that, as expected, the magnitude of glitch is decreasing with threshold voltage increase. Table I compares the original circuit and the modified one where HVTLP transistors and a substrate bias of 1.2V is used. We conclude that a total and leakage energy savings of up to 40% can be achieved at the expense of delay increase. For applications whose speed requirements are not stringent, like biomedical devices, sensor networks and radio frequency identification (RFID) [START_REF] Bol | Interests and limitations of technology scaling for subthreshold logic[END_REF], delay of the circuit can be further and further increased as long as it achieves energy savings. But this is not generally the case : given as a constraint parameter, delay has been always a major concern for digital circuit designers.

B. Proposed algorithm

We develop an algorithm that modify the netlist of a logic circuit to have a glitch free one that still maintain a given circuit delay performance. The proposed algorithm makes a bottom up analysis of the circuit : it starts at the primary inputs to finish at the output nodes. It effects a threshold voltage adjustement that may increase the overall circuit delay by a permitted amount calculated as the difference of the delay required and the critical path delay of the circuit without any modification. Effectively, this algorithm searches for an eventual glitch appearance by examining the differential path delay (DPD) of input signals of each gate in the circuit. For each glitchy gate, the highest threshold voltage achieved by using High Threshold Voltage Low Power (HVTLP) transistors with a substrate bias of 1.2V, is used. Then, if the critical path delay still maintain the delay requirements, the algorithm moves to the next gate. Else, the threshold voltage is regulated until the delay constraints are conserved. The glitch may not be totally filtered, but its magnitude will be reduced which reduce the energy consumed as shown in section IV-A . For the case of a large DPD such that the glitch will be a complete transition with a width W, increasing threshold voltage will not reduce the energy consumed until W will be 0. Even so, the proposed algorithm achieve energy savings since high threshold voltage will significantly reduce leakage currents. As using high threshold voltage reduce the leakage energy consumed, the circuit can be further optimized using the algorithm for non-glitchy gates. To apply our algorithm the arrival times of input signals should be computed for all the gates of the circuit. The proposed algorithm is given in figure 5.

C. Experimental Results

The proposed algorithm has been tested on C17 and C432 ISCAS'85 benchmark circuits for different values of delay constraints.

Figure 1 .

 1 Figure 1. Glitch filtering approach.

Figure 2 .

 2 Figure 2.Threshold voltage variation for SVTLP and HVTLP NMOS transistors as function of V sb .

Figure 3 .

 3 Figure 3. Circuit with glitch formation.

Figure 4 .

 4 Figure 4. Output voltage waveforms for different threshold voltages.

Table I ENERGY

 I REDUCTION AND DELAY PENALTY OF PROPOSED GLITCH REDUCTION METHOD APPLIED TO THE CIRCUIT IN FIGURE 3.

		original	modified
	Total-Energy(fJ)	4.76	2.85
	Leakage-Energy(fJ)	7.99e-3	4.61e-3
	Delay(ps)	74.8	107.8
	Total-Energy-reduction(%)	x	40
	Leakage-Energy-reduction(%)	x	42
	Delay-increase (%)	x	44
	Energy-Delay-Product(psxfJ)	356.04	307.23

The results show a significant energy savings after the algorithm has been applied. We observe that for C17 circuit, we can even achieve energy savings without any increase in the circuit delay.

As expected the less the delay constraints are stringent the more the energy savings we can achieve. Unlike others glitch reduction techniques, our technique do not adds an extra energy that reduce the achievable energy savings.

V. CONCLUSION

We describe an optimization technique for glitch power reduction based on threshold voltage variation. The idea is to adjust threshold voltage of specific transistors of gates with eventual glitch appearance in order to eliminate spurious transitions and hence achieve energy savings. We develop an algorithm that modifies a circuit netlist to have a glitch free one that still maintain the delay performance required. A total energy savings of 14% and a leakage energy savings of 78% at the expense of just 5% of delay increase is achieved applying our method to C17 benchmark circuit.