Louis Hawley

Rémy Rahem

Wael Suleiman

Kalman Filter Based Observer for an External Force Applied to Medium-sized Humanoid Robots

External force observer for humanoid robots has been widely studied in the literature. However, most of the proposed approaches generally rely on information from sixaxis force/torque sensors, which the small or medium-sized humanoid robots usually do not have. As a result, those approaches cannot be applied to this category of humanoid robots, which is widely used nowadays in education or research.

In this paper, we improve the external force observer in [1] to handle the case of an external force applied in any direction and at an arbitrary point of the robot structure. The new observer is based on Kalman filter formulation and it allows the estimation of the three force components. The observer is simple to implement and can easily run in real time using the embedded processor of a medium-sized humanoid robot such as Nao or Darwin-OP. Moreover, the observer does not require any change to the robot hardware as it only uses measurements from the available force-sensing resistors (FSR) inserted under the feet of the humanoid robot and from the robot inertial measurement unit (IMU).

The proposed observer was extensively validated on a Nao humanoid robot. In all conducted experiments, the observer successfully estimated the external force within a reasonable margin of error.

I. INTRODUCTION

Having a biped structure and articulated arms makes humanoid robots ideal candidates to perform a variety of manipulation and transportation tasks. However, the manipulation or transportation of an object may create instability for the humanoid robot if no measures are taken to consider the added force, which is generated by the transported object, and adjust the walking gait. Moreover, physical interaction with a human, such as collaboratively transporting an object, also generates interaction force that should be taken into account to prevent instability issues. In those cases, it is primordial to accurately estimate the force acting on the humanoid robot.

As an example, the transportation of a heavy object on a rolling cart, as depicted in [START_REF] Rioux | Humanoid navigation and heavy load transportation in a cluttered environment[END_REF], [START_REF]Autonomous SLAM based humanoid navigation in a cluttered environment while transporting a heavy load[END_REF], can be greatly improved if the external force generated by the cart and the transported object can be accurately measured or estimated. In [START_REF]Autonomous SLAM based humanoid navigation in a cluttered environment while transporting a heavy load[END_REF], a planning and navigation framework for the transportation of an object on a rolling cart have been proposed. The framework generates a feasible trajectory for the humanoid robot and the cart while considering the physical limitations of the robot. In the proposed implementation, the robot is only able to differentiate between light and heavy loads by executing a turning in place motion with the cart, and it

F ext M c ẍc M c ÿc M c zc ⌧ y ⌧ x M c g
Fig. 1. An external force applied to a Nao robot then uses the appropriate set of motion primitives in the planning phase. However, it would be useful to have a better estimation of the external force exerted on the humanoid robot for: I)-an optimal motion planning, II)-more stable locomotion by integrating the external force into the pattern generation module.

Moreover, for a cooperative transportation task of an object with two humanoid robots, as presented in [START_REF] Rioux | Cooperative SLAM-Based Object Transportation by Two Humanoid Robots in a Cluttered Environment[END_REF], [START_REF]Cooperative Vision-Based Object Transportation by Two Humanoid Robots in a Cluttered Environment[END_REF], it would be very useful to measure the force exchanged through the object in order to maintain the equilibrium of the whole system.

A. Relevant Works

There are multiple approaches to estimate external forces exerted on a humanoid robot. These approaches can be roughly separated into three main categories : sensor-based approaches, model-based approaches and alternative methods such as machine learning-based methods. Each category is summarized with some examples in the following subsections.

1) Sensor-based Approach: The sensor-based approach directly measures the external force applied to a humanoid robot by using the robot's sensors. Such a method was used in [START_REF] Harada | ZMP Analysis for Arm / Leg Coordination[END_REF], [START_REF]Pushing manipulation by humanoid considering two-kinds of ZMPs[END_REF], where a pushing manipulation task was performed by a HRP-2 robot. In that case, the external force applied to the humanoid robot can be directly measured at all time using the six-axis Force/Torque (F/T) sensors embedded in the robot's wrists. Then, by considering the external force, the authors were able to apply the required compensation on the torso position to ensure the stability of the robot while pushing the object.

A similar approach has also been proposed by the same authors to enable a humanoid robot to lift and transport a heavy object [START_REF] Harada | A humanoid robot carrying a heavy object[END_REF]. The sensor-based approach obviously provides great performance as the force can be accurately measured at all times. However, F/T sensors are not available on all robots and it would be difficult to integrate them on smaller and cheaper robots due to the size and cost of such sensors.

2) Model-based Approach: We define model-based approach as the methods that estimate the external force acting on a humanoid robot by using a dynamic model of the robot and its actual state vector. For instance, the method presented in [START_REF] Kaneko | Disturbance observer that estimates external force acting on humanoid robots[END_REF], where the magnitude of an external force is estimated using a disturbance observer. The observer uses measurements from the robot's IMU and from F/T sensors located in the ankles to estimate the magnitude of a strong force such as a kick to the robot chest. In order to estimate the disturbance, the authors used a simplified dynamic model of the humanoid where the acceleration of its Center of Mass (CoM) is a function of the forces applied to the robot, which are: the gravity, the floor reaction force and the unknown external forces. Since the acceleration of the robot's CoM is known using IMU measurements and the floor reaction force can be directly measured with the F/T sensors, their results show that a strong impact force can be estimated somewhat accurately using that model. Another interesting method that can also estimate the external moment is proposed in [START_REF] Benallegue | Model-based external force/moment estimation for humanoid robots with no torque measurement[END_REF], that method relies on a model-based estimator to fuse the data of the robot's sensors (F/T sensors and IMU) and the whole body dynamics. An advantage of this approach over the previously presented sensor-based approach is that the external force/moment do not have to be applied at a particular point in order to be estimated. However, these methods still require F/T sensors.

A similar approach has been presented in [START_REF] Kaneko | Slip observer for walking on a low friction floor[END_REF] to enable a humanoid robot to walk on low friction floor. Using the linear inverted pendulum model (LIPM) to represent the humanoid robot, the authors predicted the ground reaction force and compared it with the measurements from F/T sensors in the ankles of a HRP-2 robot. The difference between the two values was then attributed to a slip force that can be easily derived. Although a slip force is also a strong force of short duration, which is very different from the force encountered in interaction tasks between several robots or a robot and a human, these works are interesting as they show that it is possible to get a good estimation of an external force by only using a limited number of sensors and a simple dynamic model of the robot. Such an approach would be interesting for medium-sized humanoid, such as Nao humanoid robot, as they usually have a small number of sensors and limited computational power. Nao [START_REF] Gouaillier | Mechatronic design of NAO humanoid[END_REF] is a medium-sized humanoid robot manufactured by SoftBank Robotics.

In [START_REF] Mattioli | Interaction Force Reconstruction for Humanoid Robots[END_REF], the Nao robot's limbs are modeled as kinematic chains attached to the pelvis of the robot, and a relation between a force acting on each link of the robot and the joint torques is established, a method commonly used for serial manipulators. Since Nao's joints do not possess torque sensors nor do they provide a torque estimation, the authors had to define an actuator model to estimate the motor torque from the steady-state error. Then, the problem is transformed into isolating the torque caused by gravity to obtain the portion of the torque that is the result of the external force. Although this approach can be quite powerful as the magnitude of the force along with the contact point can be found, there are multiple drawbacks in the proposed implementation. For instance, it uses the motors steadystate error to estimate the actuator torque. That means the proposed algorithm is mainly able to estimate quite a big external force.

Rotella et al. [START_REF] Rotella | Humanoid Momentum Estimation Using Sensed Contact Wrenches[END_REF] presented some formulations to estimate the CoM position of a robot and to correct inaccuracies in the robot links model. With regard to our work, the most interesting contribution is an external wrench estimator that uses a linear and angular momentum model of the robot instead of the simplified LIPM dynamic. However, the approach is limited to humanoids equipped with F/T sensors at the end-effectors.

Another model-based method to estimate the magnitude of an external force applied to a humanoid robot is presented in [START_REF] Hawley | External Force Observer for Mediumsized Humanoid Robots[END_REF]. Its main advantage is that it does not require using F/T sensors but instead it uses measurements from the robot's force-sensing resistors (FSR) sensors. The approach has been successfully validated on a Nao robot and pointed out that an external force applied to the robot can be efficiently estimated with an acceptable precision. However, the main limitation of that method is that the formulations are only valid for an external force applied in the horizontal plane.

3) Alternative Methods: Machine learning based approaches have also proven to be useful in estimating an external force applied to medium-sized/small humanoid robots.

In [START_REF] Berger | Inferring guidance information in cooperative human-robot tasks[END_REF], perturbations transmitted to a Nao robot are detected by comparing the actual sensor measurements with the predicted ones obtained with a probabilistic model. Using this method, the authors were able to successfully detect and provide an estimation of the perturbation magnitude. However, the external force that generates the perturbation has never been quantified as the algorithm can either provide the difference between the sensors (IMU and FSR) measurements and the predicted values or an equivalent footstep length. Although it is certainly possible to formulate an estimation of the force from these results, such analysis has not been conducted.

In [START_REF] Carpentier | Center-of-Mass Estimation for a Polyarticulated System in Contact -A Spectral Approach[END_REF], the authors present a complimentary filter formulation to estimate the CoM of a robot or human. The proposed filter uses a spectral approach to analyze the different methods used to estimate the position of CoM (i.e. motion capture, contact forces and moments). Then, the authors propose to fuse the data by optimally using each sensor. Obviously, this approach requires a diverse set of sensors to span the whole frequency range, a condition that is often not satisfied in the case of small humanoid robots.

The main contribution of this paper is the design of a model-based observer that is able to estimate the three components of an external force applied to a humanoid robot. The force is projected onto the CoM frame, and therefore it can be directly used, for instance, in a pattern generator. The observer design and implementation architecture are given in Section IV.

This observer was also validated on a Nao humanoid robot through several experimental scenarios. The results of these experiments are analyzed in Section VI.

II. ZMP WITH AN EXTERNAL FORCE

To study the impact of an external force on the Zero-Moment Point (ZMP) formulation, we consider the two following cases:

A. An External Force Directly Applied to the Robot CoM This case study is illustrated in Fig. 2. The ZMP can be found by considering the forces acting on the simplified LIPM model, which are the gravity, the inertial force and the external force. Let the ZMP on the horizontal ground be given by the following vector 1 :

! " F ext M c ẍc M c ÿc M c zc ⌧ y ⌧ x M c g
p = p x p y T (1)
where p x and p y are respectively the projections of the ZMP on the x and y axis.

To compute p, one can use the following formula:

p = N n × τ o (f o |n) (2)
where the operators × and (.|.) design respectively the cross and scalar products, and 1 As a general rule, bold parameters are vectors or matrices

• N is a constant matrix

N = 1 0 0 0 1 0 (3)
.

• the vector n is the normal vector on the horizontal ground n = 0 0 1 T .

• The vector f o is the result of the gravity, inertial and external force:

f o = M c g -M c Ẍc + F ext (4)
where g denotes the gravity acceleration (g = -gn), M c is the total mass of the humanoid robot, and Ẍc = [ẍ c ÿc zc ] T is the acceleration of the CoM.

F ext = [F x F y F z ]
T is the external force applied to the humanoid robot.

The normal force is given by the following expression:

(f o |n) = f o n = -M c g -M c zc + F z (5)
• τ o denotes the moment of the external force and the moment of the force f o about the origin of the fixed world frame, it is given by the following expression:

τ o =X c × f o (6)
where X c = x c y c z c T is the Cartesian position of CoM. Thus, we obtain:

p =   p x p y   =    x c + Mc zc f o n ẍc -zc f o n F x y c + Mc zc f o n ÿc -zc f o n F y    (7)
B. An External Force Applied to a Point of the Robot

! " F ext ⌧ ext y ⌧ ext z ⌧ y ⌧ x M c g ⌧ ext x
Fig. 3. An external moment is applied to the CoM of LIPM model as a result of an external force applied to a point of Nao robot, the accelerations are omitted for clarity An external force applied to a point of the robot structure generates an additional external moment applied to the CoM of LIPM model as shown in Fig. 3. For instance, this force could be the result of the robot pushing, pulling or lifting an object or interacting/collaborating with a human.

Similarly to the previous case, the resultant force applied to the LIPM model is:

f o = M c g -M c Ẍc + F ext (8)
However, the moment τ o in this case becomes:

τ o =X c × f o + τ ext (9)
By applying the formula in (2), we obtain:

p =   p x p y   =    x c + Mc zc f o n ẍc -zc f o n F x -1 f o n τ ext y y c + Mc zc f o n ÿc -zc f o n F y + 1 f o n τ ext x    =     x c + Mc zc f o n ẍc -zc f o n F x + τ ext y zc y c + Mc zc f o n ÿc -zc f o n F y - τ ext x zc     (10)
It is clear from [START_REF] Benallegue | Model-based external force/moment estimation for humanoid robots with no torque measurement[END_REF] that, without additional sensors such as a six-axis F/T sensor, the components of external force and moment cannot be estimated. However, by defining the following virtual force:

F ext =   Fx Fy Fz      F x + τ ext y zc F y - τ ext x zc F z    (11)
One can figure out that, from ZMP point of view, the impact of an external force F ext applied to a point of the robot is the same as a virtual force F ext applied directly to the CoM of LIPM model. Therefore, we only consider the latter case in the sequel of the paper.

III. MEASUREMENTS

Using the robot available sensors (FSR, IMU and joint encoders), the ZMP, CoM position and acceleration, and ground reaction force can be estimated. Those variables will play the role of output (observation) signals for the Kalman filter based observer. Note that, in this work, all measurements are unfiltered (raw) signals.

A. Ground Reaction Force

The ground reaction force can be directly estimated using the FSR sensors. In the case of Nao robot, for instance, each foot possesses 4 FSR sensors with a force range between 0 and 25N, and the ground reaction force f o n is given by:

f o n = -       4 i=1 F R i right foot + 4 i=1 F L i left foot       (12)
where F * i is the force measured at the i th FSR.

B. ZMP

Recall that the ZMP and Center of Pressure (CoP) coincide when the robot is in a stable configuration [START_REF] Vukobratović | Zero-Moment Point -Thirty Five Years of Its Life[END_REF]. Thus, the ZMP can be computed as follows:

p x = x CoP = 4 i=1 F R i x R i + 4 i=1 F L i x L i |f o n | p y = y CoP = 4 i=1 F R i y R i + 4 i=1 F L i y L i |f o n | (13)
where f o n is defined in [START_REF] Gouaillier | Mechatronic design of NAO humanoid[END_REF], and x * i and y * i are the x and y positions of the i th FSR in the foot frame.

C. CoM Acceleration

The acceleration of CoM can be directly extracted from the IMU, which is a standard part of a humanoid robot sensors. As for Nao, the IMU provides measurements from a three-axis accelerometer and a two-axis gyroscope. As in [START_REF] Hawley | External Force Observer for Mediumsized Humanoid Robots[END_REF], we assume that the CoM coincides with the IMU. Therefore, the accelerometer data can be used to approximate the acceleration of the CoM.

D. CoM Position

The CoM position in the support foot frame can be easily obtained by a simple forward kinematics.

IV. EXTERNAL FORCE OBSERVER

The first component of the external force applied to the humanoid robot that can be directly estimated from the FSR information is F z using [START_REF]Cooperative Vision-Based Object Transportation by Two Humanoid Robots in a Cluttered Environment[END_REF]. Then, the horizontal components F x and F y can be estimated using [START_REF]Pushing manipulation by humanoid considering two-kinds of ZMPs[END_REF].

A. Estimating F z

To estimate F z , Eq. ( 5) is discretized with a sampling time T to obtain the following linear dynamic system:

Z(k + 1) = A Z(k) + B u z (k) Y z (k) = C z Z(k) + v z (k) (14)
where:

• The state vector:

Z(k) = z c (kT ) żc (kT ) zc (kT ) F z (kT ) Ḟz (kT ) T •
The input signal:

u z (k) = ... z c (kT ) Fz (kT )
T where ... z c (kT ) is the acceleration derivative which is usually called the jerk.

• The output signal:

Y z (k) = z c (kT ) zc (kT ) f o n (kT ) + M c g T •
The model constant matrices:

A =       1 T T 2 2 0 0 0 1 T 0 0 0 0 1 0 0 0 0 0 1 T 0 0 0 0 1       , B =        T 3 6 0 T 2 2 0 T 0 0 T 2 2 0 T        C z =   1 0 0 0 0 0 0 1 0 0 0 0 -M c 1 0   (15)
• v z is the observation noise and considered as a Gaussian white noise with a covariance matrix R z . Thus, v z ∼ N (0, R z ).

To formulate the estimation problem as a Kalman filter, we reformulate the dynamic system in ( 14) as follows:

Z(k + 1) = A Z(k) + ω z (k) Y z (k) = C z Z(k) + v z (k) ( 16
)
where ω z (k) = B u z (k) becomes the process noise, which is supposed to be a Gaussian white noise with a covariance Q, i.e. ω z ∼ N (0, Q), which is independent of v z . The matrix Q is given by:

Q = B σ 2 ... z c 0 0 σ 2 Fz B T (17)
where σ2 ... As above-mentioned F x and F y can be estimated from [START_REF]Pushing manipulation by humanoid considering two-kinds of ZMPs[END_REF]. However, because of the similarity between the two components of ZMP (p x and p y ), we only consider the estimation of F x as F y can be estimated in an analogous way.

At first, p x in ( 7) is discretized with a sampling time T to obtain the following linear dynamic system:

X (k + 1) = A X (k) + ω x (k) Y x (k) = C x X (k) + v x (k) (19)
where:

• X (k), ω x (k), v x , A and Q are defined in the same way as in Section IV-A by replacing F z by F x and z c by x c . • The output signal:

Y x (k) = x c (kT ) ẍc (kT ) p x (kT ) T • The output matrix: C x =    1 0 0 0 0 0 0 1 0 0 1 0 Mc ẑc(kT ) f o n (kT ) -ẑc(kT ) f o n (kT ) 0    ẑc (kT ) and f o n (kT ) = -M c g -M c
ẑc (kT ) + Fz (kT ) are computed using [START_REF] Mehra | On the identification of variances and adaptive kalman filtering[END_REF]. By applying Kalman filter on [START_REF] Flacco | Contact detection and physical interaction for low cost personal robots[END_REF], an estimation of the state vector X is obtained.

The observer implementation architecture is summarized in Fig. 4.

IMU FSR

CoM Acceleration ZMP CoM position Ground reaction force Kalman filter applied on Eq ( 16)

ẑc (kT ) ẑc (kT ) Fz (kT )

Kalman filter applied on Eq [START_REF] Flacco | Contact detection and physical interaction for low cost personal robots[END_REF] Fx (kT )

Kalman filter

Fy (kT )

Analogous way Joint encoders Regarding the observer parameters, the best results have been obtained with:

σ 2 ... x = σ 2 ... y = σ 2 ... z = 10 3 σ 2 Fx = σ 2 Fy = σ 2 Fz = 10 3 R z =   0.01 0 0 0 1 0 0 0 1   , R x = R y =   0.01 0 0 0 1 0 0 0 0.01   (20)
Moreover, the experimental data as well as C++/ROS implementation and the MATLAB code that has been used to produce the results in this paper can be downloaded here 2 . These data and code could be helpful for the humanoid robotics community with the objective of benchmarking and comparing with other/their approaches.

The observer has been experimentally validated on the Nao robot through three scenarios, and it has been also used to accomplish human-robot interaction tasks. In this experiment, a force in the xz plane is applied through a spring-based force meter as shown in Fig. 5(a). As the torso frame is close to the CoM position of the robot, the magnitude of resulting moment applied to the LIPM model is relatively small. Therefore, applying a force to the robot torso is practically equivalent to applying it to its CoM.

The F x and F y components of the force can be easily obtained by the following formula:

F x = F ext cos(θ) F z = F ext sin(θ) (21)
where F ext is the force magnitude measured with the force meter, and θ is the angle between the force meter and the

x axis, the positive direction of rotation is defined to be clockwise.

The experiment consists of applying a force with a constant magnitude of approximately 2N in the xz plane but with different values of the angle θ.

Fig. 5(b) shows the overall performance of the observer, it points out that the observer reasonably succeeds in estimating the three components of the external force. However, note that the main difficulty in this experiment was keeping a force with constant magnitude by only visually monitoring the force indicator. Moreover, the angle (θ) was also visually monitored using a protractor.

We can also notice that the observer detects non-zero force along y axis, this is mainly because of the difficulty of keeping the force meter in the vertical plane. The setup for the second and third scenarios are respectively presented in Fig. 6(a) and Fig. 7(a), and the attached mass m o =2.5 pounds (1.13 kg). In both experiments, a force of around 11 N has been applied to the robot right arm, first along the y axis in Fig. 6(a) and then along the z axis in Fig. 7 By applying inverse dynamics to the kinematic chain from the point at which the force is applied to the robot CoM, the wrench vectors in the second and third experiments as well as the virtual forces, as defined in [START_REF] Kaneko | Slip observer for walking on a low friction floor[END_REF], are respectively: Those experiments are very similar to the experiments conducted in [START_REF] Mattioli | Interaction Force Reconstruction for Humanoid Robots[END_REF], they are used as a benchmark to compare with their results. Fig. 6(b) and Fig. 7(b) show that the proposed observer succeeds in estimating the corresponding virtual external force, this is because the estimated force is reasonably close to the expected one.

W 2 = F ext 2 τ ext 2 ≈         0.0 -11.12 0.0 -0.05 0.0 -0.08         =⇒ F ext 2 ≈   0.0 -10.96 0   W 3 = F ext 3 τ ext 3 ≈         0
The authors in [START_REF] Mattioli | Interaction Force Reconstruction for Humanoid Robots[END_REF] point out that the estimated force, by their method, has an important magnitude instead of a small one when the mass is removed. They argue that their main sources of error are the unmodeled static friction and the lowcost hardware of Nao robot, for instance, plastic joint gears. On the other hand, the estimated force by our observer after removing the mass has a small magnitude of around 0.5 N, which is an order of magnitude smaller than the results in [START_REF] Mattioli | Interaction Force Reconstruction for Humanoid Robots[END_REF]. Furthermore, those results also show the robustness of our observer to the sources of error identified in [START_REF] Mattioli | Interaction Force Reconstruction for Humanoid Robots[END_REF].

C. Compliant Human-Robot Interaction

As an application of our observer, we have implemented an interactive walking based on the interaction between the robot and a human as shown in Fig. 8. In this scenario, the robot walks in the direction of the detected force (forward or backward) if the force magnitude exceeds a predefined threshold, and stops if the magnitude of detected force is less than the threshold. A second application is a compliant control of the robot arms as a function of the detected force in the sagittal plane (x axis) as follows:

x h = x 0 h + 1 k F x (23)
where x h and x 0 h are respectively the desired position of the hand and its initial position along the x axis. k is a virtual spring constant and F x is the detected force along x axis. A snapshot of the experiment is given in Fig. 9.

VII. CONCLUSION

In this paper, we introduced a Kalman filter based observer to estimate the three components of an external force applied to a medium-sized/small humanoid robot. In contrast to conventional methods that are mainly based on the use of expensive 6-axis force/torque sensors, the proposed observer uses measurement from force-sensing resistors located under the feet of the robot and the information from the robot's IMU sensor. This observer was extensively tested on a Nao robot and the results revealed that it overcomes the limitations of the related approaches in the literature, those limitations are mainly caused by the low-cost hardware components of the robot.

However, there is still room for improvement, for instance, the covariance matrices R and Q in the observer implementation have been experimentally identified, however, they can be estimated using some existing methods in the literature, such as [START_REF] Mehra | On the identification of variances and adaptive kalman filtering[END_REF]. Moreover, even though the objective of our observer is to estimate an equivalent force, called virtual force, that is applied to the CoM frame, a combination with methods for identifying the contact point, such as [START_REF] Flacco | Contact detection and physical interaction for low cost personal robots[END_REF] or [START_REF] Flacco | Residual-based contacts estimation for humanoid robots[END_REF], would provide the external force components at the contact point and in the world frame by a simple inverse dynamics calculation.

Future work will focus on integrating the external force observer into the pattern generator module and humanhumanoid robot collaborative tasks.

Fig. 2 .

 2 Fig. 2. An external external force applied to the CoM of a Nao robot

z c and σ 2 Fz

 2 are respectively the covariance of ... z c and Fz . By applying Kalman filter on (16), an estimation of the state vector is obtained: Ẑ(k) = ẑc (kT ) ẑc (kT ) ẑc (kT ) Fz (kT ) Fz (kT ) T (18) B. Estimating F x and F y

Fig. 4 .

 4 Fig. 4. Block diagram of the external force observer F ext = [ Fx Fy Fz] T

  Estimating the three components of Fext, the expected force is in green dashed line.

Fig. 5 . 1 :

 51 Fig. 5. Experiment 1: Applying a fixed-magnitude force in the x -z plane.

B

  . Applying a Force to a Point of the Robot Arm

  Estimating the three components of Fext, the expected force is in green dashed line.

Fig. 6 . 2 :

 62 Fig. 6. Experiment 2: Applying a force along y axis

  (a).

Fig. 7 .

 7 Fig. 7. Experiment 3: Applying an arbitrary external force

1 2 Fig. 8 . 2 Fig. 9 .

 12829 Fig. 8. Interactive walking by applying a sagittal guidance force to the robot hands

https://goo.gl/VkhejY

ACKNOWLEDGMENT

This research is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) under the grant RGPIN-2012-419406.


Denis Migdal
email: denis.migdal@ensicaen.fr

Christophe Rosenberger
email: christophe.rosenberger@ensicaen.fr

Analysis of Keystroke Dynamics For the Generation of Synthetic Datasets

Keywords: Keystroke, Synthetic dataset, Data Analysis

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Keystroke dynamics (KD) [START_REF] Giot | Keystroke dynamics overview[END_REF] is a behavioral biometric modality that allows the authentication of individuals through their way of typing a password or a free text on a keyboard. It is a behavioral biometrics which has the advantage of not requiring additional sensor than the keyboard. This biometric modality also allows continuous authentication through time [START_REF] Mondal | A study on continuous authentication using a combination of keystroke and mouse biometrics[END_REF][START_REF] Li | Enhanced free-text keystroke continuous authentication based on dynamics of wrist motion[END_REF].

User authentication with keystroke dynamics is generally done in real time (i.e., online) in a real world system. Scientists working on keystroke dynamics do not analyze the performance of their sys-tem in an online way (i.e, by asking users to authenticate themselves in real time and to impersonate other users). Indeed, they work in an offline context by using samples previously collected by other researchers, and stored in a benchmark dataset. A complete list of available keystroke dynamics datasets has been made in [START_REF] Monaco | Public keystroke dynamics datasets[END_REF][START_REF] Giot | A review on the public benchmark databases for static keystroke dynamics[END_REF]. As it can be seen, most of datasets have less than 200 individuals and few samples for each user. The collection of such datasets is very time consuming, this is the main reason why there is not more very large datasets like for the face modality as for example [START_REF] Learned-Miller | Labeled faces in the wild: A survey[END_REF].

In this paper, we present a study whose objective is to model real KD data in order to generate synthetic KD datasets. This approach has been used for the digital fingerprint modality with the SFINGE software [START_REF] Cappelli | Sfinge: an approach to synthetic fingerprint generation[END_REF] as their collection and distribution are regulated in many countries. We believe the KD model will be able to help the research community to create a new dataset of higher quality than the existing ones. We think this work is important, because it is known that KD studies are not fair as (i) acquisition protocols are different between studies [START_REF] Giot | Unconstrained keystroke dynamics authentication with shared secret[END_REF]; (ii) there is not always a comparative study [START_REF] Killourhy | Should security researchers experiment more and draw more inferences?[END_REF] when authors propose new algorithms; and (iii) there are not always a valuable statistical evaluation [START_REF] Killourhy | Should security researchers experiment more and draw more inferences?[END_REF]. Our work contributes to solve these problems.

The paper is organized as follows. Section II is dedicated to some background information on Keystroke dynamics and existing studies for this biometric 1 modality. We present in section III the definitions and the components of the KD model we propose in this study. Section IV concerns the evaluation of the proposed KD model on two real KD datasets and 4 matching methods. We analyze in section V the keystroke dynamics data of these two datasets to set the parameters of the proposed KD model. Section VI is dedicated to the validation process of the KD model showing the capability to generate similar keystroke dynamics data. Last, section VII concludes this work and gives some perspectives.

Related works

A keystroke dynamic system (KDS) is composed of two main modules: the enrollment and the verification modules. Each user must enroll himself/herself in the KDS which computes a biometric reference given multiple samples (i.e., several inputs of the password) acquired during the enrollment step. For each input, a sequence of timing information is captured (i.e., time when each key is pressed or released) from which some features are extracted (i.e., latencies and durations) and used to learn the model which characterizes each user. During a verification request, the claimant types his/her password. The system extracts the features and compares them to the biometric reference of the claimant. If the obtained distance is below a threshold, the user is accepted, otherwise he/she is rejected.

First works on KD have been done in the eighties [START_REF] Gaines | Authentication by keystroke timing: some preliminary results[END_REF], although the idea of using a keyboard to automatically identify individuals has first been presented in 1975 [START_REF] Spillane | Keyboard apparatus for personal identification[END_REF]. In the preliminary report of Gaines et al. [START_REF] Gaines | Authentication by keystroke timing: some preliminary results[END_REF], seven secretaries typed several paragraphs of text and researchers showed that it is possible to differentiate users with their typing patterns. Since then, several studies have been done, allowing to decrease the quantity of information needed to build the biometric reference, while improving the performances [START_REF] Giot | Unconstrained keystroke dynamics authentication with shared secret[END_REF][START_REF] Umphress | Identity verification through keyboard characteristics[END_REF][13][START_REF] Revett | A machine learning approach to keystroke dynamics based user authentication[END_REF][START_REF] Lee | Retraining a keystroke dynamics-based authenticator with impostor patterns[END_REF]. However, most studies are not comparable because they use different datasets or protocols [START_REF] Giot | Unconstrained keystroke dynamics authentication with shared secret[END_REF][START_REF] Killourhy | Should security researchers experiment more and draw more inferences?[END_REF].

Keystroke dynamics generative model

First, we define many terms to build the proposed model:

• Digraph: D = [C 0 , C 1 ]
, array of two characters.

• DigraphTime: We propose in this paper a generative keystroke dynamics model. We explain its different components (see also Figure 1 and2):

DT D = [d 0 ,
• d 0 = d 2 -d 4 ; • d 0 = d 1 -d 3 ; • d 1 = d 2 -d 5 ; • d 3 = d 4 -d 5 ; • d 0 ≥ 0 • d 1 ≥ 0 • d 5 ≥ 0 • Text: T n = {D i } i∈ 0,n , an array of n Di- graphs D i . A text T n is said consistent if ∀i ∈ 0, n , D i-1 [1] = D i [0].
• DigraphGen: DG D () = DT , generates a Di-graphTime for a given Digraph.

• User: U (T n ) = K, generates a Keystroke dynamics from a given Text. A User is composed of a set of DigraphGen.

• DigraphGen 2 : DG 2 D () = DG D , generates a DigraphGen for a given Digraph.

• UserGen: U G() = U , generates a User. A UserGen is composed of a set of DigraphGen 2 .

• DigraphGen 3 : DG 3 (D) = DG 2 D , generates a DigraphGen 2 for a given digraph.

Evaluation

In the scope of this study, two datasets are used : Data1 (GREYC-KeyStroke Dataset [START_REF] Giot | Greyc keystroke: a benchmark for keystroke dynamics biometric systems[END_REF]) & Data2 (GREYC-Web based KeyStroke dynamics dataset [START_REF]Web-Based Benchmark for Keystroke Dynamics Biometric Systems: A Statistical Analysis[END_REF]). These datasets are listed in [START_REF] Monaco | Public keystroke dynamics datasets[END_REF] and are enough representative of real world keystroke dynamics data. From these two datasets, the following sub-datasets are created:

• Data{1,2}D: DigraphTimes for each User and Digraph extracted from Data{1,2}.

• Data{1,2}K: a fixed Text's Keystrokes for each User, extracted from Data{1,2}. Two fixed Texts being present in Data2, extracted subdatasets are thus Data2K1 and Data2K2.

The comparison between keystroke dynamics samples is evaluated in this paper through 4 distance functions:

We suppose to compute a distance between two templates K A and K B , these functions are defined as follows:

• Blesha [START_REF] Bleha | Computer-access security systems using keystroke dynamics[END_REF]: We suppose that the template K A is associated by µ the average value of biometric samples:

ST AT 1 = (K B -µ) t (K B -µ) ||K B ||.||µ|| (1)
• Hocquet [START_REF] Hocquet | User classification for keystroke dynamics authentication[END_REF]: We suppose that the template K A is associated by µ and σ the average value of biometric samples and the standard deviation:

ST AT 2 = 1 - 1 n n i=1 e - |K B (i)-µ i | σ i (2)
• Monrose [START_REF] Monrose | Authentication via keystroke dynamics[END_REF]: The function is given as follows:

ST AT 3 = n i=1 (K B (i) -K A (i)) 2 (3)
• BioHashing: This algorithm is a template protection scheme [START_REF] Teoh | Biohashing: two factor authentication featuring fingerprint data and tokenised random number[END_REF] where the biometric template is projected given a key and quantized. The comparison is realized with the Hamming distance. We apply this protection scheme and compare the templates in the transformed domain.

In the scope of this paper, the BioHashing and Monrose distances between a template (sample) and a set of templates (references) are computed as the minimal distance of the sample with each template in the reference gallery.

Analysis of keystroke dynamics data

As previously seen, generating a keystroke dynamics template from a given Text T n consists in generating an array of DigraphTime, i.e. generating 6 * n durations. To be able to generate a keystroke dynamics similar to that one user could type, these 6 * n durations have to be transformed into a set of assumed independent variables which laws and parameters can then be estimated for a user, enabling their random generation, and thus the computation of a synthetic keystroke dynamics corresponding to a given user. In the scope of this paper, only the linear (in)dependency of variable is considered.

Variables (in)dependency

Linearly correlated variables can be transformed into a set of non-linearly correlated variables, through PCA (Principal component analysis), first introduced by Pearson in 1901 [START_REF] Kpfrs | On lines and planes of closest fit to systems of points in space[END_REF]. However, we show that durations are not strongly correlated between them, and thus, in the scope of this article, assume them to be independent. Even if the usage of PCA is irrelevant in such a case, its first step enables the computation of the inter-correlations of two variables by the computation of a correlation matrix. In a correlation matrix C = {C i,j } {i,j}∈ 0,n 2 , C i,j is the linear correlation between the variables i and j. A correlation matrix C = {C i,j } {i,j}∈ 0,n 2 , with C i,j the linear correlation between the variables i and j, is computed as follows:

1. Given a matrix M = {M k } k∈ 0,K of K entries M k = {M k,i } i∈ 0,n , with M k,i the realization of the variable i for the entry k.

M = {

M k,i -µi σi } i∈ 0,n ,k∈ 0,K where µ i is the mean of {M k,i } k∈ 0,K , and σ i , its standard deviation.

3. C = 1/K * M T * M
In the scope of this paper, we arbitrary consider that sets with less than a certain number of elements cannot provide pertinent results, and are thus discarded. We used 23 as arbitrary value in this paper (cf section 5.2).

To qualify presence of specific correlations between two variables i, j inside m subsets of entries, m correlations matrix C l , l ∈ 0, m are computed from such subsets. Each element C i,j of the final correlation matrix C is then computed as the mean of each C l i,j

: C i,j = 1 m Σ m-1 l=0 C l i,j .
If each subset corresponds to, e.g. a User, M will be said, in this paper, "splitted by User", and C will qualify the presence of Userspecifics correlations across all Users.

To qualify presence of the same correlations between two sets of variables {i x } x∈ 0,m , {j x } x∈ 0,m , of length m, entries are splitted in m sub-entries M m * k+x = {M k,ox } o∈{i,j} . The correlation matrix C is then computed from M'. If each x corresponds to, e.g. a Digraph, M will be said, in this paper, "merged by Digraph", and C will qualify the presence of non-Digraphs-specifics correlations across all Digraphs.

Laws followed by Variables

Once the variables assumed independent, or transformed in such a way, laws followed by each variable are searched through the following process:

1. Given the realizations of a variable X, and a law law p with unknown parameters p;

2. Estimate p from the median, mean, min, max, or/and standard deviation of X;

3. Compute the χ 2 (X, law, p) score qualifying the capacity of X to match the values that would be expected if X follows law p through a χ 2 test.

The χ 2 test qualifies the capacity of a set of observed values to match a set of expected values. The χ 2 test returns χ 2 (X, law, p) = 1 -α, in which α is the p-value, i.e. the probability to obtain the same 1 -α score if X follows law p . If the p-value is below an arbitrary threshold (s.a. 0.05), the hypothesis "X follows law p " can then be rejected.

However, in the scope of this paper, we do not aim to reject the hypothesis, but to select laws that seem to best represent X. The χ 2 (X, law, p) score can then be seen as a score of distance between observed values of X, and the expected values. We compute χ 2 (X, law, p) as follows:

1. Let Card(S) be the cardinal of S;

2. Let a%b be the rest of the division of a by b;

3. IR is divided in n = Card(X)/5 subspaces E i , i ∈ 0, n , each expected to contain 5 el- ements of X. E n-1 is expected to contain Card(X)%5 elements of X if 5 Card(X); 4. Let X i = X ∩ E i ; 5. Let Card(E i ) = 5, and
Card(E n-1 ) = Card(X)%5 if 5 Card(X); 6. Let Sum = Σ n-1 i=0 (Card(E i ) - Card(X i )) 2 /Card(E i ).
7. Let cdf f be the cumulative distribution function of the law χ 2 of freedom f; 8. χ 2 (X, law, p) = cdf n-Card(p)-1 (Sum).

In the scope of this paper, we arbitrary consider that sets with less than a certain number of elements cannot provide pertinent results, and are thus discarded. We used 23 as arbitrary value in this paper, thus having at least 5 subspaces and 100 users.

From tested laws, the best 3 are selected i.e. the 3 laws that minimize the score of distance χ 2 (X, law, p). These laws are tested with and without exclusion of aberrant values (here, X values that differ from ±3σ from the median value of X). When the parameters p have different estimations, only the one that minimizes the score of distance is kept. A set of 19 laws have been tested in this paper:

• arcsine

• beta • betaprime • chi • chisquare • cosine • erlang • exponential • f • gamma • gumbel • laplace • logistic • lognormal • uniform • normal • rayleigh • t • triangular
To qualify the capacity of n subsets of X, X i , i ∈ 0, n , to follow a same law law, but each with different parameters p i , the score of distance χ 2 (X, law) is computed as the mean of the χ 2 test applied on each X i : 1 n Σ n-1 i=0 χ 2 (X i , law, p i ).

Experimental observations

In this section, we analyze the statistics of real keystroke dynamics from the datasets presented in section IV.

Durations correlations

First, diagonals of correlation matrix are discarded. Correlations between two durations DT Di [START_REF] Giot | A review on the public benchmark databases for static keystroke dynamics[END_REF], and

DT Dj [0] are discarded if j = i + 1 in Data{1,2}K datasets, or if D 1 [1] = D 2 [0]
in Data{1,2}D datasets, as they are in fact, or might be, the same duration. Digraph are considered equal:

• in Data{1,2}K datasets, if their positions in the Keystroke are equals.

• in Data{1,2}D datasets, if their content are equals.

As shown in Table I, no strong stable correlation has been found between durations of DigraphTime from different Digraph, (Out: Data{1,2}K, Out U: Data{1,2}K splitted by User). DigraphTime will be thus assumed independent. Also, no strong stable correlation implying durations d 0 and d 5 of a same DigraphTime has been found (05 K: Data{1,2}K splitted by User, 05 D: Data{1,2}D merged and splitted by Digraph 05: Data{1,2}D merged by Digraph).

Stable correlations have been detected between durations d 1 , d 2 , d 3 , d 4 of a same DigraphTime (05: Data{1,2}D merged by Digraph). This may be due to the fact each of theses durations can be written as

d x = d 3 + k x * d 0 + l x * d 5 with l x ∈ {0, 1}, k x ∈ {0, 1}, and σ(d 3 ) ≈ 3 * σ(d 0 + d 5 ).
In the scope of this paper, DigraphTime is assumed to be computable from 3 independent durations (ideally d 0 , d 5 , d 3 ).

Durations laws

For the 6 DigraphTime durations DT D [i], i ∈ 0, 6 , the 5 best laws that minimize χ 2 (DT D [i], law), with parameters depending on the Digraph and User, are presented in Table 2. DigraphTime durations will then be assumed to best follow either a gumbel, normal, or logistic law, which parameters depends on the User and Digraph. In order to reduce the number of possible combinations, if a DigraphTime duration is generated with a given law, all other durations will be generated by the same law, but with different parameters. As seen in the previous sections, DigraphTime durations follow either a gumbel, a normal, or a logistic law which parameters can be estimated for each known User and Digraph. For a given User and Digraph, a DigraphGen can be then implemented as a set of 6 random engines generating the 6 Digraph-Time durations with the chosen law and estimated parameters.

We can see clearly in

We propose 10 consistency strategies, 1 for inconsistent DigraphTime, in which all durations are randomly generated (u), and 9 for partially-consistent DigraphTime, in which 3 durations are computed from the 3 others.

The durations to compute can be chosen among the 8 following lists, and be used for all Digraph and User, or be randomly chosen (n) for each new DigraphTime to generate:

• 0:

d 3 d 4 d 5 • 1: d 2 d 3 d 5 • 2: d 2 d 3 d 4 • 3: d 1 d 4 d 5 • 4: d 1 d 3 d 4 • 5: d 1 d 2 d 5 • 6: d 1 d 2 d 4 • 7: d 2 d 1 d 3
Once the DigraphGen created for a given User, the keystroke dynamics of a given Text T n is generated through the following process:

1. K[1] = T n 2. ∀i ∈ 0, n , K[0][i] = DT Tn[i] = DG Tn[i] ().
Before the consistency strategy application, and if Keystroke is expected to be consistent (or partially consistent), the DigraphTime first duration

K[0][i][0] is settled, if exists (i.e. if i > 0), to the last duration of the previous DigraphTime K[0][i -1][5].

Results analysis protocol

Synthetic datasets, SData{1,2}K L CS, are generated, for each law L, and consistency strategy CS, from each dataset Data{1,2}K. For each User of Data{1,2}K, the same number of entries (Keystroke dynamics) it has in Data{1,2} are generated (as seen in the previous section) and inserted into SData{1,2}K L CS. For each User, the first 10 entries are used as reference templates, the others as samples. User with less than 25 entries are discarded (i.e. with less than 15 samples).

For each dataset SData{1,2}K L CS, and distance function DistFct, 3 datasets are computed:

• DataSU: to qualify the capacity of synthetic Keystroke dynamics to be indistinguishable from real Keystroke dynamics;

• DataU: to qualify the KDS performance with real Keystroke dynamics data;

• DataS: to qualify, in comparison with DataU, the capacity of synthetic datasets to match the KDS performance that would be expected with real Keystroke dynamics data.

These datasets are composed of legitimate and impostor scores, computed with the distance function DistFct. Legitimate scores are obtained by comparing the reference template with samples from the same user. Impostors scores are obtained by comparing the reference template of users with samples from other users. DataU is computed from Data{1,2}K, and DataS, from SData{1,2}K L CS. In DataSU, legitimate scores are legitimate scores of DataU, and impostors scores are the distance, for each User, between real user templates, and its synthetic samples.

We consider the False Acceptance Rate (FAR) describing the ratio of accepted impostor data, the False Rejection Rate (FRR) describing the ratio of falsely rejected legitimate users. The Equal Error Rate (EER) corresponds to configuration of the biometric system when FAR equals FRR. The FAR/FRR points, and the EER are computed from EvaBio [START_REF] Mahier | Computation evabio: A tool for performance evaluation in biometrics[END_REF] considering 25,000 thresholds. Each final indicator is computed by averaging the computed indicators of each Data{1,2}K datasets.

Usurpation of keystroke dynamics

The EER value computed from DataSU is used to qualify the capacity of synthetic Keystroke dynamics data to be indistinguishable from real Keystroke dynamics data. An EER of x% means that it is not possible to choose a threshold, such as rejecting less than x% genuine users, without accepting less than x% EER impostors. Thus, an EER of 50% means that, for this threshold, the choice to accept or reject a user is not better than random. An EER > 50% means that, for this threshold, more impostors will be accepted than genuine users.

As shown in table 3, synthetic Keystroke dynamics data are indistinguishable from real one (EER ≈ 50%), when using the chosen distance functions, and are even more, for some configurations, more accepted than real Keystrokes. The fact that synthetic data are more accepted than real data, can be explained either by a lesser intra-score variance, or by a lesser intra-score mean, for synthetic data compared to real data.

Consistency strategy seems to matter more than the used law. Strategies that computes d 5 (0, 1, 3, 5) instead of randomly generating it have a worst EER value than the strategy which randomly choose the durations to be computed, which is worst than other strategies. Results can also be divided into 4 groups in which the laws giving the greater EER values are, in order of superiority, gumbel, normal, and logistic laws.

As shown by the symmetric of the FAR/FRR curves in Figure 3, our proposed Keystroke generation method is thus able to produce synthetic samples that enable identity usurpation of a known user, by imitating its keystroke dynamics.

Scores estimations

A χ 2 test is performed, in a similar way than previously, between the legitimate and the impostor scores of DataSU, to qualify the capacity of synthetic samples of having the same scores than real samples. However, as it can be seen in Figure 4, synthetic sam- We thus demonstrate that even being close to real samples scores, our Keystroke generation method is not able yet to produce synthetic samples which scores would match real samples.

EER and FAR/FRR estimations

The absolute and relative distances between the EER values, computed from the synthetic (DataS) and real (DataU) entries, qualify the capacity of the synthetic datasets to estimate the EER value of the real one.

As shown in Table 4, the EER value can be estimated with a variable accuracy (errors from 0.34 to 28.43). The logic found in previous sections is not respected, and even inverted. Indeed, consistency strategies 0, 1, 3, 5 and n, best estimate the EER value than strategies 7, 6, 4, and 2, the worst strategy being u.

As shown in Figures 5 and6, our proposed Keystroke generation is thus able to estimate EER values of real dataset. However, as shown in Table 4, the best consistency strategies to estimate EER value are the worst to generate synthetic samples scores close to real samples scores.

Conclusion

In this paper, we presented a method that enables the generation of synthetic keystroke dynamics data from known Users, to either usurp real user KD, or to estimate EER value of a KDS. These methods have been tested on fixed text, but could be as well applied to free text.

This work constitutes a first step towards the generation of large synthetic Keystroke dynamics datasets. The following step would be the generation of keystroke dynamics data for an unknown user. Such large synthetic Keystroke dynamics datasets could then be used to fairly compare KDS performances, as well as to improve learning-based KDS' performances.

We also aim to improve the current method with a better estimation of laws parameters, either through an non-User specific PCA on DigraphTime durations, a better estimation of statistics, or the usage of several laws (e.g. cosine for d 0 , d 5 , and gumbel for d 3 ).

•

  Keystroke dynamics: K = [{DT i } i∈ 0,n , T n ], an array of n DigraphTime DT i associated to the Digraph T n [i]. Keystroke is said consistent (or partially consistent) if T n , and all DT i are consistent (or partially consistent), and if ∀i ∈ 0, n , DT i-1 [5] = DT i [0].

Figure 1 :

 1 Figure 1: KD Generative model

Figure 2 :

 2 Figure 2: DigraphTime

Figure 3

 3 Figure 3: FAR/FRR curves of synthetic samples against real samples compared to real references with Hocquet distance (DataSU from SData2K1 gumbel 6) Figure 4: Hocquet distances of real (grey) and synthetic (black) Keystrokes samples with real templates (DataSU from SData2K2 gumbel 6)

Figure 5

 5 Figure 5: FAR/FRR curves for a real dataset with BioHashing distances (DataU from Data2K1) Figure 6: FAR/FRR curves for a synthetic dataset with BioHashing distances (DataS from SData2K1 logistic 0)

Table 2

 2

	that the estimated

Table 1 :

 1 Correlations found in Data1D/Data2D, and in Data1K/Data2K1/Data2K2 datasets.

	Datasets Rank	d 0	χ 2	d 1	χ 2	d 2	χ 2	d 3	χ 2	d 4	χ 2	d 5	χ 2
		1 cosine (3σ) 0.964 gumbel (3σ) 0.808 gumbel (3σ) 0.795 gumbel (3σ) 0.820 gumbel (3σ) 0.820 cosine (3σ)	0.965
		2 normal (3σ) 0.966 normal (3σ) 0.822 normal (3σ) 0.809 normal (3σ) 0.834 normal (3σ) 0.833 cosine		0.968
	Data1D	3 cosine	0.967 logistic (3σ) 0.831 logistic (3σ) 0.812 cosine (3σ)	0.842 logistic (3σ) 0.841 normal (3σ) 0.969
		4 logistic	0.967 gamma (3σ)	0.836 betaprime (3σ) 0.820 logistic (3σ) 0.842 cosine (3σ)	0.844 normal	0.969
		5 normal	0.967 cosine (3σ)	0.839 cosine (3σ)	0.820 gamma (3σ)	0.863 gamma (3σ)	0.863 logistic	0.970
		1 normal (3σ) 0.869 gumbel (3σ) 0.810 gumbel (3σ) 0.781 gumbel (3σ) 0.835 gumbel (3σ) 0.82 normal (3σ) 0.874
		2 cosine (3σ) 0.872 normal (3σ) 0.820 normal (3σ) 0.794 normal (3σ) 0.838 normal (3σ) 0.832 logistic (3σ) 0.881
	Data2D	3 logistic (3σ) 0.880 logistic (3σ) 0.880 betaprime (3σ) 0.800 logistic (3σ) 0.848 betaprime (3σ) 0.835 cosine (3σ)	0.884
		4 gamma (3σ) 0.893 gamma (3σ)	0.834 gamma (3σ)	0.804 cosine (3σ)	0.858 gamma (3σ)	0.836 gamma (3σ)	0.900
		5 gumbel (3σ) 0.893 betaprime (3σ) 0.834 logistic (3σ) 0.807 betaprime (3σ) 0.890 logistic (3σ) 0.839 betaprime (3σ) 0.879

Table 2 :

 2 Top 5 results of χ 2 test with 19 laws, with (3σ) and without exclusion of abherant values

	7 Keystroke dynamics genera-
	tive model
	7.1 Principles

Table 3 :

 3 Mean of DataSU EER (Top 5)


Loubna Ghammam
email: loubna.ghammam@unicaen.fr

Morgan Barbier

Christophe Rosenberger
email: christophe.rosenberger@ensicaen.fr

Enhancing the Security of Transformation Based Biometric Template Protection Schemes

Keywords: Template protection, biometric authentication, attack, performance evaluation

Template protection is a crucial issue in biometrics. Many algorithms have been proposed in the literature among secure computing approaches, crypto-biometric algorithm and feature transformation schemes. The BioHashing algorithm belongs to this last category and has very interesting properties. Among them, we can cite its genericity since it could be applied on any biometric modality, the possible cancelability of the generated BioCode and its efficiency when the secret is not stolen by an impostor. Its main drawback is its weakness face to a combined attack (zero effort with the stolen secret scenario). In this paper, we propose a transformation-based biometric template protection scheme as an improvement of the BioHashing algorithm where the projection matrix is generated by combining the secret and the biometric data. Experimental results on two biometric modalities, namely digital fingerprint and finger knuckle print images, show the benefits of the proposed method face to attacks while keeping a good efficiency.

I. INTRODUCTION

Biometrics is an emerging technology for authentication applications. Many biometric modalities are well known and used (such as fingerprints), the design of intelligent sensors is advanced (liveness detection) and algorithms provide very good results. Privacy issues concerning this particular personal information still limit its operational use. The General Data Protection Regulation (GDPR) 2016/679 is a regulation in EU law on data protection and privacy for all individuals within the European Union. Biometric data is of course perceived as personal and sensitive data. In Europe, as for example, the central storage of biometric data is forbidden or limited to a small amount of users. In order to solve this problem, new biometric systems have been proposed in the last decade based on the "privacy by design" paradigm. These biometric template protection schemes have as objective to guarantee the security and privacy of users to face attacks such as identity theft (e-government applications, border control, etc.) [START_REF] Jain | Biometric template security[END_REF].

Three main approaches can be distinguished dealing with template protection in biometrics. First, biometric cryptosystems or secure sketches, such as those presented in [START_REF] Chabanne | Optimal iris fuzzy sketches[END_REF], [START_REF] Juels | A fuzzy vault scheme[END_REF], [START_REF] Juels | A fuzzy commitment scheme[END_REF], resort to cryptography. Second, secure computing methods aim at computing the comparison of two biometric templates by an untrusted party [START_REF] Bringer | Gshade: faster privacy-preserving distance computation and biometric identification[END_REF], [START_REF] Chatterjee | Secure biometricbased authentication scheme using chebyshev chaotic map for multiserver environment[END_REF]. Last, we find feature transformations approaches for template protection. The BioHashing algorithm is one of the most popular technique and is based on biometric data salting. It has been developed for different biometric modalities such as those presented in [START_REF] Belguechi | Biohashing for securing minutiae template[END_REF], [START_REF] Saini | Soft biometrics in conjunction with optics based biohashing[END_REF], [START_REF] Teoh | Biohashing: two factor authentication featuring fingerprint data and tokenised random number[END_REF].

These last systems are called cancelable since the result generated from a biometric template, namely BioCode, can be revoked in case of interception or loss. This BioCode cannot be used as a cryptographic key as the generated BioCode is not exactly the same for each biometric capture. These particular biometric systems must of course address classical issues such as a high level of performance (i.e., minimizing the Equal Error Rate (EER) or Area Under the Curve (AUC) value of the system) but also new constraints concerning privacy. In the literature, many papers have been published dealing with the definition of new schemes for the protection of biometric templates (such as those presented in [START_REF] Belguechi | Biohashing for securing minutiae template[END_REF], [START_REF] Jin | Ranking-based locality sensitive hashing-enabled cancelable biometrics: Index-of-max hashing[END_REF], [START_REF] Ratha | Generating cancelable fingerprint templates[END_REF], [START_REF] Teoh | Cancellable biometrics realization through biophasoring[END_REF]).

Most of such protection schemes lack of robustness considering the stolen token scenario. In this case, an attacker knowing the secret of the transformation protection scheme, has a big advantage to impersonate an user as presented in [START_REF] Belguechi | Comparative study on texture features for fingerprint recognition: application to the biohashing template protection scheme[END_REF], [START_REF] Nagar | Biometric template transformation: A security analysis[END_REF]. This is due to the fact that the used projection matrix is computed only given the secret. The main contributions of the paper are twofolds. First, we propose a new transformation for the protection of biometric data. Its main benefit is that the used projection matrix is not only related to a secret but also embed information computed from the biometric data which limits some attacks. Second, we analyze the behavior of this new transformation in comparison with two other methods by using a recent analysis methodology of such schemes [START_REF] Rosenberger | Evaluation of biometric template protection schemes based on a transformation[END_REF].

The paper is organized as follows. Section II gives the background on template protection schemes based on a transformation. Section III is dedicated to a literature review on template protection schemes based on a transformation. Then, we present the proposed methodology in Section IV. Section V illustrates the benefits of the proposed methodology through experimental results to be compared with those of two cancelable biometric systems. Finally, in section VI, we conclude and we give some perspectives.

II. BACKGROUND

The general principle of template protection schemes based on a transformation is illustrated on Figure 1. These schemes consist in generating a binary output called BioCode given a biometric template and a secret.

We propose to keep the notations of [START_REF] Nagar | Biometric template transformation: A security analysis[END_REF]. Let T z and Tz represent respectively the template and query biometric features of user z. Let f be the feature transformation function. Let K z be a set of transformation parameters corresponding to the user z. A feature transformation is a non-invertible function using a user related key K z (i.e typically a random seed or a strong password), applied to the used biometric template T z . The BioCode f (T z , K z ) is stored in a database or in a personal device. It is generally considered that, given the transformed template f (T z , K z ) and the key K z , it is possible to recover the original template T z (or a close approximation) as presented in [START_REF] Nagar | Biometric template transformation: A security analysis[END_REF]. Thus, it is requested to store this key in a second support, even if the reconstruction of the original template strongly depends on the used biometric modality.

Let n denotes the dimension of the f (T z , K z ) BioCode for the user z. Let D T denotes a distance function between the biometric features in the untransformed (original) domain. However, the comparison is computed in the transformed domain, thus we need D T a distance function in the transformed one. The cancelable biometric system outputs a verification decision denoted R z if the distance between the reference BioCode and query BioCode is less than a decision threshold denoted as :

R z = 1 {D T (f (Tz,Kz),f ( Tz,Kz))≤ } (1)
The performance of the authentication system is generally estimated with FRR (False Rejection Rate)/FAR (False Ac-ceptance Rate) rates and the feature transformation should not decline the performance of the system. In fact, this approach tends to improve the performance of the biometric system without any protection even if the key K z is necessary for the user z to authenticate herself/himself.

III. RELATED WORKS

The concept of privacy protection of biometric data has been defined in 2001 in a seminal paper [START_REF] Nalini | Enhancing security and privacy in biometrics-based authentication systems[END_REF]. Since then, many methods have been proposed among random projections approaches [START_REF] Jaishanker K Pillai | Sectored random projections for cancelable iris biometrics[END_REF], BioHashing methods [START_REF] Teoh | Biohashing: two factor authentication featuring fingerprint data and tokenised random number[END_REF], Bloom filters [START_REF] Rathgeb | On application of bloom filters to iris biometrics[END_REF]. . . A complete review of cancelable biometric systems can be found in [START_REF] Vishal | Cancelable biometrics: A review[END_REF]. Very recently, Teoh et al. [START_REF] Jin | Ranking-based locality sensitive hashing-enabled cancelable biometrics: Index-of-max hashing[END_REF] proposed a new two-factor scheme to protect the biometric template by transformation. Compared with previous works, this method is based on localized random projection and on the rank correlation. Moreover, the obtained results show that this system is strongly resistant against the main attacks. These good results are the consequence of their technical called Index-Of-Max which can be viewed as a machine learning on the plain database. For this previous constraint, we do not compare to this method where the BioSystem is tuned for a particular basis. More generally, we can find a security analysis of the biometric system protecting the biometric template based on transformations [START_REF] Rosenberger | Evaluation of biometric template protection schemes based on a transformation[END_REF].

In this paper, We detail particularly two popular template protection schemes: BioHashing [START_REF] Teoh | Biohashing: two factor authentication featuring fingerprint data and tokenised random number[END_REF] and BioPhasor [START_REF] Teoh | Cancellable biometrics realization through biophasoring[END_REF] and in this context, we detail each algorithm.

The BioHashing algorithm is applied on biometric templates that are represented by real-valued vectors of fixed length (so the metric used to evaluate the similarity between two biometric features is the Euclidean distance). It generates binary templates of length lower or equal to the original length (here, the metric used to evaluate the similarity between two transformed templates is the Hamming distance). This algorithm has been originally proposed for face and fingerprints by Teoh et al. in [START_REF] Teoh | Biohashing: two factor authentication featuring fingerprint data and tokenised random number[END_REF], where the fingerprint features are, in a first time, transformed in a real-valued vector of fixed length to generate the biometric template (this step is not useful and not described in this paper). The BioHashing algorithm transforms the biometric template T = (T 1 , . . . T n ) in a binary template B = (B 1 , . . . B m ), with m ≤ n, as following:

Algorithm 1 BioHashing 1: Inputs 2: T = (T 1 , . . . , T n ): biometric template, 3: Seed: secret 4: Output B = (B 1 , . . . , B m ): BioCode 5: Generation with the Seed of m pseudorandom vectors V 1 , . . . , V m of length n, 6:
Orthogonalize vectors with the Gram-Schmidt algorithm, 7: for i = 1, . . . , m do compute x i =< T, V i >. 8: Compute BioCode:

B i = 0 if x i < τ 1 if x i ≥ τ,
where τ is a given threshold, generally equal to 0.

The specificity of the BioHashing algorithm is that it uses a one way function and a random seed of m bits. It is important to note that every enrolled biometric feature uses a different seed in order to create a specific BioCode. The performance of this algorithm is ensured by the scalar products with the orthonormal vectors. The quantization process of the last step ensures the non-invertibility of the data (even if n = m, because each coordinate of the input T is a real value, whereas the coordinates of the output B is a single bit). Finally, the random seed guarantees the diversity and revocability properties.

BioPhasor was proposed by Teoh et al. in [START_REF] Teoh | Cancellable biometrics realization through biophasoring[END_REF] and it was introduced as a form of cancellable biometrics which is based on iterated mixing between the user-specific pseudo-random number and the biometric feature. The BioPhasor algorithm is supposed to be an improvement of the BioHashing one. It is described in Algorithm 2. The step 8 is added in this paper in order to generate a binary output.

An interesting component of these schemes is that no learning phase is required. Nevertheless, some weaknesses have been reported in the former approach in [START_REF] Kong | An analysis of biohashing and its variants[END_REF], [START_REF] Lumini | Empirical tests on biohashing[END_REF], [START_REF] Simoens | Privacy weaknesses in biometric sketches[END_REF]. A main reason is that the projection matrix is only related to the secret. If an impostor obtains the Key (known as stolen token attack), the attack is quite easy especially by combining it with the zero effort one (use of a biometric data belonging to another user). Indeed, with the knowledge of the Key, the impostor will use the same projection basis that the legitimate

B i = 0 if h i < τ 1 if h i ≥ τ,
where τ is a given threshold, generally equal to 0.

user; thus it increases the success probability of its attack. In the sequel, we propose a new transformation whose objective is to limit this aspect.

IV. PROPOSED METHOD

We intend in this paper to enhance the security of the BioHashing algorithm by limiting the impact of attacks based on the stolen token. We modified the projection matrix in order to be computed from the biometric template and not only from the secret. We suppose in this method to have some statistics concerning the features (i.e. average and standard deviation). This hypothesis is not a high restriction of the operational use of the method and keeps its possible use for any biometric modality. We generate a projection matrix given the seed as secret and the biometric template. In the proposed method, we use:

• A constant number C that permits us to define a trade-off between performance without any attack and robustness to attacks.

• A polynomial P k such as P k (X) = 3 i=0 a i * X i , k = 1 : m, a i ∈ [-C.σ C.σ], a 3 = 0.
The choice of a polynomial of degree 3 is well justified. In general, a k-wise independent function family is presented by a random polynomials of degree k -1. In this context, Hoory et al. in [START_REF] Hoory | Simple permutations mix well[END_REF] have presented an educated conjecture that 4-wise independence suffices to achieve pseudorandomness. That is why, in our case, we use a 4-function family, so a polynomial of degree 3. In the following figure IV, we detail the general principle of the proposition of the BioHashing improvement and the construction of the projection matrix computed from the seed and from biometric data. The proposed algorithm of our method is given in detail in Algorithm 3.

The key ingredient of our proposed method is the new construction of the projection matrix P . Firstly, we choose randomly polynomials of degree 3 initialized by the Seed. On the other hands, we construct the projection matrix P

B i = 0 if X i < τ 1 if X i ≥ τ,
where τ is a given threshold, generally equal to 0.

by the evaluation of these polynomials at each point given by the biometric template; which produces the first method for biometric template protection by a transformation which depends both on a secret and of a biometric template itself.

V. VALIDATION

In this section, we present some experimental results demonstrating the benefit of the proposed transformation.

A. Dataset

In order to study the performance and robustness of the proposed transformation, we selected two biometric datasets.

The first one is well known as the FVC2002 DB1 database composed of digital fingerprint images. 100 individuals provided 8 samples of fingerprints. Figure 3 presents some examples of fingerprints. We compute well known Gabor features for each image, we obtain in this case 512 real valued features for each fingerprint. The second dataset is the PolyHK which is composed of images of knuckle prints [START_REF] Zhang | Finger-knuckle-print verification based on band-limited phase-only correlation[END_REF] (see figure 4). The database has been acquired on 4 fingers of 165 volunteers, leading to 660 different classes. Each class contains 12 images acquired during 2 sessions. We compute Gabor features for each image, we obtain 256 real valued features for each finger knuckle print.

B. Evaluation

Few works have been dedicated to the evaluation of such biometric systems in the literature [START_REF] Adler | Biometric system security. Handbook of biometrics[END_REF], [START_REF] Nagar | Biometric template transformation: A security analysis[END_REF], [START_REF] Zhou | Feature correlation attack on biometric privacy protection schemes[END_REF]. The ISO/IEC 24745 "Information technology Security techniques Biometric information protection" defines the security properties of a biometric system, we use in this paper the same terms. Cancelable systems must fulfill several properties as also mentioned in [START_REF] Maltoni | Handbook of Fingerprint Recognition[END_REF], [START_REF] Rosenberger | Evaluation of biometric template protection schemes based on a transformation[END_REF]:

• Revocability/Renewability: It should be possible to revoke a biometric template and to generate a new one from the original biometric data. • Performance:

The template protection shall not deteriorate the performance of the original biometric system. • Non-invertibility or Irreversibility:

From the transformed data, it should not be possible to obtain enough information on the original biometric data, to prevent any attack consisting in forging a stolen biometric template (as for example, it is possible to generate an eligible fingerprint given minutiae [START_REF] Maio | Fingerprint image reconstruction from standard templates[END_REF]). This property is essential for security purposes. For any attack, an impostor provides an information in order to be authenticated as the legitimate user. The success of the attack is given by:

F AR A ( ) = P (D T (f (T z , K z ), A z ) ≤ ) (2)
Where F AR A is the probability of a successful attack by the impostor for a decision threshold set to . The A z BioCode is computed by the impostor by taking into account as much information as possible within different contexts.

• Diversity or Unlinkability:

It should be possible to generate different BioCodes for multiple applications, and no information should be deduced from the comparison or the correlation of different realizations.

• Indistinguishability:

It should be infeasible to cross correlate two protected templates. Based on some of the early works [START_REF] Ratha | Enhancing security and privacy in biometrics-based authentication system[END_REF], [START_REF] Bolle | Biometric perils and patches[END_REF] which identified weak links in each subsystem of a generic authentication system, some papers considered the possible attacks in cancelable biometric systems (such as those presented in [START_REF] Jain | Biometric template security[END_REF], [START_REF] Nagar | Biometric template transformation: A security analysis[END_REF], [START_REF] Saini | Soft biometrics in conjunction with optics based biohashing[END_REF], [START_REF] Teoh | Cancellable biometrics and annotations on biohash[END_REF]). We follow the Shannon's maxim ("The enemy knows the system"), we so assume that the impostor has all necessary information on the process used to generate the BioCode (feature generation method, BioCode size...).

Based on the principle of each attack, we generate many fake attempts A z of the genuine user in an authentication case:

• Zero effort attack: an impostor user x provides its biometric feature Tx and parameter K x to be authenticated as the user z :

A z = f ( Tx , K x )
• Brute force attack:

An impostor tries to be authenticated by trying different random values of A: A z = A • Stolen token attack:

An impostor has obtained the token K z of the genuine user z and tries different random values T to generate:

A z = f (T, K z ) • Stolen biometric data attack:
An impostor knows Tz (directly or after computation of the feature on a biometric raw data) and tries different random numbers K to generate: A z = f ( Tz , K) • Worst case attack:

An impostor user x provides its biometric template Tx and parameter K x to be authenticated as the user z (zero effort attack) and has also obtained the token K z of the genuine user z to generate:

A z = f (T x , K z ) • Listening attacks:
An impostor must not be able to extract any information from different BioCodes issued from the same user. Since BioCodes can be revoked, an impostor can intercept Q of them and issue a new one by predicting an admissible value (as for example by setting each bit to the most probable value). TThese attacks consist in the following process:

-Generation of Q BioCodes for user z: These attacks allow us to quantify the robustness of cancelable biometric verification systems based on feature transformation.

B z = {f (T z , K z 1 ), .., f (T z , K z Q )} -Prediction of

C. Experimental results

In this section, we present the experimental results we obtained with the proposed transformation face to the BioHashing and BioPhasor protection schemes. Before considering their behaviors, we computed the performance of biometric systems by using Gabor features with the Euclidean distance without any transformation. For the fingerprint dataset, we obtained an Equal Error Rate (ERR) of 28.3% and 23.8% for the Finger knuckle print one. We can see clearly that these results are poor and such a biometric system could not been used in real conditions. Transformation protection schemes are known to increase performance if the secret is unknown by impostors. For the two datasets and then for the two biometric modalities, the BioHashing and BioPhasor algorithms provide a perfect performance without any attack i.e. when the secret is unknown by the impostor. For the proposed transformation, it is not always the case, in fact, it depends on the value of C). When C is set to 10, the EER value equals to 2%. The main reason of this decrease of performance is due to the fact that the projection matrix is computed by considering biometric data and not only the secret. As the performance when using raw features (i.e. without any transformation) is low, it is normal to have such an impact. We will see later that we gain of course in term of robustness. We will also show later a study of the impact of the C value on performance and robustness of the proposed transformation scheme.

We then analyze the robustness of the protection schemes related to tested attacks. Of course, the brute force attack (curve in blue) is the less successful attack and all schemes are able to avoid it. The BioHashing and BioPhasor algorithms have a similar behavior for all attacks even if the BioPhasor algorithm is more sensitive to the stolen token attack (curve in green). Only the worst case attack is completely possible (F AR A = 100%) for these two protection schemes when the decision threshold value is set to have the behavior of the biometric systems at the EER value. The proposed transformation scheme is slightly more robust to attacks (even in the worst case one) thanks to the use of information related to the biometric data in its computation of the projection matrix. In order to better identify the robustness behavior of these three protection schemes, we propose to plot the FRR value (as a performance indicator) face to the F AR A one of the worst case (as a robustness indicator in the worst case) in Figure 8. There is a compromise to find between performance and robustness in such transformation schemes. This curve demonstrates that for a given performance (here the FRR value), the robustness evaluated by the F AR A value is lower for the proposed method.

Finally, we present in Table I the results of the study concerning the impact of the C value in the proposed scheme for the fingerprint dataset. We can see clearly that a higher value of C permits to obtain a better performance but also a lower robustness in the wort case attack. As for example, for a value of C equals to 50, the performance is very good with an EER value of 0.4% and a F AR A value of 99%. With the BioHashing and BioPhasor algorithms, the F AR A value was 100%. Setting the C value can adjust the compromise between performance and robustness. We have to recall at this point that the use of raw features lead to an EER value of more that 22%, we could expect even better results of the proposed transformation with more efficient features.

VI. CONCLUSION AND PERSPECTIVES

Protecting biometric templates is actually crucial due to new law regulations on data protection (such as the GDPR in Europe) and the possible use of biometric data for authentication in cloud services. We believe feature-based template protection schemes could have a high impact on security services in the near future. In this context, features, that can be computed without any learning with unprotected templates from other users, should be protected by such cancelable schemes. The proposed transformation scheme has the advantage to better combine the secret with the biometric data in order to reduce the impact of the stolen token attack. We obtain with this scheme a configurable transformation by adjusting the compromise between the expected performance by assuming the secret is unknown and a better robustness otherwise with a single value (C).

Perspectives of this work are multiple. First, it concerns the proposal of other combination techniques of the biometric data and the secret (step 10 of the proposed algorithm). Second, we could consider its extension to multi-biometrics where a biometric modality could be used to generate the projection matrix given the secret and the features of the second modality could be projected with it. Last, this scheme can be used in real biometric authentication applications in industry as the computation is very fast while keeping a good privacy protection.

Fig. 1 .

 1 Fig. 1. General principle of template protection schemes based on a transformation.

Algorithm 2

 2 BioPhasor 1: Inputs 2: b = (b 1 , . . . , b n ): biometric template 3: Seed: secret 4: Output B = (B 1 , . . . , B m ): BioCode 5: Generation with K of m pseudorandom vectors V 1 , . . . , V m of length n, 6: Orthogonalize vectors with the Gram-Schmidt algorithm, 7: for i = 1, . . . , m do compute h i = 1/n n j arctan(b j 2 /V j i ). 8: Compute BioCode:

Fig. 2 .

 2 Fig. 2. General principle of the proposed method.

Fig. 3 .

 3 Fig. 3. Some examples of fingerprint images in the FVC2002DB1 dataset.

Fig. 4 .

 4 Fig. 4. Some examples of finger knuckle print images in the PolyHK dataset.

  a possible BioCode value by setting the most probable value of each bit given B z , ⇒ computing the F AR A value for Q = 3 and Q = 11

Figures 5 ,

 5 Figures 5,[START_REF] Chabanne | Optimal iris fuzzy sketches[END_REF] and 7 present at the same time the performance of cancelable biometric systems without any attack and the robustness to attacks described in the previous section. The yellow curve presents the evolution of the false rejection rate (FRR) depending on the decision threshold ( in equation 2). Other curves correspond to the F AR A value for all considered attacks.

Fig. 5 .

 5 Fig. 5. Security analysis of the BioHashing algorithm on the fingerprint database (a) and the finger knuckle print (b).

Fig. 6 .

 6 Fig. 6. Security analysis of the BioPhasor algorithm on the fingerprint database (a) and the finger knuckle print (b).

Fig. 7 .

 7 Fig. 7. Security analysis of the proposed algorithm on the fingerprint database (a) and the finger knuckle print (b). Experiments have been done with polynomial function of degree 3 and C = 10.

Fig. 8 .

 8 Fig. 8. Comparison between performance and robustness face to attacks for the three tested template protection schemes. The FRR value is seen as a performance indicator while the F AR A value gives the probability of successful attack in the worst case.

TABLE I STUDY

 I OF THE IMPACT OF THE C VALUE ON PERFORMANCE AND ROBUSTNESS IN THE PROPOSED TRANSFORMATION SCHEME.C value EER (without attack) zero effort brute force stolen token stolen biometric worst case listening N=3 listening N=11

	5	7%	6.9%	0%	4.0%	6.7%	94.1%	5.4%	4.1%
	20	0.6%	0.8%	0%	0%	5.7%	99.3%	0.3%	4.1%
	50	0.4%	0.3%	0%	0%	0.9%	99%	0.1%	0.1%

Acknowledgements. The authors thank Kevin Atighehchi for comments on this paper.