Michel Agoyan
email: michel.agoyan@cea.fr

Jean-Max Dutertre
email: dutertre@emse.fr

David Naccache
email: david.naccache@ens.fr

Bruno Robisson
email: bruno.robisson@cea.fr

Assia Tria
email: assia.tria@cea.fr

When Clocks Fail -On Critical Paths And Clock Faults

Whilst clock fault-attacks are known to be a serious security threat, an in-depth explanation of such faults still seems in order.

This work provides a theoretical analysis, backed by practical experiments, explaining when and how clock faults occur. Understanding and modeling the chain of events following a transient clock alteration allows to accurately predict faulty circuit behavior. A prediction fully confirmed by injecting variable-duration faults at predetermined clock cycles.

We illustrate the process by successfully attacking an fpga aes implementation using a dll-based fpga platform (one-bit fault attack).

Introduction

Fault attacks consist in modifying an electronic circuit's behavior to achieve malicious goals [START_REF] Boneth | On the importance of checking cryptographic protocols for faults[END_REF][START_REF] Biham | Differential fault analysis of secret key cryptosystems[END_REF]. Fault attacks exist in numerous variants ranging from a simple alteration of a round counter during symmetric encryption [START_REF] Choukri | Round reduction using faults[END_REF] to mathematical Differential Fault Attacks (dfa) where secret information is obtained by comparing (differentiating) correct and faulty encryption results [START_REF] Dusart | Differential fault analysis on aes[END_REF][START_REF] Piret | A differential fault attack technique against spn structures, with application to the aes and khazad[END_REF].

Faults can be caused by a variety of intrusive and non-intrusive means [START_REF] Choukri | The sorcerer's apprentice guide to fault attacks[END_REF] such as lasers [START_REF] Skorobogatov | Optical fault induction attacks[END_REF], electromagnetic perturbations [START_REF] Quisquater | Eddy current for magnetic analysis with active sensor[END_REF][START_REF] Hutter | Optical and em fault-attacks on crt-based rsa: Concrete results[END_REF], voltage variations [START_REF] Guilley | Silicon-level solutions to counteract passive and active attacks[END_REF] or clock glitches [START_REF] Fukunaga | Practical fault attack on a cryptographic lsi with iso/iec 18033-3 block ciphers[END_REF].

In this work we present a new clock alteration technique for scenarii in which the attacker is given access to the target's clock. We start by explaining and modeling the chain of events causing the faulty behavior. This theoretical analysis perfectly reflects experimental observations and allowed the injecting of precise single-bit faults into a chip running the aes algorithm (an actual implementation of the attack described in [START_REF] Giraud | DFA on aes[END_REF]).

After introducing the model, we will overview the fault injector's design, the target chip's structure and the way in which the injector was used to extract aes keys.

Why Clock Faults Occur?

We inject faults by violating synchrony, a basic assumption under which traditional digital ICs operate. In essence, most 4 ICs execute calculations by processing data by combinatorial logic blocks separated by D flip-flop register banks sharing the same clock (figure 1). Data is usually latched by registers at raising clock edges. Between two such edges, the computed data travels between registers and gets modified by the intermediate combinatorial logic blocks. The time needed to propagate data through combinatorial logic is called propagation delay. The propagation delay and a second delay element, inherent to the use of D flip-flop, called set-up time, define the circuit's maximal operating frequency (nominal circuit period). Indeed, to ensure proper circuit operation, the clock period must be strictly greater than the maximal propagation delay in concerned circuit (this maximal propagation delay is called critical path) plus the registers' set-up time. In other words:

T clock > t critical + t set-up (1)
As a matter of fact any data bit entering a register is the result of a combinatorial calculation involving several previous register output bits. The transformation of the previous registers' output into the next register's input bit takes a determined delay. This delay depends on the the logic performed as well as on the data transiting through the logic. In addition, propagation time varies with circuit temperature and power supply voltage.

Overclocking

Overclocking consists in decreasing the clock period (or, put differently, increasing clock frequency). If setup delays are not respected, the D flip-flop's input is not given sufficient time to reach the latch. Causing faulty data to be latched instead. This led several authors to use overclocking as fault injection means [START_REF] Guilley | Silicon-level solutions to counteract passive and active attacks[END_REF][START_REF] Selmane | Practical setup time violation attacks on aes[END_REF].

A decreased clock period can potentially affect logical paths whose propagation times exceed the decreased clock period minus the set-up time. From the attacker's perspective the ability to control precisely the clock period is crucial for inducing faults with precision. Note that temperature and power supply changes may also be used to exert such control.

Fault attacks consist in injecting faults at precise moments. To avoid injecting faults continuously, overclocking must be brief and transient. The fault injection technique described in the next section allows doing so.

Injecting Clock Delay Faults

An attacker needs to control two parameters: the precise moment at which the fault occurs and the clock anomaly's duration. The latter must be controlled with high resolution, typically a few of tens of picoseconds. The two waveforms differ only between the delimiters positioned at 20ns and 30ns. During that interval, the faulty clk's period is reduced by ∆ ns. The ∆ time decrement causes a set-up time violation fault. Note that the extension by ∆ of the preceding clock cycle's low state has no adverse effects.

Generating faulty clk with sufficient accuracy and precision is a challenging task. To do so, we use the embedded Delay Locked Loop (dll) of a recent fpga family (Xilinx Virtex 5). Two clocks (clk delayed i) with programmable skews are generated from clk. The skews of the clk delayed i signals, denoted δ i , are programmable. faulty clk is obtained by switching between the clk delayed i using a trigger signal. Figure 3 The device assembles faulty clk by combining clk delayed 2's raising edge and clk delayed 1's falling edge. This is controlled by the signal trigger that positions the perturbation in time. The resolution at which faulty clk's shape can be controlled (in our setting 35ps) depends on δ t , the smallest elementary delay that the dll is able to provide. We will refer to faulty clk's altered signal chunk as the faulting period (in Figure 3, the interval between 24ns and 30ns).

An oscilloscope screen-shot of faulty clk is shown on Figure 4 (uppermost signal). Here, clk's period (10ns) was reduced to [t 1 , t 2 ] = 49δ t ≃ 8.2ns in faulty clk. The lowermost signal's high level indicates the aes' start. The implementation completes an encryption round in one cycle. Hence, the diagram shows a fault injection at the ninth round (cf. section 3.2).

As we write this lines, a Xilinx Virtex 5 development board costs less than 1000.

Clock Fault dfa on aes

We tested the attack setup on a concrete aes implementation. The following sections describe the target chip, the attack's theoretical principle ( [START_REF] Giraud | DFA on aes[END_REF]) and report practical experiment results.

The Test Chip

The test chip (Xilinx Spartan 3AN fpga) implements a hardware 128-bit aes [START_REF]Announcing the Advanced Encryption Standard (AES)[END_REF][START_REF] Daemen | Aes proposal: Rijndael[END_REF] written in vhdl. The design consists of three main blocks: a communication and control module (ccm), a key expansion module (kem), and an encryption module (enm).

The ccm manages the serial link through which plaintext and key material are input. The start signal triggering the encryption and the resulting ciphertext also transit through the ccm. In addition, the ccm controls the kem and the enm's operations during encryption.

The implementation uses a 128 bit data path and runs the kem and the enm in parallel. Consequently, an encryption round is completed in one clock cycle and the entire aes computation takes 11 clock cycles.

The kem generates the round keys 'on the fly'. At each clock cycle, a new round key is transferred from the kem to the enm. We will not overview the kem in detail as it is of little relevance for our attack. We nonetheless underline that the kem's critical delay path is much smaller than the enm's one -this is essential for the attack to work.

The enm architecture is depicted on Figure 5. The enm breaks-down into five submodules: AddRoundKey, SubBytes, ShiftRows, MixColumns, and Mux. As their names suggest, the first four correspond to the standard aes transformations. They are assembled with the multiplexer module, Mux, to form a closed loop, implementing a complete aes round. The Mux module opens the loop for plaintext acquisition during the initial round and closes it afterwards. The AddRoundKey module has a dedicated bus (ciphertext) through which ciphertext is output after the final round. The MixColumns module is bypassed during the final round. SubBytes is the only clocked module (all the others being purely combinatorial blocks). This allows, as mentioned before to complete an encryption round in one clock cycle. This loop architecture features a long data propagation path. Consequently, the design's critical delay path is located in the enm. The nominal clock frequency of this aes implementation is 100 MHz.

Giraud's One-Bit Attack

This section recalls Giraud's dfa [START_REF] Giraud | DFA on aes[END_REF] (hereafter "Giraud's one-bit attack"). The attack is based on the restrictive assumption that the opponent can inject one-bit faults. We use the following notations:

M i the algorithm's state at the end of round i K i i-th round key number C a correct ciphertext D a faulty ciphertext
The attacker injects a single bit fault into one byte of state M 9 just before the final round's SubBytes operation. This allows to retrieve the last round key K 10 .

Consider the i-th state byte just before the final round M 9

i . The corresponding byte index in the ciphertext is ShiftRows(i). As per the aes' specifications, for all i ∈ {0, . . . , 15}:

C ShiftRows(i) = SubBytes(M 9 i ) ⊕ K 10 ShiftRows(i) (2)
If a one-bit fault e,5 is injected into the j-th byte of M 9 , we obtain at index j:

D ShiftRows(j) = SubBytes(M 9 j ⊕ e) ⊕ K 10
ShiftRows(j) (3)
and for index i ∈ {0, . . . , 15}\{j}:

D ShiftRows(i) = SubBytes(M 9 i ) ⊕ K 10 ShiftRows(i) (4)
Hence, a comparison (differentiation) between the correct and the faulty ciphertexts leaks information both on the fault's position and on the aes key. For i ∈ {0, . . . , 15}\{j}, equations 2 and 4 yield:

C ShiftRows(i) = D ShiftRows(i) (5)
This allows to identify the faulty byte's index j because the only index for which C ⊕ D is nonzero is ShiftRows(j). Moreover, at index j, equations 2 and 3 yield:

C ShiftRows(j) ⊕ D ShiftRows(j) = SubBytes(M 9 j ) ⊕ SubBytes(M 9 j ⊕ e) (6)
Equation 6 is then solved for the eight possible values of e (the only hypothesis made on e is that e = 2 i for some i). This provides a set of candidates for M 9 j . At this point, new one-bit fault injections targeting the same byte are required to reduce the solution set to one item, namely M 9 j . The probability to find M 9 j in three attempts (i.e. three different faults) is ≃ 99%. K 10 ShiftRows(j) is then calculated from M 9 j , using equation 2:

K 10 ShiftRows(j) = C ShiftRows(i) ⊕ SubBytes(M 9 j ) ( 7 )
Equation 7 shows that the attack works independently on each round key byte. Indeed, no MixColumns transformation follows the fault injection and MixColumns is the only transformation capable of propagating faults between bytes. Consequently, the whole round key K 10 can be progressively retrieved by the attacker. Finally, knowing K 10 , the secret key K is found by reversing the key expansion algorithm.

The Attack Process

The experimental set-up (Figure 6) consists in a test chip board (tcb), a clock fault generator (cfg), and a computer. The tcb embeds the aes implementation described in the previous section. Its clock is provided by the cfg. The secret key, the plaintext and the encryption start signal are transmitted to the tcb from the computer via a serial link. As encryption ends, the ciphertext is offloaded via the serial link. The tcb provides a trigger signal to the cfg indicating the exact beginning of the encryption process. This is done to ease synchrony between fault injection and aes rounds. Note that the trigger could be replaced by the inspection of power consumption (spa) to precisely locate the aes rounds in time.

The cfg generates a 100 MHz clock signal continuously fed into the tcb. When the trigger signal indicates that encryption was started, a countdown is started. As the countdown ends a faulting period is produced during round nine, as required by Giraud's one-bit attack. The serial link between the computer and the cfg is used to define the exact value of the faulting period decrement ∆.

The computer runs a test campaign as described in Algorithm 1.

The generation of a faulting period is automatic and hence not explicitly mentioned in the pseudo-code. Indeed, the trigger signal sent from the tcb to the cfg indicates the encryption's launch and causes a faulting period generation during the ninth round as shown in Figure 4. As the faulting period gradually decreases, more set-up time violation faults appear in the calculation process. The resulting faulty ciphertexts are ordered by decreasing faulting periods. Then, Algorithm 1 Test Campaign Pseudo-Code send the key K and the plaintext M to the test chip. ∆ ← 0. while (clock period > ∆) do encrypt and retrieve the ciphertext ∆ ← ∆ + δt end while faulty ciphertexts are successively compared with the correct ciphertext. This allows identifying the faulty bytes (Eq. 6 and Eq. 7). In addition, for each ciphertext byte, a list of induced faults is built in order of appearance. Assuming that the first injected fault stems from a one-bit fault induced just before the last SubBytes, we build the corresponding set of guessed bytes for K 10 from equations 6 and 7. For the test campaign described in Algorithm 1, we obtain a set of guesses for every byte of K 10 . To reduce progressively these sets to singletons we inject different one-bit faults repeating the test campaigns with the same key but with different plaintexts. Indeed, each data bit arriving to the SubBytes's registers possesses its own logic path and propagation time (section 2). This propagation time highly depends on the data handled during encryption. Consequently, plaintext changes modify all propagation times. As propagation times vary, the injected one-bit faults differ with a 7/8 probability at the byte level.

As a result, one needs at least three (and sometimes four) test campaigns (same key and different plaintexts) to retrieve the entire round key. Finally, K is obtained by reversing the key expansion process.

Experimental Results

A first experiment targeting the final aes round was conducted to test the cfg. We implemented successfully Giraud's one-bit attack as well as its extension to two bits.

Injecting multiple-bit faults To test the hypothesis that a decrease in the faulting period causes more internal faults, we targeted the aes' tenth round while progressively increasing ∆.

As expected, the comparison of correct and faulty ciphertexts reveals that the device progressively transits from normal operation to multi-bit faults. This is done by exhibiting none, one-bit, two-bit and multiple-bit faults as describe in Figure 7. Note that the above is not an attack but an experimental fault characterization experiment where the aes plays the role of a big propagation delay cause. Figure 7 shows the faults' timing and nature as a function of the faulting period's duration (the horizontal axis). Each horizontal bar, associated with a byte number, uses a color code to reflect the nature of faults (no fault, one-bit fault, two-bits fault, more than two bits fault) and their point of appearance in time. The first one-bit fault occurs in byte 3 for a faulting period of 7585 ps (= 10000-69×35 ps), and the last fault appears on byte 0 for a 5800 ps faulting period. With the exception of bytes 2 and 8, the first fault is always a one-bit fault. This is compatible with the theoretical fault genesis model introduced at the beginning of this paper. Figure 8 shows how the number of byte candidates is progressively reduced. Each of the four sets in the example is associated with a different one-bit fault e i (for i ∈ {1, 2, 3, 4}). Here, the correct round key byte, 0x25, is found at the intersection of sets associates with e 1 , e 3 and e 4 . The fact that the set associated with e 2 stands apart indicates that the one-bit fault assumption about e 2 was wrong. Leaving aside e 2 , the set of guesses is reduced by successive sets intersections until a singleton containing the round key byte is reached. However, taking into account set e 2 requires the more complex intersection process described in Algorithm 2.

Algorithm 2 Determining the correct key byte

i ← 2 S ← e1 while (|S| ̸ = 1) do if S ∩ ei = ∅ then S = S ∪ ei else S = S ∩ ei end if i ← i + 1
end while return the single element of S as the correct key byte.

The probability of injecting successively a one-bit and a two-bits fault when reducing the faulty period can be estimated from Figure 7, where bytes 8, 5, 2 and 0 are counterexamples. Many examples seem to indicate that this probability is greater than 70% (experimentally).

The probability to cause successively one-bit, two-bits, and three-bits faults seems to be 50% (experimentally). These statistics were obtained thanks to the very small value of the faulty-period granularity, namely 35ps. This resolution seems to alow the progressive accumulation of faulty propagation pathes as the faulting period is decreased.

Implementation of Giraud's one-bit attack We wrote a script that plays the elementary test campaign described in section 3.3. The test was run several tens of times, for at least five different plaintexts per run. We always found the correct round key. The probability of injecting one-bit faults was found to be ≥ 90%.

Extension to two-bits attack The ability to inject and identify two-bits faults has prompted us to extend the attack. Note that two-bit attacks can defeat onebit parity check countermeasures. To that end, we need to solve equation 6 with the assumption that e is a two-bits fault. The only adverse effect is an increase in the cardinality of potential solution sets. The experiment was successful: 13 to 14 round key bytes were found on average with the automated test campaign. This allowed to exhaustive-search the whole key.

Repeating the campaign with many different plaintexts, we were always able to inject two-bits fault at every byte location, even if we never succeeded in injected two-bits faults simultaneously in all bytes.

Conclusion

This paper described a new fault injection technique based on clock reshaping during a precise cycle to cause setup time violation faults. The new technique would be a good candidate to benchmark safe-error [START_REF] Yen | Checking before output may not be enough against fault-based cryptanalysis[END_REF] or differential behaviorial [START_REF] Robisson | Differential behavioral analysis[END_REF] analysis This technique is inexpensive, efficient and non-intrusive. As such, it underlines the importance of further research in countermeasure design.

Fig. 1 .

 1 Fig. 1. Synchronous Representation of Digital ICs

Figure 2 Fig. 2 .

 22 Figure 2 illustrates a correct clock signal, clk, and a modified clock meant to cause faults, faulty clk.

Fig. 3 .

 3 Fig. 3. Faulty (faulty clk) Clock Signal Generation

t1 t 1 t 2 both

 2 signals are represented on a 2.0V/div scale 8.2ns

Fig. 4 .

 4 Fig. 4. faulty clk (uppermost signal) and aes start (lowermost signal)

Fig. 5 .

 5 Fig. 5. aes Structure

Fig. 6 .

 6 Fig. 6. Experimental Set-up

Figure 7

 7 Figure 7 reports fault characterization experiments conducted with a constant key and two different plaintexts.

Fig. 7 .

 7 Fig. 7. Two fault injection experiments: same key, different plaintexts

2 Fig. 8 .

 28 Fig. 8. Reducing the number of bytes candidates

That is: e = 2 i for i = 0, . . . , 7.


Juliette Pénichoux

Thierry Moreau

Aurélien Latouche
email: aurelien.latouche@cnam.fr

Simulating recurrent events that mimic actual data: a review of the literature

come

Introduction

Simulation studies that are performed for assessing statistical methods used in recurrent events analysis often rely on standardized event generation processes that cannot take into account some features of the data [START_REF] Metcalfe | The importance of varying the event generation process in simulation studies of statistical methods for recurrent events[END_REF]. Indeed, it is well recognized in the literature on recurrent events that in most medical research applications, the underlying process is likely to be complex, involving both heterogeneity in the population, and eventdependence [START_REF] Box-Steffensmeier | Repeated events survival models: the conditional frailty model[END_REF]. Heterogeneity arises if some subjects are more prone than others to experience events, owing to unkown or unmeasured factors, e.g. associated with lifestyle or genetic predisposition. This leads to intrasubject correlations between the waiting times separating consecutive events. Event-dependence refers to the case when the occurrence of one event may render the occurrence of subsequent ones more or less likely, e.g. as a result of a learning or immunity process, of a weakening of the body after an event or of an intervention subsequent to each event. The aim of this paper is to review the simulation methods that are used in the literature, focusing on those which generate data according to the features of actual data that are encountered in practice, so that they could be recommended for wider use. While, in many articles, heterogeneity is included in the model generating the simulated data, such is not the case of event-dependence, as will be seen. This work aims to describe the useful generating methods in detail and to underline the relevance of each of them for considering data of different types. In section 2, the literature is reviewed in order to assess how recurrent events processes are generated, and how they match the events dynamic of the biomedical data motivating these articles. The simulations performed are compared with the clinical studies considered. Section 3 describes methods for generating data according to the scenario of interest. A discussion follows in which the benefits of the methods are compared and some guidelines are suggested.

2 Results of the literature search

Introduction to the intensity models used

The simulation studies included in the following review mostly rely on modeling intensity or include processes that can be written as intensity-based processes. We focus on multiplicative proportional hazards intensities as they are the most frequent in the literature. The intensity process may be defined in one of two common timescales: the calendar-time and the gap-time. In calendar-time, the time denotes the time elapsed since the origin. In gap-time, the time is reset to zero when an event occurs, so that the time corresponds to the time elapsed since the previous event [START_REF] Kelly | Survival analysis for recurrent event data: an application to childhood infectious diseases[END_REF].

In calendar-time, for a vector of external covariates X, the common multiplicative proportional intensity λ(.) can be written:

λ(t|X = x) = λ 0 (t) exp(β ′ x). (1)
In gap-time, the common multiplicative proportional intensity becomes:

λ(t|X = x) = h 0 (t -T N (t -) ) exp(β ′ x). (2)
The models are described in more detail in section 3.

How recurrent events processes are generated in the literature

We conducted a literature research in six journals: Biometrics, Biometrika, Computational Statistics and Data Analysis, Journal of the Americal Statistical Association, Lifetime Data Analysis and Statistics in Medicine. This search was performed using the journals' webpages on October 10 th , 2011, focusing on articles published between 1996 and this date. In total, we found 135 articles that featured the term 'recurrent events', 'repeated events' or 'repeated failures', in the title and/or in the key words. Among these, 102 featured simulation studies of recurrent events. The generating intensity processes could be classified according to timescale, consideration of heterogeneity, and consideration of event-dependence. In the simple case where the baseline intensity is constant, the process can be expressed equivalently in either timescale (see Section 3), so Tables 1 and2 include a category "constant baseline" in addition to the categories "gap-time" and "calendar-time". A total of 92 articles could be classified in this way. In the sequel, the focus will be on these 92 articles. The 10 remaining articles were not classified, for example because the focus was not on time of events but only on the number of events over given intervals.

The results are summarized in Table 1 in which simulation studies are classified according to timescale and type of population (homogeneous or heteregeneous), and in Table 2, in which simulation studies are classified according to timescale and whether or not the generation processes included event-dependence. Note that one article can appear in more than one cell in each of these tables if several scenarios are considered in the simulation study.

Among the 92 articles, 54 (59%) featured exlusively processes with constant baseline. Note that a constant baseline does not exclude time-dependent covariates. Only 13 (14%) considered processes with a baseline that was a function of time since the previous event (gap-time), and 28 (30%) featured processes with a non-constant calendar baseline.

Among the articles with a baseline that varied with time, 3 articles [START_REF] Agustin | Goodness-of-fit of the distribution of time-to-first-occurrence in recurrent event models[END_REF][START_REF] Gonzalez | Modelling intervention effects after cancer relapses[END_REF][START_REF] Han | Parametric latent class joint model for a longitudinal biomarker and recurrent events[END_REF]] also considered a timescale that was in-between the gap-time and the calendar-time: some events caused a "perfect repair" that led the baseline to restart at 0, while for other events the baseline continued as in a calendar timescale.

Heterogeneity among subjects was widely considered: 80 articles (87%) included in the generation process a subject-specific random effect that was usually gamma-distributed.

Event-dependence was more rarely considered : only 14 (15%) of the simulation studies feature this scenario. It was modeled through an internal covariate, an event-specific baseline or an event-specific random effect.

In addition, Table 1 suggests that the situation of a heterogeneous population is seldom considered when the baseline hazard is time-varying and formulated in gap-time (in 7 articles among the 13 featuring processes generated according to this timescale). Table 2 suggests that the situation of event-dependence is less widely considered in simulation studies using a time-varying baseline hazard in a calendar timescale (in 2 articles [START_REF] Gonzalez | Modelling intervention effects after cancer relapses[END_REF][START_REF] Kvist | Repeated events and total time on test[END_REF] among the 28 using this timescale).

The list of papers included in this review are available as supplementary material.

Comparison of generating process and data of interest in two examples

In most of the simulation studies considered, the statistical methods investigated were also used in the analysis of biomedical data. In this section, we focus on two examples of applications that are commonly encountered in the recurrent events analysis setting, and investigate whether the processes generated in the simulation study have similar properties to those expected in the data of interest.

Recurrent events associated with asthma

Asthma is the disease of interest in several studies involving recurrent events. Three articles in our literature search analyze data from asthma studies [START_REF] Cai | Marginal means/rates models for multiple type recurrent event data[END_REF][START_REF] Chen | The analysis of multivariate recurrent events with partially missing event types[END_REF][START_REF] Cook | Robust tests for treatment effects based on censored recurrent event data observed over multiple periods[END_REF]. The events of interest are not always the same. The focus was on hospitalisations or physicians' visits attributable to asthma in the study by Cai and Schaubel [START_REF] Cai | Marginal means/rates models for multiple type recurrent event data[END_REF], on exacerbations in the study by Chen and Cook [START_REF] Chen | The analysis of multivariate recurrent events with partially missing event types[END_REF], and on coughing episodes in the study by Cook et al [START_REF] Cook | Robust tests for treatment effects based on censored recurrent event data observed over multiple periods[END_REF]. However, these different types of events are expected to reflect the general evolution of the disease.

Asthma is expected to be a stable disease that does not evolve much with time, at least in adults [START_REF] Duchateau | Evolution of recurrent asthma event rate over time in frailty models[END_REF], which is the population of interest in two of the studies mentioned [START_REF] Chen | The analysis of multivariate recurrent events with partially missing event types[END_REF][START_REF] Cook | Robust tests for treatment effects based on censored recurrent event data observed over multiple periods[END_REF]. However, in young children, Cai and Schaubel [START_REF] Cai | Marginal means/rates models for multiple type recurrent event data[END_REF] found that the risk of event was increased strongly between 0 and 2 years old and was decreased after 3 years old. Heterogeneity is to be expected as it is well-known that there are various degrees of severity of asthma [START_REF] Wildfire | Development and validation of the composite asthma severity index-an outcome measure for use in children and adolescents[END_REF]. The fact that asthma is considered as a stable disease suggests that at least in adults there is no event-dependence and no time-dependence.

The three simulation studies considered [START_REF] Cai | Marginal means/rates models for multiple type recurrent event data[END_REF][START_REF] Chen | The analysis of multivariate recurrent events with partially missing event types[END_REF][START_REF] Cook | Robust tests for treatment effects based on censored recurrent event data observed over multiple periods[END_REF] all relied on simulated processes that had a constant baseline and featured no event-dependence but were heterogeneous across the sample. This matches the expected events dynamic of the data processes, except maybe for the asthma study by Cai and Schaubel [START_REF] Cai | Marginal means/rates models for multiple type recurrent event data[END_REF] which focuses on children, for whom time-dependence is to be expected.

Recurrent infections in patients with chronic granulomatous disease

Among the articles considered in the review of simulation studies, seven [START_REF] Box-Steffensmeier | Repeated events survival models: the conditional frailty model[END_REF][START_REF] Lin | Semiparametric regression for the mean and rate functions of recurrent events[END_REF][START_REF] Sun | A class of accelerated means regression models for recurrent event data[END_REF][START_REF] Sun | A class of Box-Cox transformation models for recurrent event data[END_REF][START_REF] Sun | Marginal regression models with time-varying coefficients for recurrent event data[END_REF][START_REF] Tong | Variable selection for recurrent event data via nonconcave penalized estimating function[END_REF][START_REF] Zeng | Semiparametric transformation models with random effects for recurrent events[END_REF] investigated an application to a dataset from a clinical trial that was originally introduced by Flemming and Harrington [START_REF] Fleming | Counting processes and survival analysis[END_REF], which studies the effect of gamma interferon in patients with chronic granulomatous disease (CGD). This double-blind trial focused on the occurrence of recurrent infections among 128 patients randomised between placebo and gamma interferon.

The expected clinical course of CGD is not well-known since it is a relatively rare condition (about one in 250000 individuals) [START_REF] Van Den Berg | Chronic granulomatous disease: The european experience[END_REF]. It is therefore not obvious whether one should expect time-varying hazards, but in the absence of evidence this possibility should not be discarded. Moreover, one study [START_REF] Martire | Clinical features, long-term follow-up and outcome of a large cohort of patients with Chronic Granulomatous Disease: an Italian multicenter study[END_REF] pointed out that the incidence for some particular types of infections increases during follow-up. Similarly, event-dependence has received little attention, but the analysis of the dataset of the gamma interferon trial by Box-Steffensmeier and De Boef [START_REF] Box-Steffensmeier | Repeated events survival models: the conditional frailty model[END_REF] provides an indication of event dependence, specifically that the occurrence of an infection makes recurrence more likely. Additionally, it is well-known that there is much clinical heterogeneity among patients with CGD [START_REF] Dinauer | Chronic granulomatous disease and other disorders of phagocyte function[END_REF][START_REF] Winkelstein | Chronic granulomatous disease. Report on a national registry of 368 patients[END_REF].

Among the seven simulation studies presented in articles featuring an application to the CGD data, two considered only homogeneous processes [START_REF] Sun | A class of Box-Cox transformation models for recurrent event data[END_REF][START_REF] Sun | Marginal regression models with time-varying coefficients for recurrent event data[END_REF]. A time-varying baseline in a calendar timescale was considered in four [START_REF] Sun | A class of Box-Cox transformation models for recurrent event data[END_REF][START_REF] Sun | A class of accelerated means regression models for recurrent event data[END_REF][START_REF] Tong | Variable selection for recurrent event data via nonconcave penalized estimating function[END_REF][START_REF] Zeng | Semiparametric transformation models with random effects for recurrent events[END_REF] of the seven studies, which is more than expected given the total proportion found in the review (30%). Event-dependence, however, was considered only in one simulation study [START_REF] Box-Steffensmeier | Repeated events survival models: the conditional frailty model[END_REF].

Simulations based on intensity-based models

The simulation studies included in the review in Section 2 were based on the modeling intensity or could be written as intensity-based processes. This intensity could be additive or more commonly multiplicative. Accelerated failure time models, which correspond to non-proportional intensities, are seldom used in this setting. In this section, we present examples of intensity processes for each scenario defined in Section 2. We focus on multiplicative proportional hazards intensities as they are the most frequent in the literature.

For a given suject, let T j denote the time elapsed from the origin to the j th event, T 0 = 0, and let W j = T j -T j-1 denote the j th waiting time between two consecutive events. The W j are referred to as gap times.

Whichever intensity process is considered, a first step is to generate a censoring time C, and external covariates X if necessary. Then, the event times T j (or the gap times W j ) are generated until the censoring time is reached. The censoring time is usually generated according to an exponential or uniform distribution.

Standard model in each timescale

If the baseline is constant, the recurrent event process can be considered equivalently as an homogeneous Poisson process or a sequence of exponentially distributed gaps, which makes the generation process straightforward. A constant baseline is the most widely considered scenario. However, it may be beneficial to consider different generation processes that involve time-varying baselines using the two possible timescales.

Calendar timescale

When the relevant timescale is the time elapsed since inclusion, the baseline is specified in a calendar timescale. In particular, calendar-time is the proper timescale in the case of a disease that evolves in the long run. For example, the calendar timescale is used when analyzing exacerbations in patients with chronic bronchitis, which often arise as a consequence of impaired lung function due to the underlying disease [24, p233].

The standard model is a non-homogeneous Poisson process, with an intensity given by (1). The most common choice is a Weibull intensity, with λ 0 (t) = λνt ν-1 where λ is a scale parameter and ν a shape parameter. The Weibull distribution allows the baseline to increase (if ν > 1) or decrease (if ν < 1) over time.

There are several ways to simulate a non-homogeneous Poisson process.

The first is the inversion method. The j th gap time W j of a non-homogeneous Poisson process has a cumulative distribution function F j (.) equal to:

F j (w) = 1 -P (W j > w) = 1 -P {N (T j-1 + w) -N (T j-1 ) = 0}. (3)
The number of events observed for a non-homogeneous Poisson process with intensity λ(.) between two times s and t follows a Poisson distribution with parameter t s λ(u)du, which leads to:

N (T j-1 + w) -N (T j-1 ) ∼ P T j-1 +w T j-1 λ(u)du . (4)
And (3) becomes :

F j (w) = 1 -exp - T j-1 +w T j-1 λ(u)du . (5)
The successive event times are then simulated according to the following algorithm:

1. Let T 0 = 0.

2. For the j th event, simulate

V j ∼ U [0, 1]. 3. Let W j = F -1 j (V j ) and T j = T j-1 + W j .
4. If T j < C, return to step 2 for the next event. Otherwise, the j th event time is censored: set

T j = C and W j = C -T j-1 .
In the Weibull intensity in model ( 1), the event times simulated are given by:

T j = T j-1 + W j = - log(1 -V j ) λ exp(β ′ x) + T ν j-1 1 ν . Note that if V j ∼ U [0, 1], then 1 -V j ∼ U [0, 1],
and log(1 -V j ) can be replaced by log(V j ) in this expression.

A limitation of this method is that it requires F j (.) to be invertible. If not, numerical methods may be used to obtain W j such that F j (W j ) = V j .

Another common method is referred to as thinning [START_REF] Lewis | Simulation of nonhomogenous poisson processes by thinning[END_REF]. First, a value λ is chosen such that ∀t, λ(t) ≤ λ. If times T * 1 , T * 2 ... are generated from a homogeneous Poisson process with rate λ, the process in which each time T * j has a probability λ(T * j ) λ of being an actual event time, and the complementary probability 1 -

λ(T * j ) λ
to be discarded, is a process with intensity λ(.). The corresponding generating algorithm is as follows:

1. Let T 0 = 0, T * = 0, j = 1.

2. Simulate E ∼ Exp(λ), and let T * = T * + E.

Simulate

V ∼ U [0, 1]: if V ≤ λ(T * j )
λ , T * = T j is the j th event time. Then set j = j + 1, and return to

step 2. If V > λ(T * j )
λ , T * is not an event time, go back directly to step 2.

The assumption that the intensity λ(.) is an upper-bounded by a value λ is not as restrictive as it might seem, since processes are usually generated over a finite time interval [START_REF] Cook | The Statistical Analysis of Recurrent Events[END_REF]. This second method is often slower than the inversion method and is usually used when the inverse of the cumulative distribution functions F j (.) is not easily calculated.

Gap timescale

When the relevant timescale is the time elapsed since the previous event, the baseline is specified in gaptime. The gap-time is relevant for example when generating a risk of relapse after a complete response to treatment [START_REF] Gonzalez | Modelling intervention effects after cancer relapses[END_REF]. The basic model is a renewal process, with an intensity given by [START_REF] Box-Steffensmeier | Repeated events survival models: the conditional frailty model[END_REF].

The Weibull distribution is again a natural choice, with h 0 (w) = λνw ν-1 allowing for a risk of event that increases or decreases starting from the occurrence of the previous event.

The successive i.i.d. gap times can be simulated in the same way survival times are simulated when the hazard function is time-varying [START_REF] Bender | Generating survival times to simulate Cox proportional hazards models[END_REF]. The cumulative distribution function for the j th gap time W j is :

F j (w) = 1 -exp{H 0 (w) exp(β ′ x)}. Since F j (W j ) ∼ U [0, 1], exp{-H 0 (w) exp(β ′ x)} ∼ U [0, 1].
The gap time W j can therefore be simulated by generating V j ∼ U [0, 1], and setting

W j = H -1 0 - log(V j ) exp(β ′ x) .
In the case of the Weibull distribution, it leads to

W j = - log(V j ) λ exp(β ′ x) 1 ν .
New gap times are simulated until the sum of the gap times exceeds the censoring time: the last gap time is then censored.

When the cumulative hazard is not easily invertible, one may use numerical methods to solve H 0 (W j ) exp(β ′ x) ≈ -log(V j ).

Again, an alternative method is to use an acceptance-rejection type algorithm [27, p170] adapted from the thinning method proposed by Lewis and Shedler for non-homogeneous Poisson processes [START_REF] Lewis | Simulation of nonhomogenous poisson processes by thinning[END_REF].

Denote h an upper-bound of the hazard function h(.) over the period considered. Each gap time W j is simulated according to the following algorithm:

1. Set W 0 = 0 and W * = 0.

Generate E ∼ Exp(h), and let W

* = W * + E. 3. Generate V j ∼ U [0, 1]. If V j ≤ h(W * )
h , W * is the j th simulated gap time W j . Otherwise return to step 2.

Again, new independent gap times are simulated until the sum of the gap times exceeds the censoring time: the last gap time is then censored.

Heterogeneity

Heterogeneity in the population is often realistic in biomedical studies and is widely considered in simulation studies designed to evaluate the properties of statistical methods. In both timescales, heterogeneity can be easily introduced by generating a subject-specific random effect. This random effect is usually assumed to follow gamma or lognormal distribution, but other distributions are possible. For instance, one may use a binary random effect to investigate model properties when binary covariates are omitted.

In the calendar timescale, the resulting process is referred to as a mixed Poisson process, i.e. the individual processes are Poisson processes with subject-specific intensities. The usual multiplicative proportional intensity can then be written :

λ(t|U = u, X = x) = λ 0 (t)u exp(β ′ x), ( 6
)
where u is the subject-specific random effect. A simple simulation scheme consists in first generating independent individual frailties according to the chosen distribution. The individual intensities are obtained from [START_REF] Han | Parametric latent class joint model for a longitudinal biomarker and recurrent events[END_REF]. The individual processes can then be simulated as Poisson processes as described in Section 3.1.1.

In gap-time, the result is a shared frailty model, with the corresponding intensity :

λ(t|U = u, X = x) = h 0 (t -T N (t -) )u exp(β ′ x).
Similarly, independent individual frailties are first generated. The gap times for each subject can then be simulated as described in Section 3.1.2.

Generating recurrent event with event-dependence

Event-dependence is less widely considered in simulation studies than heterogeneity. However, for many recurrent events data, it is realistic to consider that the occurrence of an event may alter the risk of occurrence of new events. For this reason, it is useful to investigate the properties of statistical methods under eventdependence, e.g. when the occurrence of an event increases the risk of occurrence of subsequent ones.

Event-dependence in gap-time

It should be noted that a baseline that is a function of time in gap-time constitutes a form of event-dependence, since the occurrence of an event causes the hazard to reset back to its original level. However, additional event-dependence may be appropriate. A common assumption is to consider that the occurrence of an event increases the risk of further ones (e.g. by making subjects more frail and therefore prone to experiencing events), or on the contrary decreases it (by making subjects immune or more vigilant). A natural way to generate such data is to take an event-specific baseline or an internal time-dependent covariate that induces dependence [START_REF] Kalbfleisch | The statistical analysis of failure time data[END_REF], as did Cook et al [START_REF] Cook | Robust estimation of mean function and treatment effects for recurrent events under event-dependent censoring and termination: Application to skeletal complications in cancer metastatic to bone[END_REF] in their simulation studies.

A simple example is to multiply the baseline by α j-1 with j -1 the number of previous events. Setting α < 1 will cause the risk of event to decrease with the number of past events while setting α > 1 will cause the risk to increase with the number of events. This is similar to including the number of previous events as an internal covariate in a Cox-type model with associated parameter ln(α).

However, if the risk increases exponentially after each event, the event process may "explode", with a number of events that tends to infinity for a finite time [START_REF] Gjessing | Recurrent events and the exploding Cox model[END_REF]. A common solution is to consider the risk to be stabilised after a pre-specified number of events, e.g. 4 events. Considering a Weibull baseline and a risk that is stable after 4 events, the data generation model without covariates becomes :

h 0j (w) = λνw ν-1 α j-1 j = 1, 2, 3, 4 λνw ν-1 α 4 j > 4.
Usually, when considering a gap-timescale, the internal covariate that induces the dependence between events is constant for each gap-time and includes the number or timing of previous events. The successive event times W j can then be simulated as described in Section 3.1.2, with a different cumulative distribution function for each event.

However, more general processes can be considered, with internal time-dependent covariates that can be non-constant over the at-risk interval for a given event. The methods described in the next section for calendar-time are general and valid for such processes.

Event dependence in calendar-time: generating dynamic covariates

Inducing event-dependence in a calendar timescale is usually done by considering an intensity with an internal (or dynamic) covariate that depends on the past of the process. An obvious choice is N (t -) to allow the intensity to increase or decrease with the number of past events. Such a generating process was considered by Gonzalez et al [START_REF] Gonzalez | Modelling intervention effects after cancer relapses[END_REF]. However, just as in gap-time, the risk is that the number of events will increase enormously when considering N (t -) as a covariate. Other possible covariates are N (t -)/t, in order to simulate processes in which each event leads to an increased intensity, but where the effect of this increase fades with time, min(N (t -), K 0 ) to consider the case where the risk becomes stable after K 0 events, or N (t -)-N (t-u) u where u ≤ t, in order to focus on recent history (on a period of length u) as suggested by Elgmati et al [START_REF] Elgmati | A martingale residual diagnostic for longitudinal and recurrent event data[END_REF].

The inversion method described in Section 3.1.1 for non-homogeneous Poisson processes can be used for more general processes, with an intensity that depends on the past. Indeed, Cook and Lawless [24, p30] showed that:

P (W j > w|T j-1 = t j-1 , H(t j-1 )) = exp - t j-1 +w t j-1 λ(u|H(u))du ,
where H(t) is the history of the event process up to time t (not included). The intensity is therefore considered deterministic on the interval ]t j-1 , t j-1 + w]. This makes it possible to simulate a process defined by an intensity conditional on the past by the inversion method in the same way as a non-homogeneous Poisson process, by generating F j (W j ) ∼ U [0, 1] with F j (w) = P (W j ≤ w|T j-1 , H(T j-1 )).

For example, the process with intensity:

λ(t) = λ 0 (t) exp(β ′ x + ϕN (t -)),
has a cumulative distribution function for the j th gap time:

F j (w) = 1 -exp - t j-1 +w t j-1 λ 0 (u) exp(β ′ x + ϕ(j -1))du .
When the baseline intensity is Weibull, generating 1 -F j (W j ) = V j ∼ U [0, 1] leads to:

T j = T j-1 + W j = - log(V j ) λ exp(β ′ x + ϕ(j -1)) + T ν j-1 .
When considering the covariate min(N (t -), K 0 ) with a Weibull baseline intensity, the event times become:

T j = - log(V j ) λ exp(β ′ x + ϕ(min(j -1, K 0 ))) + T ν j-1 1 ν
.

More generally, an intensity process:

λ(t) = exp(g 0 (t) + I(N (t -) > 0)g 1 (t -N (t -)) + g 2 (N (t -)) + β ′ x),
which is a function of the time in a calendar timescale, the time in gap-time and the number of past events, leads to a cumulative distribution function of the j th gap time :

F j (w) = exp - w 0 exp{g 0 (t j-1 + u) + I(j -1 > 0)g 1 (u) + g 2 (j -1) + β ′ x}du . (7)
F j (W j ) ∼ U [0, 1] is then generated and W j is obtained from [START_REF] Kvist | Repeated events and total time on test[END_REF] analytically or numerically. In the case of the internal covariates N (t -)/t or N (t -)-N (t-u) u suggested before, a closed-form expression for W j cannot easily be obtained.

Another possibility is to use the thinning method described in Section 3.1.1, which was introduced by Lewis and Shedler [START_REF] Lewis | Simulation of nonhomogenous poisson processes by thinning[END_REF] for non-homogeneous Poisson processes, but was later proved to be valid when the intensity depends on the past [START_REF] Ogata | On Lewis' simulation method for point processes[END_REF]. In the algorithm described in Section 3.1.1, λ(T * ) can therefore depend on the past.

A final possibility that does not require an upper bound for the intensity or a numeric inversion of the cumulative distribution function is based on a discrete approximation. A similar method was used for example by [START_REF] Aalen | Survival and event history analysis : a process point of view[END_REF]. A partition of the time period is defined by intervals of length ∆t: [k∆t, (k + 1)∆t], k = 0, 1, 2, .... The length ∆t is chosen to be small, so that one can assume that a subject cannot have more than one event in each interval [k∆t, (k + 1)∆t]. The occurrence of an event during the interval [k∆t, (k + 1)∆t] is generated with an event probability λ(k∆t)∆t, which approximates the probability of an event occurrence corresponding to (4). The intensity is then updated for the next interval. For example, if ∆t corresponds to a day, the occurrence or not of an event is simulated every day. According to this approximation, simulating a process with an intensity λ 0 (t) exp(β ′ x + ϕ N i (t -) t ) where λ 0 (.) is expressed as a function of time in days leads to consider that the subject has an event on the k th day with a probability λ 0 (k) exp(β ′ x + ϕ N i (k-1) k-1 ). This method of simulating general event processes is simple, but can also be slow.

Discussion

This review complements the extensive simulation study conducted by Metcalfe and Thompson [START_REF] Metcalfe | The importance of varying the event generation process in simulation studies of statistical methods for recurrent events[END_REF]. Their study exemplified the fact that no model is adequate for all settings, and that it is necessary to investigate several situations in simulation studies, including heterogeneity and event-dependence, in order to have satisfactory knowledge of the properties of a given statistical method.

In this review, we show that heterogeneity has received much attention. By contrast, event-dependence is seldom considered in simulation studies, especially when the event-generation process is expressed in calendar-time. This is unfortunate since event-dependence is known to have a critical effect on statistical models, in terms of properties and interpretation of the covariates effects [START_REF] Cheung | Estimation of interention effects using first or multiple episodes in clinical trials: the Andersen-Gill model re-examined[END_REF]. One reason for neglecting this possibility may be that the generation of non-Poisson processes is not obvious. In this paper, the general relation between the intensity process and the successive gap times is recalled and examples of event-dependent processes are provided, in both timescales. Additionally, as noted by Bender et al [START_REF] Bender | Generating survival times to simulate Cox proportional hazards models[END_REF] for univariate survival times, most simulation studies use only constant baseline hazards corresponding to exponential gap times. However, it was pointed out by Duchateau et al [START_REF] Duchateau | Evolution of recurrent asthma event rate over time in frailty models[END_REF] that the choice of the timescale in the statistical model may considerably impact the results. This suggests that a given model may demonstrate different properties in simulation studies according to the timescale of the generation process, and that it would be wise to consider generation processes with time-varying baselines in both timescales.

Extensive simulation studies that investigate the properties of a given statistical method should therefore include, in our opinion, datasets generated by using a mixed Poisson process, processes with constant baseline, Weibull calendar-time, and Weibull gap-time baseline, and processes with event-dependence e.g. with the number of previous events included as a covariate in the intensity generation process. Such studies would assist the applied statistician in choosing the appropriate method in a particular biomedical context. When a statistical method is suggested for a biomedical context, the scenario considered for the simulation studies should match the expected dynamics of the event-process in this context. If a long-term evolution of the disease cannot be ruled out, calendar-time-varying intensity processes should be included. This would be the case for exacerbations in patients with chronic brochitis [24, p233]. A baseline varying with time in gap-time should be considered in the case when an event is followed by an intervention that leads to a kind of renewal of the subject. This was considered in particular in the context of cancer relapse where a complete response was possible [START_REF] Gonzalez | Modelling intervention effects after cancer relapses[END_REF]. Heterogeneity should in general be included since there is no biomedical context in which one can be sure that all risk factors are known and always measured. Event-dependence should be considered in contexts in which an event could alter the risk of subsequent events. For example, this was considered in a simulation study by Kvist et al [START_REF] Kvist | Repeated events and total time on test[END_REF], in the context of affective disorder, in which an important theory is the sensitisation theory, i.e. that the occurrence of an event may increase the risk of subsequent events.

Table 1 :

 1 Results of literature search : generated processes classified according to timescale and type of population (homogeneous/heterogeneous) (cells are non-exclusive)

		Homogeneity Heterogeneity	Total
	Constant baseline hazard	28 (30%)	62 (67%)	70 (76%)
	Baseline hazard in gap-time	8 (9%)	7 (8%)	13 (14%)
	Baseline hazard in calendar-time	12 (13%)	25 (27%)	28 (30%)
	Total	38 (41%)	80 (87%)	92

Table 2 :

 2 Results of literature search : generated processes classified according to timescale and inclusion of event-dependence (cells are non-exclusive)

		No event-dependence Event-dependence	Total
	Constant baseline hazard	64 (70%)	13 (14%)	70 (76%)
	Baseline hazard in gap-time	12 (13%)	3 (3%)	13 (14%)
	Baseline hazard in calendar-time	27 (30%)	2 (2%)	28 (30%)
	Total	83 (90%)	14 (15%)	92
		13


Baptiste Bouillot

Duyen Le Quang

Quang-Du Le

Jean-Michel Herri

Philippe Glenat

Pierre Duchet-Suchaux

Experimental procedure and results to measure the composition of gas hydrate, during crystallization and at equilibrium, from N2-CO2-C2H6-C3H8-C4H10 gas mixtures

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

EXPERIMENTAL PROCEDURE AND RESULTS TO MEASURE THE COMPOSITION OF GAS HYDRATE, DURING CRYSTALLIZATION AND

AT EQUILIBRIUM, FROM N 2 -CO 2 -C 2 H 6 -C 3 H 8 -C 4 H 10 GAS MIXTURES (1) BOUILLOT Baptiste, (1) LE QUANG Duyen, (1) LE QUANG Du, (1) HERRI Jean-Michel , (2) GLENAT, Philippe, (2) DUCHET-SUCHAUX, Pierre (1) Gas Hydrate Dynamics Centre, Ecole Nationale Supérieure des Mines de Saint-Etienne, France (2) TOTAL S.A, France

Context

-Study, from an experimental point of view, the gas hydrate formation of various hydrocarbon gas mixtures in an instrumented batch reactor.

-Determine the equilibrium points (P, T) and Composition of each phase (gas x G , liquid x L , hydrate x H ). -The molar composition at equilibrium for gas mixtures are given and present uncertainties inferior to 4% (for both gas phase and hydrate phase).

-The modeling part is presented in [6] and discuss the relevancy of thermodynamics in this case The offshore extraction of petroleum in deeper conditions increases day after day and favors the conditions of hydrate formation (low temperature and high pressure), representing a great concern in flow assurance. Actually, another parameter comes into play, because as the oil/gas fields matures, the fraction of water increases.

Experimental results

G x

Others experimental results on gas mixtures (example)

-Other data available in the full text.

-Most of the experimental data were obtain from quick crystallization process -One experiment from slow crystallization to detect a possible kinetic effect (see [6])

•--

  Δx i /x i = 3% (mol. Frac.), • ΔT = 0.2°C, • ΔP = 0.1 bar, • Δm w = 0.1g (of mass of injected water), • Δm w (liq. sample) = 0.1g, • ΔV R = 3% (Reactor Volume), • Δρ w = 10 kg.m -3 , • Δρ H β = 40 kg.m -3 . • Ρ w = 1000 kg/m 3 . • ρ H β = 790 kg/m 3 . Gas hydrates equilibrium experiments of hydrocarbon gas mixtures (N 2 -CO 2 -CH 4 -C 2 H 6 -C 3 H 8 -C 4 H 10 ) were performed in an instrumented batch reactor. Pressure and Temperature of equilibrium for quick crystallization (high crystallization rate) are given and are in accordance with literature data.

Figure 3 .

 3 Figure 3. Experimental procedure

Figure 2 .

 2 Figure2. Calculation procedure (see also[1])

Figure 4 .

 4 Figure 4. Experimental and literature [2,3,4] equilibrium data of pure CO 2 gas hydrates.

Figure 5 .

 5 Figure 5. Experimental and literature [3,4,5] equilibrium data of pure CH 4 gas hydrates.

Table 1 .

 1 Molar compositions of the studied gas mixtures (std. dev. ≈ 3%)

	Gas	CO 2	N 2	CH 4	C 2 H 6	C 3 H 8 n-C 4 H 10
	mixtures
	Gas 1	-	-	1	-	-	-
	Gas 2	1	-	-	-	-	-
	Gas 3	0.24	-	0.76	-	-	-
	Gas 4	0.225	-	0.775	-	-	-
	Gas 5	0.821	0.179	-	-	-	-
	Gas 6	0.057	-	0.915	0.028	-	-
	Gas 7	0.053	-	0.919	0.028	-	-
	Gas 8	0.052	-	0.919	0.029	-	-
	Gas 9	-	-	0.953	0.028	0.019	-
	Gas 10	-	-	0.864	0.057	0.053	0.026
	Gas 11	-	-	0.848	0.05	0.076	0.026
	Gas 12	0.055	-	0.917	0.028	-	-

Figure 1. Experimental set-up

Validation on others experimental results Pure CO 2 and Pure CH 4 :

  Experimental data of CH 4 -C 2 H 6 -C 3 H 8 gas mixture… … and numerical results in the table 2 below.

	Figure 6.
	Xgas (C1)	Xgas (C2)	Xgas (C3)	Xhyd (C1)	Xhyd (C2)	Xhyd (C3)	N Hyd .
	99.18 0.57 0.24	88.18 7.08 4.73	6.27
	98.96 0.77 0.27	86.60 7.95 5.45	6.46
	98.63 1.01 0.37	85.10 8.76 6.14	7.32
	98.18 1.33 0.49	84.43 8.94 6.63	6.65
	97.52 1.87 0.61	84.31 8.85 6.84	8.05
	97.16 2.08 0.77	83.90 8.74 7.36	8.90
	96.75 2.32 0.93	82.44 8.31 9.25	7.59
	96.39 2.48 1.13	81.96 8.12 9.93	8.38
	96.02 2.60 1.38	80.69 7.74 11.57 10.05


I Terreros

I Iordanoff

J L Charles

D Coupard

S Tcherniaieff

DISCRETE ELEMENT METHOD, A TOOL TO INVESTIGATE COMPLEX THERMO MECHANICAL BEHAVIOUR: APPLICATION TO FRICTION STIR WELDING

Keywords:

In material forming or cutting, the contact zone between the tool and the working piece is often very difficult to analyse because of diverse local phenomena. In the special case of FSW, a specific difficulty is the study of material mixing. In this work, the Discrete Element Method is applied as a tool to understand/propose/confirm physical scenarii involved in FSW process. A simple micro behaviour law that takes into account thermal, mechanical and material aspects is proposed and a two dimensional simulation based on this law is presented. The first results show qualitatively good agreement with real nuggets observations. One of the main advantages of DEM approach is the simplicity of the local input law compared with laws proposed at the continuous media mechanics level. Nevertheless, numerical aspect must be improved to carry out three dimensional simulations with reasonable calculation times.

INTRODUCTION

Friction Stir Welding (FSW) is a process invented at The Welding Institute (TWI) in the early 90's. This process allows to join metals without fusion or filler materials. The welds are created by the combined action of frictional heating and mechanical deformation due to a rotating tool.

At the present time the approach used to simulate FSW process are based on the Finite Elements Method (FEM) or on the Finite Differences Method (FDM) [START_REF] Ch | Thermomechanical and Figure 8: Proposed method microstructural modelling of the Friction Stir Welding process[END_REF][START_REF] Mccune | Modelling Friction Stir Welding with the Finite Element Method -a comparative study[END_REF][START_REF] Palm | Improved verification of FSW-process modelling relating to the origin of material plasticity[END_REF]. With those methods it is possible to simulate temperature fields [START_REF] Ch | Thermomechanical and Figure 8: Proposed method microstructural modelling of the Friction Stir Welding process[END_REF][START_REF] Mccune | Modelling Friction Stir Welding with the Finite Element Method -a comparative study[END_REF][START_REF] Palm | Improved verification of FSW-process modelling relating to the origin of material plasticity[END_REF], strain and residual strain fields [START_REF] Ch | Thermomechanical and Figure 8: Proposed method microstructural modelling of the Friction Stir Welding process[END_REF][START_REF] Mccune | Modelling Friction Stir Welding with the Finite Element Method -a comparative study[END_REF], heat source distribution, plastically-deformed zone and material flow while looking at deformation velocities [START_REF] Palm | Improved verification of FSW-process modelling relating to the origin of material plasticity[END_REF]. We can also carry out simulations of the material flow using a fluidmechanics based approach [START_REF] Heurtier | Mechanical and thermal modelling of Friction Stir Welding[END_REF].

Despite the efforts done, to achieve an accurate simulation of material flow when it is between a solid state and a "slurry" state is still a delicate problem to manage with classic FEM. A possible suggestion is to use the DEM for the modelling, because of the exceptional skills of this method to carry out calculations of solid flows. There are works that demonstrate the validity of DEM to simulate and understand wear problems (in this field the difficulty to simulate the flux has been always a big challenge) like [START_REF] Iordanoff | Modélisation du Comportement Tribologique des 3e Corps Solides et Gazeux[END_REF][START_REF] Richard | Thermal study of the dry sliding contact with third body presence[END_REF][START_REF] Seve | Using Discrete Models to Simulate Solid Third Bodies: Influence of the Inter-granule Forces on the Macroscopical Behaviour[END_REF][START_REF] Iordanoff | A Discrete Element Model to Investigate Sub Surface Damage due to surface polishing[END_REF], among others. * Postal address: Pio XII, 1, 5C, 01400, Laudio, Araba, Spain. Phone: +34 670 749 016. e-mail: iterrerosgalindez@gmail.com

NUMERICAL TOOL

GENERAL DESCRIPTION

The method described is based on smooth particle dynamic method. The material is considered as a set of discrete spheres that move under a force field. The forces are contact forces to simulate multi-body interactions. Their magnitudes depend on a thermomechanical law which will be briefly described in next paragraphs. Particle movements are calculated using an explicit algorithm to integrate the dynamic Newton law. It is vital to remember that each sphere doesn't represent any kind of particle, but all the spheres as a whole represent the behaviour of the material.

THERMOMECHANICAL LAW

As it is said above, the magnitude of contact forces is given by a thermomechanical law. First of all mechanical aspect is described, then thermal aspect, and finally a description of how mechanical and thermal aspects are linked will be given.

Mechanical part

The particles' (spheres in this study) interaction laws define micro mechanical properties of the material. Every particle is subjected to contact force. This force acts only when two particles geometrically interact and it is divided into three parts: repulsion, adhesion (both are energy conservative) and energy dissipation. Geometrical interaction δ allows the force calculation (Figure 1). • Repulsion is represented by a linear spring, whose stiffness is K. Repulsive force F r is:

F r = K • δ (1)
• Adhesion has been simplified to a parameter γ which depends on particles' temperature and accumulated displacement, γ = γ(T, U ). Relationship between γ, T and U is given in the section 2.2.3.

F a = γ = γ(T, U ) (2)
• The energy dissipation in the contact is due to the damping force written as:

F d = 2α K • M ij • δ (3)
Where α id the damping coefficient (< 1), δ the impact speed and M ij the equivalent mass of the contact.

The sum of the interaction forces is:

F contact = -K • δ + γ -2α K • M ij • δ (4)

Thermal part

For the thermal part, model proposed by Richard et al. [START_REF] Richard | Thermal study of the dry sliding contact with third body presence[END_REF] has been implemented. However, some modifications to reach the purpose of modelling FSW process have been carried out. The thermal analysis is divided in two parts, heat creation and heat diffusion.

It is supposed that all the dissipated mechanical energy is transformed into thermal energy. Following the model of Richard, created thermal power is described like:

Q ij = 2α K • M ij • δ2 (5)
Assuming that heat is equally distributed between particles i and j, the temperature rise of each particle because of thermal energy creation in ∆t time interval is described by the following equation:

dT i = 0, 5 • Q ij ∆t ρ i c pi V i (6)
At the moment of analysing heat diffusion, one assumption that differs from reference [START_REF] Richard | Thermal study of the dry sliding contact with third body presence[END_REF] is taken:

• The transmission surface between two spheres is supposed to be 4R 2 m (a 2R m × 2R m square), where R m is the average radius of spheres, instead of the contact surface between both of them.

This assumption is closer from a continuum medium like FSW than hypotheses from [START_REF] Richard | Thermal study of the dry sliding contact with third body presence[END_REF], which try to describe a real granular field. Then, the temperature rise of a particle is:

T + i = T - i + j dT + ij ( 7
)
Where T + i is the temperature of a particle i at the moment t + ∆t, T - i is the temperature of the same particle at time t and dT + ij is the thermal contribution of each j particle to i during ∆t period. From the viewpoint of discrete elements dT + ij is defined as:

dT + ij = 4λ ij R m ∆t(T - j -T - i ) ρ ij c pij V ij (8)
Where ρ ij , c pij , V ij and λ ij are the equivalent density, equivalent heat capacity, equivalent volume and equivalent conductivity of particle i respectively. Mixing equations 7 and 8, the equation that explains the temperature rise of particle i due to contact with the j particles around is:

T + i = T - i + ∆t j 4λ ij ρ ij c pij V ij • R m (T + j -T - i ) (9)

Coupling of mechanical and thermal aspects

Once described thermal and mechanical parts, the micromechanical law which links both of them will be explained. The following analysis is based on Aluminium alloy Al 2017-T451. This analysis micro-structural changes due to temperature and displacement. When material's temperature rises, the precipitates that harden it disappear and mechanical properties of material become weaker. Moreover, the hardest material is located near the trajectory of the pin of the tool (nugget), that is, the most displaced material which has the smallest grain-size. So, a relationship among material's temperature cycle, grain size, displacement and hardness can be established.

Having in mind those ideas, a surface that describes the cohesion between particles depending on the cycles mentioned above has been created. The curve A from Figure 3 describes the first approach to that surface. In axes x and y temperature and displacement are shown respectively. The third axe represents the hardness. The adhesion force of the contact between two particles is taken from this graphic and it is supposed proportional to hardness.

There are two surfaces apart A, which describes the law mentioned before. The surface marked as B represents the adhesion force between two particles of different plates and the other surface, marked as C, represents the adhesion between two particles of raw material of the same plate. Initially all the contacts between particles are classed in one of those two surfaces. Cohesion in the contact of particles of different plates will be null if they are cold, independently of the displacement. Nevertheless, cohesion will rise with the temperature until this curve crosses with A. From this moment the contact between the two particles won't be described by B curve but by A. The same happens with C. Equal material's particles' cohesion will be inversely proportional to the temperature until the temperature reach a limit. Here again we pass to the other curve to modell the adhesion force. The main idea behind that behaviour law is that adhesion decreases with temperature and increases with displacement.

The results of the modelling of FSW with the proposed micromechanical law show qualitatively good agreement with real nuggets observations, as it will be shown in section 3.

NUMERICAL INTEGRATION

A Verlet velocity algorithm is chosen to integrate the velocities and new positions of spheres over a small time step ∆t while accelerations have been calculated thanks to Newton laws described above [START_REF] Richard | Thermal study of the dry sliding contact with third body presence[END_REF][START_REF] Iordanoff | A Discrete Element Model to Investigate Sub Surface Damage due to surface polishing[END_REF].

GEOMETRY AND BOUNDARY CONDI-TIONS

For this work simulations with a two dimensional model have been carried out. The two walls surrounding the field are fixed and they are at ambient temperature.

RESULTS

MECHANICAL PROPERTIES

Figure 4 shows experimental hardness of Al 2017-T451 across the weld at the top and a cross section taken from the results of the simulations carried on at the bottom. The meaning of coloured parts will be explained in the next paragraph.

THERMAL PROPERTIES

Simulated temperature contours have at the beginning good agreement with works like [START_REF] Mccune | Modelling Friction Stir Welding with the Finite Element Method -a comparative study[END_REF] and [START_REF] Zhu | Numerical simulation of transient temperature and residual stresses in friction stir welding of 304l stainless steel[END_REF]. Nevertheless, the current model has problems to dissipate heat, and temperature fields have tendency to have parallel isotherm lines, as shown in Figure 6. This problem may be solved with the solutions proposed in section 4.

MATERIAL MIXING

The Discrete Element Method allows also to see how materials of both plates are mixed while tooling. Figure 7 is an example of that. The field before tooling is shown in the first picture and second one represents the field after the process.

CONCLUSIONS AND PERSPECTIVES

Looking at the results one can conclude that:

• Discrete Element Method models qualitatively good FSW process as similarities between real and experimental hardness, weld characteristics and temperature fields have been obtained.

• It is necessary to develop thermal boundary conditions because in the present work heat is evacuated only through two walls and this is not realistic.

• The proposed micromechanical law must be validated.

• Calculation time must be reduced.

In order to reduce the calculation time, a combination among Finite Element Method, Dynamic Discrete Element Method and Static Discrete Element Method is proposed. Effectively, FEM allows us to calculate rapidly thermal effects in areas where the movement of particles is always insignificant; the current dynamic DEM in areas where particles' movement is being important can be used; finally a static DEM which doesn't search contacts among particles in areas where movements are negliable can be used. In fact, 80% of calculation time is spend in searching contacts, and may be reduced using this method.

In Figure 8 we can see a schematic draw of this method.

Figure 1 :

 1 Figure 1: Geometrical interaction

Figure 2 :

 2 Figure 2: Grain-size of raw material and nugget

Figure 3 :

 3 Figure 3: Adhesion forces

Figure 4 :Figure 5

 45 Figure 4: Comparison of hardness & cohesion

Figure 5 :

 5 Figure 5: Cohesion Numerical simulation is qualitatively in accordance with experiences. The center of the weld is harder than raw material an there are two bands composed by weaker material surrounding it.

Figure 6 :

 6 Figure 6: Temperature fields

Figure 7 :

 7 Figure 7: Material mixing


Saman Noorzadeh

Bertrand Rivet

Pierre-Yves Guméry

Enhancing Fetal ECG Using Gaussian Process

This paper investigates the problem of the extraction of fetal ECG from the noisy abdominal channel. The method which is presented here is a multi-channel one which uses two data channels: one is the abdominal channel and the other is the reference channel. The reference channel, which is recorded from the chest of the mother, can well define the covariance function to model the maternal ECG contribution in the abdominal channel. The maternal ECG can then be extracted using this model which proposes a non-linear extraction. The reference signal is also proposed to be replaced with a 1-bit thoracic signal which can be cheaply recorded and has less time complexity to be processed.

I. INTRODUCTION

During the pregnancy the measurement of electrocardiogram (ECG) from the maternal abdomen can provide information about the health and physiological status of the fetus. However, the recorded abdominal ECG does not contain only the fetal ECG (fECG), but is contaminated by several sources of noise and disturbances, and the strongest source of disturbance is the maternal ECG (mECG). Therefore, the separation of fetal ECG from the maternal ECG has become a challenging task.

Several methods have already been applied to tackle this problem. One can cite Kalman filtering [START_REF] Sameni | A nonlinear bayesian filtering framework for ECG denoising[END_REF], and non-parametric modeling [START_REF] Niknazar | Fetal ECG extraction from a single sensor by a non-parametric modeling[END_REF], which are based on a single channel of data. While Kalman filter considers a prior information directly on the data, non-parametric method models the signal using only its second order statistics.

There are also multi-channel methods which have addressed this problem like blind source separation [START_REF] Zarzoso | Maternal and foetal ECG separation using blind source separation methods[END_REF], singular value decomposition [START_REF] Callaerts | Comparison of SVD methods to extract the foetal electrocardiogram from cutaneous electrode signals[END_REF], or adaptive filtering [START_REF] Widrow | Adaptive noise cancelling: Principles and applications[END_REF]. The most efficient method is the blind source separation based on independent component analysis (ICA) [START_REF] Zarzoso | Maternal and foetal ECG separation using blind source separation methods[END_REF]. However, this method need a lot of sensors (at least the sum of the dimensions of the source sub-spaces, typically three for each ECG) and this thus increases the cost of the device. Consequently, adaptive filtering [START_REF] Ferrara | Fetal electrocardiogram enhancement by time-sequenced adaptive filtering[END_REF] which can require only two sensors (one for recording the noisy signal and one for the recording of the noise reference) is an effective solution to reduce the cost of the device while keeping good fECG extraction performance.

The proposed method is a multi-channel one which benefits from a non-parametric method to model the ECG signal by Gaussian process (GP). The mECG signal is modeled in this study using another data channel, and it can then be estimated and extracted from the noisy abdominal channel in order to obtain the fECG through the rest of the signal. Although this method uses another reference channel like adaptive filters, it provides a non-linear solution for fECG extraction despite adaptive filters. Moreover it is proposed to use a 1-bit reference channel instead of the full range (e.g. 16 bits) one, which can reduce the further cost of the recording device and is more memory efficient.

The rest of the paper is organized as follow: the method and its principals are described in section II and the results of fetal ECG extraction using the proposed method is shown and compared to that of adaptive filtering in III. The conclusions are finally presented in section IV.

II. EXTRACTION OF FETAL ECG

As proposed in the previous studies an ECG signal can be modeled by GP using its second order statistics [START_REF] Niknazar | Fetal ECG extraction from a single sensor by a non-parametric modeling[END_REF], [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. Therefore, an ECG signal as a quasi-periodic signal, noted as x(t), can be modeled using a GP by its mean and covariance function as

x(t) = GP m(t), k(t, t ) , (1)
where m(t) is the mean function that is set to 0 in this study, and k(t, t ) is the covariance function which has to be well define the similarity between data points. Different definitions of the covarince function has been proposed before [START_REF] Niknazar | Nonparametric modelling of ECG: applications to denoising and to single sensor fetal ECG extraction[END_REF], [START_REF] Noorzadeh | Modeling quasi-periodic signals by a non-parametric model: Application on fetal ECG extraction[END_REF], which require the knowledge on R-peaks. Here we propose to define a covariance function to model the maternal ECG using a reference ECG which is recorded from the chest.

Consider the noisy abdominal ECG channel which is a mixture of maternal ECG, s m (t), fetal ECG, s f (t), and other noises like the electromyography of the mother's and/or the fetus' or the environmental noise , n(t):

s(t) = s m (t) + s f (t) + n(t). (2)
Another maternal reference channel, r(t), is also recorded from the maternal chest which mainly contains the maternal ECG, r m (t), and the noise, r n (t):

r(t) = r m (t) + r n (t). (3)
The two signals s m (t) and r m (t), which are attributed to the maternal ECG, are assumed to be correlated; however, they are decorrelated with the fetal ECG, s f (t), and with the noises n(t) and r n (t). Let's define y N (t) for any signal y(t) as

y N (t) = [y(t), • • • , y(t -N + 1)] T .
Now using the mentioned assumption the maternal ECG, s m (t), is modeled to depend on a window of the reference ECG, r N (t), by a Gaussian process:

s m (t) = f r N (t) + (t), (4)
with f r N (t) ∼ GP 0, k r N (t), r N (t ) , (5)
where we will then propose the following covariance function

k r N (t), r N (t ) = σ 2 exp -r N (t) -r N (t ) T r N (t) -r N (t ) 2l 2 , ( 6
)
where σ is used to model the amplitude of the signal while l is the length-scale and models the smoothness of he ECG. N is also the length of the window of the input of the covariance function.

The advantage of using such a covariance function is that there is no need for a prior knowledge on the beats and the quasi-periodicity of the ECG signal, but the points' relations of mECG are assumed to be well defined by the covariance function depending on a reliable reference signal. It has to be noted that the GP model does not assume a linear relationship between the input ( r N (t)) and the output (s m (t)) of the Gaussian process [START_REF] Pérez-Cruz | Gaussian processes for nonlinear signal processing: An overview of recent advances[END_REF] in equation ( 4). The rest of the abdominal signal can then be modeled by a GP using the covariance function

k n (t, t ) = σ 2 n δ(t -t ), (7)
where σ n is the amplitude and δ(.) is the delta Dirac function.

The proposed method can then extract the maternal ECG that is modeled with the GP:

ŝm (t * ) = k(t * )K -1 s. (8)
Here ŝm (t * ) is the maternal ECG which is estimated according to the method at t * points, and s is the vector of the abdominal observation defined as [s(t 1 ) • • • s(t τ )] T considering τ as the length of the signal. k is the maternal covariance which is defined as

k(t * ) = [k r N (t * ), r N (t 1 ) , • • • , k r N (t * ), r N (t τ ) ].
K is the covariance matrix of the abdominal signal whose i-jth element is defined as

K(i, j) = k r N (t i ), r N (t j ) + k n (t i , t j ).
Finally, the estimated maternal ECG can be omitted from the abdominal recording to leave an estimation of fetal ECG:

ŝf (t * ) = s(t * ) -ŝm (t * ). (9)
Apparently this estimation results in a noisy fetal ECG, since the abdominal signal contains not only the maternal ECG noise but also other kinds of noises (equation ( 2)); however, the main noise source, which is the maternal ECG, is eliminated; and the task of maternal and fetal ECG separation is accomplished. Now we propose to replace the reference signal (r(t)) by a 1-bit signal recorded from the maternal chest noted as r q (t). Since r and r q are highly correlated the assumption of the correlation between s m and r q is also true; according to this fact, the Gaussian process model for the maternal ECG can use the 1-bit signal as the input of the covariance function, which is less costly to be recorded using cheap sensors. Another remark that should be noted is that the the Gaussian process is a nonlinear solution to separate maternal and fetal signals [START_REF] Pérez-Cruz | Gaussian processes for nonlinear signal processing: An overview of recent advances[END_REF].

III. RESULTS

The proposed method is applied to the data we have recorded and the results are compared to the results obtained from the adaptive filter. The signals that are used here are the thoracic and abdominal signals recorded from a women in the 8 th month of pregnancy, using a 16-bit ADC resolution device with the sampling frequency of 400Hz. These recordings are shown in Fig. 1 which are the abdominal observation and thoracic ECGs noted as s(t) and r(t) respectively, in equations ( 2) and ( 3). It can be visually seen that the abdominal signal contains the maternal and fetal ECGs. Considering the large peaks as the R-peaks of maternal ECG, it is clear that the amplitude of the maternal signal is much higher than fetal ECG, whose lower amplitude R-peaks are visually seen in this data channel.

The Gaussian process model is applied to model the maternal ECG using both the 16-bit and the 1-bit digitized reference thoracic signals. Fig. 2 shows the extraction of maternal ECG from the abdominal channel which was depicted in Fig. 1, using the 16-bit thoracic ECG as a reference for the Gaussian process model, while the results of Fig. 3 are obtained by considering the 1-bit thoracic ECG as the reference. We have computed the 1-bit signal of r(t) as sgn(r). The figures show that the proposed GP method is capable of separating maternal and fetal ECGs. Furthermore, the use of 1-bit thoracic ECG as the reference does not worsen the performance. These results are then compared to the ones obtained by an adaptive filter designed by the LMS (Least Mean Square) algorithm that uses the same reference signals (full range and 1-bit) to separate maternal and fetal ECGs from the same abdominal signal. These results are shown in Fig. 4 and Fig. 5. While the efficiency of the Gaussian process method is almost the same for the 1-bit and the 16-bit reference signals, it is clear that adaptive filter cannot separate the fetal ECG using the 1-bit reference and this is due to the fact that adaptive filter is a linear filter, however the non-parametric model does a nonlinear separation. The performance of the adaptive filter using the full range reference is almost equivalent to that of GP; However, the extraction using a 1-bit reference is only possible with the GP method.

IV. CONCLUSIONS

This paper has investigated the problem of fetal ECG extraction by considering a thoracic reference signal for the maternal ECG. The maternal ECG is modeled by its second order statistics and is subtracted from the noisy abdominal channel to obtain the fetal ECG. The covariance function which is used can well define a model for maternal ECG using the thoracic signal as the input of the covariance function. The use of another reference channel would let to define a covariance function without considering any previous knowledge on the shape or quasi-periodicity of the ECG. This can be an advantage in the presence of arrhythmia, since the model does not need to define any prior information on the beats or the shape of the signal. On the other hand, non-parametric model does not consider a linear relationship between the input and the output of GP. This study has also proposed to use a 1-bit thoracic signal instead of a full range one, and the results have verified that using a 1-bit reference the method is also efficient. Although adaptive filtering method can also extract the fetal ECG from the abdominal channel using a full range thoracic ECG, this method does not work considering a 1-bit reference. The adaptive filter is a linear filter and the 1-bit reference will not provide sufficient information for extraction of maternal ECG from the abdominal channel. The advantage of using a 1-bit signal is that it is less costly to record by using cheap sensors or 1-bit ADCs, and is also less memory consuming and can be processed faster than the high-costly full range signals which are normally recorded with more that 10-bit ADCs. Fig. 2: Extraction of fetal ECG using a 16-bit reference with Gasussian process: r(t) is the thoracic reference, ŝm (t) is the maternal estimation, and ŝf (t) is the fetal estimation. Fig. 3: Extraction of fetal ECG using a 1-bit reference with Gasussian process: r(t) is the thoracic reference, ŝm (t) is the maternal estimation, and ŝf (t) is the fetal estimation. Fig. 4: Extraction of fetal ECG using a 16-bit reference with the adaptive filter: r(t) is the thoracic reference, ŝm (t) is the maternal estimation, and ŝf (t) is the fetal estimation. Fig. 5: Extraction of fetal ECG using a 1-bit reference with the adaptive filter: r(t) is the thoracic reference, ŝm (t) is the maternal estimation, and ŝf (t) is the fetal estimation.

V. ACKNOWLEDGMENT

Fig. 1 :

 1 Fig. 1: The abdominal noisy observation channel (on top), and the 16-bit thoracic ECG channel (on bottom)

This work has been supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025) and ERC AdG-2012-320864-CHESS.