Estimating the Signal-to-Noise ratio under repeated sampling of the same centered signal : applications to side-channel attacks on a cryptoprocessor Gilles R. Ducharme, IMAG, Univ Montpellier, CNRS, Montpellier, France Philippe Maurine, LIRMM, Univ Montpellier, CNRS, Montpellier, France

Abstract-This paper introduces an estimator of the signalto-noise ratio in the framework where a noisy source emits the same signal a number n of times. The estimator has the structure of a U -statistic from which derives many desirable properties : it is unbiased, consistent and, being a Rao-Blackwellisation of existing proposals, is closer to optimal variance-wise. However, its variance is numerically difficult to evaluate and two approximations are obtained to facilitate its use in practice. These allow to quantify the improvement in variance, which is found to be substantial as the estimator needs roughly one third of the data previously required to perform similarly. Moreover, a simulation shows that the estimator is approximately normally distributed for n as small as 10, which allows for accurate inference. The estimator is then applied to data arising in a cryptanalysis, where the numerical security of a cryptoprocessor is tested against a side-channel attack. This problem is a representative of situations where the signal-to-noise ratio must be precisely estimated for small n. We derive a rigorous data-driven approach that is shown to much enhance the efficiency of standard side-channel attacks.

Index Terms-Asymptotic normality, Cryptanalysis, Sidechannel attack, Signal-to-noise ratio, U -statistic.

I. INTRODUCTION

We assume the following framework: a received signal M = (M 1 , ..., M T ) (the prime denotes transposition) is observed at various time points t in the range {1, ..., T }. This received signal consists of an unknown deterministic signal τ = (τ 1 , ..., τ T ) perturbed by a random noise η = (η 1 , . . . , η T ) . Hence the model is

M = τ + η. (1)
We assume that the signal can be centered, i.e. t τ t = 0, and that η obeys the T -dimensional multinormal distribution, η ∼ N T (0, σ 2 η I) with expectation 0 and covariance matrix σ 2 η I, where I denotes the T -dimensional identity matrix. The signal-to-noise ratio is defined as

SN R = τ τ /(T σ 2 η ). (2)
In our context, we further assume the same signal τ of fixed length T is emitted n times, each time perturbed by a different random noise η. Also T > 2 for reasons to be explained in Section ??. This context is encountered in particular in cryptanalysis and more precisely in side-channel attacks (SCA) on cryptoprocessors (Diop et al. [?]), which is our motivating application. In SCA, a plaintext is sent for encryption to a cryptoprocessor, which mingles it with a secret key and produces the signal τ . An attacker tries to spy on this signal but gets to observe a version M of τ perturbed with random noise η. The experiment can be replicated n times, thus generating M i = τ + η i , i = 1, ..., n where nothing can be assumed about τ beyond its centering. For more details, see Section ??.

Estimating the SN R is an important and difficult problem in signal processing. Here, we introduce the estimator

SN R = 2 n(n -1) n i=1 n j=i+1 h(M i , M j ), (3)
where

h(M i , M j ) = (T -2) M i + M j 2 2T M i -M j 2 - 1 2 , (4)
and • 2 is the squared Euclidean norm. We establish its stochastic behavior and show that SN R has theoretical and practical advantages over competitors, e.g. those of Diop et al.

[?], Simon and Dolinar [?] and Coppola et al. [?]. In particular, we show it is unbiased and approximately normal for n as small as 10. We also show that (??) derives from an application of the Rao-Blackwell theorem so its variance is improved with respect to the competitors. This variance is however numerically difficult to evaluate and two simplified expressions are obtained. A comparison shows that these are close to the exact values and that our estimator requires roughly three times less data to perform similarly as these competitors. These advantages allow to establish accurate confidence intervals and derive powerful statistical tests pertaining to the SN R. Then, SN R is applied to data arising in a cryptanalysis, where the numerical security of a cryptoprocessor is tested against a sidechannel attack. This problem is a representative of situations where the SN R must be precisely estimated for small n. We derive a rigorous and principled data-driven approach that is shown to much enhance the efficiency of standard side-channel attacks.

The paper is organized as follows. Section ?? recalls some useful statistical results. Section ?? derives the theoretical properties of SN R and presents the two approximations. Section ?? explains how accurate statistical inference can be produced about the SN R. Section ?? investigates the variance improvement of our estimator. Side-channel attacks (SCA) on cryptoprocessors, the motivating application of the present work, are explained in Section ??. A new data-driven SN Rbased SCA is presented to illustrates the usefulness of the methods of the paper. A conclusion closes the paper.

II. PRELIMINARY STATISTICAL RESULTS

A. Stochastic behavior of a U -statistic

An expression of the form (??) is a U -statistic. The function h(•, •) is its kernel and is of order 2. Many theoretical results about the stochastic behavior of U -statistics are listed in Chapter 5 of Serfling [?]. In particular, let U n be a U -statistic of order 2 whose variance exists and such that (hereafter E(X) and V(X) denote the expectation and variance of the random quantity X)

E(h(M i , M j )) = θ. (5)
Then U n is an unbiased estimator of θ, i.e.

E(U n ) = θ. Moreover, setting h 1 (m 1 ) = E [h(m 1 , M 2 )] (here m 1 is a fixed value of M 1 ) with h 1 (m 1 ) = h 1 (m 1 ) -θ, we have V(U n ) = 2 n(n -1) 2(n -2)ζ 2 1 + ζ 2 2 , (6)
where

ζ 2 1 = E( h 2 1 (M 1 )) and ζ 2 2 = V(h(M 1 , M 2 )). Impor- tant for statistical inference is the fact that if ζ 1 > 0, U n is approximately normally distributed with expectation θ and variance V(U n ), whose dominant term is 4ζ 2 1 /n.
B. Stochastic behavior of quadratic forms in normally distributed random vectors.

We recall that if Z ∼ N T (0,

σ 2 I), then Z + µ 2 /σ 2
has the non-central chi-square distribution χ 2 T (µ µ/σ 2 ) with T degrees of freedom. We also recall for convenience the following well-known theorem.

Theorem 1. Let W 1 ∼χ 2 d1 (δ 2 ) and W 2 ∼ χ 2 d2 (0) be in- dependent random variables. Then F = (d 2 W 1 )/(d 1 W 2 ) ∼ F d1 d2 (δ 2 )
, the non-central Fisher distribution with degrees of freedom (d 1 , d 2 ) and non-centrality parameter δ 2 . Moreover

E(F ) = d 2 (d 1 + δ 2 ) d 1 (d 2 -2) , if d 2 > 2, (7)
V(F ) = 2d 2 2 (d 1 + δ 2 ) 2 + (d 1 + 2δ 2 )(d 2 -2) d 2 1 (d 2 -4)(d 2 -2) 2 if d 2 > 4.
(8) Finally E(F ) and V(F ) exists when d 2 > 2 and d 4 > 4 respectively.

III. STOCHASTIC BEHAVIOR OF THE ESTIMATOR

A. Bias

Write V (η i , η j ) = η i + η j + 2τ 2 / η i -η j 2 and rewrite h(M i , M j ) = h(η i , η j ) = (T -2) 2T V (η 1 , η 2 ) -1 2 , h 1 (M i ) = h 1 (η i ), etc. It is easy to see that η i -η j η i + η j + 2τ ∼ N 2T 0 2τ , 2σ 2 η I 0 0 I . (9)
From Section ??,

η i -η j 2 is distributed as a 2σ 2 η χ 2 T (0) while η i + η j + 2τ 2 has for distribution a 2σ 2 η χ 2 T (2ρ), upon introducing the notation ρ = τ τ /σ 2 η = T ×SN R. Thus, V (η i , η j ) is distributed as χ 2 T (2ρ) /χ 2 T (0)
, where the two χ 2 are independent. From Theorem ??, this is the stochastic representation of a non-central Fisher distribution with both degrees of freedom equal to T and non-centrality parameter 2ρ. Using (??) shows that SN R is unbiased when T > 2.

B. The term ζ 2 2

From (??), when T > 4,

ζ 2 2 = T 2 + 2ρ(ρ -2) + T (4ρ -1) T (T -4) . ( 10
) C. The term ζ 2 1
Let η 1 denote the fixed value of the random η 1 corresponding to m 1 and set x =τ η 1 , y = η 1 η 1 . After completing the square and using the notation Ψ d (•; δ 2 ) for the cumulative distribution function of a χ 2 d (δ 2 ), we get

P [V (η 1 , η 2 ) ≤ v] = P η 2 + η1 + 2τ 2 ≤ v η 2 -η1 2 = Ψ T b(v); δ 2 (v) if v ≥ 1 1 -Ψ T b(v); δ 2 (v) if v < 1, (11)
where

δ 2 (v) = 4 (1 -v) 2 ρ + (1 + v)x + (1 + v) 2 y/4 , (12) b
(v) = 4v (1 -v) 2 (ρ + 2x + y) , (13)
after noticing that the term σ 2 η appears solely through τ /σ η and η 1 /σ n , so that scale invariance allows to set σ 2 η = 1 without loss of generality.

Write V (η 1 ) = E (V (η 1 , η 2 )). The positiveness of V (•, •) entails V (η 1 ) = ´∞ 0 P [V (η 1 , η 2 ) ≤ v] dv which can in prin- ciple be computed via (??). Moreover because E (h 1 (η 1 )) = SN R, we have E (V (η 1 )) = (2ρ + T )/(T -2), so that ζ 2 1 = T -2 2T 2 × E V 2 (η 1 ) -((2ρ + T )/(T -2)) 2 . ( 14
)
It is shown in the Appendix that

E V 2 (η 1 ) = ∞ 0 √ ρy -√ ρy    1 0 (1 -Ψ T b(v); δ 2 (v) dv + ∞ 1 Ψ T b(v); δ 2 (v) dv    2 × Γ( T 2 ) √ πΓ( T -1 2 ) (1 -x 2 /ρy) (T -3)/2 / √ ρy × ψ T (y; 0)dxdy. (15)
This triple integral is difficult to evaluate numerically. A first approximation is based on the following proposition whose proof is in the Appendix.

Proposition 1.

V (η 1 ) = T + 4x + y + 4ρ T + y + O(T -1/2 ). (16)
Injecting into E(h 2 1 (η 1 )) leads to the explicit expression :

ζ 2 1 = (T -2) 2 T 3 ρ [2(ρ -1) -2e T /2 (T (ρ -1) -ρ)E T /2 (T /2)- T ρe T /2 E T /2 (T /2) 2 ] + O(T -1/2 ), (17)
where E n (z) denotes the exponential integral function. This approximation is only first order accurate and, for small T , is degraded by the asymmetry of a non-central χ 2 . As shown in the Appendix, a Cornish-Fisher expansion leads to :

Proposition 2. Let c = -2 3 /(6 √ T ) . Then V (η 1 ) 1 (T + y) 2 -2c 2 (T + 2y) c 2 (4y + 8x -2T ) + (T + y)(T + 4x + y + 4ρ) + 4 c 2 (T 2 (2x + y + ρ)(1 + y(4x + y + 4ρ)) + 2T 3x 2 + y 2 + 3yρ + ρ 2 + 4x(y + ρ) -2c 2 2yρ -2x 2 + T (2x + y + ρ) 1/2 . ( 18
)
Injecting into E(h 2 1 (η 1 ))) does not yield an explicit expression but turns out to be easy to evaluate numerically.

We numerically evaluate the accuracy of these two approximations for ζ 1 . First the value of the triple integral defining ζ 1 is evaluated by 2000 Monte Carlo replications for various values of T and SN R. Then the two approximations are computed. The results appear in Table ??, where MC refers to the Monte Carlo values of ζ 1 .

As expected, the 1 st order approximation is less accurate, but can be useful when T ≥ 25 . The 2 nd order approximation offers much improvement in all investigated cases.

IV. INFERENCE FOR SNR

Summing up, SN R is an unbiased and consistent estimator of SN R approximately distributed as

N SN R, σ2 SN R = 2 n(n -1) 2(n -2) ζ2 1 + ζ2 2 , (19)
where ζ1 , ζ2 are obtained by plugging in the value of SN R in (??) and (??). Statistical inference about SN R derives from (??). Let z α denote the α-th quantile of the N (0, 1) distribution. A twosided confidence interval of approximate level

1 -α is SN R ± z 1-α/2 σSNR . ( 20
)
Similarly the null hypothesis H 0 : SN R ≤ θ 0 can be rejected at approximate level α in favor of H 1 : SN R > θ 0 when

Z n = ( SN R -θ 0 ) σSNR > z 1-α/2 . ( 21
)
Adaptation to other inference problems is obvious. The only point that remains to clarify regards the accuracy of the risks associated with these inferential procedures.

To investigate this, a simulation study was performed. For the pairs (SN R, T ) listed in Table ??, samples of size n = 10 and 20 where taken according to model (??). 95% confidence intervals (??) were computed using the 1 st and 2 nd order approximation and their coverage probabilities obtained from 5000 replications. Results appear in Table ?? and show that the true coverage probabilities are close to nominal for n as small as 10 when SN R ≥ 1 and the 2 nd order approximation is used. When SN R<1, these approximations are less precise because the distribution of SN R is bounded to the left. The 1 st order approximation performs close to nominal when T ≥ 25. In view of these results, we recommend the 2 nd order approximation while the explicit and convenient 1 st order approximation can be used when T ≥ 25.

V. IMPROVEMENT IN VARIANCE

Another unbiased estimator of SN R is

SN R i,i+1 = 2 n n i=1 i:odd h(M i , M i+1 ), (22)
which is a variant of proposals in e.g. Computationally, SN R i,i+1 and its other linear (e.g. a single sum of terms) variants are easier to work with, but inferior to SN R . First, they are not invariant to the pairs taken in h(•, •), a serious flaw, and their variance is larger than that of SN R. To see this, note that the set of all M i , ranked by their first component, is a sufficient statistic. From general results on U -statistic (see Serfling [?], Theorem, p.176), SN R in (??) is the expectation of h(M i , M j ) conditioned on this

V( SN R) ≈ 4ζ 2 1 n ≤ 2ζ 2 2 n = V( SN R i,i+1 ). (23)
The ratio √ 2ζ 1 /ζ 2 compares the length of confidence interval (??) to that based on SN R i,i+1 , It serves as a measure of the relative efficiency of statistical inferences (confidence interval, hypotheses test) based on SN R. Table ?? shows the values of this ratio for the pairs (SN R, T ) in Table ??, with a ratio of 0.5 indicating that a statistical procedure based on a linear estimator requires a sample of size 4n to perform similarly to one based on n received signal using SN R. The average of the entries in this table is approximately 0.55 showing that, overall, statistical inference based on SN R requires roughly over three times less observations to perform similarly.

Note that other estimators of the SN R in the present context have been proposed, notably some based on correlations (Bershad and Rockmore [?]). Their stochastic behavior is largely unknown and will not investigate here.

VI. APPLICATION TO SIDE CHANNEL ATTACKS

Cryptoprocessors perform the calculations ciphering sensitive information. They apply algorithms that use secret keys and if an opponent obtains the key, the ciphered information can easily be transformed back into plaintext. When the cryptoprocessor is ciphering, the plaintext, mingled with the secret key, gets converted into an analog signal τ . This conversion can leak information about τ to be exploited by an opponent to recover the key. This is referred to as a side-channel attack (SCA) and constitutes one form of cryptanalysis. Among the most successful SCA are those based on the electromagnetic (EM) emanations as the cryptoprocessor manipulates the text/key combination. EM-based SCA are easy to perform cryptanalyses because the attacker can send the plaintext for encryption a number n of times and observe, through an inexpensive sensor, the emanated EM radiations at various time points. The resulting EM curves M i , i = 1, ..., n, are modeled as in (??), where the η i are random noises arising from various causes. The process is calibrated to provide a centered τ .

If the attacker can separate the noise from τ , parts of the key can be recovered. To give an idea of how this can be done, Figure ?? shows a sample of n = 10 EM curves M i (in blue) while a cryptoprocessor was ciphering the same plaintext with the AES algorithm. The black curve is the (unknown to an attacker) signal τ . The AES algorithm involves so-called "rounds" and the oscillations of the black curve are related to the calculations performed during the last round using a part of the secret key. The knowledge of these oscillations, coupled with some other available informations, can reveal the calculations being performed and thus that part of the key. This indicates that a sophisticated SCA on the blue curves could be successful in exploiting the leakage to recover the key. The whole SCA process is however rather involved; we refer the interested reader to e.g. Brier et al. [?] for details.

When dealing with noisy measurements, filtering is a preprocessing technique that seeks to reduce the noise by focussing on the frequencies that best carry the signal. In general, a fruitful filtering requires some a priori knowledge about the informative frequencies, otherwise these may get filtered out. This is especially true in EM-based SCA where a small change in the position of the sensor can drastically alter the EM curves. In the present SCA context, as in many other applications, such a priori knowledge is unfortunately often unavailable. Tiran et al. [?] have devised an approach to detect the informative frequencies in a EM-based SCA. They theorize a quantity they refer to as the Leakage-to-Noise Ratio (LN R) = τ * τ * /(T σ 2 η ), where τ * is the part of τ that is informative about the key. They argue that, in general, LN R ≤ SN R. Then they translate by Fast Fourier (FF) transform the time domain received signal M i of length T into the frequency domain over an interval [F M IN , F M AX ] whose bounds are chosen in relation with T . For the frequency interval of width δf centered on f , they introduce the quantity LN Rp(f ) = SN R(f )/f , where SN R(f ) is the SN R computed from the data at frequency f and the division by f penalizes the large frequencies. Thus large values of LN Rp(f ) are indicative of interesting frequencies. They evaluate from the FF transformed M i a crude approximation of SN R(f ) and retain only those frequencies where the corresponding LN Rp(f ) is greater than an ad hoc threshold. They set all others frequencies to zero and translate back the results into the time domain to get filtered signal that are then processed by a standard SCA to extract the leaked information about the secret key.

Here we use SN R to supplement their approach with a better estimate of LN Rp(f ) and exploit its stochastic properties to derive a rigorous and data-driven way of identifying the candidate frequencies to be retained. Our goal is to illustrate that the efficiency of standard SCA can be enhanced by the methods of the paper.

For this, we develop a multiple testing procedure based on (??) to test all pairs of hypotheses H 0 (f ) :

LN Rp(f ) ≤ λ 0 vs H 1 (f ) : LN Rp(f ) > λ 0 in the grid of values f ∈ [F M IN , F M AX ].
The f values for which the tests reject H 0 (f ) are the candidate frequencies carrying leakage information.

Setting LN Rp(f ) = SN R(f )/f , the test statistic Z * n (f ) = f × ( LN Rp(f ) -λ 0 )/σ SN R is under H 0 (f ) approximately N (0, 1
) and H 0 (f ) gets rejected when Z * n (f ) > z α * where α * is taken to control the risk of at least one false rejection (a family-wise Type 1 error) in this multiple testing situation.

To show the efficiency of the above data-driven filtering approach, an experiment was conducted in the following way. Recall from Section ?? that with n = 10, asymptotic normality offers a good approximation to the distribution of SN R. In a preliminary study, two sets of n = 10 EM curves of length T = 3000 were collected above an AES mapped into a cryptoprocessor-emulating FPGA (field-programmable gate array) device. The first set was collected with the EM sensor positioned far from the FPGA surface in order to obtain noisy data while the second set was collected with the sensor much closer to get small noise measurements. The SN R(f ) values were then computed for the 499 frequencies in [0.8M Hz, 400M Hz] by steps of 0.8M Hz with δf = 0.8M Hz. Here SN R(f ) varies from 1.21 to 14781, so it is reasonable to use the explicit first order approximation for ζ 2 1 given by (??). Also we took the standard α = 5% with α * = 1 -α/499, derived from the Bonferroni inequality. Subject-matter considerations, coupled with the attractive equilibrium between signal and noise, lead to the choice λ 0 = 1.

Panel (a) of Figure ?? shows the evolutions of LN Rp(f ) (the square root is used for increased clarity) with f for the noisy measurements while Panel (b) shows the same for the small noise data. The black dots correspond to the frequencies f where H 0 (f ) has been rejected at level α * . They suggest that measurements should be filtered so as to retain the frequencies in [0.8; 54M Hz] ∪ [130M Hz, 160M Hz] ∪ [240M Hz, 260M Hz] for noisy measurements, while for small noise, measurements should be filtered so as to keep frequencies [0.8, 104M Hz] ∪ [120M Hz, 260M Hz] ∪ [300M Hz, 360M Hz], if one aims at preserving the leakage that can be exploited by the subsequent SCA.

With this knowledge gained from the filtering step, a fullfledged SCA called a correlation power analysis (Brier et al. [?]) was applied to 5000 new EM-curves observed from the same apparatus and on their corresponding 5000 filtered measurements as described above. This allowed to compute the global Guessing Entropy (gGE), a figure of the merit of an attack, which is used to compare the efficiency of the SCA with and without filtering. The Guessing Entropy (GE) gives the average position of the correct sub-key of 8 bits in the ranking of the 256 possible sub-keys provided by the SCA. If the rank of a sub-key is one, the SCA successfully discloses the sub-key. The gGE is the average of the 16 GE values obtained for each of the 16 sub-keys constituting the whole 128 bits AES key. When gGE = 1, the whole key has been correctly disclosed by the SCA. Note that in SCA, the security of a ciphering device is typically evaluated using information theoretic metrics such as the entropy, while actual attacks are evaluated by empirical security metrics such as the above guessing entropy. See Standaert et al. [?] for a complete description of this and other information metrics and their relationships in the context of SCA.

Figure ?? shows the evolutions of the gGE with the number of processed curves for noisy and small noise data, with and without application of our data-driven filtering method. With small noise data, our method leads to a slight improvement in the convergence rate toward success (e.g. gGE = 1.0 with 3200 measurements with filtering versus 4000 without) of the SCA. For noisy data, the improvement is more significant.

VII. CONCLUSION

This paper proposes a new unbiased estimator of the SN R. Being a Rao-Blackwellisation, the estimator improved linear competitors variance-wise. THis operator is idempotent so that using it on SN R gives back the same estimator. However, one intriguing possibility arises from the realization that one could extend (??) to a third (fourth etc.) order unbiased U -statistic that could yield further improved estimators of the SN R, perhaps initiating a convergence process toward optimality.

More research is needed to evaluate the cost-benefit of this idea. for some constant c near 0 to be determined. Now, the two terms Cornish-Fisher expansion of a χ 2 T (δ 2 (v)) at 1/2 is -κ/6 + O(T -1 ), where κ = 3 √ 2(T + 3δ 2 (v))/ 3 T + 2δ 2 (v). Analysis of κ shows that it is rather stable over a large range of (SN R, T ) values and that c = -2 3 /(6 √ T ) provides a good compromise .

Fig. 1 .

 1 Fig. 1. EM curves (n = 10) from an AES

Fig. 2 .Fig. 3 . 1 - 1 -

 2311 Fig. 2. (a) LN Rp(f ) for noisy measurements (b) LN Rp(f ) for small noise measurements. The gray line is the threshold value λ 0 = 1. The dots correspond to values of f where H 0 (f ) are rejected at a family-wise error rate of 5%.

TABLE I VALUES

 I OF ζ 1 BY MONTE-CARLO (MC) SIMULATIONS (2000 REPLICATIONS) ALONG WITH THOSE OBTAINED FROM THE 1 st AND 2 nd ORDER APPROXIMATIONS, FOR VARIOUS VALUES OF T AND SN R.

	SN R	Approx.	0.1	1.0	5.0	10.0	25.0	50.0 100.0 200.0
		MC	0.132 0.551 1.877 3.660 8.636 16.78 33.38 65.86
	T = 5	2 nd 1 st	0.149 0.555 1.883 3.476 8.221 16.12 32.90 63.47 0.083 0.311 1.078 2.000 4.752 9.332 18.49 36.80
		MC	0.106 0.403 1.379 2.555 6.041 11.84 23.61 46.79
	T = 10	2 nd 1 st	0.105 0.394 1.340 2.476 5.860 11.49 22.75 45.26 0.080 0.303 1.042 1.933 4.592 9.017 17.86 35.56
		MC	0.066 0.248 0.851 1.585 3.747 7.316 14.64 28.99
	T = 25	2 nd 1 st	0.065 0.247 0.843 1.560 3.696 7.251 14.36 28.58 0.059 0.224 0.766 1.420 3.372 6.619 13.11 26.10
		MC	0.046 0.174 0.595 1.103 2.612 5.119 10.23 20.32
	T = 50	2 nd 1 st	0.046 0.174 0.594 1.100 2.607 5.115 10.13 20.16 0.043 0.166 0.567 1.050 2.492 4.892 9.690 19.29

TABLE II ACTUAL

 II COVERAGE PROBABILITIES (BASED ON 5000 REPLICATIONS) OF 95% CONFIDENCE INTERVAL (??) USING THE 1 st AND 2 nd ORDER APPROXIMATION FOR ζ 1 .

							SN R
	T	n Approx. 0.1	1	5	10	25	50	100 200
	5	10 20	1 st 2 nd 1 st 2 nd	94.6 90.2 90.9 98.9 96.6 94.7 90.7 85.2 87.8 98.2 95.8 93.9	91.8 95.1 86.4 94.5	91.2 91.3 90.8 91.6 94.5 95.4 95.2 95.1 86.1 85.8 85.6 85.2 94.2 93.9 94.8 94.1
	10	10 20	1 st 2 nd 1 st 2 nd	89.9 89.4 91.2 98.1 96.8 95.8 88.2 90.7 90.9 98.9 96.3 93.4	91.7 94.9 91.7 94.1	90.1 91.9 91.3 90.8 94.2 95.0 94.4 94.3 91.3 91.7 91.4 90.6 94.1 94.3 94.2 94.1
	25	10 20	1 st 2 nd 1 st 2 nd	92.2 93.6 94.4 98.0 96.4 94.9 92.5 93.3 92.2 98.8 96.4 94.1	94.1 94.5 92.5 94.1	93.4 93.5 92.9 92.8 94.9 95.1 94.7 95.0 93.2 92.9 92.9 92.6 94.4 94.2 94.8 94.9
	50	10 20	1 st 2 nd 1 st 2 nd	94.1 93.8 93.6 97.7 96.4 94.5 93.9 94.9 94.4 98.3 96.9 94.9	94.1 94.7 93.6 94.4	94.3 94.0 94.4 94.4 95.0 94.8 95.4 95.2 94.4 94.2 94.3 94.4 95.2 95.1 95.6 95.0
	sufficient statistic. The Rao-Blackwell theorem ensures that
	this conditional expectation has a smaller variance than linear
	competitors and thus improves on them. We can quantify
	this improvement in variance by exploiting expression (i) of
	Lemma A in Serfling [?], p. 183, which yields :

TABLE III RATIO

 III OF LENGTHS OF CONFIDENCE INTERVALS BASED ON SN R AND SN R i,i+1 .

ACKNOWLEDGEMENTS

The authors would like to thank Prof. Serge B. Provost from The Western University, Ontario (Canada), for his helpful suggestion regarding the distribution of V (η 1 ). They also acknowledge the contribution of two anonymous reviewers and the associate editors for comments that led to much improvements.

APPENDIX

We prove (??). In principle, the evaluation of E V 2 (η 1 ) requires a T -fold integration over the density of η 1 . But δ 2 (v) and b(v) being functions of x, y solely we need only the joint density f (X,Y ) (x, y), where X = τ η 1 , Y = η 1 η 1 and, by scale invariance, η 1 ∼ N T (0, I). Now the marginal density of Y is easily seen to be that of a χ 2 T (0) :

y T /2-1 e -y/2 , y > 0.

As for the conditional density of (X | Y = y), we use the following argument. The random vector η 1 has the stochastic representation η 1 = √ Y U, where U = η 1 / η 1 and Y and U are independent with U being a directional vector uniformly distributed over the T -dimensional unit sphere. Hence if t denotes the unit vector τ / τ , the event


Y Aoustin

S Bellavoir

G Branchu

C Chevallereau

A Formal'sky

Ph P Lemoine

Molina

A Semi Quadruped Walking Robot: First Experimental Results

SemiQuad is walking robot prototype with a platform and two two-link legs. Thus, it is five-link mechanism. Front leg models identical motions of two quadruped's front legs, rear leg models identical motions of two quadruped's rear legs. We call this robot also as "virtual" quadruped. The mechanism is driven by four DC motors. Its control system comprises a computer, hardware servo-systems and power amplifiers. The locomotion of the prototype is planar curvet gait. In the double support our prototype is statically stable mechanism. In the single support it is unstable and under actuated system. There is no fly phase. We describe here the scheme of the mechanism, the characteristics of the drives and the control strategy.

The dynamic model of the planar walking is recalled for the double and single support phases.

The designed control strategy overcomes the difficulties appeared due to unstable and under actuated motion in the single support. In the first experiments, the designed algorithm was successfully implemented.

INTRODUCTION

It is possible to distinguish three main families of the walking robots. Firstly, the static stable robots, their motion is such that the projection of the mass center is always inside the support polygon i.e, polygon formed by the stance feet. Only slow motions are achieved in the static stability mode. Secondly, the semi dynamic stable robots, their motion is such that the Zero Momentum Point (ZMP) is inside the support surface during some time-intervals [START_REF] Vukobratovic | Biped Locomotion[END_REF].

The feet have a surface contact with the ground. The third family of the walking robots is composed of bipeds, or quadrupeds walking robots, which have sometime less than three point feet on the ground. These walking robots are under actuated mechanisms for some phases of the motion. A free rotation of the robot around the point contact or around the axe defined by the two point contacts is possible. Many papers are devoted to this kind of walking robots, for example to the bipeds, which are adapted to the human environment, see for example [START_REF] Aoustin | Stability of a cyclic biped gait and hastening of the convergence to it[END_REF]. However the quadrupeds are also interesting objects of study, see [START_REF] Furusho | Realization of Bounce Gait in a Quadruped Robot with Articular-Joint-Type Legs[END_REF]. There are several possible kinds of gait for a quadruped, for example, amble, trot, curvet [START_REF] Sutherland | Machines that walk[END_REF][START_REF] Perrin | Modelling and Control of a Quadrped for a Dynamically Stable Gait[END_REF][START_REF] Formal | On Ballistic Walking Locomotion of a Quadruped[END_REF][START_REF] Muraro | Optimal Trajectories for a Quadruped Robot with Trot, Amble, and Curvet Gaits for two Energetic Criteria[END_REF].

Usually two following half steps of a quadruped are not symmetric in the saggital plane.

The dynamical quadruped walking robots present an interesting challenge. They are useful due to their possibilities to transport payloads with fast displacement. For these reasons we decided to build a prototype, which is a semi quadruped or virtual quadruped. It has two legs, which are linked by a platform. In the saggital plane with a walking curvet gait it can be viewed as a quadruped robot using the virtual leg principle [START_REF] Sutherland | Machines that walk[END_REF], if the motions of the front legs (respectively the motions of the rear legs) are synchronized.

The article is organized as follows: Section 2 describes our virtual quadruped mechanism. Section 3 gives a description of the experimental environment. A dynamical model of the mechanism is defined in Section 4. The control strategy is explained in Section 5. Section 6 presents the first experimental results. Section 7 contains our conclusion.

DESCRIPTION OF THE MECHANISM

A photo of SemiQuad is presented in Figure 1. The prototype is composed of a platform and two identical two-link legs with knees. The legs have passive (uncontrolled) feet extending in the frontal plane. Thus, the mechanism can fulfill planar walking only. So, only 2-D motion in the sagittal plane is considering. There are four electrical DC motors with gear reducers to actuate the joints between the platform and the thighs and the knee-joints. These four actuated joints of the robot are each equipped with one encoder to measure angular position. The angular velocities must be calculated, using the signals about the angular positions. The four encoders have 2000 pts/rev and are attached directly to the motor shaft.

The absolute orientation of the support leg shin is given by potentiometer sensor, fixed to the lateral stabilisation device. It does not appear in Figure 1. The speed of response or bandwidth of each axis of the robot is determined by the transfer function of the mechanical power train (motors, gears) and the motor amplifiers that drive each motor. In the case of SemiQuad we have approximately a 53 Hz bandwidth in the mechanical portion of the axes and approximately 1.7 kHz for the amplifiers. Using dynamic simulations, we have chosen parameters of the prototype (the sizes, masses, inertia moments of the links), the convenient actuators. The lengths, masses and inertia moments of each link of SemiQuad and the parameters of the four actuators with their gear reducer are specified in All pieces of SemiQuad are designed using a parametric DAO software. Then it is easier to make some improving modifications for the prototype or to build in future a version 2 of our prototype eventually.

SemiQuad is a two-legged mechanism as an anthropomorphic biped. But unlike a biped its body is a horizontal platform and the legs are linked to the ends of this platform. Unlike a biped the front leg is always in front of the SemiQuad and the rear leg is always behind. There is no a swap of the legs as during the biped walking.

EXPERIMENTAL ENVIRONMENT

To get the measurement data and to implement the control algorithm, a dSPACE system was selected as the real-time control platform. The low-level computation, digital-to-analog and analog-to-digital conversion, as well as the user interface are provided by the dSPACE package. The control computations were performed with a sample period of 5 ms (200 Hz).

The software is developed in C language.

DYNAMIC MODEL

The models for single support and double support motions are presented in this Section.

Single Support Phase

Letthevector  with four components denote the applied torques in the haunch and knee joints. Let

* 1 2 3 4
( , , , )

      j j x j y
R ( R ,R ) be the force applied to the stance leg tip j -the ground reaction. Using the second Lagrange method, motion equations of the SemiQuad in the swing phase are obtained and they have the following well-known form for j=1or 2: j A( q )X H( q,q ) D D

( q )R       j  (1)
* j j D ( q )X H ( q,q ) 0    (2)
The two Cartesian coordinates of the platform mass center and the five orientation angles

  * 1 2 3 4
, , , , q       of the legs and the platform define vector X. Here is the symmetric, positive definite 7 matrix of kinetic energy; is the 7 vector of the centrifugal, Coriolis and gravity forces; D is a 7 A( q ) 1 7  H( q,q )   4  fixed matrix, consisting of zeros and units; is a 7 matrix; is the transposed matrix

2 j D ( q ) 2  * j D ( q ) 7  of , is a matrix. Setting null acceleration condition j D ( q ) j H ( q,q )  2 7 
(2) of a contact between the supporting leg tip and the ground implies that both abscise and ordinate of that leg tip do not change. Thus, in the single support the number of freedom degrees is five, but there are only four torques. This means that SemiQuad is under actuated mechanism in the single support motion.

Also it is unstable during this motion.

Double Support Phase

In the double support, SemiQuad has both legs on the ground. Instead of matrix equations ( 1) and ( 2) we have the following matrix equations for this phase:

(3)

1 1 2 A( q )X H( q,q ) D D ( q )R D ( q )R        2 * 1 1 * 2 2 H ( q,q ) D ( q ) X 0 H ( q,q ) D ( q )                  (4)
In the double support, SemiQuad has three degrees of freedom and four torques. In this phase, our walking robot is over-actuated. However it is static stable during double support motion.

CONTROL STRATEGY

By the analogy to the animal walk, the walk gait modeled here consists of alternating phases of single and double support. There is no a fly. SemiQuad must jump with front or rear leg to realize this walking gait. It is not possible to adopt another way without leg sliding.

(Let us note, we can define a motion of SemiQuad, using only double support phase with sliding of the leg tips along the bearing surface.) Then we have to deal with the underactuated problem in the single support. We tested this strategy in simulation [START_REF] Aoustin | Control Strategy for Dynamic Locomotion Design of Walking Quadruped[END_REF] to study its feasibility. Let us consider a current half step n, which is composed of a double support and a single support on the rear stance leg. The next half step n+1 is composed of another double support and a single support on the front stance leg. At the beginning of the current half step n, SemiQuad has to prepare its front leg tip to take off. Therefore, initially in the double support, SemiQuad transfers its mass center backward tracking suitable reference trajectories for the inter-link angles. After, also in the double support, the front leg bends in the knee and then it sharply unbends. As maximal as possible torque is applied to the front leg knee to unbend it. (At this time the front leg pushes the ground.) After we invert this torque and the front leg bends in the knee once more. At this time front leg tip looses its contact with the ground and the single support phase starts. In the single support of half step n, four torques are applied to the joints to track the reference trajectories. The reference trajectories are such that the distance between the leg tips becomes lager to the end of the single support. Due to this the front leg tip moves forward. The reference trajectory for each inter-link angle is a polynomial function of time, but before the impact, the desired inter-link angle becomes a constant. That way we are sure to get the desired final configuration of our prototype before the impact even the duration of the single support phase is not exactly the expected one. The single support phase finishes at the instant when the front leg touches the ground (with an impact).

Just after the impact of the front leg, the following half step n+1 begins. At the beginning of this half step, SemiQuad has to prepare its rear leg tip to take off. Therefore, initially in the double support of half step n+1, SemiQuad transfers its mass center forward.

Then also in the double support, the rear leg bends in the knee, after it sharply unbends and then bends once more. The rear leg tip looses its contact with the ground and the new single support phase starts with stance front leg. In the single support of half step n+1, four torques are applied to the joints to track the reference trajectories. The reference trajectories are designed such that the distance between the leg tips becomes smaller to the end of this next single support. Due to this the rear leg tip moves forward.

To track the reference trajectories we use PD -controller. However, sometime we apply constant torques as it is explained above. In Figure 3, we can see the detailed graph of the torque, which is applied to the front leg knee. This graph corresponds to only one half step of SemiQuad. The experimental robot has performed almost thirty five half steps.

Figure 1 .

 1 Figure 1. Semi Quadruped (virtual quadruped).

Figure 2 .

 2 Figure 2. SemiQuad's Scheme.

Table 1 .

 1 The gear ratio is 50. The maximal value of the torque in the output axe of each motor gear is 40 N.m.

		platform	thigh	shin	DC motor in	DC motor in
					Haunch	knee
	length in m	0.37	0.15	0.15
	mass in kg	0.57	0.47	0.40	2.82	2.74
	inertia in kg.m 2 0.32	0.84	0.06	3.25 10 -5	2.26 10 -5

Table 1 .

 1 Parameters of SemiQuad.

CONCLUSION

In this paper, we present a walking robot, which is a statically stable in the double support and unstable in the single support. Its motion is planar to illustrate a curvet gait of a quadruped robot. This robot has performed first steps. The control strategy mixes the definition of a reference trajectory and its control to the actuated joints. The main difficulties of this control are to define the switching time for each control action because the control algorithm consists of several different parts (for double supports, for single supports).

Currently we are improving our mechanism to obtain a more robust device. We are going also to construct a sensor to detect the instant when the swing leg touches the ground. In the future we hope to include also the force sensors to measure the ground reaction in the stance leg.