Saeed Mian Qaisar
email: saeed.mian-qaisar@imag.fr

Laurent Fesquet
email: laurent.fesquet@imag.fr

Marc Renaudin
email: marc.renaudin@tiempo-ic.com

Effective Resolution of an Adaptive Rate ADC

Most real life signals are of non-stationary nature. An efficient acquisition of such signals can be achieved by adapting the acquisition rate according to the input signal local characteristics. In this context, an ARADC (Adaptive Rate Analog to Digital Converter), based on the level crossing sampling is proposed. The ADC effective resolution is a classical parameter to characterize its performance. In this context, a novel method is devised to measure the ARADC resolution. A criterion for properly choosing the different system parameters in the aim of acquiring the desired effective resolution is also described.

Introduction

The signal acquisition process dictates the performance of the complete signal processing chain. A smarter signal acquisition results into an efficient system and vice versa [START_REF] Bilinskis | Digital alias free signal processing[END_REF]. The classical ADCs are time-invariant, hence are parameterized by taking into account the worst possible case for the considered application. They capture the input signal at a constant rate. It causes an increased number of samples to be processed, especially in the case of low activity sporadic signals [START_REF] Mark | A nonuniform sampling approach to data compression[END_REF][START_REF] Sayiner | A Level-Crossing Sampling Scheme for A/D Conversion[END_REF][START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF][START_REF] Qaisar | Computationally efficient adaptive rate sampling and filtering[END_REF][START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF][START_REF] Qaisar | Computationally efficient adaptive resolution short-time Fourier transform[END_REF]. This shortcoming is resolved up to a certain extent by employing the ARADC in various applications [START_REF] Qaisar | Computationally efficient adaptive rate sampling and filtering[END_REF][START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF][START_REF] Qaisar | Computationally efficient adaptive resolution short-time Fourier transform[END_REF]. The motivation behind the ARADC is to achieve a smart A/D conversion. The idea is to acquire only the relevant signal parts and to adapt the sampling rate according to the input signal local variations. An efficient solution is proposed by smartly combining the features of both the uniform and the non-uniform signal processing tools.

The ARADC

The block diagram of the ARADC is shown in Figure 1. In the ARADC, a band-limited input signal is acquired by a LCADC (Level Crossing based ADC) [START_REF] Sayiner | A Level-Crossing Sampling Scheme for A/D Conversion[END_REF][START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF].

According to [START_REF] Mark | A nonuniform sampling approach to data compression[END_REF], the sampling instants of a nonuniformly sampled signal obtained with the LCADC are defined by Equation 1. Where t n is the current sampling instant, t n-1 is the previous one and dt n is the time delay between the current and the previous sampling instants.

n n n dt t t + = -1
Eq. 1

LCADCs drastically reduce the activity of the post processing chain, because they only capture the relevant information [START_REF] Sayiner | A Level-Crossing Sampling Scheme for A/D Conversion[END_REF][START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF][START_REF] Qaisar | Computationally efficient adaptive rate sampling and filtering[END_REF][START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF][START_REF] Qaisar | Computationally efficient adaptive resolution short-time Fourier transform[END_REF]. Let ∆V in and ∆x(t) be the LCADC and the input signal x(t) amplitude dynamics respectively. In order to avail the complete LCADC resolution in the studied case, ∆x(t) is always adapted to match ∆V in . For a M-bit LCADC, the maximum and the minimum sampling frequencies are defined by Equations 2 and 3 respectively [START_REF] Qaisar | Computationally efficient adaptive rate sampling and filtering[END_REF]. Where, Fs max and Fs min are the maximum and the minimum sampling frequencies. f max is the highest and f min is the lowest frequency of x(t).

( )

1 2 . . 2 max max - = M f Fs Eq. 2 ( ) 1 2 . . 2 min min - = M f Fs Eq. 3
LCADCs deliver a non-uniform time repartitioned output. The non-uniformly sampled signal obtained with LCADCs can be used for further non-uniform digital processing [START_REF] Qaisar | Computationally efficient adaptive rate sampling and filtering[END_REF]. However in the studied case, the nonuniformity of the sampling process, which yields information on the signal local features, is employed to select only the relevant signal parts. Furthermore, the characteristics of each signal selected part are analyzed and are employed later on to adapt the system parameters accordingly. This selection and local-features extraction process is performed by the ASA (Activity Selection Algorithm) or the EASA (Enhanced Activity Selection Algorithm). They split the signal into a series of active windows. A detailed description of these algorithms is given in [START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF][START_REF] Qaisar | Computationally efficient adaptive resolution short-time Fourier transform[END_REF]. The ASA and the EASA display interesting features with the LCSS, which are not available in the classical case. They select only the active parts of the non-uniformly sampled signal, obtained with the LCADC. Moreover, they correlate the length of the selected window with the signal local characteristics. In addition, they provide an efficient reduction of the spectral leakage phenomenon [START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF][START_REF] Qaisar | Computationally efficient adaptive resolution short-time Fourier transform[END_REF]. Finally, the selected signal obtained with the activity selection algorithm is resampled uniformly. The resampling frequency Frs i of the i th selected window W i can be specific. According to [START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF][START_REF] Qaisar | Computationally efficient adaptive resolution short-time Fourier transform[END_REF], Frs i can be computed by employing the following Equations.

i i i t t L min max - = Eq. 4 i i i L N Frs = Eq. 5
Where, tmax i and tmin i are the final and the initial times of W i . L i is the length in seconds and N i is the number of samples laying in W i . The upper and the lower bounds on Frs i are Fs max and Fs min respectively.

The ARADC Resolution

Each stage has its impact on the overall effective resolution. In order to quantify the impact of each stage, the error sources of each stage are discussed and a method to compute the SNR (Signal to Noise Ratio) at each step is devised in the following subsections.

The LCADC SNR

Classically during an ideal A/D conversion process, the sampling instants are exactly known, where as samples amplitudes are quantized at the ADC resolution [START_REF] Kester | Data conversion handbook[END_REF][START_REF] Walden | Analog-to-Digital converter survey and analysis[END_REF], which is defined by the ADC number of bits. This error is characterized by the SNR, which can be expressed by Equation 6.

76

. 1 . 02 . 6 ) ( + = M dB SNR Eq. 6
Here, M is the ADC number of bits. It follows that the SNR of an ideal ADC depends only upon M and it can be improved by 6.02 dB for each increment in M. The A/D conversion process, which occurs in LCADCs, is dual in nature. Ideally in this case, samples amplitudes are exactly known, while the sampling instants are quantized at the timer resolution T timer . According to [START_REF] Sayiner | A Level-Crossing Sampling Scheme for A/D Conversion[END_REF][START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF], the SNR in this case is given by Equation 7.

(

)

timer sig T f dB SNR . log 20 19 . 11 ) ( - - = Eq. 7
Here, f sig is the input signal frequency. It shows that in this case, the SNR does not depend on M any more, but on x(t) characteristics and T timer . An improvement of 6.02 dB in the SNR can be achieved by simply halving T timer . Equations 6 and 7 respectively give the theoretical SNR of the classical and the level crossing converters. Practically, the A/D converters introduce further errors like the time jitter, the comparator ambiguity, etc in addition to the quantization phenomenon [START_REF] Walden | Analog-to-Digital converter survey and analysis[END_REF]. Usually, the converter practical SNR real is measured from the spectrum of a windowed sequence of the ADC output samples. A method of calculating the converter practical SNR is detailed in [START_REF] Baker | What does the ADC SNR mean[END_REF]. By knowing the SNR real of an ADC, its ENOB (Effective Number of Bits) can be calculated by employing the following Equation [START_REF] Kester | Data conversion handbook[END_REF][START_REF] Walden | Analog-to-Digital converter survey and analysis[END_REF]. Theoretically, the LCADC SNR can be improved as far as it is required by reducing T timer . But practically there is a limit, which is imposed by the analog blocks accuracy [START_REF] Sayiner | A Level-Crossing Sampling Scheme for A/D Conversion[END_REF][START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF]. In fact, the analog blocks determine the threshold levels precision. If these levels are known with uncertainty a, then this error must be taken into account and it will result into SNR degradation. Usually, the ADC SNR real is computed by employing the spectral analysis [START_REF] Baker | What does the ADC SNR mean[END_REF]. The LCADC output is non-uniformly distributed in time. Hence, its spectrum can not be properly computed with the classical tools. Several methods have been developed for the spectral analysis of the non-uniformly sampled. In [START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF], performances of the GDFT (General Discrete Fourier Transform) and the Lomb's algorithm are studied for the case of level crossing sampled signal. It is shown that these methods are erroneous because of the presence of wideband spectral noise. Hence, they can not provide a proper calculation of the LCADC SNR real . In context of the above discussion, a novel approach is proposed for the LCADC SNR real measurement. It does not require frequency domain transformation and calculates the SNR directly in time domain. The practical ADC is characterized by employing a monotone sinusoid [START_REF] Kester | Data conversion handbook[END_REF][START_REF] Walden | Analog-to-Digital converter survey and analysis[END_REF]. Therefore a similar signal given by Equation 9is employed in this case.

( )

θ π + = t f A t x sig . . . 2 sin . ) ( Eq. 9
Here, A, f sig and φ are the amplitude, the frequency and the initial phase. For the ease of process understanding ϴ = 0 is considered in the studied case. In the case of a mono harmonic signal it is possible to analytically calculate the level crossing instants [START_REF] Gretains | Time-frequency representation based chirp like signal analysis using multiple level crossings[END_REF]. Thus, in this case t n can be calculated by employing Equation 10.

      = A level f t m sig n arcsin . . . 2 1 π Eq. 10
Here, level m is the m th level crossing threshold. Amplitude of the n th level crossing sample x n can be calculated as: x n =level m . Hence, by employing this method, firstly an ideal LCSS is implemented for x(t), which provides an exact knowledge of both timeamplitude values of the level crossing samples. The only error occurs in the ideal LCADC is the time quantization [START_REF] Sayiner | A Level-Crossing Sampling Scheme for A/D Conversion[END_REF][START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF]. By assuming that the time error δt is uncorrelated to the input signal, it is modeled as a white noise. If δt n is the time quantization occurs for t n , then it can randomly takes a value between 0 and T timer . Thus, tq n (the quantized version of t n ) can be obtained by employing Equation 11. The time quantization also affects the amplitude value of the corresponding level crossing sample. The erroneous sample amplitude value is calculated by using Equation 12. In the case of a real LCADC, there also exists error due to the threshold levels ambiguity [START_REF] Sayiner | A Level-Crossing Sampling Scheme for A/D Conversion[END_REF][START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF]

The Activity Selection SNR

For a monotone sinusoid, the activity selection algorithm (ASA/EASA) parameters can be easily adjusted to avoid the signal truncation. In this case, the windowing is performed with the adaptive length rectangular function [START_REF] Qaisar | Computationally efficient adaptive resolution short-time Fourier transform[END_REF]. Since, the employed window shape is rectangular, it has no impact on the ARADC output resolution. It just selects the relevant parts of the LCADC output and passes them to the resampler block (cf. Figure 1).

The Resampler SNR

The resampling process requires interpolation, which changes properties of the resampled signal compared to the original one [START_REF] De Waele | Time domain error measures for resampled irregular data[END_REF]. For the practical LCADC, there exist uncertainties in the time-amplitude pairs of the levelcrossing samples (cf. Section 3.1). These uncertainties accumulate in the interpolation process and deliver the overall error at the ARADC output. If (tr n , xr n ) represents the time-amplitude pair of the n th interpolated sample. Then the n th reference sample amplitude xo n , which should be obtained by sampling x(t) at tr n , can be calculated by employing the following Equation.

(

)

n sig n tr f A xo . . . 2 sin . π = Eq. 16
The resampling error per interpolated observation Ie n is given by the absolute difference between xo n and xr n . The RMS (Ie) for N resampled observations can be calculated by employing a similar relation to Equation 13. Finally, the ARADC SNR real can be computed as a ration between the RMS (signal) and the RMS (Ie).

The Simulation Results

In order to illustrate the proposed method, a simulation is performed. In this case, x(t)=V max .sin(2.π.2300.t) is employed as input. V max = 0.9 v is chosen. Equation 7shows that for a fixed f sig , the ideal LCADC SNR varies as a function of T timer . In order to demonstrate this statement, the simulations are performed for fixed M=3 and by varying T timer between [2 0 ; 2 - Although the LCADC SNR is independent of M, yet an appropriate value of M should be chosen in order to ensure a proper reconstruction of the acquired signal [START_REF] Sayiner | A Level-Crossing Sampling Scheme for A/D Conversion[END_REF][START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF][START_REF] Qaisar | Computationally efficient adaptive rate sampling and filtering[END_REF][START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF][START_REF] Qaisar | Computationally efficient adaptive resolution short-time Fourier transform[END_REF].

In the case of a practical LCADC, the threshold levels ambiguity error a also occurs along with the time quantization error δt. The modeling of a is not straight forward and it depends upon the circuit architecture and technology employed for its implementation. A study on a for different LCADC implementations is out of the scope of this article. Here, the example of the AADC [START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF], is taken into account. In this case, the threshold levels are generated with the DAC (D/A Converter). Hence, a 3bits DAC is implemented in the Cadence circuit design tool using the STMicroelectronics 0.13-µm CMOS technology.

a mainly occurs because of the process and the mismatch variations, introduced during the circuit fabrication. The effect of the process and the mismatch variations on v out (threshold level) is studied for different input combinations by employing Monte Carlo simulations. It is found that due to the effect of process variations v out varies within the range of ± 0.21% of v out . Similarly, due to the mismatch variations v out varies within the range of ± 0.11% of v out . Finally, the variation of v out due to the combined effect of process and mismatch variations is calculated and it is ± 0.23% of v out . Following this, a is chosen equal to ± 0.23% of x n and xe n is computed by employing Equation 15. The practical LCADC SNR real is computed by employing the method discussed in Section 3.1. The simulation is performed for the same parameters, employed in the case of ideal LCADC. The obtained results are summarized in Table 2 Here, ENOB values are calculated by employing Equation 8. In a practical LCADC, the conversion error mainly consists of δt and a. Table 2 demonstrates that how a is limiting the ENOB. In the studied case, for higher T timer values [2 0 ; 2 -2 ] µs, the major error occurs because of δt and the employed value of a has minor impact on the SNR real . Contrary, with a further reduction of T timer the error occurs because of a is getting significant compared to the error introduced by δt. Hence, for lower T timer , a is the main limiting factor on the SNR real improvement. For the employed a, the limit on the achievable SNR real is around 59 dB for T timer =2 -5 µs. Further reduction in T timer will not introduce noteworthy gain in the SNR real , except by achieving an appropriate reduction in a.

The selected signal obtained at the ASA/EASA output is resampled uniformly (cf. Figure 1). A large range of interpolation functions is available. Computationally efficient solutions, such as the NNR (Nearest Neighbor Resampling) and the linear interpolations are employed for the resampling purpose.

For a fixed T timer , the resampled data SNR increases with the increase in M. The reason behind is that for any kind of employed interpolation, the upper bound on Ie n is imposed by q [START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF]. Here, q is the LCADC quantum and is given as: q=2V max /(2 M -1). Here, V max is the LCADC half amplitude range. It follows that an increase in M causes a reduction in q, which consequently results into a reduced Ie n .

The SNR of the uniformly sampled data, obtained in cases of the NNR and the linear interpolations is calculated by employing the method discussed in Section 3.3. It is performed by varying M and T timer between [3; 8] and [2 2 ; 2 -5 ] µ seconds respectively. The ARADC conversion error mainly consists of δt, a and Ie. Once the threshold levels are established, a remains constant [START_REF] Sayiner | A Level-Crossing Sampling Scheme for A/D Conversion[END_REF]. For given a once T timer is decided, then the next step is an appropriate choice of M and the interpolation order. From simulation results it is found that while employing the linear interpolation, M=8 is sufficient to approach the upper SNR real bound of 59 dB. On the other hand, in the NNR interpolation case, M=10 bits is required to achieve the upper SNR real bound. It follows that for certain δt and a, the upper achievable SNR real bound can be obtained for lower M with the increase in interpolation order. As an example, the same SNR real bound can be achieved for M=4 and T timer =2 -5 µ seconds, while performing the resampling with a fourth order interpolator. Note that higher order interpolators provide better results at the cost of an increased computation per resampled observation. Therefore, an appropriate interpolation order should be employed, which keeps the system computationally efficient, while not much affects the SNR real for the chosen parameters.

In order to compare the ARADC performance with the classical ADC, their SNR curves are plotted on Figure 2. The SNR values for the classical case are obtained by employing Equation 6. In Figure 2, SNR ARADC-NNR and the SNR ARADC-Linear , represent SNR ARADC obtained in the case of the NNR and the linear interpolations respectively. Figure 2 shows that in the studied case, for each value of M the SNR ARADC-NNR and the SNR ARADC-Linear remain higher than the corresponding classical one. It shows that for an appropriate choice of T timer , a and the interpolation order, a higher ENOB can be achieved for a given M, in the case of ARADC compared to the classical ADC.

Conclusion

The ARADC is well suited for low activity sporadic signals. For such signals, it leads towards a drastic computational gain compared to the counter classical approaches [START_REF] Qaisar | Computationally efficient adaptive rate sampling and filtering[END_REF][START_REF] Qaisar | Spectral Analysis of a signal Driven Sampling Scheme[END_REF][START_REF] Qaisar | Computationally efficient adaptive resolution short-time Fourier transform[END_REF]. A novel method to compute the ARADC SNR has been devised. It is shown that results obtained with the proposed method are in coherence with the theoretical ones, which verifies the proposed approach correctness. The ARADC SNR depends on M, T timer , a and the interpolation order. For a targeted application, an appropriate set of these parameters should be found, which provides an attractive trade off between the system computational complexity and the delivered output quality, while ensuring the proper signal reconstruction.

Figure 1 :

 1 Figure 1: The block diagram of the ARADC.

  LCADC conversion error per sample pointCq n is given by the absolute difference between x n and xq n . The RMS (root mean square) value of Cq for N level crossing samples can be calculated by employing Equation 13. Finally the SNR of an ideal LCADC can be calculated by employing Equation 14.

Figure 2 :

 2 Figure 2: The SNR curves for the ADC and the ARADC.

  Ce n is given by the absolute difference between x n and xe n . The RMS (Ce) for N level crossing samples can be calculated by employing a similar relation, shown by Equation13. Finally, the LCADC SNR real can be computed as a ratio between the RMS (signal) and the RMS (Ce).

	n xe	=	xq	n	±	a ∆	n	Eq. 15

. Let a n is the error introduced due to the quantization levels ambiguity into xq n . Then, the n th erroneous level crossing sample amplitude xe n , contains effect of both δt n and a n and it can be calculated by employing Equation

15

. The real LCADC conversion error per sample point

Table 1 :

 1 2 ] µs. The ideal LCADC SNR is measured by employing Equation 14. The obtained results are summarized in Table 1. The ideal LCADC SNR for fixed M=3 and varying T timer .

	T timer (µs) 2 0	SNR LCADC-TH (dB) 41.58	SNR (dB) 41.21
	2 -1	47.60	47.44
	2 -2	53.62	53.44
	In

Table 1 ,

 1 SNR LCADC-Th , represents the LCADC theoretical SNR, computed for the given parameters by employing Equation7. These results show accordance between SNR LCADC-Th and the obtained SNR, which verifies the authenticity of the proposed LCADC SNR measurement method.

Table 2 :

 2 . The real LCADC SNR for fixed M=3 and varying T timer .

	T timer	SNR LCADC-TH	SNR real	ENOB
	(µs)	(dB)	(dB)	(dB)
	2 0	41.58	41.16	6.54
	2 -1	47.60	47.05	7.52
	2 -2	53.62	52.88	8.49
	2 -3	59.64	56.71	9.13
	2 -4	65.66	58.93	9.50
	2 -5	71.66	59.45	9.58


Boyan Valtchanov
email: boyan.valtchanov@univ-st-etienne.fr

Viktor Fischer
email: fischer@univ-st-etienne.fr

Alain Aubert
email: alain.aubert@univ-st-etienne.fr

A coherent sampling-based method for estimating the jitter used as entropy source for True Random Number Generators

The paper presents a method, which can be employed to measure the timing jitter present in periodic clock signals that are used as entropy source in true random number generators aimed at cryptographic applications in reconfigurable hardware. The method uses the principle of a coherent sampling and can be easily implemented inside the chip in order to test online the jitter source. The method was carefully validated in various simulations that have shown that the measured jitter size corresponds perfectly to that of the jitter injected to the model. While the primary aim of the proposed measuring technique was the evaluation of the quality of jitter as an entropy source in random number generators, we believe that the same principle can be used in order to characterize the jitter in fast communication links as well.

Introduction

In the global communication era, more and more recent industrial applications need to secure data and communications. Many cryptographic primitives and protocols that are used to ensure confidentiality, integrity and authenticity use random number generators in order to generate confidential keys, initial vectors, nonces, padding values, etc. While random bit-stream generators can be easily implemented in analog or mixed-signal devices, the generation of random bit-streams is a challenging task when the generator should be implemented in a logic device like FPGAs (Field-Programmable Gate Arrays). Clearly, logic devices are well suited for algorithmic (pseudo) random number generators, but the true-random number generators need sources of randomness that are difficult to find and explore in logic devices. A mathematical model of the true random number generator (TRNG) is also a crucial element of the cryptographic application design since the final entropy of the generated random bit-stream could be characterized and thus certified if one is able to characterize the physical phenomenon that is used as the entropy source. If the model does not exist, there would be no guarantee that the final entropy of the output stream is true-random, pseudo-random or perhaps a mixture of random and pseudo random phenomena. Characterizing and monitor the entropy source (the jitter) and proposing a mathematical model is the main motivation of the paper.

Jitter as an entropy source for TRNGs

Many of the TRNGs known up to date [START_REF] Fischer | High performance True Random Number Generator in Altera Stratix FPLDs[END_REF], [START_REF] Kohlbrenner | An embedded true random number generator for FPGAs[END_REF], [START_REF] Sunar | A Provably Secure True Random Number Generator with Built-In Tolerance to Active Attacks[END_REF], use the jitter present in clock signals (generated using ring oscillators, phase-locked loops or delay-locked loops) as a source of entropy. The quality of the generated random bits is related to the parameters of the clock jitter. In order to avoid jitter manipulations and attacks, it is important to measure these parameters on-line and, if possible, inside the device. The jitter can be defined as a short-term variation of an event from its ideal position [START_REF] Technologies | Synchronous Optical Network (SONET) Transport Systems: Common Generic Criteria[END_REF]. In general, it is expressed as the variation in time of the zero crossing (rising or falling edge) of the clock signal. The jitter can be a good candidate for randomness generation, since its behavior is closely related to the thermal noise inside semiconductors [START_REF] Hajimiri | A general theory of phase noise in electrical oscillators[END_REF]. The advantage of the thermal noise employed as a source of randomness is that it is relatively difficult to manipulate it in order to realize an attack on the TRNG. The method presented in this paper considers only a truerandom (Gaussian) jitter component and it does not take into account the deterministic behavior of the jitter at this stage of our research. For a deeper understanding of the jitter behavior we recommend to read [START_REF] Sw Wedge | Predicting random jitter-Exploring the current simulation techniques for predicting the noise in oscillator, clock, and timing circuits[END_REF]. The proposed method allows to accurately quantify the random component of the jitter present in clock signals generated inside logic devices. Although the technique can be used to measure the jitter, it has been developed not for measurement or testing purposes, but rather for modeling a TRNG that uses the jitter as a source of randomness. The proposed measurement technique (see Figure 1) is based on a coherent sampling: the sampling of a periodic signal by another periodic signal featuring similar frequency [START_REF] Huang | An On-Chip Short-Time Interval Measurement Technique for Testing High-Speed Communication Links[END_REF]. The signal on the output of the sampler is called a beat signal and it is a low-frequency signal depending on the frequency difference ∆ between the two clock signals T ro 1 and T ro 2 .

Principle and theoretical background

Figure 2 shows the principle of the coherent sampling using two (clock) signals having similar frequencies and the resulting beat signal T Beat , representing the image of T ro1 . An example of this T Beat signal captured on oscilloscope is given in Figure 3. Using the infinite persistence of the oscilloscope, we can clearly see the variations of the period of the beat signal. These variations are the consequence of the jitter present in T ro 1 and T ro 2 signals.

Because of the coherent relationship between the two frequencies, each "half-period" of the beat signal is an integer number of the clock period T ro 2 . A counter clocked with this clock signal can thus be used in order to represent these variations. In the next section, we will discuss how we can compute the jitter present in T ro 1 by observing the variations in a population of several T Beat periods.

If the proposed technique would be used to measure precisely the jitter of the internal clock signal, one should use an accurate external low phase-noise VCO (Voltage Controlled Oscillator) as a sampling clock and accurately tune its period in relationship to the internal clock period in order to obtain a small ∆. Instead, in order to model the TRNG behavior and to measure the jitter inside the device, we have used two ring oscillators, implemented in the same FPGA. Both oscillators have the same number of inverters. In order to guarantee a small difference between clock periods (∆), the placement and routing have to be done manually. The final period difference is thus caused mainly by the different delays of the routing scheme selected by the placement and routing tool. Next, we will analyze the case, when only random (Gaussian) jitter component is present in the generated clock signals.

Measurement of the true-random jitter component

Let us assume that the two clock signals are derived from two internal ring oscillators, and let T ro1 Ideal and T ro 2 Ideal be the two ideal jitter-free periods. We need to achieve a small time period difference between T ro 1 Ideal and T ro 2 Ideal , namely:

T ro 2 Ideal = T ro 1 Ideal + ∆ Ideal . ( 1
)
This difference comes from the fact that even with the same number of delay elements the two ring oscillators differs due to process variations during manufacturing. With a careful placement, one can obtain ∆ of several tens of picoseconds. However the ∆ wont be reproducible from one chip to another. If a random jitter would be included in the previous equations, we obtain:

T ro 1 = T ro 1 Ideal + N (0, σ 1 ) = N (T ro 1 Ideal , σ 1 ) (2)
T ro2 = T ro2 Ideal + N (0, σ 2 ) = N (T ro2 Ideal , σ 2 ) (3)
Where N (0, σ) denote a zero-mean Normal distribution with standard deviation σ.

We can then express the difference ∆ by:

∆ = N (T ro 2 Ideal , σ 2 ) -N (T ro 1 Ideal , σ 1 ) (4)
∆ = N (∆ Ideal , σ 2 1 + σ 2 2 ) (5)
If σ 1 is the same as σ 2 , what is the case when the two signals are derived from internal ring oscillators, we get

∆ = N (∆ Ideal , √ 2σ) (6)
Otherwise one should make precise characterization of the VCO used to match the frequencies in order to measure the σ V CO .

According to [START_REF] Valtchanov | Modeling and observing the jitter in ring oscillators implemented in FPGAs[END_REF], we can express the length of T Beat as:

T Beat ∆ Ideal = N ( T ro 1 Ideal ∆ Ideal , T ro 1 Ideal ∆ Ideal σ 2 1 + σ 2 2 ) (7)
which, if σ 1 equals σ 2 , simplifies to:

T Beat ∆ Ideal = N ( T ro1 Ideal ∆ , T ro1 Ideal ∆ Ideal √ 2σ) (8)
The length of the resulting beat signal, T Beat can be then expressed as a normal process:

T Beat ∆ Ideal = N (µ T Beat , σ T Beat ) (9)
with the mean and standard deviation:

µ T Beat = T ro 1 Ideal ∆ , σ T Beat = T Ro Ideal ∆ Ideal √ 2σ (10)
In consequence, if we measure the µ T B eat and σ T Beat using the principle presented in Figure 1, which is based on the use of an 8-bit counter, we can precisely calculate the amount of the random jitter, expressed in 1σ ps, i.e. the RMS jitter (root mean square) present in the two clock signals using equation (11).

σ = σ T Beat ∆ Ideal T Ro Ideal ∆ Ideal √ 2 (11)

Simulation results

In order to validate equation (11), we have used a simulation model presented in [START_REF] Valtchanov | Modeling and observing the jitter in ring oscillators implemented in FPGAs[END_REF] and depicted in Figure 4. The random jitter is generated in text files using Matlab and then injected in VHDL simulation using the Textio package. We have injected different amounts of random jitter (RMS) to the clock signals and analyzed the obtained values of the counter. The T ro1 Ideal was set to 5 ns (200Mhz) and ∆ to 40 ps. The results of the simulations and recalculated jitter values using equation 11 are presented in Table 4. As it can be seen, the measurement precision that can be achieved is close to 1 ps RMS. Figure 5 present the case for 7 ps RMS jitter present in both T ro 1 and T ro 2 signals.

Injected 1σ

Measured

Discussion and conclusions

We have proposed a jitter measurement technique that can be embedded in FPGA devices for evaluating and monitoring of the source of randomness employed in true random number generators. The measurement technique can be used as well to characterize the jitter present in high-speed clock signals, if an external VCO (Voltage Controlled Oscillator) is used. The use of an external and precise clock source is necessary in order to closely match the period of the signal under test to the period of the reference clock signal. We have shown by simulation that the measurement error of the proposed method is less than 1 ps RMS of the random component of the jitter. However, in real world situations and especially inside FPGAs, the jitter can exhibit a non negligible deterministic component due to various factors (power supply variations, cross-talks, R-F interference, etc...). In this case, equation (11) cannot be used for random component jitter quantification and the deterministic jitter has to be considered, too. However, we believe that it is possible to integrate this deterministic behavior of the jitter in the proposed model. This integration is the objective of our current research.

Figure 1 :

 1 Figure 1: Random jitter component measurement based on the coherent sampling.

Figure 2 :

 2 Figure 2: Principle of the coherent sampling.

Figure 3 :

 3 Figure 3: Experimental T Beat signal example.

Figure 4 :

 4 Figure 4: Simulation setup.

Figure 5 :

 5 Figure 5: Histogram of the simulated T Beat .

Table 1 :

 1 Simulation results of the random jitter quantification.

			Measured	Calculated
	RMS jitter [ps]	µ T beat	σ T beat	1σ RMS [ps]
	10	121.93	4.03	10.19
	9	121.98	3.64	9.20
	8	121.97	3.24	8.19
	7	121.97	2.81	7.10
	6	121.98	2.47	6.24

This paper was partially supported by the Rhône-Alpes Region and Saint-Etienne Métropole,


Serge Samper -Pierre-Antoine

Adragna -Hugues Favreliere
email: hugues.favreliere@univ-savoie.fr

Maurice Pillet
email: maurice.pillet@univ-savoie.fr

Modeling of 2D and 3D assemblies taking into account form errors of plane surfaces

Keywords: 3D assembly, form errors, modal analysis, positioning force, statistical analysis, non conformity rate

The tolerancing process links the virtual and the real worlds. From the former, tolerances define a variational geometrical language (geometric parameters). From the latter, there are values limiting those parameters. The beginning of a tolerancing process is in this duality. As high precisions assemblies cannot be analyzed with the assumption that form errors are negligible, we propose to apply this process to assemblies with form errors through a new way allowing to parameterize forms and solve their assemblies. The assembly process is calculated through a method allowing to solve the 3D assemblies of pairs of surfaces having form errors using a static equilibrium. We have built a geometrical model based on the modal shapes of the ideal surface. We compute the completely deterministic contact points between this pair of shapes according to a given assembly process. The solution gives an accurate evaluation of the assembly performance. Then we compare the results with or without taking into account form errors. When we analyze a batch of assemblies, the problem is to compute the Non-Conformity Rate (NCR) of a pilot production according to the functional requirements. We input probable errors of surfaces (position, orientation, and form) in our calculus and we evaluate the quality of the results compared to the functional requirements. The pilot production can then be validated or not.

Introduction and purpose

Variational geometry is the first step of any tolerancing process even if most of the time this step is not conscious. Assumptions are made in order to simplify the "real" features. There are many ways to give geometry to some parameters. For example, we could say that a sphere could be given by one scalar if we do not care about its form errors. We could add to this parameter several errors like sphericity, cylindricity, conicity, and define deeply the geometry. In an extreme way, we could give the whole set of coordinates of the measured points as parameters. How many parameters do we need to define a given geometry? We propose a method allowing to define several levels to adjust the set of form parameters based on a geometric filtering.

A first answer is to start from the function from a top-down analysis (from the need to the solution). Then we could say that a minimum set of parameters could be sufficient, and the difference between the measured objects and the ideal has to be constrained within allowable values. We could want to know more about those deviations in order to better understand the influence of form errors on the function (finer analysis). The second answer is given by a bottom-up analysis where the manufacturer would like to know the forms of his production according to the customer's needs.

When the set of parameters is fixed (geometric model), the functional requirements (most of the time, precisions and assembly conditions) fix their limits. Tolerancing leads to quantifying the allowable limits fixed to those parameters. There are two ways to solve the compatibility of limits, the worst case and the statistical analyses. Simulations and industrial applications showed that statistical cases authorize larger tolerances. Thus, in this paper, we will show how we analyze a batch of assemblies by taking into account form errors. The aim is to validate a pilot production of assemblies.

In order to take into account form errors, the matting of surfaces has to be solved. If there are no defined criteria, the matting is undefined. We propose a method allowing to solve the matting of a pair of surfaces from the skin model by solving the static equilibrium of the pair depending on the matting force and the assembly device.

The solved values can come from a batch of measured surfaces of a pilot production that gives a batch of forms. In this paper, we use simulations that have the same statistic properties. Then we compute the mean and standard deviation of a probable batch of assemblies, which is checked with the functional requirements. This process gives NCR of a supplier assemblies according to the functional requirements of his customer taking into account full geometric deviations (dimension, position, orientation and form) of a pilot production in order to optimize the risks and costs of the function.

Assembly analysis model of form errors

Literature review

Parts assembly simulation is a heterogeneous research domain. Some papers are related to the stack of tolerances analysis but very few of them propose to take into account form errors.

Models for geometric description of features are the first step for a precision analysis. The complexity of the description is associated to the corresponding solving. From the most simple model based on a one-dimensional description assembly to a several dimensional one (discrete decomposition of features), we can find increasing data (and process) models. The one-dimensional models ([1], [START_REF] Graves | Tolerance Analysis Taylored to your Organization[END_REF]) are used to solve size variations. With those scalar models the analyses can be done with the worst cases or the statistic solving.

In order to take into account orientations and positions of associated surfaces, literature on the subject gives several models using a set of parameters. Requicha [START_REF] Requicha | Toward a theory of Geometric Tolerancing[END_REF] uses a geometric model with possible offsets of surfaces. Srinivasan [START_REF] Srinivasan | Role of sweeps in tolerancing semantics[END_REF] modified this model changing the offset by a sweeping sphere. Wirtz [START_REF] Wirtz | Vectorial tolerancing[END_REF] introduced the vectorial parameter concept which Gaunet [START_REF] Gaunet | Vectorial tolerancing model[END_REF] added the TTRS [START_REF] Clement | The TTRS: 13 constraints for dimensioning and tolerancing[END_REF] (Technologically and Topologically Related Surfaces) parameters of Clement. Sacks and Joskowicz [START_REF] Sacks | Parametric kinematic tolerance analysis of general planar systems[END_REF] developed a method for mobile planar mechanisms. The Jacobian Torsor of Desrochers and Laperriere [START_REF] Laperrière | Modeling Assembly Quality Requirements Using Jacobian or Screw Transforms : A Comparison[END_REF] is a kinematic model associated to an interval method. The T-MAP ® (Tolerance Map) model of Davidson and Shah [START_REF] Davidson | A New Mathematical Model for Geometric Tolerances as Applied to Round Faces[END_REF] allows to define those parameters in a constrain space. The Small Displacement Torsor (SDT) model of Bourdet and Clement [START_REF] Bourdet | A study of optimal-criteria identification based-on the small-displacement screw model[END_REF] is a kinematic based model that unifies the three-translation and the three-rotation parameters. Giordano and Duret [START_REF] Duret | Clearance Space and Deviation Space[END_REF] [START_REF] Giordano | Tolerance Analysis and synthesis by means of clearance and deviation spaces[END_REF] have developed the domain model in order to limit the SDT components to respect the zone concept of GPS standards [14].

The skin model [START_REF] Dantan | Geometrical product specifications -model for product life cycle[END_REF] is a measurable surface, which must be parameterized (filtered) in order to take into account errors. The filtering can be made by several methods. Most of them use periodic functions that can build simple to complex shapes. Fourier series is the first used decomposition method used in order to define periods on circular ( [START_REF] Cho | Quantitative circularity tolerance analysis and design for 2D precision assemblies[END_REF]) and linear shapes. Discrete harmonic transforms are well known in signal processing. 2D fourier is used by Capello et al. ([17]- [START_REF] Capello | The harmonic fitting method for the assessment of the substitute geometry estimate error. Part I : 2D and 3D theory[END_REF]) used to define form parameters of rectangular shapes. Huang et al. [START_REF] Huang | Mode-based Decomposition of Part Form Error by Discrete-Cosine-Transform with Implementation to Assembly and Stamping System with Compliant Parts[END_REF] developed a mode-base method for form error decomposition and characterization by using Discrete Cosine Transform (DCT). In the case of disk shapes, the periodic method of Zernike Polynomials [START_REF] Wyant | Basic wavefront aberration theory for optical metrology[END_REF] can be very useful especially in optics.

A Tchebychef Fourier Series model is used by Henke ([21]) in order to describe the forms of a cylinder and identify specific types of error shapes. Some models based on fractals [START_REF] Zahouani | Fractal models of surface topography and contact mechanics[END_REF] and wavelets [START_REF] Lee | Morphological characterisation of engineered surfaces by wavelet transform[END_REF] give solutions on local simple geometry (roughness to waviness on rectangular shapes). In the case of cylindrical shapes, Fourier Tchebychef models used by Gouskov [START_REF] Gouskov | Influence of the machine-tool defects on the formation of the residual defects on a surface: application to turning[END_REF] defines specific form errors such as eccentric or conical.

We proposed in [START_REF] Samper | Form Defects Tolerancing by Natural Modes Analysis[END_REF] a new way to define form error parameters based on the eigen-shapes of natural vibrations of surfaces. The originality is that the set of form parameters can be computed for any kind of shape. The modal shapes have vectorial space properties and are defined by using a discretization of the ideal surface. This method is resumed in section 2.3 .

In more general cases, some authors propose to define a basis built on the shape parameters of the CAD system as Cubeles ( [START_REF] Cubeles-Valade | Nominal and actual geometry explicit declaration. Application to dimensional inspection[END_REF]) and Pottman ( [START_REF] Pottmann | A concept for parametric surface fitting which avoids the parameterization problem[END_REF]) who propose to define a complex shape by a mechanism built on the set of control points of NURBS. Those methods based on polynomials can be used for a complex shape but they do not define systematically a basis with its properties and are difficult to use for a set of shapes. They are mainly used to define simple form error shapes.

Most of the models are top-down type ("a priori" shapes), they can be used to simulate shapes without knowing the production process. But some form error models uses bottom-up methods ("a posteriori" shapes) in order to use informations from measurements. A model of eigen-shapes derived from a Principal Component Analysis (PCA) is proposed by Summerhays and Henke in [START_REF] Henke | Methods for evaluation of systematic geometric deviations in machined parts and their relationships to process variables[END_REF] and [START_REF] Summerhays | Optimizing discrete point sample patterns and measurement data analysis on internal cylindrical surfaces with systematic form deviations[END_REF]. Camelio and Hu [START_REF] Camelio | Compliant Assembly Variation Analysis Using Component Geometric Covariance[END_REF] use a PCA of a production data to study sheet metal assemblies. In order to analyze 3D medical images, Nastar, in [START_REF] Nastar | Frequency-based nonrigid motion analysis : application to four dimensional medical images[END_REF] uses the modal shapes to define a set of instances (over time) of a human organ. In image processing [START_REF] Cootes | Timeline of developments in algorithms for finding correspondences across sets of shapes and images[END_REF] the PCA eigen-shapes is often used and gives interesting research directions.

In order to assess assemblies, some authors ( [START_REF] Ameta | Influence of form on Tolerance-Map-generated frequency distributions for 1D clearance in design[END_REF], [START_REF] Ameta | Using Tolerance-Maps to Generate Frequency Distributions of Clearance and Allocate Tolerances for Pin-Hole Assemblies[END_REF]) uses a form error zone which is added to a position and orientation zone. This type of method does not allow to mix position form and orientation and can increase costs. In order to analyze assemblies with a better accuracy, some works give matting solutions for circular [START_REF] Cho | Quantitative circularity tolerance analysis and design for 2D precision assemblies[END_REF] or prismatic joints [START_REF] Pierce | A method for integrating form errors into tolerance analysis[END_REF]. The solution of matting between two surfaces is computed by minimization of a geometric criterion [16] [35]. Huang and Kong [START_REF] Huang | Simulation and integration of geometric and rigid body kinematics errors for assembly variation analysis[END_REF] use Monte Carlo method in order to define variations of assembled surfaces described by DCT. An other way is to study the static equilibrium of assemblies (with an assembly force) with form errors [START_REF] Adragna | Analysis of Shape Deviations of Measured Geometries with a Modal Basis[END_REF].

Position and orientation models: SDT domain model

The skin model is reduced to its associated surfaces [START_REF] Duret | Clearance Space and Deviation Space[END_REF]. A frame is attached to each surface. Their possible deviations are size, position and orientation. In order to model those last two types of deviations, we can use a kinematic based model (displacements parameters). Bourdet and Clement [START_REF] Bourdet | A study of optimal-criteria identification based-on the small-displacement screw model[END_REF] defined the Small Displacement Torsor (SDT) that gives to positions and orientations of frames (of associated surfaces) a unified model. The SDT model allows calculating deviations without taking into account the observation frame. We can make simple operations on SDT (addition, transportation…). The model of SDT domains of Giordano [START_REF] Duret | Clearance Space and Deviation Space[END_REF] is based on the duality of contact conditions and SDT components inequalities. It builds a mathematical representation of those inequalities by a convex hull in a space of 6 dimensions (or less as in the case of planar or axisymmetric mechanisms with 3D domains). Each specification or functional requirement is expressed as a domain. The functional requirement (FR) verifications consist in checking that the resulting SDT of all the involved deviations are within the FR domain (inclusion test in a 6D space).

2.3

Form errors parameters: modal model

Calculation of the modal model

We propose to define geometry variation by using a geometrical basis built on the mode shapes. They are calculated by using the dynamic equations of either continuum or discrete mechanics. We summarize those equations in the following.

⎩ ⎪ ⎨ ⎪ ⎧ ρ → Γ = → div([σ]) at each point of the volume [σ]{n} → = f n → dS at each point of the outer surface (1)
Equations ( 1) are the local dynamic equations and the associated boundary conditions, where ρ is the density, → Γ the acceleration vector, → f v is the field forces vector (nil here) [σ] is the stress tensor, n → the normal vector to the dS local surface and f n → the applied local pressure on dS. The solving of the set of equations (1) for a continuum domain leads to obtaining the displacement solutions. We can use this method for some simple geometry (as line, square or disks) cases [START_REF] Leissa | Vibrations of plates[END_REF]. In most of cases, there are no analytic solution. We have to use a numerical method in most cases. As we compute modal form bases for any type of shape, we use the Finite Element Analysis (FEA) [START_REF] Zienkiewicz | The Finite Element Method for Solid and Structural Mechanics[END_REF] in order to determine the modal shapes. The discrete analysis is resumed in the following equations [START_REF] Graves | Tolerance Analysis Taylored to your Organization[END_REF][START_REF] Requicha | Toward a theory of Geometric Tolerancing[END_REF][START_REF] Srinivasan | Role of sweeps in tolerancing semantics[END_REF][START_REF] Wirtz | Vectorial tolerancing[END_REF][START_REF] Gaunet | Vectorial tolerancing model[END_REF][START_REF] Clement | The TTRS: 13 constraints for dimensioning and tolerancing[END_REF][START_REF] Sacks | Parametric kinematic tolerance analysis of general planar systems[END_REF].

The form error modes are the eigen-shapes of the perfect surface. If we use a FEA, the nodes of the meshing will be the corresponding measured points. The calculus of the natural mode shapes is made by the solving of the dynamic conservative equilibrium of the discrete geometry. The following equations ( 2) and the boundary conditions (which are described in displacements here) will give the eigen-shapes of the geometry.

M q ¨ + K q = 0 (2)
M is the mass matrix and K the stiffness matrix of the structure and q the nodal displacements vector. The solutions of (2) are written as follows:

q i =Q i Cos(ω i t) ( 3
)
Q i is the amplitudes shape vector and ω i the corresponding pulsation. The linear system (4) gives the n natural pulsations ω i and modes shapes

Q i . ⎝ ⎛ ⎠ ⎞ M -1 K - 1 ω i ² I d Q i = 0 (4)
If there are no boundary conditions we obtain as rigid body modes as allowed global displacements. Those modes would be used, in our model, to define global position and orientation of the associated surface.

Projection operator

We suppose that we have measured (by a CMM for example) a given geometry in order to obtain a set of errors (displacements of measured points). This set of Euclidian displacement vectors becomes the V vector of the measured feature. If we do measure the position of nodes, we have to make an interpolation in order to set V. As we have defined the modal basis Q, we can obtain the projections of V in this basis by calculating the projection equation ( 5) in a nonorthonormal basis:

V = ∑ i=1 n λ i Q i = Q λ (5)
Here, Q is the matrix of eigen-vectors Q i . λ is the vector of modal coefficients λ i . We choose to give Q i an infinity norm (||Qi||∝ = 1) in order to give λ i a metric sense. We obtain the vector λ by solving equations [START_REF] Gaunet | Vectorial tolerancing model[END_REF].

(Q t Q) -1 Q t V = λ (6)
The λ i coefficient can be showed as a histogram chart called modal spectrum. This representation is understood for metrologists and designers. We compute the vector residue of the projection in the reduced basis by

R(m) = V -∑ i=1 m λ i Q i with m<n (7)
The corresponding scalar value is given by equation [START_REF] Sacks | Parametric kinematic tolerance analysis of general planar systems[END_REF].

r(m) = ||R(m)|| (8)
As the natural shapes are defined by modal masses and modal stiffness their complexity are globally growing with i (growing order of the frequency). This property is very useful in order to filter shapes, by keeping lowest values of λ i . Most of measured shapes have got a decreasing modal spectrum.

Properties of modal bases

In Fig. 1 we show the first seventeen modes of a square plate with free boundary conditions. Solving can be made by a Finite Element Analysis (FEA) or with the help of analytical solutions [START_REF] Leissa | Vibrations of plates[END_REF]. The six first are rigid bodies that will reveal the position and orientation errors. We will convert them into SDT components of the associated surface. From the seventh, we can see form error modes. At a higher level of mode number, we will observe undulation and further, roughness. Form, undulation and roughness of a surface are distinguished by the level of periods of a flexure behavior.

In addition, we could solve, size variations (depending on the size parameters) by using membrane behavior. We can see with this elementary example the properties of the modal bases. The mode shapes are sorted by the value of the associated natural frequency. As we do not care about this value, we keep this sorting operator that naturally gives mode shapes in increasing complexity order.

⎩ ⎨ ⎧ ⎭ ⎬ ⎫ λ R1 λ R2 λ R3 λ R4 λ R5 λ R6 = ⎣ ⎢ ⎡ ⎦ ⎥ ⎤
α 61 α 62 … α 66 ⎩ ⎨ ⎧ ⎭ ⎬ ⎫ Tx A Ty A Tz A Rx A Ry A Rz A (9)
The geometric parameterization built on mode shapes has the following properties.

Geometric Basis:

Q i and Q j are independent if i ≠ j Uniqueness:
To a set of coefficients, λ i corresponds only one form given by V Stability:

Continuum parameters (no conditional test, bifurcations, …).

Growing complexity:

Q j is a more complex form than Q i if j>>i.

Exhaustiveness:

Any variation of a shape can be decomposed in the modal basis Metric of shape variations: λ i represents the metric coefficient of the mode shape Q i in the shape V

Analytic or numeric method?

We are making software that decomposes form errors of a measured shape. We can use analytic solutions of (1) in the cases of simple geometries. This way is easy to program, but in most of cases, we have to solve an FEM of eigenshapes of a given geometry with the help of (4). Those shapes are input in a matrix Q of vector parameters. This one can be an input of the software. Then the measured (or simulated) shape V (thus numeric) decomposition is obtained by the same solving for the two types of eigen-shapes (if eigen-shapes are analytic we discretize them). If possible, we measure displacements of points that are taken from the FEM meshing (points are nodes). If not, we use interpolations.

Separation of geometric errors

Dimension, position and orientation errors are easy (as concepts) to separate to form, waviness and roughness. Form errors are not easy to define, compared to waviness and roughness. Those geometric errors can be defined by the same type of mathematic model but do not lead to the same mechanic behavior. Surfaces sizes can be very different compromising a size separation criterion. We propose to define form errors as the first mode shapes having the main influence on assembly results (i.e. final position of surfaces in contact). If the contact force is small and the surfaces are filtered to remove waviness and roughness, the two surfaces are assembled with a number of points equal to the possible components of the contact torsor. In the case of planes, we can apply one force and two moments thus there must be three contact points according to a small contact force and two form error surfaces. The waviness is defined as intermediate mode shapes of the surface that are involved in a lot of small contact zones. They should remain smooth and could be deformed (elastically or plastically) with small forces. Roughness is defined as numerous very small contact zones (spots) that should (in most of cases) have a very small influence on the assembly position. In this paper we focus on position, orientation and form errors of surfaces.

2.4

Single assembly of two form errors

Contact determination

First we make a filtering of each shape by keeping the firsts (i<m) coefficients λ i . We obtain the deviation surface by subtracting the two matting shapes. We keep the m first coefficients considering the residue r(m). We can then build a meshing adapted to m (the Nyquist-Shannon evaluation is used in order to define a minimum number of FEM nodes). By this way, the number of possible matting facets is lower and their sizes are bigger. Then we define a difference surface (surface #1 becomes perfect and #2 is the difference surface). In order to select all the possible matting facets, we compute the convex hull of this surface. The solving of the intersection of the matting force whit this convex hull gives the matting facet thus the contact configuration.

Position parameters and FR verification

Assembly is obtained by knowing contact points and their displacements. The rigid coefficients λ Ri are deduced and equations [START_REF] Laperrière | Modeling Assembly Quality Requirements Using Jacobian or Screw Transforms : A Comparison[END_REF] gives the corresponding SDT. As explained in 2.2 , we can verify by an inclusion test if those SDT components are inside the FR domain or not. The assembly of form errors can be compared with the assembly of associated surfaces (model without form errors) in order to know the influence of form errors on the functional requirement.

Multiple assemblies verification and NCR

Our proposal is to evaluate the assembly of two batches having form errors. The statistical characteristics of the batches are supposed to be known (from measurement of the pilot production batches for instance). Then, the assembly deviation or the NCR can be evaluated through simulations. The following part 2.5.2 presents a solution to simulate a batch of form error based on its statistical characteristics. As this paper deals with no measure, the part 2.5.1 presents our solution to simulate a virtual batch characteristics (mean and covariance matrix).

Generation of statistical characteristics of a virtual batch

We have developed routines in order to make virtual skin models. We built virtual assemblies that have the properties of observed batches in order to explain our method even if we can study any kind of batch. The statistical description of a batch is based on a mean value of position, orientation and form and a covariance matrix of these parameters. The statistical characteristics can be obtained from modal descriptions of a measured or simulated batch. The modal characteristics of the simulated pilot production (mean vector µ λ and the covariance matrix C λ of modal coefficients) are obtained with the following method:

-A mean shape is identified, called the mother shape of the batch. Based on observed modal analyses, modal coefficients λ i follows a decreasing distribution law with a µ 0 /i type evolution (in Fig. 2, µ 0 =0.2 mm). The number of modal coefficients is limited to m that could be associated to a measurement filtering. -Based on this mother shape, a small-simulated batch is randomly drawn. The λ i modal coefficient are normally drawn with a σ 0 /i standard deviation (in Fig. 2, σ 0 =0.01 mm). The size n of the simulated batch directly affects the covariance matrix. Because random modal coefficients are independently obtained, a batch of a large size of shapes will have a quite non-covariant covariance matrix. Then, the smaller n is (n≥2), the higher the covariance factors are. -The mean modal signature µ λ and covariance matrix C λ (Fig. 2) of the simulated batch are computed to create the simulated batch of n shapes.

Generation of a virtual batch based on statistical characteristics

As the aim is to check that the pilot productions respect the functional requirement, simulations are used. The NCR evaluation is based on the statistical characteristics of the pilot production. Our approach is to simulate two batches with a large number of parts based on the two statistical descriptions of the two pilot productions. One batch example is presented in Fig. 3. The generation of one virtual batch is made as follows:

-A set of N random modal signature Λ r (of m modal coefficients) is drawn. Each modal coefficient follows a standard distribution law (for instance, normal distribution with a mean µ λi = 0 and a standard deviation σ λi = 1).

-We compute a diagonal matrix C λ diag (that contains the variance on the main direction) and the transfer matrix P (from the diagonal variance matrix to the batch covariance matrix) from the covariance matrix C λ of the pilot production batch. A modal signature Λ v of a virtual form error of the pilot production is obtained from the random modal signature Λ r with the following relation:

Λ v = P . C λ diag . Λ r +µ λ ( 10
)
Mean modal signature

Modal coefficient standard deviation

The mean modal signature : µλ

The modal covariance matrix : Cλ

Fig. 3 Generation of a virtual batch based on statistical characteristics

The assemblies can be analyzed as presented in section 2.4 . The NCR is computed by an inclusion test of the assembly SDT components inside of the FR domain.

Application on an assembly of two parts

The following example is a classic chain of parameters closed by a functional requirement (FR B ) of the B pair of surfaces. The pair of surfaces A (A 1 form part 1 and A 2 form part 2) presented in Fig. 4 has to keep a given distance. The location specification t gives the value of the functional requirement (between surfaces B). We assume that the pair of faces of the functional requirement B i has no geometric errors. We will compute how deviations on the contact pair A i influence the relative positions of B i . The positioning device is supposed to be perfect and the assembly is rigid.

Fig. 4 Functional Requirement of the assembly

If the supplier is asked to give interchangeable parts 1 and 2, tolerancing for each part should be made and the corresponding NCR should be higher according to the same production than with the non-interchangeability one because some compatible form errors assemblies would not be allowed. In this paper, we focus on the taking into account of form errors on the assembly NCR without interchangeability condition.

In this example, the contact surface The question for the producer is "are my products compatible with those specifications and those functions?". In order to help this relation between the customer and the supplier, we will make in the following some simulations of virtual assemblies.

We will analyze a simple assembly (Fig. 5) in order to show the method. The two matting surfaces A 1 and A 2 , may have position, orientation and form errors. The matting position of a pair of parts is given by the position set composed of an external position device and a matting force. If there is no mobility of the assembly, this set stops six DOF. In the example shown here, the external device is composed of three rods that could give three contact points because the plane joint is dominating (3 fixed DOF). The position of the resultant force F determines the three other DOF. We use this principle in order to compute the relative position SDT of the parts according to the matting set and their form errors. We simulate probable form errors of A 1 and A 2 with decreasing values of their λ 1i and λ 2i coefficient according to their i mode number (classic evolution of shapes).

a) 2D view of the assembly b) 3D view of the assembly Fig. 5 Assembly

2D assembly

We assume here that geometric errors are nil in the x direction. In order to model position, orientation and form errors, we build the associated structural model by a free-free beam. As the measuring is on the z direction, the FEM DOF are Ty and Rz (translation about y and rotation around z). The six first modes are presented in fig. 6. The modes #1 and 2 are rigid in order to identify the position and rotation parameters. All the modes have an infinite norm; the maximum displacement of each one is equal to ±1. For example, mode #2 is a rigid body rotation and the maximum value of a node's displacement is 1. The modes represent possible forms of each surface.

Fig. 6 of a free-free beam

If we measure a form ("A 1 " surface of Fig. 7) in order to extract its modal parameters, we obtain a vector of values shown here as a histogram. The most significant mode is the first mode (-0.01 mm amplitude of position), the second one is the third mode (-0.007 mm amplitude "banana" form), the third significant is the second (0.003 mm rotation amplitude) and the fourth significant is the ninth (0.003 mm amplitude). If we filter by using the lowest modes, we obtain the filtered shape of the focus in Fig. 7.b. This method uses the same scheme as Fourier and is well understood by geometric quality users.

Contact facet

We assume that the contact between a pair of surfaces is given by a static solution. Thus we have to analyze the static equilibrium of the two parts linked by the matting surface. This one is ideally plane (in this example) and has three kinematics DOF in 3D analysis. The DOF are all removed by external links (fig. 5). A joint force F creates the matting pair.

Our approach consists on identifying all the possible contact points of the two faces regarding their form deviations and a given pre-positioning mechanism. A positioning force identifies the contact points that define the stable positioning. The first step of the method is the identification of the possible contact points. We define the difference surface A linking surfaces A 1 and A 2 . This surface corresponds to the difference of the corresponding form deviations. The set λ A12 of modal parameters λ 1A2i of the equivalent surface A is obtained by the difference (eq. 11, Fig. 8) between the modal parameters λ A2 (surface A 2 ) from λ A1 (surface A 1 ). . By this way, we could consider that A 1 is perfect and A 2 has equivalent geometric errors. The second step of our approach is the identification of all (two here) the possible contact points between the two surfaces. The proposed solution is the use of the convex hull that identifies all the possible contact facets (segments here), hence all the possible contact points. Fig. 9.b shows the convex hull (polygon in 2D) of the difference surface (curve) and the potential contact points (small circles). In Fig. 10, we can observe assemblies with (Fig. 10.b) or without (Fig. 10.a) form errors. In Fig. 10.a we can see interpenetrations and gaps when associated surfaces are assembled (mean-square criterion). If we don't analyze form error contacts, we can either assemble mean-square surfaces or shifted associated surfaces, where the shift is two times the flatness value (one for each surface). Those two models enclose the values founded by our method.

λ 1A2 = λ A2 -λ A1 (11)

Form to position parameters

By using equations ( 9), we translate the rigid λ i coefficients in the corresponding SDT. In

Functional requirement verification of a single assembly

In Fig. 12, we present the resulting SDT of the assembly of associated surfaces (without taking into account form errors) of A 1 and A 2 . E R1A2 = E RA1 -E RA2 is the rigid SDT of the assembly (without taking into account form errors). E 1A2 is the assembly SDT with form errors. We have added E F1A2 , the difference between E 1A2 and E R1A2 in order to evaluate the taken into account form errors. As those torsors are calculated in the centre of surface A and the functional FR B requirement is calculated in the centre of surface B, we move the E 1A2 SDT in order to be compared to the FR B domain. This domain is built by the four inequalities of each point of B 2 from B 1 (zone condition). If each (four) summit S B2i of B 2 perimeter must remain within the segment of t range. The translation (along y axis) of each summit (called Ty SB2i ) respects the following inequalities then FR B is verified.

-t/2<Ty SB2i <t/2 [START_REF] Duret | Clearance Space and Deviation Space[END_REF] Those eight inequalities [START_REF] Duret | Clearance Space and Deviation Space[END_REF] are coupled so that four are sufficient. The corresponding FR B domain is shown Fig. The two SDT components verify the FR but we can see in this example that results are quite opposite.

Functional requirements verification of a simulated pilot production

The supplier has made a pilot production in order to check the NCR. He produces two independent batches of one hundred parts each (for part 1 and part 2). If we assume that any part 1 can be assembled with any part 2, we will have to test all the hundred assemblies. We created simulated a set of random shape having decreasing values of modal coefficients λ i (an hyperbolic law was used) and covariances such as a batch production should have. This batch is usefull to assess our method but we can analyse other set of shapes. Cvetko et al. [START_REF] Cvetko | New metrics for evaluating Monte Carlo tolerance analysis of assemblies[END_REF] gives an estimation of the NCR dispersion by the simulation dimension. In our example of 100 assemblies, the standard deviation of the NCR is about 0.5%. Here, we fixed t=0.1mm (L=20mm). As shown in Fig. 13, the conformity of this batch is better without taking into account form error in the assemblies. The risk is then to deliver some (9% here) non-conform assemblies without knowing the NCR. As those assemblies (Fig. 13.a) are made by using the associated surfaces (position, orientation but no form errors), the corresponding deviations are lower than those with form errors (Fig. 13.b). In addition, the results give the mean STD values (vertical cross in the centre of the ellipse) and the corresponding six-sigma ellipse. The supplier can identify the component involved in the NCR (here the mean of the batch is shifted on the translation direction at -0.025 mm). Here, the effective (with form errors) batch deviation could be inside the FRB if the mean were nil.

3D assemblies

The same assembly is presented with the assumption that the form error shape is a 3D one. The DOF are Ty (translation upon Y axis), Rx and Rz (rotations upon X and Z axis).

Surfaces form errors

We have built one hundred different form errors for A 1 and A 2 . The modal basis is computed by using the natural mode shapes of a free BC's (Fig. 1) of a square surface. We create probable (decreasing coefficients) form errors (Fig. 14

FCR verification of an assembly

We build a difference surface (Fig. 15.a) with the same method as the 2D one. A convex hull computation of this surface gives the set of possible contact facets. The intersection of the F axis with the convex hull gives the contact facet (thick triangle). Then we identify the three contact points belonging to this facet. Their z values give the three "rigid" coefficients of A 2 (λ 1 A2 , λ 2 A2 , λ 3 A2 ). By moving A 2 from A 1 with those three coefficients, we obtain the matting conditions. In Fig. 15, the contact facet is defined by three points. In order to see both surfaces and the contact facet, we can see, in Fig. 17 the topographic top view of the two surfaces (A 2 in filled lines and A 1 in dashed lines) and the contact facet i.e. a triangle defined by the three summits (P 1 ,P 2 ,P 3 ). Surface A 2 is repositioned with the three (λ 1A2 , λ 2A2 , λ 3A2 ) coefficients that are solution of the contact solving. In the two zooms, we can observe that in P 2 , the topographic lines have the same color (or gray) value (contact is effective) and in P 4 , different values (gap).

Functional requirements verification of a pilot production

After the virtual assemblies of one hundred pair of surfaces, we can compute the Non Conformity Rate of this batch. In Fig. 19 the two analyses of assemblies with or without taking into account form errors are shown (with mean SDT components and six sigma ellipsoid). In the case of associated surfaces, the whole batch is acceptable (within the FR B domain). The NCR of assemblies with form errors gives 14% out of FR B assemblies. -In this article, we show several graphics in order to explain our method, but we do not need any of them for making our solving.

-As the modal description of form errors is easy to define for any shape (with the help of a FEA), this method can be extended to other matting assemblies. The FEA has to be made only one time for each type of geometry. The set of modal shapes is stored in a library of modal shapes used by form deviations software. We can then solve assemblies in order to obtain SDT of a set of assemblies and computing NCR according to FR's.

Conclusion

Assembly simulations are needed to analyze deviation influences on the functional requirements. We have presented how this goal can be reached taking into account form errors on a static mechanism. Surface form errors are defined by the modal method. The natural mode shapes of the theoretical feature automatically build a geometric basis with several interesting properties as automaticity, exhaustiveness and natural sorting of shape complexity. Any form error can be decomposed in this geometric basis. The modal coefficients are used as shape parameters of 2D or 3D features. Keeping the most important coefficients makes a filtering operation. In order to compute the assembly of two surfaces, we define a filtered equivalent form error. Contact points must be on their convex hull. The set of corresponding facets define all the possible contact configurations. Assembly determination is solved by a static method. The use of the force assembly axis intersecting the convex hull of the filtered equivalent form error selects the contact facet. Then SDT model is used to transport this contact deviation to the functional requirement surface. By the domain method, this result is compared to the functional requirement domain. A simulation of a pilot production shows Non Conformity Rate with or without taking into account form errors. The supplier can use this approach in order to better answer the customer needs (NCR mastering at a good cost). We are making routines in order to analyse several type of contact surfaces. Libraries of modal shapes are made in this aim.

Fig. 1

 1 Fig. 1 Mode shapes of a free BC's square In Fig. 1, a Finite Element Model (FEM) is shown, the square is discretized by 441 nodes (21*21). We could solve 2646 (441*6) mode shapes. Sorting them by increasing frequency give increasing complexity. Since the first six mode shapes are rigid body type, we can write a linear operator [α ij ] in order to represent the corresponding modal coefficients λ Ri into the equivalent SDT components (9). As there are 6 rigid body mode shapes for a free structure and 6 kinematic Degrees Of Freedom (DOF), [α ij ] is 6*6 constants. The [α ij ] matrix is non-singular and its opposite gives SDT by knowing the λ Ri coefficients. It is then easy to move the position and orientation of frames.

Fig. 2

 2 Fig. 2 Generation of the modal statistical characteristics of a virtual batch

  A i is a square which side has a length of 40 mm (L AX =L z ). The functional requirement surfaces B i is a rectangle of 20*40 mm (L BX * L Z ) side lengths and the intermediate surface is a 40*40 mm square. The other dimensions of the parts do not influence our analysis. The two parts are not supposed to be the same and have the same type of deviations. We propose to use the tolerances of Fig. 4.b. Then the question for the designer is "what minimum values can we input for the tolerance t ?"

Fig. 8

 8 Fig. 8 Modal coefficients of the difference surface On Fig. 9.a surfaces A 1 and A 2 in their measured frames are shown. Fig. 9.b presents the difference surface (A 2 -A 1 ). By this way, we could consider that A 1 is perfect and A 2 has equivalent geometric errors. The second step of our approach is the identification of all (two here) the possible contact points between the two surfaces. The proposed solution is the use of the convex hull that identifies all the possible contact facets (segments here), hence all the possible contact points. Fig.9.b shows the convex hull (polygon in 2D) of the difference surface (curve) and the potential contact points (small circles).

Fig. 10

 10 Fig. 10 Comparison between assembly methods

  Fig. 11 the SDT of associated surfaces A 1 (E RA1 ) and A 2 (E RA2 ) are relative deviations from the nominal surface. The SDT components of E R1A2 (Error Rigid SDT) of A 2 from A 1 are shown in the 2D (T yA ,R zA ) space of SDT coordinates (Fig. 12.a). Here, we can observe the influence of the form error on the position of this plane joint.

Fig. 11 Deviations

 11 Fig. 11 Deviations SDT of associated surfaces

  12.c. The two assembly SDT (with or without taking into account form errors) are first given in the local frame attached in the centre of the surface A (Fig. 12.a). The set of torsors E X1A2 are transported in the same frame as the FR B domain (Fig. 12.b). The functional verification is made by an inclusion test showed in Fig. 12.c. In the presented case, deviations are compatible with the functional requirements. a) Assembly of A1 and A2 observed in the centre of surface A b) Assembly of A1 and A2 observed in the centre of surface B c) FR checking observed in the centre of surface B Fig. 12 Resulting SDT and Functional Requirements of one assembly

Fig. 13

 13 Fig. 13 Conformity of a set of 100 assemblies

  .a,b,c) for A 1 and A 2 . The three first λ i are associated to translations and rotations of surfaces. We solve the assembly of the associated surfaces (Fig. 14.d) in which form errors remains. Some interferences between A 1 and A 2 are shown in Fig. 14.d. a) 35 firts λ i A1 coefficients of A and contact facet of the mean assembly c) A 2 surface d) Assembly of A 1 and A 2 associated surfaces Fig. 14 Surfaces A1, A2 and assembly of associated surfaces.

  Fig. 16 3D view of the contact

Fig. 17

 17 Fig. 17 Assembled surfaces and contact points, a 2D viewWe calculate the SDT of relative position of A 2 from A 1 corresponding to those displacement values. Then we transport this torsor in the centre of the specification B. The functional requirements are respected if the SDT components are within the FR B domain (Fig18). The assembly of form error surfaces is closer to the limits than the one with associated surfaces. It is easy to add SDT to several joints and check the FR B . The contact conditions of the FR B are given the by the eight inequalities of[START_REF] Duret | Clearance Space and Deviation Space[END_REF]. Here, as there are three DOF parameters, those eight inequalities (half spaces) define an octahedron.

Fig. 18 FR

 18 Fig. 18 FR B domain and assembly verifications

Fig. 19

 19 Fig. 19 Functional requirement (FC B ) 3D Domain and SDT of assembliesRemark:-In this article, we show several graphics in order to explain our method, but we do not need any of them for making our solving.-As the modal description of form errors is easy to define for any shape (with the help of a FEA), this method can be extended to other matting assemblies. The FEA has to be made only one time for each type of geometry. The set of modal shapes is stored in a library of modal shapes used by form deviations software. We can then solve assemblies in order to obtain SDT of a set of assemblies and computing NCR according to FR's.

Acknowledgement

The authors acknowledge the financial supports of the French Ministry of Research, the European Union, the "Assemblée des Pays de Savoie" and the Competitiveness Pole "Arve Industrie Haute Savoie Mont-Blanc".