Guillaume Bonfante

Matthieu Kaczmarek

Jean-Yves Marion

Architecture of a Morphological Malware Detector

à 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.

Introduction

The identification of malicious behavior is a difficult task. Until now, no technologies have been able to automatically prevent the spread of malware. Several approaches have been considered but neither syntactic analysis nor behavioral consideration were really effective. Presently, human analysis of malware seems to be the best strategy, next malware detectors based on string signature remains the most reliable solution. From this point of view, we have tried to easier the task which consist in finding a good signature within a malicious programs. Our technique has been inspired from the article [START_REF] Bruschi | Detecting self-mutating malware using control-flow graph matching[END_REF] where control flow graphs are used to detect the different instances of the computer virus MetaPHOR.

Generally speaking, detection strategies based on string signatures uses a database of regular expressions and a string matching engine to scan files and to detect infected ones. Each regular expression of the database is designed to identify a known malicious program. There are at least three difficulties tied to this approach. First, the identification of a malware signature requires a human expert and the time to forge a reliable signature is long compared to the time related to a malware attack. Second, string signature approach can be easily bypassed by obfuscation methods. Among recent work treating this subject, we propose to see for example [START_REF] Beaucamps | On the possibility of practically obfuscating programs towards a unified perspective of code protection[END_REF][START_REF] Christodorescu | Testing malware detectors[END_REF][START_REF] Filiol | Malware pattern scanning schemes secure against black-box analysis[END_REF]. Third, as the quantity of malware increase, the ratio of false positives becomes a crucial issue. And removing old malware signatures would open doors for outbreaks of reengineered malware. Thus, a current trend in the community proposes to design next generation of malware detectors over semantic aspects. [START_REF] Dalla Preda | A Semantics-Based Approach to Malware Detection[END_REF][START_REF] Christodorescu | Semantics-aware malware detection[END_REF][START_REF] Walenstein | Normalizing metamorphic malware using term rewriting[END_REF]. However, most of semantic properties are difficult to decide and even heuristics can be very complex as it is illustrated in the field of computer safety. For those reasons, in [START_REF] Bonfante | Control Flow Graphs as Malware Signatures[END_REF] we try to propose and to construct a morphological analysis in order to detect malware. The idea is to recognize the shape of a malicious program. That is, unlike string signature detection, we are not only considering a program as a flat text, but rather as a semantics object, so adding in some sense a dimension to the analysis. Our approach tries to combine several features: (a) to associate syntactic and semantic analysis, (b) to be efficient and (c) to be as automatic as possible.

Our morphological detector uses a set of CFG which plays the role of a malware signature database. Next, the detection consists in scanning files in order to recognize the shape of a malware. This design is closed to a string signature based detector and so we think that both approaches may be combined in a near future. Moreover, it is important to notice that this framework make the signature extraction easier. Indeed, either the extraction is fully automatic when the malware CFG is relevant or the task of signature makers is facilitated since they can work on an abstract representation of malicious programs.

This detection strategy is close to the ones presented in [START_REF] Christodorescu | Semantics-aware malware detection[END_REF][START_REF] Bruschi | Detecting self-mutating malware using control-flow graph matching[END_REF] but we put our strengths to optimize the efficiency of algorithms. For that sake, we use tree automata, a generalization to trees of finite state automata over strings [START_REF] Comon | Tree automata techniques and applications[END_REF]. Intuitively, we transform CFG into trees with pointers in order to represent back edges and cross edges. Then, the collection of malware signatures is a finite set of trees and so a regular tree language. Thanks to the construction of Myhill-Nerode, the minimal automaton gives us a compact and efficient database. Notice that the construction of the database is iterative and it is easy to add the CFG of a newly discovered malicious program.

Another issue of malware detections is the soundness with respect to classic mutation techniques. Here, we detect isomorphic CFG and so several comon obfuscation methods are canonically removed. Moreover, we add a rewriting engine which normalizes CFG in order to have a robust representation of the control flow with respect to mutations. Related works are [START_REF] Bruschi | Detecting self-mutating malware using control-flow graph matching[END_REF][START_REF] Christodorescu | Software transformations to improve malware detection[END_REF][START_REF] Walenstein | Normalizing metamorphic malware using term rewriting[END_REF] where data flow of programs is also considered.

The design of this complete chain of process is summarized by Figure 1.

We also provide large scale experiments, with a collection of 10156 malicious programs and 2653 sane programs. Those results are promising, with a completely automatic method for the signature extraction we have obtained a false positive ratio of 0.1%.

This study is organized as follows. First we expose the principles of CFG extraction and normalization. Then, we present a matching engine for CFG that is based on tree automata. Finally we

CFG in x86 languages

Road-map. Since we focus on practical aspects we choose to work on a concrete assembly language. This language is close to the x86 assembly language. We detail how to extract CFG from programs, we underline the difficulties that can be encountered and we outline how they can be overcome with classic methods. Finally, we study the problem of CFG mutations. We propose to normalize the extracted CFG according to rewriting rules in order to remove common mutations.

An x86 assembly language. We present the grammar of the studied programming language. The computation domain is the integers and we use a subset of the commands of the x86 assembly language. The important feature is that we consider the same flow instructions as in x86 achitectures, as a result the method that we develop can be directly applied to concrete programs.

Addresses N Offsets Z Registers R Expressions E ::= Z | N | R | [N] | [R] Flow instructions I f ::= jmp E | call E | ret | jcc Z Sequential instructions I d ::= mov E E | comp E E | . . . Programs P ::= I d | I f | P ; P
Next, a program is a sequence of instructions p = i 0 ; . . . ; i n-1 . The address of the instruction i k is k. In order to ease the reading and without loss of generality, we suppose that i 0 is the first instruction to be executed, the address 0 is the so called entry point of the program. We observe that the control flow of programs is driven by only four kinds of flow instructions. Given an instruction i k ∈ I f , the possible transfers of control are the following.

• If i k is an unconditional jump jmp e. The control is transferred to the address given by the value of the expression e.

• If i k is a conditional jump jcc x. If its associated condition is true, the control is transferred to the address k + x. Otherwise, the control is transferred to the address k + 1.

• If i k is a function call call e. The address k + 1 is pushed on the stack and the control is transferred to the the value of the expression e.

• If i k is a function return ret. An address is popped from the stack and the control is transferred to this address.

Prerequisites. The extraction of the CFG from a program is tied to several difficulties. First, we need access to the instructions of the program. As a result packing and encryption techniques can thwart the extraction. This problem is part of the folklore, indeed classical string signature detectors also have to face those techniques. Many solutions such as sand-boxes and generic unpackers have been developed to overcome this difficulty. The presentation of those solutions exceeds the scope of the current study then we refer to the textbooks [START_REF] Filiol | Advanced viral techniques: mathematical and algorithmic aspects[END_REF][START_REF] Filiol | Computer Viruses: from Theory to Applications[END_REF][START_REF] Ször | The Art of Computer Virus Research and Defense[END_REF]. Second, we are confronted to obscure sequenceq of instructions such as push a; ret which have the same behavior as the instruction jmp a. This is also part of the folklore and we will suppose that such sequence of instructions are normalized during the disassembly phase of the extraction.

Third, the target addresses of jumps and function calls have to be dynamically computed. For example, when we encounter the instruction jmp eax we need the value of the register eax in order to follow the control flow. In such cases we rely on an heuristic (|e| ) which provides the value of the expression e if it can be statically computed. If the value cannot be computed then (|e| ) = ⊥. Such an heuristic can be based on partial evaluation, emulation or any other static analysis technique.

The extraction procedure. The control flow consists in the different paths that might be traversed through the program during its execution. It is frequently represented by a graph named a control flow graph (CFG). The vertices stand for addresses of instructions and the edges represent the possible paths that the control flow can follow.

We suppose that we have access to the code of programs and that we have an heuristic (| | ) to evaluate expressions. Table 1 presents a procedure to extract CFG from programs. We observe that this procedure closely follows the semantics of flow instructions. Indeed, the vertices of the CFG are labeled accordingly to the instruction at the night address and the nodes are linked according to the possible control transfers.

• The symbol inst of arity 1 labels addresses of sequential instructions. There is only one successor: the address of the next instruction.

• The symbol jmp of arity 1 labels addresses of unconditional jumps. There is only one successor: the address to jump to.

• The symbol jcc of arity 2 labels addresses of conditional jumps. There is two successors: the address to jump to when the condition is true and the address of the next instruction where the control is transfered when the condition is false.

• The symbol call of arity 2 labels addresses of function calls. There is two successors: the address of the function to call and the return address, that is the address of the next instruction.

• The symbol end of arity 0 labels addresses of function returns and undefined instructions.

There is no successor.

The entry point of the program correspond to the root of the CFG.

Instruction Graph

i n ∈ I d i n = jmp e (|e| ) = k i n = call e (|e| ) = k i n = jcc x

Otherwise Table 1: Control flow graph extraction

Normalizing mutations. Our CFG representation is a rough abstraction of programs. Indeed we do not make any distinction between the different kinds of sequential instruction, all of them are represented by nodes labeled with the symbol inst. This first abstraction level makes the CFG sound with respect to mutations which substitutes instructions with the same behavior. For example the replacement of the instruction mov eax 0 by the instruction xor eax eax does not impact our CFG representation.

However, the soundness with respect to classic mutations techniques remains an important issue. Indeed, some well know mutation techniques can alter the CFG of malicious programs. In order to recover a sound representation of the control flow we apply reductions on CFG. A reduction is defined by a graph rewriting rule. As a case study, we consider three reductions associated to classic mutation techniques. Of course several other reductions can be defined in order to handle more mutations techniques. We use the following reductions

• Concatenate consecutive instructions into blocks of instructions.

• Realign code removing superfluous unconditional jumps.

• Merge consecutive conditional jumps.

Those abstractions can be defined through the graph rewriting rules of Table 2. Figure 2 presents an assembly program and its reduced CFG. We remark that each rewriting rule impose a diminution of the size of the rewritten graph then the reduction clearly terminates. Moreover, since there is no critical pair we have no problem of confluence. Nevertheless, normalizing mutations through rewriting rules is a generic principle that could be applied on sophisticated cases. Then, the issues of termination and confluence shall be carefully considered.

Table 3

Efficient database

Road-map. Morphological detection is based on a set of malware CFG which plays the role of malware signatures. This collection of CFG is compiled into a tree automaton thanks to a term representation. Since tree automata fulfill a minimization property, we obtain an efficient representation of the database. Next, we apply this framework for the sub-CFG isomorphism problem in order to detect malware infections.

From graphs to terms. A path is a word over {1, 2} * , we write ǫ the empty path. We define the

ρ1 < ρ2 ρ < ρi ρ < τ ⇒ ρρ ′ < τ τ ′
A tree domain is a set d ⊂ {1, 2} * such that for any path ρ ∈ {1, 2} * and any integer i ∈ {1, 2} we have

ρi ∈ d ⇒ ρ ∈ d
A tree over a set of symbols F is a pair t = (d(t), t) where d(t) is a tree domain and t is a function from d(t) to F.

We consider the set of symbols F = {inst, jmp, call, jcc, ret} ∪ {1, 2} * and the trees overs this set. In such trees, a nodes labeled by a path ρ = {1, 2} * is thought as pointers to the night node of the tree. Then, a tree have two kinds of nodes: the inner nodes labeled by symbols of {inst, jmp, call, jcc, ret} and the pointer nodes labeled by path in {1, 2} ρ . In the following we write d(t) the set of inner nodes of the tree t, that is

d(t) = ρ ρ ∈ d(t) t(ρ) ∈ {inst, jmp, call, jcc, ret}
Next a tree t is well formed if for any paths ρ, τ ∈ d(t)

t(ρ) = τ ⇒ τ ∈ d(t) and ρ ≤ τ
We observe that any CFG can be represented by a unique well formed tree.

Tree automata. A finite tree automaton is a tuple A = (Q, F, Q f , ∆), where Q is a finite set of states, F is a set of symbols, Q f ⊂ Q is a set of final states and ∆ is a finite set of transition rules of the type a(q 1 . . . q i ) → q with a ∈ F has arity i and q, q 1 , . . . , q i ∈ Q.

A run of an automaton on a tree t starts at the leaves and moves upward, associating a state with each sub-tree. Any symbol a of arity 0 is labeled by q if a → q ∈ ∆. Next, if the direct sub-trees t 1 , . . . , t n of a tree t = a(t 1 , . . . , t n ) are respectively labeled by states q 1 , . . . , q n then the tree t is labeled by the state q if a(q 1 , . . . , q n ) → q ∈ ∆. A tree t is accepted by the automaton if the run labels t with a final state. We observe that a run on a tree t can be computed in linear time, that is O(n) where n is the size of t, that is the number of its nodes.

For any automaton A, we write L(A) the set of trees accepted by A. A language of trees L is recognizable if there is a tree automaton A such that L = L(A). We define the size |A| of an automaton A as the number of its rules.

Tree automata have interesting properties. First, it is easy to build an automaton which recognize a given finite set of trees. This operation can be done in linear time, that is O(n) where n is the sum of the sizes of the trees in the language. Second, we can add new trees to the language recognized by an automaton computing a union of automata, see [START_REF] Comon | Tree automata techniques and applications[END_REF]. Given an automaton A, the union of A with an automaton A ′ can be computed in linear time, that is O(|A ′ |).

Finally, for a given recognizable tree language, there exists a unique minimal automaton in the number of states which recognizes this language. This property ensures that the minimal automaton is the best representation of the tree language.

Theorem 1 (From [START_REF] Comon | Tree automata techniques and applications[END_REF]). For any tree automaton A which recognize a tree language L we can compute in quadratic time (O(|A| 2 )) a tree automaton A which is the minimum tree automaton recognizing L up to a renaming of the states.

Building the database. We explain how this framework can be used to detect malware infections. Suppose that we have a set {t 1 , . . . , t n } of malware CFG represented by trees. Since this set is finite, there is a tree automaton A which recognizes it.

Next, consider the tree representation t of a given program. Computing a run of A on t, we can decide in linear time if this tree is one of the the trees obtained from malware CFG. This means that that we can efficiently decide if a program has the same CFG as a known malware.

Finally, we can speed-up the detection computing the minimal automaton which recognize the language {t 1 , . . . , t n }. From a practical point of view this is the most efficient representation of the malware CFG database.

Detecting infections. Actually, when a malicious program infects an other program, it includes its own code within the program of its host. Then, we can reasonably suppose that the CFG of the malicious program appears as a subgraph of the global CFG of the infected program. As a result, we can detect such an infection by deciding the subgraph isomorphism problem within the context of CFG.

First we have to observe that we are not confronted with the general sub-graph isomorphism since CFG are graphs with strong constraints. In particular the edge labeling property implies that a CFG composed of n nodes accepts at most n subgraphs. As a result, the sub-CFG isomorphism problem is not NP-complete. Then to detect sub-CFG it is sufficient to run the automaton on the tree representations of any sub-CFG.

Experiments

Road-map. We consider the win32 binaries of VX Heavens malware collection [2]. This collection is composed of 10156 malicious programs. Then, we have collected 2653 win32 binaries from a fresh installation of Windows Vista TM . This second collection is considered as sane programs.

Using those samples we experiments with our implementation of the morphological detector. We focus our attention on false positive ratios in order to validate the our method. Indeed, we have to know if it is possible to discriminate sane programs from malicious ones only considering their CFG. The following experimental results agree with this hypothesis.

CFG extraction in practice. To overcome the difficulties of CFG extraction we have chosen the following solutions

• In order to deal with crypted and packed samples, we use the unpacking capabilities of ClamAV TM [3].

• We have implemented a dynamic disassembler based on the disassembler library Udis86 [1].

Our module is able to follow the control flow and it keeps track of the stack in order to remove push a; ret sequences.

• The evaluation heuristic (|e| ) proceed as follows.

When we encounter a dynamic flow instruction, we emulate the preceding block of sequential instructions in order to recover the value of the expression e. Our emulation technology is also build over Udis86. It is limited to a subset of x86 assembly instruction, interruptions and system calls are not taken into account.

• We reduce the obtained CFG according to the rules of Table 2.

Figure 3 gives the sizes of the reduced CFG extracted from the programs of those collections. On the X axis we have the upper bound on the size of CFG and on the Y axis we have the percentage of CFG whose size is lower than the bound. We obeserve that about 5% of the database are programs with a non valid PE header, they produce an empty graph. Then we are able to extract a CFG of more that 15 nodes from about 65% of the samples. The remaining 30% produces a CFG which have between 1 and 15. We think that those graphs are too small to be relevant. We are curruntly working on this part of the samples to improve our extraction procedure.

Building the database. The size of malware control flow graphs clearly impact the accuracy of the control flow detector. We have observed that the graphs extracted from some malware were too small to be relevant and the resulting detector makes many false alerts because of a few such graphs. As a result, we impose a lower bound on the size of the graphs that we include in the database. Next, we have done several tests using different lower bounds.

Let N ∈ N be the lower bound on the size of CFG. We build the minimized automaton A N M which recognizes the set of tree representations of reduced malware CFG that are composed of more than N nodes. We define the morphological detector D N M as a predicate such that for any program p ∈ P we have D N M (p) = 1 if a malware CFG appears as a subgraph of the CFG of p and D N M (p) = 0 otherwise. We have seen in the previous sections that D N M can be decided using A N M . This design has several advantages. First, when a new malicious program is discovered, one can easily add its CFG to the database using the union of tree automata and a new compilation to obtain a minimal tree automaton.

The computation of the 'not minimal' automata takes about 25 minutes. The minimization takes several hours but this delay is not so important. Indeed, within the context of an update of the malware database, during the minimization we can release the 'not minimal' automaton. Indeed, even if this is not the best automaton it still recognize the malware database and it could be used until the minimization is terminated.

Evaluation. We are interested in false positives, that is sane programs detected as malicious. For that, we have collected 2653 programs from a fresh installation of Windows Vista We do not evaluate false negatives, that is undetected malicious programs. Indeed, by construction all malicious programs of our malware collection are detected by the morphological detector. Nevertheless, this methods seems promising for this aspect. Indeed, the study [START_REF] Bruschi | Detecting self-mutating malware using control-flow graph matching[END_REF] has shown that a CFG based detection allows to detect the high-obfuscating computer virus MetaPHOR with no false negative.

However, the presence of the lower bound on the size of malware CFG that enter in the database implies that some malware are undetected. Those can be considered as false negatives, even if they more related to the technical problem of CFG extraction than to the methods of morphological detection. The results of the experiments mention this ratio of undetected malware.

Experimental results. We have built tree automata from the malware samples considering different lower bounds N on the size of CFG. According to the previous section we obtain the morphological detectors D N M . We have tested those detectors on the collection of saneware in order to evaluate the false positives. It takes about 5 h 30 min to analyze the collection of saneware, this represents the analysis of 2 ′ 319 ′ 294 sub-CFG. Table 4 presents the results. The first column indicates the considered the lower bound N . The second column indicates the ratio of false positives that is the number of sane programs detected as malicious out of the toalt number of sane programs. The third column indicates the ratio of undetected malicious programs that is the number of malware samples with a CFG size lower than the bound out of the total number of malware samples.

Analysis.

As expected, we observe that the false positives decrease with the lower bound on the size of CFG. Over 15 nodes, the CFG seems to be a relevant criterium to discriminate malware.

Concerning the remaining false positives. The libraries ir41 qc.dll and ir41 qcx.dll, and the malicious program Trojan.Win32.Sechole have the same CFG composed of more than 1 ′ 000 nodes. We have tested those programs with commercial antivirus software and the libraries ir41 qc.dll and ir41 qcx.dll are not detected whereas the program Trojan.Win32.Sechole is detected as malicious. The malicious programs seems to be based on the dynamic library and the extraction algorithm was not able to extract the CFG related to the malicious program.

Concerning the ratio of undetected malware, The only way to improve the detector is to implement a better heuristic for control flow graph extraction. In its current version our prototype only use a few heuristics.

For comparison, statistical methods used in [START_REF] Kephart | Automatic Extraction of Computer Virus Signatures[END_REF] induce false negatives ratios between 36 % and 48 % and false positive ratios between 0.5 % and 34 %.

A detector based on artificial neural networks developed at IBM [START_REF] Gj Tesauro | Neural networks for computer virus recognition[END_REF] presents false negatives ratios between 15 % and 20 % and false positive ratios lower than 1 %. The data mining methods surveyed in [START_REF] Schultz | Data Mining Methods for Detection of New Malicious Executables[END_REF] present false negatives ratios between 2.3 % and 64.4 % and false positive ratios between 2.2 % and 47.5 %. Heuristics methods from antivirus industry tested in [START_REF] Gryaznov | Scanners of the Year[END_REF] present false negatives ratios between 20.0 % and 48.6 % and false positive ratios lower than 0.2 %.

Figure 1 :

 1 Figure 1: Design of the control flow detector

0:

  Figure 2: A program and its CFG

  presents mutations of the program of Figure 2. All of them have the same reduced CFG as the original program.

  0: cmp eax 0 1: jne +7 2: mov ecx eax 3: sub ecx 1 4: mul eax ecx 5: cmp ecx 1 6: jne -3 7: jmp +2 8: mov eax 1 9: ret 0: cmp eax 0 1: jne +8 2: push eax 3: pop ecx 4: dec ecx 5: mul eax ecx 6: cmp ecx 1 7: jne -3 8: jmp +2 9: inc ecx 10: ret 0: cmp eax 0 1: jne +7 2: mov ecx eax 3: dec ecx 4: mul eax ecx 5: cmp ecx 1 6: jne -3 9: ret 8: inc ecx 9: jmp -2 0: cmp eax 0 1: jne +9 2: mov ecx eax 3: dec ecx 4: mul eax ecx 5: cmp ecx 2 6: ja -3 7: cmp ecx 1 8: jne -5 9: jmp +2 10: inc ecx 11: ret 0: cmp eax 0 1: je +2 2: jmp +10 2: push eax 3: pop ecx 4: sub ecx 1 5: mul eax ecx 6: cmp ecx 2 7: ja -3 8: cmp ecx 1 9: jne -5 10: ret 11: mov ecx 1 12: jmp -2

Figure 3 :

 3 Figure 3: Sizes of control flow graphs

  TM . Let us note S this set of programs. Let N ∈ N be a lower bound on the size of malware CFG, we consider the following approximation of the false positives of the detector D N M False positives {p | D N M (p) = 1 and p ∈ S}

Table 3 :

 3 Control flow graph mutations path order for any path ρ, τ ∈ {1, 2} * and any integer i ∈ {1, 2} as follows.

Table 4 :

 4 Results of the experiments


domain: spi phys math info

text:
Helene Piet-Lahanier

Nonlinear Estimation of Sensor Faults With Unknown Dynamics for a Fixed Wing Unmanned Aerial Vehicle

Enzo Iglesis 1 Nadjim Horri 2 Karim Dahia 3 James Brusey 4 and Helene Piet-Lahanier 5

Abstract-In this paper, the estimation of additive inertial navigation sensor faults with unknown dynamics is considered with application to the longitudinal navigation and control of a fixed wing unmanned aerial vehicle. The faulty measurement is on the pitch angle. A jump Markov regularized particle filter is proposed for fault and state estimation of the nonlinear aircraft dynamics, with a Markovian jump strategy to manage the probabilistic transitions between the fault free and faulty modes. The jump strategy uses a small number of sentinel particles to continue testing the alternate hypothesis under both fault free and faulty modes. The proposed filter is shown to outperform the regularized particle filter for this application in terms of fault estimation accuracy and convergence time for scenarios involving both abrupt and incipient faults, without prior knowledge of the fault models. The state estimation is also more accurate and robust to faults using the proposed approach. The root-mean-square error for the altitude is reduced by 77 % using the jump Markov regularized particle filter under a pitch sensor fault amplitude of up to 10 degrees. Performance enhancement compared to the regularized particle filter was found to be more pronounced when fault amplitudes increase.

I. INTRODUCTION

Sensor fault detection and diagnosis systems are increasingly important to aircraft mission integrity. Inertial navigation malfunctions in particular have often been the cause of flight incidents and crashes as in the case of the Qantas F72 and Croatia Boeing 737-200 flights [START_REF] Wang | Fault detection and exclusion for tightly coupled gnss/ins system considering fault in state prediction[END_REF]. Small unmanned aerial vehicles are also becoming increasingly autonomous, with a need to further develop their ability to detect, estimate and recover from sensor faults. The certification of autonomous unmanned aerial vehicle (UAV) systems is also subject to the development of systems to monitor sensor data and fault flags.

Fault detection often consists of statistical tests used to detect a change on residuals. This is either performed by comparison against a threshold or by assuming a residual distribution, such as the Student's t-test, the cumulative sum France. helene.piet-lahanier@onera.fr (CUSUM) or a likelihood ratio test such as the generalized likelihood ratio test (GLRT) [START_REF] Marzat | Modelbased fault diagnosis for aerospace systems: a survey[END_REF].

A higher level of fault diagnosis is fault estimation. Indeed, fault detection does not provide quantitative information on the impact of fault on the system. Hence, fault estimation is required to obtain valuable information that can be exploited e. g. to compensate the fault. For linear Gaussian models, the Kalman filter (KF) is known to be the optimal estimator and the approach has been extended to fault estimation. For nonlinear models, the extended Kalman filter (EKF) has been used for both fault detection [START_REF] Eykeren | Sensor fault detection and isolation using adaptive extended kalman filter[END_REF] and fault estimation [START_REF] Alcalay | An adaptive extended kalman filter for monitoring and estimating key aircraft flight parameters[END_REF]. When a variety of faults is considered, interacting multiple model (IMM) approaches are often used, and they rely on the use of several filters, all of which are associated to a specific fault mode. IMMs usually consist of interacting KFs or EKFs depending on the linearity of the system. They have been proven efficient for various applications as in [START_REF] Rago | Failure detection and identification and fault tolerant control using the immkf with applications to the eagle-eye uav[END_REF] but present the major drawback of requiring detailed knowledge of the models of the various faults.

Particle filters have proved their efficiency in estimation problems for non-linear dynamics and non-Gaussian noises [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking[END_REF]. In the context of fault estimation it has been combined with a jump Markov transition model that enables to switch from a nominal to faulty model using the transition probability matrix that accounts for abrupt changes of modes [START_REF] Doucet | Particle filters for state estimation of jump markov linear systems[END_REF], [START_REF] Tafazoli | Hybrid system state tracking and fault detection using particle filters[END_REF]. This approach is combined in this paper with a regularized particle filter (RPF) [START_REF] Musso | Improving Regularised Particle Filters[END_REF] and a Kalman correction step to handle unexpected fault dynamics. It is called a jump Markov regularized particle filter (JMRPF).

The dynamical model of fixed wing UAVs is well known [START_REF] Beard | Small unmanned aircraft: Theory and practice[END_REF], a model-based approach [START_REF] Ding | Model-based fault diagnosis techniques: design schemes, algorithms, and tools[END_REF]- [START_REF] Castillo | Robust model-based fault detection and isolation for nonlinear processes using sliding modes[END_REF] is therefore used in this paper. With no prior knowledge of the fault dynamics, a zero-order constant fault model is used by the JMRPF, which is initialized by default in a fault free mode. The actual incipient sensor faults applied to the system do not match this zero-order model and the process and regularization noises are not set to handle abrupt fault amplitudes, but the JMRPF is shown in Section VI to accurately estimate sensor faults and robustly estimate the UAV states despite this model mismatch. This is due to the fact that piecewise constant approximations of the faults can be rapidly tracked using the JMRPF that switches more swiftly when a mode change is detected.

The proposed approach is applied to a nonlinear model of longitudinal dynamics for a fixed wing UAV.

The main contributions of the paper are as follows:

• A jump strategy between the fault free and faulty sensor modes is proposed, where the a priori distribution of the fault is computed using sensor innovation terms. Under both fault free and faulty modes, a small set of sentinel particles is allowed to test the alternate mode, leading to fast mode transitions and higher fault and state estimation accuracy. • A Kalman correction is added to the JMRPF to place the particles in the most likely areas of the state space. This correction further improves state and fault estimation accuracy and was not used in reference [START_REF] Iglésis | A jump-markov regularized particle filter for the estimation of ambiguous sensor faults[END_REF].

The paper is organized as follows. In Section II, a nonlinear model of fixed wing UAV flight dynamics is presented, including the longitudinal UAV autopilot. Section III is focused on the sensor fault types and modes that the nonlinear filters have to estimate. In Section IV, a more general stochastic nonlinear Markovian jump model is presented including the effect of the fault modes Section V details the JMRPF approach for state and sensor fault estimation in a stochastic Markovian jump model framework. In Section VI, the fault and state estimation performance of the JMRPF is evaluated and compared against the RPF for scenarios including abrupt and incipient faults with no prior knowledge of the fault dynamics and amplitudes. Section VII concludes the paper.

II. LONGITUDINAL FLIGHT DYNAMICS AND CONTROL SYSTEM

A. Longitudinal UAV Dynamics

A nonlinear model of longitudinal fixed wing UAV dynamics is used in this paper. The state vector representing UAV longitudinal dynamics is denoted z = p d u w θ q . The state p d denotes the position component in the downward direction along k v (see Fig. 1), u represents the variation with respect to the trim condition on longitudinal velocity along the i b axis, w represents the vertical velocity along the k b axis (see Fig. 1) and the states θ and q respectively represent the pitch angle and pitch rate variations. The nonlinear longitudinal model of the Aerosonde UAV is

Earth center Horizon i v k v i b k b θ
Center of mass obtained from [START_REF] Beard | Small unmanned aircraft: Theory and practice[END_REF] and is given by:

                         ṗd = -sin (θ)u + cos (θ)w u = -qw + F x m ẇ = qu + F z m θ = q q = - J xz J y q 2 + 1 J y M (1a) (1b) (1c) (1d) (1e)
where F x and F z respectively represent the external forces along the i b and k b axes (see Fig. 1), m is the mass of the UAV, J xz is a product of inertia, J y is the moment of inertia about the pitch axis and M is the pitching moment. All these parameters are defined in [START_REF] Beard | Small unmanned aircraft: Theory and practice[END_REF].

The discrete time model derived from (1) can be expressed as:

z k+1 = F k (z k , u k ) (2)
where u = δ e δ t T is the control input vector and δ e and δ t respectively represent the elevator deflection and throttle input. The input δ e is bounded in the interval [-25°, 25°] and the input δ t is bounded in the interval [0, 1]. The full state is observed using an inertial navigation system (INS), which is hybridized with a global navigation satellite system (GNSS) receiver, a magnetometer and a barometer. The measurement equation is given in Section IV where a stochastic model including the effects of sensor faults is presented. The observation function is denoted H k and is given by:

H k (z k ) = -p dk u k w k θ k q k (3)

B. Longitudinal Autopilot

In the following, we assume that the desired trajectory to be followed by the UAV consists of a desired flight path angle γ c and a desired velocity norm V c that are used as reference inputs to the longitudinal control loop. The actuator inputs δ e and δ t are given by:

δek+1 = -L θ ẑk -L θi θik+1 δtk+1 = -L u ẑk -L ui ūik+1 (4a) (4b)
where the bar notation represents a variation around the trim condition 1 . The gains L θ , L u , L ui , L θi are obtained after solving the Riccatti equation of the linear quadratic regulator (LQR) problem for the desired steady trim condition with integral correction and weighting matrices Q = diag 1 0 4 0 0 , R = I 2×2 . The LQR approach is suboptimal for the nonlinear system, but the computed gains will be shown to achieve robust stability despite the nonlinearity of UAV dynamics, particularly when the initial condition does not lead to a severe coupling between rotational and translational dynamics. The integral gains are L θi = 1.00 and L ui = -1.00. Integrated state deviations θi and ūi are given by:

     θik+1 = (γ c k + A u ûk + A w ŵk -θ)dt + θik ūik+1 = (( V c k -V w ŵk ) 1 V u -ûk )dt + ūik (5a) (5b)
where the parameters V u = 1.00, V w = 0.05, A u = 0 and A w = 0.03 relate pitch and speed reference guidance commands to γ c and V c .

III. SENSOR FAULT TYPES AND MODES As in [START_REF] Isermann | Terminology in fault detection and diagnosis[END_REF], a fault is defined here as an "unpermitted deviation of at least one characteristic property of the system". In this paper, sensor outputs can be faulty and the fault states have two possible modes, a fault-free mode denoted m (0) and a faulty mode denoted m (1) .

For the UAV model under consideration, the number of sensors denoted n y is set to 5 without loss of generality but in the next section, the fault and state estimation algorithm presented here is written for a more general case where n y depends on the application. Fault estimation is implemented for n Θ sensors where n Θ ≤ n y is user-defined depending on the sensors deemed likely to be faulty. Each sensor has an index j ∈ [1, n y ], a fault variable Θ j (m j k ) and an associated fault mode m j k at time step k. A fault and state estimation algorithm is developed here for the purpose of fault tolerant navigation, with no prior knowledge of fault dynamics or amplitudes.

The fault types applied to the system are additive and of unknown amplitude. Their dynamics are also unknown as the approach tackles abrupt or incipient faults. Exponential fault models are for example considered in the numerical simulation section, but the estimation algorithms use a zeroorder fault model, which is initialized by default in mode m (0) . Abrupt faults occur suddenly (stepwise), while incipient faults occur gradually and driftwise [START_REF] Varga | Solving fault diagnosis problems[END_REF], [START_REF]Fault-diagnosis systems: an introduction from fault detection to fault tolerance[END_REF].

The estimation algorithms also have to determine when the fault is no longer active. Both abrupt and incipient faults scenarios under consideration in this paper are therefore intermittent, such that the fault is deactivated after a certain time [START_REF] Varga | Solving fault diagnosis problems[END_REF].

In the next section, a stochastic Markovian jump system model is used to represent the dependency of the nonlinear aircraft model on sensor faults.

IV. STOCHASTIC JUMP MARKOV NONLINEAR SYSTEM

MODEL A discrete-time stochastic Markovian jump system model is used to represent the dynamics of the system, including the transitions between fault modes. This generic system model is given by:

     m k+1 ∼ p (m k+1 |m k ) z k+1 = F k (z k , u k ) + η k y k = H k (z k ) + G k (Θ k (m k )) + ν k (6a) (6b) (6c)
where z k ∈ R nz represents the system state u k ∈ R nu is the control input, y k ∈ R ny is the vector of measurements.

F k (•) represents the dynamics of the system (see [START_REF] Marzat | Modelbased fault diagnosis for aerospace systems: a survey[END_REF]), H k (•) represents the measurement function (see [START_REF] Eykeren | Sensor fault detection and isolation using adaptive extended kalman filter[END_REF]) and G k (•) represents the impacts of the sensor faults on the measurements. The process and sensors noises are η k ∈ R nz and ν k ∈ R ny . They are assumed to be of zero mean and the covariance matrices are respectively defined as

E η k η k = Q k and E ν k ν k = R k . They are assumed to be independent E η k ν k = 0.
A sensor fault vector Θ k ∈ R nΘ at time step k is defined. It is a function of the fault mode vector m k , for which each element m j k where j is a sensor index can either be faulty m (1) or fault free m (0) , at time step k. The number of sensors n Θ ≤ n y for which fault diagnosis is applied is userdefined, but taken to be equal to one in Section VI, where the method is applied to a single sensor fault on the pitch measurement. To simplify notations, Θ k = Θ k (m k ) is used for the remainder of this article.

In this Markovian jump system, the mode switching between fault free and faulty modes is managed by a jump strategy, which is presented in Subsection IV-A

A. The Transition Probability Based Jump Strategy

At every time step k, the probability P m

(j) k+1 |m (i) k to switch from mode m (i) to m (j) is constant π ji .
Hence, π 10 is the probability to switch from nominal mode m (0) to a faulty mode m (1) while the probability π 01 is the probability to switch from a faulty mode m (1) to a nominal mode m (0) . The transition probability matrix Π represents the probability of switching from one mode to another. It is given by: Π = π 00 π 10 π 01 π 11 [START_REF] Doucet | Particle filters for state estimation of jump markov linear systems[END_REF] Each state of the fault vector Θ is associated with a Π matrix. Note that for the numerical analysis in Section VI, a single scalar fault is considered without loss of generality and to focus the analysis on estimation performance for abrupt and incipient fault types, but the proposed approach applies to multiple observable sensor faults. The diagonal elements of the π jj matrices represent probabilities to remain in the same mode for the given sensor.

The Markov chains can be represented by the transitions diagram shown in Fig. 2:

m (0)
m (1) π 11 π 00 π 10 π 01

Fig. 2. Markov chain of the JMRPF applied to fault estimation

The objective of the method presented here is to simultaneously detect the occurrence of a fault and to estimate its amplitude. An extended state vector x k ∈ R nx is introduced and given by:

x k = z k Θ k (8)
which extends the original system state z k to account for the fault values. The system model ( 6) can then be written as an extended state model:

x k+1 = f k (x k , u k ) + η k y k = h k (x k ) + ν k (9a) (9b)
where η k ∈ R nx now represents the process noise of the extended state space model. The dynamics of the system and measurements functions f k and h k are given by:

f k (x k , u k ) = F k (z k , u k ) T k (Θ k ) (10a) h k (x k ) = H k (z k ) + G k (Θ k ) (10b)
where T k (•) represents the dynamics of the sensor fault.

For this representation of the problem under consideration, estimation of x k requires nonlinear filtering techniques. As introduced previously, jump Markov based filters such as the JMRPF are well-suited to handle such systems.

V. JUMP MARKOV REGULARIZED PARTICLE FILTER

The proposed JMRPF combines the RPF [START_REF] Musso | Improving Regularised Particle Filters[END_REF] with a new jump Markov strategy to manage transitions between faulty and fault free modes using a transition probability matrix. As in a RPF, particles are resampled from a continuous approximation to the posterior density to reduce degeneracy without compromising the diversity of solutions. A first version of a JMRPF was introduced in [START_REF] Iglésis | A jump-markov regularized particle filter for the estimation of ambiguous sensor faults[END_REF] to recover from ambiguous abrupt sensor fault scenarios. In this paper, a Kalman correction is introduced to correct particles states by placing them in more likely state space regions. The Kalman innovation covariance matrix S k is taken into consideration in the calculation of the particles weights. This reduces the variance of the weights. The estimated state vector x k and the estimated covariance matrix Pk are obtained by the JMRPF algorithm using estimated state feedback, previous estimated state and the measurement as inputs to the filter. The total number of particles is denoted N p .

The proposed JMRPF algorithm is introduced in Algorithm 1. The algorithm has prediction, update, estimation and regularization-resampling steps. The PREDICT and UPDATE functions of Algorithm 1 are improved compared to [START_REF] Iglésis | A jump-markov regularized particle filter for the estimation of ambiguous sensor faults[END_REF] to deal with sensor fault estimation and a JUMP function is introduced within the prediction step.

A. Prediction Step

The prediction step is described in Algorithm 2. The i th state variable is propagated using the following probability transition density for the state x k :

x i k|k-1 ∼ p x k|k-1 |x i k-1 , m i k (11)
Then, one obtains a predicted cloud of particles (x 1 k|k-1 ,

x 2 k|k-1 , • • • , x Np k|k-1
). The mode of each state and particles of Θ k are updated in the prediction step. This update is performed here using a uniform law and compared to userdefined probabilities π ji to switch from one mode to another as described in Fig. 2.

Algorithm 1 jump Markov regularized particle filter k ← 0 . . . Initialization loop k ← k + 1 for each i ∈ [1, N p ] do PREDICT(x i k|k-1 , x i k-1 , m i k , u k , y k ) end for ESTIMATE(x k|k-1 , Pk|k-1 , w 1:Np k-1 , x 1:Np k|k-1 ) ŷk|k-1 ← Np i=1 w i k-1 h k x i k|k-1 S k ← Np i=1 w i k-1 y i k|k-1 -ŷk|k-1 y i k|k-1 -ŷk|k-1 + R k PXY ← Np i=1 w i k-1 x i k|k-1 -xk|k-1 y i k|k-1 -ŷk|k-1 K k ← PXY S -1 k Kalman gain for each i ∈ [1, N p ] do UPDATE(x i k , w i k , w i k-1 , x i k|k-1 , K k , S k , y k ) end for ESTIMATE(x k , Pk , w 1:Np k , x 1:Np k ) N eff ← 1 Np i=1 w i k 2 if N eff ≤ N p Γ then if true then resample MULTINOMIAL(x 1:Np k , x 1:Np k , w 1:Np k ) for each i ∈ [1, N p ] do w i k ← 1 Np Reset the weights REGULARIZE(x i k , xi k ) end for end if end loop
To simplify notations in this section, it is assumed that the maximum number of possible sensor faults n Θ is equal to the number of sensors n y , although this number can easily be reduced in practice to restrict fault management to specific sensors. The same index j is also used to denote the j th state of vector Θ k , representing fault estimates of the j th sensor of y k . The jump amplitude of the i th particle of the j th state of Θ at time step k is computed as follows:

Θ i,j k|k-1 =      β i,j k if U ≤ π j 10 and m i,j k = m (0) Θ i,j k|k-1 if U < π j 11 and m i,j k = m (1) 0 else ( 12
)
where U ∼ U (0, 1) and β i k is given by:

β i k = y k -h k x i k|k-1 (13)
A new jump strategy for sensors fault modes is proposed. It uses a priori distribution of the fault, which is computed using sensor innovation terms from [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF]. Irrespective of the current estimated mode, the alternate mode will continue to be tested using a small subset of fault state particles that will be called sentinel particles. Those particles are selected from a uniform distribution and their number are smaller Algorithm 2 Detail of the function PREDICT from Algorithm 1

function PREDICT(x i k|k-1 , x i k-1 , m i k , u k , y k ) η i k ∼ N (0, Q k ) x i k|k-1 ← f k x i k-1 , u k + η i k β i k ← y k -h k x i k|k-1
See ( 13)

for each j ∈ [1, n Θ ] do Jump step of Θ JUMP(Θ i,j k|k-1 , m i,j k , β i,j k )
See Algorithm 3 end for end function when the off-diagonal terms π 01 and π 10 of the transition probability matrix are smaller. Those terms are set based on empirical sensor false alarm and missed detection rates. This small number of sentinel particles will continuously test the probability of transition to the alternate fault modes. The use of a small number of sentinel particles to test the alternate mode enhances real time operation prospects compared to previously published particle filter jump strategies, where the total or an arbitrary number of particles is used to evaluate the probabilities of both modes before any jump between them [START_REF] Tafazoli | Hybrid system state tracking and fault detection using particle filters[END_REF].

The JUMP function used in Algorithm 2 is described in Algorithm 3.

Algorithm 3 Detail of the function

JUMP from Algorithm 2 function JUMP(Θ i,j k|k-1 , m i,j k , β i,j k ) U ∼ U (0, 1) if m i,j k = m (0) then Θ i,j k|k-1 in mode m (0) if U ≤ π 10 then Transition m (0) → m (1) Θ i,j k|k-1 ← β i,j k m i,j k ← m (1) else Transition m (0) → m (0) Θ i,j k|k-1 ← 0 end if else if m i,j k = m (1) then Θ i,j k|k-1 in mode m (1) if U ≤ π 01 then Transition m (1) → m (0) Θ i,j k|k-1 ← 0 m i,j k ← m (0) end if end if end function

B. Update Step

The update step is described in Algorithm 4. The i th particle x i k is assigned to a weight w i k that is proportional to its likelihood:

wi k = w i k-1 p y k |x i k|k-1 , m i k (14a) w i k = wi k Np i=1 wi k (14b)
In (14b) a normalization is applied to ensure that Np i=1 w i k = 1. Compared to the update step described in [START_REF] Iglésis | A jump-markov regularized particle filter for the estimation of ambiguous sensor faults[END_REF], an additional feature was introduced. Indeed, a Kalman update on the particles x i k|k-1 is applied to place the particles in more likely regions of the state space, which is shown in numerical simulations to enhance estimation performance accuracy. The Kalman update is given by:

x i k = x i k|k-1 + K k ỹi k (15)
At this step, the likelihood is assumed to follow a Gaussian distribution.

Algorithm 4 Detail of the function UPDATE from Algo

- rithm 1 function UPDATE(x i k , w i k , w i k-1 , x i k|k-1 , K k , S k , y k ) ỹi k ← y k -h k x i k|k-1 Innovation wi k ← w i k-1 N (ỹ i k ; 0, S k ) See (14a) w i k ← wi k Np i=1 wi k x i k ← x i k|k-1 + K k ỹi k See (15) end function

C. Estimation

The estimation step aims to perform a global estimation of the state vectors xk and xk|k-1 , with its associated covariance matrices Pk and Pk|k-1 respectively. This step is described in Algorithm 5.

Algorithm 5 Detail of the function ESTIMATE from Algorithm 1 function ESTIMATE(x, P, x 1:Np , w 1:Np )

x ← Np i=1 w i x i P ← Np i=1 w i x i -x x i -x end function

D. Regularization-Resampling Step

This step consists of two stages, the resampling and the regularization of the selected particles. Its purpose is to remove the particles with a low likelihood and to replace them by duplicating the particles with a high likelihood and regularizing the duplicated particles.

a) Resampling step: In the MULTINOMIAL function of algorithm 1, the particles are selected according to a multinomial law with w i k as parameter. Then the probability to choose a particle is:

P xj k = x i k = w i k (16)
b) Regularization step: The regularization step is described in Algorithm 6. The particles are randomly moved according to a regularization kernel K (x). The regularization equation is:

x i k = xi k + hD k ε i k ( 17
)
where h ∈ R + * is the bandwidth factor in the re-scaled kernel density K (•) and with P k = D k D k and ε ∼ K (x). The kernel density is a symmetric probability density function such that:

xK (x) dx = 0, x 2 K (x) dx < ∞ (18)
The optimal kernel K (•) and bandwidth factor h are chosen to minimize the mean integrated square error (MISE) between the theoretical and estimated posterior density, and is given by:

MISE (p) = E (p (x k |Y 1:k ) -p (x k |Y 1:k )) 2 dx k (19
) where p (x k |Y 1:k ) is the JMRPF approximation of the conditional density. The Epanechnikov kernel is used as a regularization kernel [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF].

K (x) = nx+2 2cn x 1 -x 2 if x < 1 0 else ( 20
)
where c nx is the volume of the unit hypersphere in R nx .

Algorithm 6 Detail of the function REGULARIZE from Algorithm 1

function REGULARIZE(x i k , xi k ) ε i k ∼ K x i k see (20) x i k ← xi k + hD k ε i k end function
The resampling step is only performed if N eff /N p is lower than a user-defined threshold Γ ∈ (0; 1), where the efficiency N eff is given by:

N eff = 1 Np i=1 w i k 2 (21)
The conditional density is then approached by:

p (x k , m k |Y 1:k ) ≈ Np i=1 w i k K h x k -x i k δ m i k (m k ) (22)
where:

K h (x k ) = 1 h nx K 1 h x k (23)
VI. SIMULATION RESULTS AND ANALYSIS The application example is the fixed-wing UAV described in Section II. A pitch sensor fault in an inertial navigation system is assumed. The UAV is initialized in straight level flight with an initial velocity of 40 m s -1 and initial altitude of 500 m. The simulation time step is 50 ms.

During all the simulation, the desired flight path angle is set to γ c = 0°and the desired norm of the velocity vector is V c = 40 m s -1 .

A. Filter Parameters

For the two filters compared, the JMRPF and the RPF the stochastic process model of the filter is given by:

     z k+1 Θ k+1 = F k (z k , u k ) Θ k + η k y k = H k (z k ) + G k (Θ k ) + ν k (24a) (24b)
where F k (•) is the longitudinal dynamic of the UAV given by ( 2), and H k (•) is the measurement function given by (3).

As previously explained, a zero-order fault model is used by the filter, but the actual dynamical model of the fault does not always match the fault model of the filter. The process noise η k and measurement noise ν k are both Gaussian and independent with a zero mean and their covariance matrices are respectively Q k and R k . The impact of the sensor fault on the state and measurements is given by:

G k (Θ k ) = 0 0 0 Θ k 0 (25)
The JMRPF and RPF parameters are2 :

• The standard deviation vector used to compute the covariance matrix P 0 = diag σ 2 0 for the extended state vector x k is given by:

σ 0 = 1 1 1 0.3 0.1 0.3 (26)
• The standard deviation vectors used to compute the covariance matrices

Q k = diag σ 2 Q and R k = diag σ 2
R are respectively given by: σ Q = 0.1 0.1 0.1 0.03 0.01 0.1

(27) σ R = 1 1 1 0.3 0.1 (28)
• The resampling threshold Γ and regularization bandwidth h are:

Γ = 0.5 (29) h = 0.2817 (30)
• The number of particles N p is set to 1000.

• The transition probability matrix for sensor faults (JMRPF only) is equal to:

Π = 0.99 0.01 0.01 0.99 (31)
Selection of the transition probability matrix Π was driven by the expected device false alarm probability.

B. Fault Scenario

The simulation scenario was chosen to illustrate the efficiency of the JMRPF and RPF algorithms to estimate intermittent pitch sensor faults with unknown amplitude and unknown dynamics. To assess those two aspects of estimation performance, the fault scenario has two phases. On the first phase, an abrupt bias type sensor fault is injected with unexpected amplitude starting at t = 10 s and ending at t = 20 s. On the second phase, the fault also has unexpected dynamics, starting from t = 30 s until t = 40 s.

1) Unexpected Amplitude: Estimation performance is assessed for faults amplitudes much larger than process and regularization noises. To ensure a fair comparison, the process and regularization noises of the fault state are taken to be the same for both filter and are given by the 6 th element of ( 27) and (30). The abrupt amplitude fault follows zero-order dynamics, as expected by both filters where the assumed fault dynamics are given by (24a). The fault estimates are initialized at zero. The amplitude of the pitch fault is taken to be 5°which is significant compared to the process and regularization noises of both filters. The fault equation for this phase is given by:

Θ (t) = 5 (32)
2) Unexpected Fault Dynamics: When the fault has unexpected dynamics, both filters continue to assume zeroorder filter dynamics. To compare the effect of unexpected fault dynamics on the two filters, an incipient fault with an exponential dynamic and a maximum amplitude of 10°is considered for this phase. The actual sensor fault equation for this phase is given by:

Θ (t) = 10 exp (t -40) (33)
The full sensor fault sequence with both phases is illustrated in Fig. 3.

C. Estimation Performance Comparison

In this section, the JMRPF and RPF are compared in terms of root-mean-square error (RMSE) and of median fault and state estimation for N MC = 100 Monte Carlo runs. The RMSE is defined for any state variable x as:

RMSE x k = NMC i=1 (x i k -x i k ) 2 N MC ( 34
)
where x i k and xi k respectively represent the true and estimated state variable at time step k for the i th Monte Carlo run and x is a generic notation that could for example represent pitch θ or altitude p d . True states depend on the Monte Carlo run because the autopilot uses estimated state feedback. The average RMSE over the total simulation time is also used as an estimation performance metric and is given by:

RMSE x = Ns k=1 RMSE x k N s (35)
where N s is the total number of time steps for the simulation. In Fig. 4, both filters estimate the abrupt fault between 10 s and 20 s and return to the fault free mode between 20 s and 30 s when the fault is no longer active. The JMRPF has higher accuracy and significantly lower convergence time compared to the RPF during both fault free and faulty phases When the transient fault phase starts at t = 30 s, the JMRPF estimates the unexpected exponential fault more and more accurately than the RPF as time approaches 40 s. At that time, the fault is no longer active and the JMRPF returns to the fault free mode almost instantaneously and with very high estimation accuracy compared to the RPF that converges much more slowly to the fault free mode. The jump strategy of the JMRPF was indeed designed to return more rapidly to the correct mode. Moreover, between 40 s and 42.55 s, the RPF does not converge. This non convergence is due to the fact that the estimated fault amplitude of approximately 6.12°at those times is high compared to the process and regularization noises. The JMRPF handles this abrupt change more efficiently. In Fig. 5, the RMSE is indeed shown to be significantly lower with the JMRPF at all except for a brief 1.35 s interval of time at 37.5 s when the exponential fault has not yet compensated the slow convergence dynamics of the RPF. The RMSE is only very briefly higher with the JMRPF at t = 20 s and t = 40 s before a sharp decrease to substantially smaller errors. This is due to the fact that the fault is suddenly activated and deactivated and during one time step, there is an error between the current fault estimate and the new fault. In Fig. 6, the pitch estimation error is also shown to be significantly smaller with the JMRPF from t = 10 s until the end of the simulation. Maximum error is below 0.74°with the JMRPF at 38.65 s and exceeds 2.23°with the RPF at 10.60 s.

Note that the true pitch signal differs between the two simulations because estimated pitch is used by the LQR feedback controller with integral action, which was presented in the control system Subsection II-B.

The JMRPF also provides significantly higher altitude estimation accuracy compared to the RPF, as shown in Fig. 7. This is particularly the case for the duration of the faults and when the faults are deactivated. True altitudes are different between the two filters because the feedback control scheme uses state estimation feedback.

The controller is designed to maintain a zero flight path angle. Fault estimation errors have a visible impact on both filters as also shown in Fig. 7 when the altitude tracking accuracy temporarily deteriorates. The same figure also shows that the altitude fluctuations are less significant with the JMRPF compared to the RPF. This is due to the fact that the JMRPF has higher estimation accuracy. The nonconvergence of the RPF is also visible between 40 s and 42.55 s when the estimated altitude does not follow the true altitude. For this scenario, the average altitude RMSE was The amplitude used for this simulation scenario has been bound by the constraints on RPF parameters tuning. The difference between the convergence characteristics of the JMRPF and the RPF is more pronounced when larger amplitudes are considered. To illustrate this, the same simulation has been run with fault amplitudes 10 times larger than the one previously introduced. The filters parameters remains unchanged.

In Fig. 8, the fault estimation performance of the JMRPF is not significantly affected by a tenfold increase in fault amplitude, without changing the JMRPF parameters. On the other hand, the RPF was unable to handle this fault amplitude and does not converge to the abrupt fault, and at 38 s the incipient fault is also too steep for it to converge.

In Fig. 9, the RMSE of the pitch sensor fault of the RPF is also shown not to converge when the amplitude of the fault is too large compared to its process and regularization noises. Indeed, with no a priori knowledge of the fault amplitude, process noise cannot be adjusted to increase when the unknown fault amplitude increases. An arbitrarily large process noise would also degrade performance in the fault free case or when the fault is small. That is why the same moderate value of process noise was adopted for both RPF and JMRPF. Using the JMRPF, the jump Markov and Kalman correction compensate for the lack of process and regularization noises when fault amplitudes are high and drive the particles towards the high-likelihood regions. This promotes the particles' diversity and reduces particle degeneracy and sample impoverishment issues, and thus enhances the ability to track high amplitude abrupt and incipient faults. For the JMRPF, this figure illustrates its capacity to handle a variety of fault amplitudes and dynamics despite them being significantly different from the zero-order fault model used by the filter. A new nonlinear filter that combines regularized particle filtering with a jump Markov strategy and a Kalman correction was used in this paper to estimate abrupt and incipient sensor faults with unexpected dynamics and amplitudes. The approach was applied to state and fault estimation for a nonlinear longitudinal fixed wing unmanned aerial vehicle (UAV) model with a pitch sensor fault in the inertial navigation system. Numerical simulations have shown that the jump Markov regularized particle filter (JMRPF) significantly outperforms the regularized particle filter (RPF) with a higher accuracy and a faster convergence to both faulty and fault free modes. State estimation accuracy is also significantly enhanced with better robustness to faults using the JMRPF, particularly when the sensor fault is active and when the system switches back to a fault free mode. The simulation results also illustrate that the JMRPF has a much higher ability to estimate larger fault amplitudes compared to the RPF which is better suited to small fault amplitudes with slower dynamics. This presents a practical benefit for a fixed wing UAV with no sensor redundancy to maintain a steady flight in the presence of sensor faults with unexpected dynamics and amplitudes.

Fig. 1 .

 1 Fig. 1. Side view of an UAV with references axis and angles

Fig. 3 .

 3 Fig. 3. Faults scenario used for the simulation

Fig. 4 .

 4 Fig. 4. Median result of the estimated pitch sensor fault

Fig. 5 .

 5 Fig. 5. RMSE of the pitch sensor fault

Fig. 6 .

 6 Fig. 6. Comparison of the median results for estimated pitch

Fig. 7 .

 7 Fig. 7. Comparison of the median results for estimated altitude

Fig. 8 .

 8 Fig. 8. Median result of the estimated pitch sensor fault with 10 times larger fault amplitudes

Fig. 9 .

 9 Fig. 9. RMSE of the pitch sensor fault with 10 times larger fault amplitudes

  with the JMRPF and 0.932 m with the RPF.This average RMSE is therefore reduced by 77 % by the JMRPF and the pitch RMSE is also reduced by 89 %, as shown in as shown in TableI.

			TABLE I
	RMSE RESULTS FOR THE SCENARIO WITH 10°FAULT AMPLITUDE
	Times (s)	10	21	30	41	RMSE
	Altitude
	RMSE JMRPF (m)	0.201 0.215 0.203 0.239 0.215
	RMSE RPF (m)	0.223 1.278 0.527 1.246 0.932
	Pitch
	RMSE JMRPF (°)	0.061 0.082 0.044 0.123 0.114
	RMSE RPF (°)	0.208 0.884 0.704 0.899 1.014
	found to be 0.215 m

In our case a straight flight at an altitude of 500 m with an air speed of

Note: The angle unit used is degree


Songlin Li
email: songlin.li@lip6.fr

Julien Denoulet
email: julien.denoulet@lip6.fr

Olivier Tsiakaka
email: olivier.tsiakaka.1@ulaval.ca

Sylvain Feruglio
email: sylvain.feruglio@lip6.fr

Multi-domain modeling and simulation of an oximeter: PVT variations impact of opto-electronic devices on the SpO 2 quantification

Keywords: Heterogeneous system, modeling, Monte-Carlo simulation, multi-domain, Near-Infrared Spectroscopy (NIRS), oxygen saturation (SpO 2 ), PhotoPlethysmoGraphy

PPG), PVT (Process, Voltage, Temperature), SystemC(-AMS).

INTRODUCTION

SpO 2 (oxygen saturation) index is measured to show the level of oxygenation in blood using an oximeter device. It is based on Near-InfraRed Spectroscopy (NIRS) and on PhotoPlethysmoGraphy (PPG) technologies [START_REF] Webster | Design of Pulse Oximeters[END_REF]. Typically, as shown in Fig. 1, two colored light sources (red and infrared) are used on body parts, usually a finger, to detect volumetric changes of blood. However, deviations always occur when performing a measurement. Therefore, device calibration is indispensable to obtain an acceptable error rate [START_REF]Guidelines for SpO2 measurement usig the Maxim MAX32664 sensor hub[END_REF]. Many well-known factors have an impact on the SpO 2 measurement, for instance motion artifacts, ambient light, skin color, etc. However, so far, influences coming from the system itself haven't been carefully studied. We propose to contribute to this field through the study of the impact of PVT (Process, Voltage, Temperature) variations in the SpO 2 monitoring system. To do this, we intend to model and simulate the embedded system in association with the monitored biological tissue.

Evaluating PVT variations is difficult to achieve through experimental methods, because it is not easy to vary certain parameters in the monitoring system (such as the peak wavelength of the red/IR LED spectrum, the spectral sensitivity of the photo-detector, etc.) to evaluate their impact on SpO 2 quantification. Parasitic signals and noise might also influence the results. Therefore, we will add these features in our modeling methodology, to give us a better knowledge of the measurement device.

We chose the SystemC and SystemC-AMS modeling tools to develop our models due to three main advantages compared to other modeling tools/languages, which are summarized in Table 1:  Multi models of computation: SystemC/ SystemC-AMS offer several Models of Computation (MoC) with Discrete Event (DE), Timed Data Flow (TDF), Linear Signal Flow (LSF) and Electrical Linear Networks (ELN). Each MoC operates at different abstraction levels and can therefore be chosen according to the desired precision;  Fast simulation: with well-chosen abstraction level and simulation time step, simulation run time can be shortened as much as possible;  Open source: SystemC and SystemC-AMS are freely distributed C++ libraries, which means we have access to these tools as well as all C++ supported functions. If necessary, it is also relatively easy to integrate other software in our model by implementing interfaces. Other advantages, like high data accuracy, lightweight installation and simplicity of usage, are also reasons why we chose this tool. We used versions 2.3.2 of SystemC and 2.1 of SystemC-AMS. Both of these libraries can be freely downloaded on Accellera's website (Accellera, n.d.).

SystemC/SystemC-AMS was proposed as an efficient tool for modeling complex systems, especially in the biomedical domain [START_REF] Pecheux | Is SystemC-AMS an Appropriate "Promoter" for the Modeling and Simulation of Bio-Compatible Systems?[END_REF], but it doesn't natively include Monte-Carlo (MC) statistical analysis method [START_REF] Menčík | Concise Reliability for Engineers[END_REF]. Indeed, this analysis is particularly relevant to observe the impact of key parameters on the system behavior. Therefore, a MC simulation method in SystemC-AMS must be developed.

Nowadays, modeling and simulation analysis methods have been widely used in the field of bio- medicine (e.g. bionic nervous systems modeling [START_REF] Cacciapuoti | Receiver Design for a Bionic Nervous System: Modeling the Dendritic Processing Power[END_REF]) and the research of medical device [START_REF] Mundt | Applications of SPICE for Modeling Miniaturized Biomedical Sensor Systems[END_REF]. The interaction modeling between medical device and the human body also exists, such as the interaction between ears and cochlear implants [START_REF] Tran | Development of HEATHER for Cochlear Implant Stimulation Using a New Modeling Workflow[END_REF], and the interaction simulation between the cardiac rhythm and pacemakers [START_REF] Greenhut | A Stochastic Network Model of the Interaction between Cardiac Rhythm and Artificial Pacemaker[END_REF].

In this context, the main objective of this article is to present our modeling and simulation approach of a highly multi-domain application with opensource tools, and to explore impacts of PVT variations in this biomedical device.

The paper is organized as follows. In section 2, we present our model structure. Then, the simulation process and the model validation method are shown in section 3. PVT variation analysis is described in section 4. Finally, we conclude in the last section.

MODEL STRUCTURE

There are mainly two parts in our model: the monitoring system, with opto-electrical components and the biological environment. Figure 2 shows the model block diagram (with the MoC used for each sub_block), with the finger as the monitored biological medium.

In this figure, Num_core block was modeled in SystemC, and a state machine was implemented to emulate a micro-controller. This block was also used to provide pulsed power supply for both LEDs by controlling a LED Driver, and to receive data from the ADC channel.

LED

A red LED (LHQ974) and an IR LED (SFH4080) manufactured by OSRAM were modeled in SystemC-AMS (note that any other type of light source can subsequently be modeled, as long as parameters describing it are available). Both LED models are similar, except for their optical and electronic parameters. The LED model is divided into five sub-blocks to simulate static and dynamic behaviors, in association with its optical characteristics. Sub-block TDF_R, simulates the static behavior described by Shockley equation [START_REF] Sze | Physics of Semiconductor Devices: Sze/Physics[END_REF], mainly. Blue sub-blocks are implemented with the ELN MoC to simulate the dynamic behavior of the LED. Figure 3 (left) shows the physical equivalent circuit of both LHQ974 and SFH4080 sub-blocks. In this equivalent circuit, C d and C j represent the diffusion and the junction capacitance, respectively. i D0 is the equivalent current source, which takes the value of the output current of TDF_R and TDF_IR. The feedback subblocks are used to transfer two feedback parameters: the current i D and the voltage V D of the LED, from LHQ974 to TDF_R or from SFH4080 to TDF_IR. Figure 3 (right) represents the equivalent circuit of LHQ974o and SFH4080o, which are used to take into account the intrinsic opto-electrical effect of the device. i equ is equal to i D . The RC pole, represented by the resistance R RC and the capacitance C RC , designates the frequency characteristic of the LED [START_REF] Bian | A Novel Spice Model of Photodetector for OEIC Design[END_REF]. Then, the P_R and P_IR subblocks are used to convert the output current i R of LHQ974o and SFH4080o into luminous intensity in Watt. The actual spectral distribution of both LEDs is usually modeled by a Gaussian curve [START_REF] Tsiakaka | Source-Detector Spectral Pairing-Related Inaccuracies in Pulse Oximetry: Evaluation of the Wavelength Shift[END_REF] and the total radiant flux is the integral of the entire spectrum.

For the oxygen rate calculation, since the molar extinction ε(λ) varies according to the wavelength, the attenuation for different wavelengths of light is calculated separately, when a beam of light passes through the biological tissue. Therefore, we should separate the luminous intensity of each wavelength in the LED model. Still, we cannot consider all wavelengths of the LED spectrum, since it would greatly increase the amount of calculation and, therefore, the simulation time. As a result, we approximated the LED spectrum by dividing it into five intervals around the peak wavelength. For the red LED (LHQ974), the peak wavelength is ideally 660nm. Thus, we chose five wavelengths from 640nm to 680nm with an interval of 10nm. For the IR LED (SFH4080), with its peak at 880nm, the chosen range was 860nm to 900nm, with the same interval.

Photo-detector

A model of the BP104S photo-detector (PD), also manufactured by OSRAM, was developed. The PD operates in reverse bias. It provides a spectral sensitivity in accordance with the spectrum of the received light to convert the light intensity into electric power. The current generated by the PD mainly consists of two parts:

 I ph : the photonic current generated by the light received from the external environment by the PD, which carries the detection information;  I dc : the parasitic dark current generated by the PD in the absence of wanted light. Sub-block BP104S in Fig. 2 is used to describe the PD spectral sensitivity and the dark current. ELN_RC and ELN_2 are made to simulate the dynamic behavior, as in Fig. 3. The PD output is sent to the ADC block, so that the pulsed current can be sampled, after amplification and voltage conversion, and data can be registered in Num_core.

Biological environment

A model of a finger was constructed to simulate the oxygen concentration variation in the blood and its detection with lights through the tissue. In the blood, there are mainly two chromophores that impact the measurement of SpO 2 : HbO 2 (oxy-hemoglobin) and HHb (deoxy-hemoglobin). The light attenuation by the biological tissue is given by the Beer-Lambert Law, shown in equation ( 1): The change of blood volume according to time is simulated by the normalized equation ( 2), where f 0 is the cardiac frequency in Hz.

d HHb HbO HHb HbO e I I ⋅ ⋅ + ⋅ - ⋅ = ]) [ ] [ ( 0 2 2 ) ( ) (
) 2 2 sin( 4

1 ) 2 sin( 2 1 0 0 t f t f A ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ = π π (2)

SIMULATION PROCESS AND MODELS VALIDATION

This section first introduces the simulation process. Then, in subsection 3.2, the accuracy of our models (at the device and at the system level) is evaluated by comparing our simulation with other experimental results in the existing literature.

Simulation process

Once all models were implemented, basic transient simulation was performed to get two PPG signals that represent the red and IR responses generated by the PD after passing through the biological medium.

Then, as shown in Figure 2, data was sampled every 1ms by the ADC block. The result is presented in Figure 4, where inputs for each LED are light pulses of 150µs every 1ms. The simulation step was set to 1µs, to be able to observe the devices transient behaviors. The heart rate was set to 10Hz (10 times higher than in reality), to reduce the simulation time.

In this figure, the pink signal is the pulsed current generated by the PD. The top envelope corresponds to the IR LED, while the bottom envelope is due to the red LED. The red LED signal amplitude is lower because the luminous intensity of red light is relatively weaker, when the same voltage is applied to both LEDs. The crosses (in the zoomed block) are samples of LED signals. On a laptop with an Intel Core i5 9 th Gen processor (2.4GHz, 8 CPUs) and 8192M of RAM, it takes 29s to simulate a 0.1s duration.

Then, after detecting the peak and valley of the red/IR signals by signal processing (filtering and extremum detection, mainly), we can calculate the R OS (Ratio-of-ratio) value with the following equation [START_REF] Webster | Design of Pulse Oximeters[END_REF]:

                = H L H L OS IR IR R R R ln ln (3)
Here, R L , R H are the valley and peak value of the red signal, respectively, and IR L , IR H correspond to the valley and peak value of the IR signal. SpO 2 value was set before simulation run-time, by setting [HbO 2 ] and [HHb] in the biological environment, according to equation (4) [START_REF] Webster | Design of Pulse Oximeters[END_REF]. We set [HbO 2 ] + [HHb] = 0.3mM, which is an approximate value for the human body [START_REF] Dash | Erratum to: Blood HbO2 and HbCO2 Dissociation Curves at Varied O2, CO2, PH, 2,3-DPG and Temperature Levels[END_REF]. In any case, in the simulation, the value of SpO 2 can be easily and dynamically changed by modifying the ratio of [HbO 2 ] and [HHb].

% 100 ] [ ] [ ] [ 2 2 2 × + = HHb HbO HbO SpO (4)
With a single period of both PPG (red and IR responses), only one R OS value can be obtained for one specific SpO 2 value. However, to obtain a quadratic SpO 2 -R OS curve [START_REF]Guidelines for SpO2 measurement usig the Maxim MAX32664 sensor hub[END_REF], we need to vary this SpO 2 value. In our simulation, we varied this parameter from 90% to 100% (a below 90% SpO 2 value means the person is in a danger state). Figure 5 shows several simulation periods. Red lines correspond to PPG due to the red source for three different SpO 2 settings (90%, 95%, and 100%). For IR responses (blue curves), a shift of the peak wavelength has also been taken into account, to illustrate PVT variation, albeit in a basic manner. We chose four values above and four values below the 880nm nominal value, each with a 5nm step increment, to constitute nine peak wavelength values in the [860nm; 900nm] interval. Combined with the three different SpO 2 parameter values, we were able to simulate twenty-seven configurations. Section 4 will present a more thorough approach to evaluate the influence of PVT variations.

Models validation

There are two steps to validate our models. First, we individually validate each component (i.e. both LEDs and PD) by comparing the simulation output data with key-parameters in the datasheets. These results are summarized in Table 2 (obtained by time domain simulation). Only the optical parameters weren't very well modeled, but with acceptable errors according to datasheets. We assume it was due to approximations made to model certain characterizations and also to the fact that we had the influence of the other parameters. In any case, we can calibrate this error rate by ourselves in the modeling, to be closer to the datasheet value.

The second step is to validate the whole system by comparing our simulation data with previous results [START_REF] Tsiakaka | Source-Detector Spectral Pairing-Related Inaccuracies in Pulse Oximetry: Evaluation of the Wavelength Shift[END_REF], where authors used six different pairs of red/IR LED to plot six SpO 2 -R OS curves. We were able to measure these LEDs optical and electrical characteristics, in order to simulate these devices with our developed model. For the PD, a model of the BP104S was used. As for the biological medium, we applied the same finger Figure 6: a) Simulation result with our models for the six different red/IR LED pairs; b) SpO 2 deviation curves between our works and [START_REF] Tsiakaka | Source-Detector Spectral Pairing-Related Inaccuracies in Pulse Oximetry: Evaluation of the Wavelength Shift[END_REF]. model. Results are presented in Figure 6.a. It shows the "110-25R OS " curve, as it is a widely used reference to basically estimate SpO 2 . We then compared our simulations to Figure 9.b in [START_REF] Tsiakaka | Source-Detector Spectral Pairing-Related Inaccuracies in Pulse Oximetry: Evaluation of the Wavelength Shift[END_REF], with the only difference being the used PD (BP104S for us and FDS100 for the previous work). Figure 6.b shows the SpO 2 deviation curves between our models and the previous results, for all six red/IR LED pairs. We can notice that the absolute value of ΔSpO 2 is always less than 1%. Then, we also calculated the SpO 2 Root-Mean-Square Deviation (RMSD) for each pair of SpO 2 -R OS curves in both works. It ranges from 0.22% to 0.77%, which proves that our simulation results correlate well with the previous results. In this way, our whole system model is validated.

PVT VARIATION ANALYSIS

The dispersion of the IR LED peak wavelength, as simulated in Figure 5, relies on a deterministic approach. For a more in-depth analysis of the impact of PVT variations, it is necessary to introduce a randomness factor. As a result, we performed a MC analysis on our system, based on the variation of selected parameters. As their impact on SpO 2 quantification can be rather unpredictable, a Corner simulation will not necessarily give us enough information. Therefore, MC analysis will allow us to get a more complete view of the system behavior.

MC simulation in SystemC-AMS

In this analysis, we varied selected parameters in a reasonable range, with a chosen Probability Density Function (PDF) of amplitude centered in a specific value. Since this feature isn't natively included in SystemC-AMS, we had to develop our own code. Firstly, we wrote a function to generate a set of numbers varying within a certain range in a specific distribution around 1. In this article, we used a Gaussian distribution. However, other PDF are possible (e.g., uniform). Then, these numbers were multiplied by the nominal value of the parameter to be varied. It constituted a set of input values. Then, we assigned these input values to a SystemC signal and connected it to the systemC-AMS input port of the target parameter. A loop was then implemented in the top file, where we assigned the random input values to the SystemC signal, one by one, to do repeated simulations. Finally, we generated a .dat file to save the output data.

Introduction of different PVT parameters in the system

PVT considers the variation of selected parameters in an electronic system: Process (P) dispersion due to manufacturing, mainly, supply Voltage (V) of the various system components, and the operating Temperature (T) of the device. The parameters we took into account are presented in Table 3, with their typical values and variation ranges. T-variation in the device is mainly due to selfheating during operation and also to the light sources during emission. Note that a variation of more than 4K is generally unacceptable in a medical device. A device power supply voltage is often affected by many factors causing instability and fluctuations (IR-drop due to current draw, crosstalk because of electromagnetic interference, etc.). Moreover, in an embedded system with battery, the voltage supply decreases with the battery usage. Such instability may impact the accuracy of our SpO 2 measure. In this study, we focused on the power supply of both LEDs and of the PD, knowing that the V REF for the ADC doesn't result in any error at all. P-variation is the deviation of component attributes during fabrication. For both LEDs, n is the ideality factor, a key parameter describing the diode junction and a solar cell's electrical behavior [START_REF] Hadj | Theoretical Models for Anomalously High Ideality Factor in a Au/SnO2-Si(n)/Al Solar Cell[END_REF]. n can slightly vary during fabrication. As shown in equation ( 5), it has an influence on the forward current I, which is directly proportional to the optical power of the light source. V F, the threshold voltage of the LED is also affected by manufacturing dispersion. Both of the parameters can be concatenated in a same Gaussian PDF. In equation ( 5), I S is the saturation current, and V T the thermal voltage. We found in the datasheets that the maximal V F for the red and IR LED is 2.5V and 1.8V, respectively. From simulations, we estimated the maximal value of n for both LEDs. It is 2.55 and 1.64, respectively, and can't be less than 1 [START_REF] Sze | Physics of Semiconductor Devices: Sze/Physics[END_REF]. So, we varied n from 1 to 2.55 for the red LED and from 1 to 1.64 for the IR LED.

Deviation of the LED peak wavelength is another relevant parameter. Even for the same type of LED, the peak wavelength of the optical spectrum is subject to deviation due to the fabrication process. We found a typical range of ±20nm for the peak wavelength of both LEDs in their datasheets.

The FWHM (Full Width at Half Maximum) parameter of the optical spectrum can also be impacted by process dispersion. Based on experimental results found in [START_REF] Filippo | LEDs: Sources and Intrinsically Bandwidth-Limited Detectors[END_REF], we chose a range of ±10% of the typical value.

Finally, the PD spectral sensitivity is also affected by an offset due to process variation. Thus, we studied the impact of this factor by shifting the PD optical response spectrum vertically and laterally.

Simulation Results

To get an idea of the influence of each of the parameters presented above, we first performed the MC simulation (consisting of 200 runs), varying only one single parameter. We were then able to calculate for this parameter the maximal SpO 2 RMSD for its two extreme values. These results are presented in the last line of Table 3. We found that the main factors that impact the SpO 2 quantification come from the deviation of both LED peak wavelengths (RMSD is 11.31% for the red LED and 1.79% for the IR LED) and from temperature (RMSD is 3.19%). The gap between both LED RMSD values can be explained by the fact that around the red light band (660nm), the slope of molar extinction curves of HbO 2 and HHb is greater than around the IR light band (880nm). The impact of other parameters is negligible.

After that, we investigated the combined influence of several parameters dispersion on the quantification of SpO 2 . Figure 7 presents the result of the MC simulation correlated with transient analysis. We varied the three key parameters (both LED peak wavelengths and T). We performed 400 runs (i.e., different configurations) to obtain 400 SpO 2 -R OS curves, as in Figure 7.a. To get these results, the simulation time was around 11h. Figure 7.b shows the distribution of SpO 2 when R OS is equal to 0.4. The RMSD of SpO 2 for the two extreme cases Table 3: PVT parameters to be varied and analyzed in the system. is 9.32% in the critical 90-100% saturation window, which is close to the dispersion value associated to the red LED peak wavelength variation. Consequently, it could be said that the red LED peak wavelength variation has the greatest impact.

We can conclude from the above results that in the oximeter manufacturing process, it is necessary to tightly control the peak wavelength deviation of the light source, to avoid an otherwise necessary calibration. At the same time, the device operating temperature influence on the SpO2 measurement accuracy cannot be ignored. For other parameters of our discussion, there is no strict requirement.

CONCLUSION

In this paper, the process to quantify SpO 2 on the finger with an oximeter is simulated with an optoelectronic model built in SystemC/SystemC-AMS. Then, the impact of PVT variations in the device on the SpO 2 quantification is explored, through a MC method combined with transient simulation, performed on the developed models. We found that the main parameters of PVT variations on the quantification of SpO 2 were the red/IR LED peak wavelength deviation and temperature. Other factors seem to have a negligible impact individually.

This modeling method in SystemC/SystemC-AMS, which associates the monitoring system, including its opto-electronic components, and the biological environment, is proved fast, accurate and flexible. This methodology can be employed for any cyber-physical system to estimate its performances, to optimize the design phase, and to help the understanding of measurement data (e.g. reproduce results close to the experimental measurements in the simulation and vary certain parameters to understand their impacts).

Figure 1 :

 1 Figure 1: The principle of oximeter.

Figure 2 :

 2 Figure 2: Model structure to emulate oximeter on finger in SystemC/SystemC-AMS.

Figure 3 :

 3 Figure 3: Equivalent circuit of LED in dynamic mode. On the left is the physical equivalent circuit. The block on the right has no physical reality, but is used to consider the device intrinsic opto-electrical effect.

  λ) the input light intensity at the wavelength λ, I(λ) the output light intensity passing through and attenuated by the biological tissue, ε HbO2 , ε HHb the molar extinction of HbO 2 and HHb, respectively, [HbO 2 ] and [HHb] the concentration of HbO 2 and HHb, and d the length of optical path.

Figure 5 :

 5 Figure 5: PPG obtained with red and IR sources for three SpO 2 value. IR LED curves (in blue) also feature variation of the peak wavelength.

Figure 4 :

 4 Figure 4: Red signal (bottom), IR signal (top) generated by photo-detector and the samples of the signals (crosses on the top and on the bottom).

Figure 7 :

 7 Figure 7: Simulation result of the variation of three key factors in the same time.

Table 1 :

 1 SystemC/SystemC-AMS versus other modeling tools.

Table 2 :

 2 Models key-parameters compared to datasheets.