Serge Romaric Tembo
email: serge.tembo@orange.com

Jean - Luc

Sandrine Vaton
email: sandrine.vaton@telecom-bretagne.eu

A 3-Layered Self-Reconfigurable Generic Model For Self-Diagnosis of Telecommunication Networks

Keywords: Self-diagnosis, Self-reconfiguration, fault propagation, Bayesian network, Probabilistic Inference, GPON, FTTH. www.conference.thesai.org

The dynamic and distributed nature of telecommunication networks makes complex the design of model-based approaches for network fault diagnosis. Most model-based approaches assume the prior existence of the model which is reduced to a static image of the network. Such models become rapidly obsolete when the network changes. We propose in this paper a 3-layered self-reconfigurable generic model of fault diagnosis in telecommunication networks. The Layer 1 of the model is an undirected graph which models the network topology. Network behavior, also called fault propagation, is modeled in Layer 2 using a set of directed acyclic graphs interconnected via the Layer 1. We handle uncertainties of fault propagation by quantifying strengths of dependencies between Layer 2 nodes with conditional probability distributions estimated from network generated data. Layer 3 is the junction tree representation of the loopy obtained Layer 2 Bayesian networks. The junction tree is the diagnosis computational layer since exact inference algorithms fail on loopy bayesian networks. This generic model embeds intelligent self-reconfiguration capabilities in order to track some changes in network topology and network behavior. These self-reconfiguration capabilities are highlighted through some example scenarios that we describe. We apply this 3-layered generic model to carry out active self-diagnosis of the GPON-FTTH access network. We present and analyze some experimental diagnosis results carried out by running a Python implementation of the generic model.

I. INTRODUCTION

Reliability, robustness, availability, accessibility are the most important requirements that telecommunication networks should guarantee. As the design of network architectures, network management also became a central issue for telecommunication operators, which have triggered significant researches about autonomic networking. The main goal of autonomic networking is to automate as much as possible numerous tedious operations of network management like fault management.

Fault diagnosis is a central aspect of network fault management [START_REF] Steinder | A survey of fault localization techniques in computer networks[END_REF]. Traditionally, fault diagnosis has been performed manually by an expert or a group of experts experienced in managing communication networks [START_REF] Steinder | Multi-layer fault localization using probabilistic inference in bipartite dependency graph[END_REF]. However, the development of telecommunication networks has increased the size and complexity of their architectures. Fault diagnosis has become too complex for humans, who can keep track of only few hypotheses in their reasonings. Humans need a long training to fully master their network segment [START_REF] Hounkonnou | Active self-diagnosis in telecommunication networks[END_REF]. Automating fault diagnosis is critical in large scale telecommunication networks. A network fault or failure is a root cause of one or multiple network anomalies observed in the form of alarms or parameters outside of the standard. Alarms or symptoms are external manifestations of failures [START_REF] Jakobson | Alarm correlation[END_REF]. A telecommunication network is naturally a distributed system, a fault occurred spreads, triggering other faults and alarms which in turn trigger further faults and alarms. The consequence of both fault and alarm propagation is that a single root cause may result in a complex and distributed pattern of subsequent failures and their corresponding alarms [START_REF] Hounkonnou | Active self-diagnosis in telecommunication networks[END_REF]. This is especially true when multiple faults propagate simultaneously.

Fault diagnosis also called alarm correlation or root cause diagnosis [START_REF] Jakobson | Alarm correlation[END_REF], [START_REF] Katzela | Schemes for fault identification in communication networks[END_REF], [START_REF] Yemini | High speed and robust event correlation[END_REF], isolates the most probable set of faults based on their external manifestations. It is a process of searching intelligible explanations to observed alarms. Fault diagnosis is a complex process due to the non-deterministic nature of fault propagation phenomenon. A single fault may generate multiples alarms and a single alarm may be triggered by several faults.

The fault diagnosis problem was addressed in the past two decades and numerous techniques have been proposed. Expert systems [START_REF] Pau | Survey of expert systems for fault detection, test generation and maintenance[END_REF], [START_REF] Scherer | A survey of expert systems for equipment maintenance and diagnostics. Fault detection and reliability: knowledge based and other approaches[END_REF] encode specialized reasonings on specific diagnosis tasks in computer applications. Model-based approaches [START_REF] Gardner | Alarm correlation and network fault resolution using the kohonen self-organising map[END_REF] develop reasonings based on an explicit representation of the network. The techniques based on machine learning algorithms like artificial neural networks [START_REF] Patton | Fault diagnosis in nonlinear dynamic systems via neural networks[END_REF], [START_REF] Goel | Towards a 'neural' architecture for abductive reasoning[END_REF], probabilistic networks [START_REF] Steinder | Multi-layer fault localization using probabilistic inference in bipartite dependency graph[END_REF], [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF] and case-based reasoning [START_REF] Lewis | A case-based reasoning approach to the resolution of faults in communication networks[END_REF] infer diagnosis based on past experiences.

A brief description and discussion of related work on fault diagnosis as well as our objectives is presented in section 2. We describe and formalize in section 3, our contribution: a generic three layered self-reconfigurable probabilistic model for telecommunication networks fault diagnosis. In section 4, we focus on self-reconfiguration capabilities of the generic model by showing how this model can be turn into an intelligent autonomous system for self-diagnosis purposes. In section 5, we present an application of the generic model to fault diagnosis of the FTTH (Fiber To The Home) access networks based on GPON (Gigabit capable Passive Optical Network). In section 6, some experimental diagnosis results are presented and analysed with respect to the GPON-FTTH network. We conclude and present future works in section 7.

II. RELATED WORKS

Early self-diagnosis approaches were called expert systems [START_REF] Pau | Survey of expert systems for fault detection, test generation and maintenance[END_REF], [START_REF] Scherer | A survey of expert systems for equipment maintenance and diagnostics. Fault detection and reliability: knowledge based and other approaches[END_REF]. An expert diagnosis system attempts to infer the cause of a problem from symptoms recognized in sensor data [START_REF] Hounkonnou | Active self-diagnosis in telecommunication networks[END_REF]. It is a problem-solving software that embodies specialized reasonings on narrow diagnosis tasks usually performed by a trained skilled human called expert. The specialized reasonings can be formalized with rules, list of facts, logic predicates, etc. The inference engines are commonly based on forwardbackward [START_REF] Madsen | Lazy propagation: A junction tree inference algorithm based on lazy evaluation[END_REF] chaining algorithms. However, expert systems mainly suffer of brittleness. The system fails when it faces a novel problem that comes out of its expertise.

Model-based diagnosis approaches [START_REF] Gruschke | Integrated event management: Event correlation using dependency graphs[END_REF], [START_REF] Houck | Towards a practical alarm correlation system[END_REF], [START_REF] Jordaan | Event correlation in heterogeneous networks using the osi management framework[END_REF], [START_REF] Kätker | A modeling framework for integrated distributed systems fault management[END_REF], [START_REF] Kätker | A generic model for fault isolation in integrated management systems[END_REF], [START_REF] Kätker | Fault isolation and event correlation for integrated fault management[END_REF] develop reasonings based on formal and explicit representation of network structure and network behavior. Network structure describes the network architecture. Network behavior describes the process of alarm propagation and alarm correlation [START_REF] Jakobson | Real-time telecommunication network management: extending event correlation with temporal constraints[END_REF]. Network structure and network behavior are then modeled [START_REF] Gardner | Alarm correlation and network fault resolution using the kohonen self-organising map[END_REF]. The obtained model is the support of reasoning algorithms which must be designed. The modelbased approach is easy to deploy and is appropriate for a largescale network if information on network resources is available [START_REF] Miyazawa | Scalable root cause analysis assisted by classified alarm information model based algorithm[END_REF]. It has the ability to deal with novel problems although its performances degrade in this case.

Model-based approach seems natural when relationships between objects are graph-like and easy to obtain [START_REF] Steinder | A survey of fault localization techniques in computer networks[END_REF]. The model can be designed in a modular or incremental fashion facilitating updates as new knowledge about the network is acquired. However, it is quite difficult to build a model close enough to the structural and functional reality of the network while maintaining a high level of abstraction to make the model independent of the various engineering techniques implemented in telecommunication networks. In addition, the model built is reduced to a static image of the network and becomes rapidly obsolete when the network changes.

In order to deal with the difficulty to obtain and selfmaintain an accurate model of fault propagation in large scale telecommunication networks, we propose in this paper a generic 3-layered self-reconfigurable probabilistic model for self-diagnosis of telecommunication networks. The model integrates two fields: a decision field and an artificial learning field. The decision field is based on Bayesian probabilistic reasoning [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF] in order to deal with uncertainties of fault propagation process. The artificial learning [START_REF] Cornuéjols | Apprentissage Artificiel, Concepts et Algorithmes[END_REF] field brings selfreconfiguration capabilities to the generic model in order to deal with the dynamic nature of telecommunication networks. The generic model has capabilities to compute diagnosis decisions and to automatically learn changes in network topology and network behaviour.

III. A LAYERED SELF-RECONFIGURABLE PROBABILISTIC MODEL

We focus in this section on the description and the formalism of the 3-layered generic model which can be applied to any distributed system for self-diagnosis purposes.

A. Description of the generic model

As in any distributed system, in telecommunication networks, faults typically propagate between related system components. There is also some cases where a fault can propagate only inside a system component without affecting the state of neighbor system components. We call this, local fault The generic model that we propose clearly considers local and distributed fault propagation by separating network topology modeling and fault propagation modeling. Note that this separation brings important properties to the generic model: modularity and extensibility which facilitates selfreconfiguration. The generic model has three layers. Layer 1 models the network topology. Layer 2 models fault propagation. Layer 3 is a simplification of layer 1 and layer 2.

The layer 1 models the network topology and distributed fault propagation between related network components. Layer 1 is an undirected graph whose nodes represent network components and edges represent bidirectional links between them (see Figure 1). Any node of this layer is a Directed Acyclic Graph (DAG) which embeds distributed dependencies between a network component and its neighbors. These distributed dependencies model fault propagation between the related network components. The layer 1 of the generic model is an undirected graph of DAGs, i.e, a layer 1 node is seen as a subset of layer 2 nodes linked via distributed dependencies. In Figure 1, layer 1 represents the topology of a network of n components E 1 , E 2 , ..., E n . The network component E 1 is linked to one neighbor E 2 . The network components E 1 and E 2 are each DAGs containing distributed dependencies that carry fault propagation between them. For example, in Figure 2, each of the components E 1 and E 3 embeds distributed dependencies with its neighboor E 2 . The component E 2 embeds distributed dependencies with its neighboors E 1 and E 3 . In Figure 1, the roles of layer 1 agents cL 1 A, rL 1 A, layer 2 agents cL 2 A, rL 2 A and layer 3 agents cL 3 A, rL 3 A will be discussed later in the paper.

The layer 2 is a set of DAGs. Each DAG models local fault propagation inside a network component, i.e. a layer 1 node. In Figure 1, a layer 2 DAG B i models local fault propagation on network component E i . The DAGs of layer 2 are interconnected via layer 1 nodes. In order to deal with uncertainties of fault propagation phenomenon, each DAG can be transformed into a Bayesian network by quantifying strength of dependencies with Conditional Probability Distributions (CPD). A CPD quantifies the degree of influence that a subset of nodes of a DAG has on their common successor node. In Figure 1, Layer 2 contains n belief networks interconnected via Layer 1. The belief network B i models local fault propagation inside the network component E i . We call Layer 2 node a variable of any belief network.

Note that the union of DAGs of layer 2 and DAG of each layer 1 node gives one large DAG which models fault propagation in the entire network. But the separation between network topology modeling and network behavior (fault propagation) modeling breaks down this large DAG in many interconnected parts. So the design of layer 1 and layer 2 implements this separation and brings us the modularity and extensibility properties necessary to provide easy selfreconfiguration capabilities to the generic model.

The distributed nature of a telecommunication network introduces mutual dependencies between related network components. This means that distributed undirected loops are unavoidable between some layer 2 nodes. In addition, undirected loops may also appear in local fault propagation inside a network component. Note that, carrying out a message passing inference algorithm like the sum-product algorithm [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF] at layer 2 will fail due to local and distributed undirected loops. In order to deal with local and distributed undirected loops between layer 2 nodes, we have built a third layer upon layer 2. The layer 3 is the junction tree representation [START_REF] Madsen | Lazy propagation: A junction tree inference algorithm based on lazy evaluation[END_REF], [START_REF] Miyazawa | Scalable root cause analysis assisted by classified alarm information model based algorithm[END_REF], [START_REF] Stefen Laurizen | Local computations with probabilities on graphical structures and their application to expert systems[END_REF] of layer 1 and layer 2.

Note that layer 1 and layer 2 are useful for learning field of the generic model since they provide modularity and extensibility properties useful for easy self-reconfiguration. Layer 3 is useful for decision field of the generic model since inference can be easily done in this layer regardless of network topology complexity at layer 1 and network behavior complexity at layer 2.

B. Formalism of the generic model

We formalize in this subsection the three layers of the generic model. We first define a local dependency and a distributed dependency. Let B i and B j =i be two distinct layer 2 Directed Acyclic Graphs (DAGs). A local dependency on B i is a tuple (u, v) such that u and v belong to the set of nodes of B i . A distributed dependency between B i and B j =i is a tuple (u, v) such that u belongs to the set of nodes of B i and v belongs to the set of nodes of B j . For example, in Figure 2,(B,D) is a local dependency which is part of the local fault propagation model on network component

E 1 . (G, D) is a distributed dependency between the related network components E 1 and E 2 . Layer 1 denoted by L 1 Graph is an undirected graph of DAGs, L 1 Graph = (L 1 N odes, L 1 Edges). L 1 N odes is the set of n layer 1 nodes or network components considered, L 1 N odes = {L 1 N 1 , L 1 N 2 , ..., L 1 N n }. L 1
Edges is the set of bidirectional links between layer 1 nodes.

Each layer 1 node L 1 N i is a DAG defined as follows:

L 1 N i = (L 1 V i , L 1 E i ). L 1 V i is the subset of layer 2 nodes embedded in layer 1 node L 1 N i . Each layer 2 node belonging to L 1 V i is part to at least one distributed dependency. L 1 E i
is the set of distributed dependencies between L 1 N i and its k i ≥ 1 layer 1 neigboor nodes. L 1 V i is defined as follows:

L 1 V i = N i ∪ ki j=1 N i j (1)
where,

N i ⊂ L 2 V i , N i j ⊂ L 2 V j , j ∈ {1, ..., n} \ {i} such that ∀u ∈ N i , ∃v ∈ N i j |u ∈ pa(v) ∨ v ∈ pa(u).
We suppose that layer 1 node L 1 N i has k i neighboors. L 2 V i is the subset of layer 2 nodes belonging to layer 2 DAG L 2 N i above layer 1 node L 1 N i (see Figure 1). The layer 2 DAG L 2 N i models local fault propagation on layer 1 node L 1 N i . For example, in Figure 2,

L 1 N 2 = (L 1 V 2 , L 1 E 2 ) such that L 1 V 2 = N 2 N 2 1 ∪ N 2 3 , where the subset N 2 = {G, I, J}, N 2 1 = {D} and N 2 3 = {M, O}. L 1 E 2 is the set of distributed dependencies between the component L 1 N 2 and its neighboors L 1 N 1 and L 1 N 3 . L 1 E 2 = {(G, D), (G, M ), (D, I), (O, I), (D, J), (O, J)} Layer 2 is a set of n DAGs. Each DAG L 2 N i = (L 2 V i , L 2 E i ) models local fault propagation on network com- ponent L 1 N i . L 2 E i is the set of local dependencies between layer 2 nodes belonging to L 2 V i .
Layer 3 is the junction tree representation of the large bayesian network G obtained by combining layer 1 and layer 2 as follows:

G = n i=1 [L 2 N i ∪ L 1 N i ].
Layer 3 is constructed in three steps. The moralization [START_REF] Naim | Réseaux Bayésiens[END_REF] of G and the triangulation [START_REF] Naim | Réseaux Bayésiens[END_REF] of the obtained moralized graph lead to a chordal graph or clique graph. The moralization consist to disorient edges of the graph G and adding a disoriented edge between each couple of parents of each node of G. The triangulation consist to add an edge to every cycle of the moralized graph whose length exceeds 3. The maximal weight spanning tree of the obtained clique graph is guaranteed to be a junction tree. The weight is the size of the intersection between adjacent cliques, i.e, the number of layer 2 nodes shared by adjacent cliques. We call layer 3 node, a clique C i of the junction tree. As every junction tree, layer 3 satisfies the running intersection property which ensures that, the intersection C i ∩C j is a subset of every clique and separator on the path between two cliques C i and C j . For example, in Figure 2, G = C GDBC ∩ C GKM L : G belongs to the clique C GM O and the separators S GO and S GM which form the path between the cliques C GDBC and C GKM L .

C. Diagnosis computations of the generic model

The layer 3 is an equivalent of layer 1 and layer 2. This means that layer 3 is sufficient to compute diagnosis decisions using for example the well known HUGIN exact inference algorithm [START_REF] Madsen | Lazy propagation: A junction tree inference algorithm based on lazy evaluation[END_REF] on a junction tree. In Figure 1, layer 3 nodes C i , C j , S ij are respectively the cliques C i and C j of the junction tree and their common separator S ij . A layer 3 node is a compound variable of some layer 2 nodes. For example, in Figure 2, layer 3 node ABC is a compound variable of Layer 2 nodes A, B, and C. Layer 3 is initialized by associating a potential to each layer 3 node. At initialization, a clique C i has the potential φ Ci and a separator S ij has the potential φ Sij as follows:

φ Ci = X∈Layer2,X∈Ci,pa(X)⊂Ci∨pa(X)=∅ P (X|pa(X)) (2) And φ Sij = 1 (3)
We note pa(X), the parent set of layer 2 node X. The potential φ Ci of a layer 3 node C i represents the joint conditional probability of layer 2 nodes that compose it (see equation 2). For example, in Figure 2: φ C ABC = P (B|A).P (C|A).P (A), φ C GDBC = P (D|B, C, G).P (G) and φ S BC = 1.

The diagnosis decisions computed at layer 3 are based on evidences propagation on a junction tree. Assume C i , C j to be neighbooring layer 3 nodes with their common separator S ij (see Figure 1). The potential φ Ci is updated when some layer 2 nodes belonging to clique C i are observed. This evidence propagates from clique C i to clique C j towards their separator S ij as follows: updating the potential of the separator S ij by the marginalization of equation 4, and updating the potential of the clique C j by the product of equation 5. The notation C i \ S ij reprensents the layer 2 nodes of the clique C i which does not belong to the separator S ij .

φ * Sij = Ci\Sij φ * Ci (4) φ * Cj = φ Cj φ * Sij φ Sij (5)
We say that the clique C j absorbs evidences from C i [START_REF] Madsen | Lazy propagation: A junction tree inference algorithm based on lazy evaluation[END_REF] or that the clique C i brings evidences to C j . Note that, update operations of clique potentials and separator potentials are done in two recursive stages. The first stage called collect is initiated by collecting evidences from observed nodes to the root nodes. The second stage is initiated by distributing evidences from the root nodes to leave nodes. Collection of evidences to a clique C i is done by collecting evidences to all the children of C i followed by absorption of evidences from each child. Similarly, distribution of evidences from a clique amounts to bring evidences to each child followed by distribution of evidences from the child.

Note that, the update operations of potentials of layer 3 nodes are made by an agent called computing Layer 3 Agent cL 3 A. The cL 3 A updates the potential of each layer 3 node when it receives some evidences (observed layer 2 nodes) from its counterpart of layer 2 called cL 2 A. When cL 2 A receives updated layer 3 node potentials, it computes the marginals (beliefs) of layer 2 nodes. These two agents may communicate using a simple mechanism like shared memory through an interface between layer 2 and layer 3 that we call L2-L3 Interface (see Figure 3). The marginal P (X|e) of a layer 2 node X consistent with evidences e, is computed from the updated potential φ * Ci of any layer 3 node C i such that X ∈ C i : P (X|e) = Ci\X φ * Ci . Note that, after updating the potential of all layer 3 nodes, the intersection property ensures coherence between updated potentials, i.e, the marginal of a layer 2 node X is the same regardless of the clique C i on which the marginalization operation is made. The most probable state of the layer 2 node X that is consistent with evidences, is one that has the highest probability.

Note that, there is an alternative approach to find the diagnosis without performing summation operations on updated potentials. With this approach, the diagnosis r * , is computed from the most probable explanation w * , of evidences as follows:

w * = C∈Layer3 w * C such that φ C (w * C ) = max w C φ C (w C )
, where w C is a configuration of layer 2 nodes belonging to layer 3 node C and w * C maximizes the potential φ C of C. The diagnosis r * is the most probable configuration of layer 2 root nodes defined by r * = w * \ i * , where i * is the most probable configuration of non root layer 2 nodes consistent with evidences.

IV. SELF-RECONFIGURATION CAPABILITIES OF THE GENERIC MODEL

In this subsection, we illustrate the potential of the generic model to self-reconfigure in order to track changes in network topology and network behavior (how faults propagate). Each of the three layers of the generic model has a reconfiguration agent which communicates by shared memory with its counterpart of the adjacent layers. The reconfiguration Layer 2 Agent rL 2 A is an artificial learning [START_REF] Cornuéjols | Apprentissage Artificiel, Concepts et Algorithmes[END_REF] system which implements knowledge discovery and data mining algorithms from a tremendous amount of network generated data. This agent must have the ability to discover new local and distributed dependencies, new states of layer 2 nodes to consider, new network statistical knowledges which require a new estimation of conditional probabilities of some layer 2 nodes. When the rL 2 A discovers some new important knowledge about network behavior, it reconfigures the layer 2 of the model and sends a reconfiguration command (reconf cmd) to the reconfiguration Layer 3 Agent rL 3 A and to the reconfiguration Layer 1 Agent rL 1 A (see Figure 4). The message reconf done is an acknowledgement of the command reconf cmd.

Note that rL 2 A will send a reconfiguration command to rL 1 A only if the new knowledge learned concerns distributed dependencies. Indeed, local dependencies are unknown at layer 1. Note also that it is the type of the new knowledge learned at layer 2 that will permit rL 3 A, after receiving a reconfiguration command from rL 2 A, to know if it should reconfigure or rebuild the layer 3. For example, if the conditional probability distribution of a layer 2 node changes, the rL 3 A needs only to raise an event notifying the computing Layer 3 Agent cL 3 A, to recompute the initial potential of all layer 3 nodes which this layer 2 node belongs to. If new local and/or distributed dependencies are learned, the rL 3 A must rebuild the layer 3. Note that any change in network topology requires the reconfiguration of the three layers of the model. The Figure 5 illustrates this reconfiguration process initiated by the layer 1 reconfiguration agent. The dynamic behaviour of the rL 1 A is modelled by a 3-states machine (see Figure 6). Initially, the model is at equilibrium, i.e, there is no running reconfiguration process. So the rL 1 A is in Data mining state, i.e, it waits for specification of a new network component by a network expert, it searchs for new knowledges about network topology from data collected by carried out active measures on network. When, the rL 1 A discovers new knowledges about network topology, e.g, distributed dependencies, it does not reconfigure the layer 1 immediately. It sends a reconf cmd to the rL 2 A and transits to the Blocked state. In Blocked state, the rL 1 A waits for the completion of upper layers reconfiguration before entering in Layer 1 Reconf state. The rL 1 A transits from Blocked state to Layer 1 Reconf state when it receives a reconf cmd from rL 2 A. In the case where the rL 2 A discovers new distributed dependencies, the rL 1 A transits from Data mining to the Layer 1 Reconf state when it receives a reconf cmd from rL 2 A, and returns to the Data mining state when it completes the layer 1 reconfiguration. The completion of the layer 1 reconfiguration yields a new equlibrium state of the model.

Layer 1 reconfiguration adds or removes new layer 1 node (network component) to/from the model. For adding a node for example, the new layer 1 node embeds distributed dependencies with its layer 1 neighboor nodes which already exist in the model. The rL 1 A sends a reconf cmd to rL 2 A (see Figure 5), waits for the rL 2 A completes the layer 2 reconfiguration and adds the distributed dependencies to the neighboor nodes of the new layer 1 node. To reconfigure the layer 2, the rL 2 A builds a layer 2 Bayesian network which models local fault propagation on the new network component. The rL 2 A notifies the computing Layer 2 Agent cL 2 A to recompute the conditional probability distribution of each layer 2 nodes having new distributed dependencies learned or specified at layer 1. Remark that only the reconfiguration agent of layer 2 takes a reconfiguration decision based on knowledges learned from network generated data. The other reconfiguration agents (rL 1 A and rL 3 A) wait for a reconf cmd from rL 2 A before starting a reconfiguration process. Therefore, the rL 2 A is the coordinator of our multi-agent based self-reconfigurable model. The Figure 7 specifies the dynamic behaviour of the layer 2 reconfiguration agent using a 3-states machine. Initially, the rL 2 A is in Data mining state, i.e, it searchs for new knowledges about network behaviour from tremendous amount of data generated by the network components. When the rL 2 A learns new local dependencies, new distributed dependencies, new states of a layer 2 node to consider or new important statistical knowledges, it transits to the Layer 2 Reconf state in order to integrate the new knowledges learned to the layer 2 of the generic model. The rL 2 A should trigger the reconfiguration of adjacents layer by sending a reconfiguration command to the reconfiguration agent of these layers. The Blocked state implements synchronisation between two successive reconfiguration processes, i.e, the rL 2 A waits for the running reconfiguration process of the model completes before triggering another process. Doing so, the completeness and consistency of the model remain guaranteed after the integration of new learned knowledges to the generic model.

V. APPLICATION OF THE GENERIC LAYERED MODEL TO GPON-FTTH SELF-DIAGNOSIS

We show in this section how the generic model previously proposed can be applied to perform self-diagnosis of FTTH (Fiber To The Home) access networks based on GPON (Gigabit capable Passive Optical Network) [START_REF] Of | [END_REF], [START_REF] Of | 1 Recommendation[END_REF]. For now, we used prior knowledges acquired from ITU-T standards [START_REF] Of | [END_REF], [START_REF] Of | 1 Recommendation[END_REF] of the GPON-FTTH network to build a 3-layered model of this network. We plan to transform this GPON-FTTH network model into an autonomous system by implementing self-reconfiguration capabilities of the generic model. The Figure 8 depicts the architecture of the GPON-FTTH network. The GPON-FTTH access network has two main network components. The Optical Line Termination (OLT) is located on operator side. The Optical Network Termination (ONT) is located on customer side. The OLT and ONT are connected through an Optical Distribution Network (ODN). The ODN is the optical infrastructure made of fibers and passive components like splitters. A Passive Optical network (PON) is a point-to-multipoint link inside the ODN. A PON has a tree-like topology which connects an OLT with a maximum of 64 ONTs (see Figure 8). Each ONT is connected to a RG (Residential Gateway) via an Ethernet link.

Since there is no interaction between PONs in ODN, and all PONs have the same behavior, we modeled one single PON. This model can be replicated to any PON of a GPON-FTTH access network. All ONTs connected to the same PON temporally share the upstream optical channel of the PON. The downstream channel of the PON is a secured broadcasting channel. The composition of Figure 9 and Figure 10 forms the application of the generic model for modeling the topology and behavior of a PON of the GPON-FTTH access network. The obtained model has two layer 1 nodes (OLT, ONT) and some layer 2 nodes are vectors in order to consider the tree-like topology of a PON.

A. Upstream modeling

The Figure 9 depicts the layer 1 OLT node and its layer 2 local dependencies. Distributed dependencies embedded in layer 1 OLT node are also depicted. We have three types of layer 2 nodes in Figure 9: faults or root causes, intermediate faults and alarms. The root causes are highlighted in Figure 9.

The transport optical fiber of the PON denoted by F iber T (OK, AT, BR) can take three states. The state OK means that there is no transmissions anomaly on this fiber. The states AT and BR mean respectively that fiber experiences high attenuation or that fiber is broken. The temperature of OLT denoted by T c OLT , is a continous variable that we discretize. The power supply of OLT denoted by Alt OLT which can be faulty or not. Now, suppose we need to extend this model by adding another layer 1 node: a Residential Gateway RG node for example (see Figure 8). A RG is a home network component that provides services to customers. It is not a GPON-FTTH network component, but adding it to the model will bring to the model the ability to correlate faults and alarms of the GPON-FTTH network with malfunctions of customer services. For example, in Figure 10, we can make such correlation with the distributed dependency N 1 , N 2 , LOF ds , LCGD ds -→ IP Connect, where IP Connect denotes the Internet access service provided by a RG. To add the layer 1 RG node, we will need only to specify and quantify uncertainties of its local dependencies and distributed dependencies with the layer 1 ONT node, which already exists in the model. We start the evaluation of the designed probabilistic fault propagation model for GPON-FTTH self-diagnosis by ckecking if inference on this model returns no fault when no alarm is observed on the PON and when counters and scalar parameters of OLT and of the two ONTs are nominal. The Table I presents the result of this test and shows that the positive state in bold, i.e, the good working state of each layer 2 root cause node is the most probable state. The Table II shows that the ON T 1 is diagnosed to be faulty when it loses upstream and downstream communication with OLT although the optical signal between them is not degraded (no alarm related to signal quality is observed). Note the missing values of counters and scalar parameters of ON T 1 since the GPON-FTTH network management system has no way to get these values from MIB (Management Information Base) of ON T 1 . These missing values useful to compute the most probable diagnosis are inferred by the probabilistic model by using its conditional probability distributions. The

VI. EXPERIMENTAL RESULTS OF THE GPON-FTTH

VII. CONCLUSION

We have presented in this paper a generic 3-layered probabilistic model of fault diagnosis of telecommunication networks. The generic model embeds two fields: a decision field and an artificial learning field. The generic model has capabilities to compute diagnosis decisions and to automatically learn changes in network topology and network behaviour.

We have shown how this generic model can be applied to design a probabilistic model-based approach of fault diagnosis of GPON-FTTH access network. We have presented and analyzed some network fault diagnosis results carried out by running a Python implementation of the decision field of the generic model. The results obtained are very consistent with experiments carried out on a GPON-FTTH network testbed.

We plan to transform this GPON-FTTH network faults diagnosis model into an autonomous system by implementing self-reconfiguration capabilities of the generic model. Another perspective is to bring more intelligence to the reconfiguration agent of layer 3 by designing a novel paradigm allowing it to take into account new dependencies learned to the lower layers without rebuilding all the layer 3. Doing so, scalability and reconfiguration efficiency of the generic model will increase.

Fig. 1 :

 1 Fig. 1: The 3-Layered Generic Model.

Fig. 2 :

 2 Fig. 2: A simple example of the 3-Layered Generic Model.

Fig. 3 :

 3 Fig. 3: Communication between layer 2 and layer 3 for beliefs updating of layer 2 nodes.

Fig. 4 :

 4 Fig. 4: Self-reconfiguration process initiated at layer 2.

Fig. 5 :

 5 Fig. 5: Self-Reconfiguration process initiated at layer 1.

Fig. 6 :

 6 Fig. 6: State machine of the layer 1 reconfiguration agent.

Fig. 7 :

 7 Fig. 7: State machine of the layer 2 reconfiguration agent.

Fig. 8 :

 8 Fig. 8: GPON-FTTH network architecture.

  The node F aultyON T denotes an ONT which transmits upstream signal outside of its granted time slot, which may conflict with data sent by other ONTs on the PON and cause data disruption for a random set of ONTs, making the PON unusable. A FaultyONT can cause a Drift of Windows DOW . OLT raises a DOW [i] alarm when an ON T [i] transmits signal beyond the slot time allocated to it. See ITU-T G984.3[START_REF] Of | [END_REF],[START_REF] Of | 1 Recommendation[END_REF] for more details.The Software Version SW V [i] alarm means that there is an incompatibility between the Image Operating System (IOS) of ON T [i] and those of the OLT. The node ON T Conf M is (ONT Configuration Mistake) denotes a configuration error during ONT provisioning.The OLT transmitted power T x OLT is regulated by the bias current I OLT . This leads to the local dependency I OLT -→ T x OLT . The OLT received power Rx OLT [i] from an ON T i depends of the OLT voltage V OLT and the state of the transport fiber F iber T .The OLT received power of ON T i also depends of the state of the branching fiber denoted by F iberDB[i] and the transmitted power of this ON T i denoted by T x ON T [i]. Note two local dependencies (F iber T -→ Rx OLT , V OLT -→ Rx OLT ) and two distributed dependencies (F iberDB -→ Rx OLT , T x ON T -→ Rx OLT ). Remark in Figure9that the distributed dependencies are part of edges of the layer 1 OLT node which is a DAG as designed in the generic model.

Fig. 9 :

 9 Fig. 9: Faults and alarms propagation raised by an OLT.

Fig. 10 :

 10 Fig. 10: Faults and alarms propagation raised by an ONT.

  F aultyON T1 [+f to, -fto] [0.007, 0.993] F aultyON T2 [+f to, -fto] [0.014, 0.986] IOS1 [OK, DF ] [0.980, 0.020] IOS2 [OK, DF ] [0.980, 0.020] ON T Conf M is1 [+ocm, -ocm] [0.002, 0.998] ON T Conf M is2 [+ocm, -ocm] [0.002, 0.998] F iberDB1 [OK, AT, BR] [0.91, 8.e-02, 9.e-07] F iberDB2 [OK, AT, BR] [0.91, 8.2-02, 9.e-07] AltON T1 [OK, DF ] [0.99, 1.e-07] AltON T2 [OK, DF ] [0.99, 1.e-07]

  NETWORK SELF-DIAGNOSISWe present and analyse in this section the network fault diagnosis results carried out by our application and implementation of the generic model to perform active self-diagnosis of GPON-FTTH access network. The experiments are performed on a GPON-FTTH testbed on which we have considered a PON with two ONTs named ON T 1 and ON T 2 . This experimental testbed allows us to emulate faults and to collect alarms raised by the GPON-FTTH network components. We also read, if availaible, the values of counters BIP us , BIP ds and scalar parameters Rx OLT , Rx ON T , T x OLT , T x ON T , voltages V OLT , V ON T , bias current I OLT , I ON T , temperatures T c OLT , T c ON T of OLT and of the two ONTs connected to the PON considered. Each table from TableIto Table VII shows the beliefs of Layer 2 root cause nodes computed by our Python implementation of the decision field of the generic model based on evidences observed and decribed in the title of each table.

TABLE I :

 I No Fault: no observed alarm, counters and scalar parameters are known and are nominal for both ON T 1 and ON T 2 .

		.
	Faults	States	Beliefs
	F iber

T [OK, AT, BR] [0.99, 6.e-03, 7.e-12] AltOLT [OK, DF ] [0.99, 4.e-07]

TABLE VII :

 VII Loss downstream communication between OLT and ON T 1 : LOSds 1 alarm is observed, counters and scalar parameters of ON T 1 are unknown, but known and nominal for ON T 2 . system does not read the values of counters and scalar parameters of ON T 1 . Note that, this is a one way loss of communication. Contrary to a bidirectional loss of communication as illustrated in TableII, ON T 1 is not diagnosed to be faulty. The most probable diagnosis computed by the model is the branching fiber downstream attenuation of ON T 1 , i.e, only the downstream optical channel between OLT and ON T 1 experiences attenuation.

		.
	Faults	States	Beliefs
	F iber T	[OK, AT, BR]	[0.92, 7.e-02, 3.e-11]
	F aultyON T1	[+f to, -f to]	[0.094, 0.906]
	F aultyON T2	[+f to, -f to]	[0.014, 0.986]
	F iberDB1	[OK, AT, BR]	[8.e-02, 0.91, 4.e-06]
	F iberDB2	[OK, AT, BR]	[0.92, 7.e-02, 1.e-06]

management

www.conference.thesai.org

loss of upstream and downstream communication between

OLT and ON T 1 may also be due to the cut of the branching fiber of ON T 1 . But this is a very rare event in reality since infrastructures of telecommunication operators always protect optical fibers from possible cuts. That is why the most probable explanation of a bidirectional loss of communication is a faulty ONT. However, other observations in addition to LOSus 1 and LOSds 1 alarms may change this explanation. In Table III, the observation of alarms notifying optical signal degradation, brings additional information to the model and changes the decision of the inference algorithm. The ON T 1 is no longer diagnosed to be Faulty but its branching fiber experiences important attenuation leading to communication loss with OLT. Remark in Table III that, contrary to the Table II, the belief that ON T 1 is faulty decreases in favor of increasing the belief that the branching fiber of ON T 1 experiences important attenuation.

In Table IV, the most probable diagnosis is the attenuation of the branching fiber of ON T 1 when SD (Signal Degraded) alarms are raised for ON T 1 with its counters and scalar parameters known and non-nominal (upstream and downstream BIP variation, low level upstream and downstream receive power). If for some reasons, the management system of GPON-FTTH network fails to get the values of counters and scalars of ON T 1 , the model is always able to take the best decision although in this case, the belief in the branching fiber attenuation of ON T 1 decreases slightly (from 0.92 to 0.76, see Table V). Now, if we observe signal degradation between OLT and ON T 1 , and between OLT and ON T 2 (see Table VI), the probabilistic model more likely believes that it is the transport fiber of the PON which experiences attenuation rather than each of the two branching fibers. Remark in Table VI that the probability of the transport fiber attenuation is 0.58. The probability of attenuation of each branching fiber is 0.50. The difference between these two probabilities is not very significant. Therefore, we may tend to believe that the model has trouble to take the best decision between the transport fiber attenuation or branching fibers attenuation. But, if we compute the joint probability of the two branching fibers attenuation (0.5 × 0.5 = 0.25), we remark that the difference between this joint probability and the probability of the transport fiber attenuation is now very important. This result enforces the belief that it is the transport fiber of the PON which experiences attenuation. Note that the joint probability of the two branching fibers is the product of the probabilities of the branching fibers since they are independent. In fact the Markov Blanket M B(F iberDB i ) of the variable F iberDB i is defined by the union of the following sets: the parent set P A i of F iberDB i , the children set CH i of F iberDB i and the set OP i of other parents of children of F iberDB i . M B i = P A i ∪ CH i ∪ OP i , where P A i = ∅ (the empty set) since F iberDB i is a layer 2 root cause node. The set CH i = {Rx OLTi , Rx ON Ti } and the set OP i = {T x ON Ti , V ON Ti , V OLT , T x OLT , F iber T }. The branching fibers F iberDB 1 and F iberDB 2 are independent each other because: F iberDB 1 / ∈ M B(F iberDB 2 ) and F iberDB 2 / ∈ M B(F iberDB 1 ).

The Table VII illustrates the case when a downstream loss of communication between OLT and ON T 1 is observed, i.e, the alarm LOSds 1 is raised and the GPON-FTTH network www.conference.thesai.org


Holali Kwassi

Denis Degue

Abderrahman Efimov

Iggidr

K H Degue

D Efimov

Interval observer for sequestered erythrocytes in malaria

For malaria patients, a usual observation problem consists in estimation of sequestered parasites Plasmodium falciparum from measurements of circulating ones. The model of an infected patient is rather uncertain, and for all rates (death, transition, recruitment and infection) in the model it is assumed that only intervals of admissible values are given. In addition, the measurements of the concentration of circulating parasites are subjected by a bounded noise, while some parameters, like the rate of infection of blood cells by merozoites, are completely unknown and highly timevarying. In order to evaluate the concentration of sequestered parasites, an interval observer is designed, which provides intervals of admissible value for that concentration, with the interval width proportional to the model uncertainty. Stability of the proposed observer can be verified by a solution of LMI. The efficiency of the observer is demonstrated in simulation for a model and for real measured data.

I. INTRODUCTION

Malaria is an important disease causing at least one million deaths around the world each year (with ninety percents among children in Africa), and it is spreed by the Plasmodium parasite. The most dangerous type of malaria is summoned by the most virulent species called Plasmodium falciparum. Sequestration is one of the characteristics of Plasmodium falciparum, which is related with the Plasmodium life cycle. The cycle begins when a parasite enters the human body through the bite of an infected mosquito, after which it migrates to the liver and starts to multiply within. The free forms resulting from this multiplication (called merozoites) are able to invade the red blood cells (erythrocytes). The infected erythrocytes are matured during the erythrocytic cycle. At roughly the middle stage of trophozoite development (in 24 hours), molecules on the surface of infected erythrocytes can link to receptors of endothelial cells. This bind has the effect of holding infected erythrocytes within vessels of organs (such as the brain), where they remain until the rupture of the erythrocyte and the release of merozoites. This period of attachment is called sequestration and during it, the infected erythrocytes are not detectable in the blood flow, they are "sequestered". Also it is widely accepted that antimalarial drugs act differently depending on the stage of parasite development [START_REF] Gravenor | Estimating sequestered parasite population dynamics in cerebral malaria[END_REF], [START_REF] Gravenor | A model for estimating total parasite load in falciparum malaria patients[END_REF].

In practice, to know the stage of infection for a patient, the total parasite concentration n i=1 y i in the bloodstream is needed, where y i represent population of parasites of certain age, from the youngest y 1 till the oldest y n , n < 1 determines the grid of age differentiation. However, only the presence of peripheral infected erythrocytes, i.e. the young parasites y 1 + y 2 + ...y k for some k < n, also called circulating, can be detected on peripheral blood smears and the other ones (sequestered y k+1 , ...y n ), which are hidden in some organs like brain and heart, cannot be observed. Nowadays, there is no clinical method of measuring the sequestered infected cells directly, and the measurements of circulating parasites are rather costly and corrupted by noise.

That is why the estimation of sequestered parasite population is an important challenge, with many authors having studied this problem [START_REF] Gravenor | Estimating sequestered parasite population dynamics in cerebral malaria[END_REF], [START_REF] Gravenor | A model for estimating total parasite load in falciparum malaria patients[END_REF], [START_REF] Ochola | Estimating sequestered parasite load in severe malaria patients using both host and parasite markers[END_REF], [START_REF] Bichara | On the estimation of sequestred infected erythrocytes in plasmodium falciparum malaria patients[END_REF]. In this work an interval observer is designed in order to estimate the interval for sequestered parasite population, which is admissible for given uncertainty in patient model. In the presence of uncertainty, which has an important impact in this application, design of a conventional estimator, converging to the ideal value of the state, cannot be realized. In this case an interval estimation becomes more feasible: an observer can be constructed that, using input-output information, evaluates the set of admissible values (interval) for the state at each instant of time. The interval length is proportional to the size of the model uncertainty (it has to be minimized by tuning the observer parameters). There are several approaches to design interval/set-membership estimators [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous time systems[END_REF], [START_REF] Kieffer | Guaranteed nonlinear state estimator for cooperative systems[END_REF], [START_REF] Olivier | Closed loop observers bundle for uncertain biotechnological models[END_REF]. This work is devoted to interval observers, which form a subclass of set-membership estimators and whose design is based on the monotone systems theory [START_REF] Olivier | Closed loop observers bundle for uncertain biotechnological models[END_REF], [START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF], [START_REF] Raïssi | Interval observers design for consistency checks of nonlinear continuous-time systems[END_REF], [START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF], [START_REF] Efimov | Interval estimation for LPV systems applying high order sliding mode techniques[END_REF]. Since the estimated variables take positive values (all elements of the state vector for infected patient are concentrations), then the obtained estimates also should take positive values only, which poses an additional constraint to satisfy in this rather complex estimation problem.

The outline of this paper is as follows. After preliminaries in Section II, the problem statement is given in Section III. The interval observer design is presented in Section IV. Numerical experiments with really measured data and selected model are described in Section V.

II. PRELIMINARIES

The real numbers are denoted by R, R + = {τ ∈ R : τ ≥ 0}. Euclidean norm for a vector x ∈ R n will be denoted as |x|, and for a measurable and locally essentially bounded input u : R + → R the symbol ||u|| [t0,t1] denotes its L ∞ norm:

||u|| [t0,t1] = ess sup t∈[t0,t1] |u(t)|,
if t 1 = +∞ then we will simply write ||u|| ∞ . We will denote as L ∞ the set of all inputs u with the property ||u|| ∞ < ∞. The symbols I n , E n×m and E p denote the identity matrix with dimension n × n, the matrix with all elements equal 1 with dimensions n × m and p × 1, respectively.

A. Interval relations

For two vectors x 1 , x 2 ∈ R n or matrices A 1 , A 2 ∈ R n×n , the relations x 1 ≤ x 2 and A 1 ≤ A 2 are understood elementwise. The relation P ≺ 0 (P 0) means that the matrix P ∈ R n×n is negative (positive) definite. Given a matrix A ∈ R m×n , define A + = max{0, A}, A -= A + -A (similarly for vectors) and denote the matrix of absolute values of all elements by

|A| = A + + A -. Lemma 1. [6] Let x ∈ R n be a vector variable, x ≤ x ≤ x for some x, x ∈ R n . (1) If A ∈ R m×n is a constant matrix, then A + x -A -x ≤ Ax ≤ A + x -A -x. (1)
(2) If A ∈ R m×n is a matrix variable and A ≤ A ≤ A for some A, A ∈ R m×n , then

A + x + -A + x --A -x + + A -x -≤ Ax (2)
≤ A + x + -A + x --A -x + + A -x -. Furthermore, if -A = A ≤ 0 ≤ A, then the inequality (2) can be simplified: -A(x + + x -) ≤ Ax ≤ A(x + + x -).

B. Nonnegative continuous-time linear systems

A matrix A ∈ R n×n is called Hurwitz if all its eigenvalues have negative real parts, and it is called Metzler if all its elements outside the main diagonal are nonnegative. Any solution of the linear system

ẋ = Ax + Bω(t), ω : R + → R q + , ω ∈ L q ∞ , (3)
y = Cx + Dω(t), with x ∈ R n , y ∈ R p and a Metzler matrix A ∈ R n×n , is elementwise nonnegative for all t ≥ 0 provided that x(0) ≥ 0 and B ∈ R n×q + [9], [20]. The output solution y(t) is nonnegative if C ∈ R p×n + and D ∈ R p×q + .
Such dynamical systems are called cooperative (monotone) or nonnegative if only initial conditions in R n + are considered [START_REF] Farina | Positive Linear Systems: Theory and Applications[END_REF], [START_REF] Smith | Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems[END_REF]. For a Metzler matrix A ∈ R n×n its stability can be checked verifying a Linear Programming (LP) problem

A T λ < 0 for some λ ∈ R n + \ {0}, or Lyapunov matrix equation A T P + P A ≺ 0
for a diagonal matrix P ∈ R n×n , P > 0 (in general case the matrix P should not be diagonal). The L 1 and L ∞ gains for nonnegative systems (3), i.e. the gains of transfer function from input to output in different norms, have been studied in [START_REF] Briat | Robust stability analysis of uncertain linear positive systems via integral linear constraints: l 1 -and l∞-gain characterizations[END_REF], [START_REF] Ebihara | L1 gain analysis of linear positive systems and its application[END_REF], for this kind of systems these gains are interrelated.

The conventional results and definitions on L 2 /L ∞ stability for linear systems can be found in [START_REF] Hassan | Nonlinear Systems[END_REF].

C. A non-homogeneous sliding mode differentiator

Let ỹ(t) = y(t) + ν(t) be a measured signal, where y : R + → R is a signal to be differentiated and ν ∈ L ∞ is a bounded measurement noise, then a differentiator can be proposed [START_REF] Efimov | A hybrid robust non-homogeneous finite-time differentiator[END_REF]:

ẋ1 = -α |x 1 -ỹ(t)|sign(x 1 -ỹ(t)) + x 2 , (4) ẋ2 = -sign(x 1 -ỹ(t)) -χsign(x 2 ) -x 2 , x 1 (0) = ỹ(0), x 2 (0) = 0,
where x 1 , x 2 ∈ R are the state variables of the system (4), α, and χ are the tuning parameters with α > 0 and > χ ≥ 0. The variable x 1 (t) serves as an estimate of the function y(t) and x 2 (t) is an estimate of ẏ(t), i.e. it provides the derivative estimate. Therefore, the system (4) has ỹ(t) as the input and ŷ(t) = x 2 (t) as the output. Lemma 2. [START_REF] Efimov | A hybrid robust non-homogeneous finite-time differentiator[END_REF] Let ẏ, ÿ, ν ∈ L ∞ , then there exist α > 0 and > χ ≥ 0 such that x 1 , x 2 ∈ L ∞ and there exist T 0 > 0, c 1 > 0 and c 2 > 0:

|x 2 (t) -ẏ(t)| ≤ c 1 ||ν|| ∞ + c 2 ||ν|| ∞ ∀t ≥ T 0 .
Estimates on T 0 > 0, c 1 > 0, c 2 > 0 and guidelines for tuning α, , χ can also be found in [START_REF] Efimov | A hybrid robust non-homogeneous finite-time differentiator[END_REF].

III. ESTIMATION OF THE HIDDEN ERYTHROCYTES

The dynamics of the parasitized erythrocytes is described by the following model [START_REF] Bichara | On the estimation of sequestred infected erythrocytes in plasmodium falciparum malaria patients[END_REF]:

. z(t) = A(t)z(t) + Eβ(t)x(t)m(t) + e 1 Λ(t) ∀t ≥ 0, Y (t) = Cz(t) + v(t), (5)
where z = (x, y 1 , . . . , y 5 , m) T ∈ R 7 + is the state vector and Y ∈ R + is the measured output, v ∈ L ∞ is the measurement noise, ||v|| ∞ ≤ V for some known V > 0; y 1 and y 2 correspond to the concentrations of free circulating parasitized erythrocytes and y 3 , y 4 , y 5 correspond to the sequestered ones; x is the concentration of healthy cells, and m is the concentration of merozoites; Λ(t) ∈ R + , Λ ∈ L ∞ represents recruitment of the healthy Red Blood Cells (RBC) and β(t) ∈ R + , β ∈ L ∞ is the rate of infection of RBC by merozoites. The variables β(t) and Λ(t) serve as exogenous uncertain inputs in [START_REF] Ebihara | L1 gain analysis of linear positive systems and its application[END_REF]. The time-varying matrix A and the constant matrices C, E, e 1 are defined as follows:

C = [ 0 1 1 0 0 0 0 ], E = [ -1 1 0 0 0 0 -1 ] T , e 1 = [ 1 0 0 0 0 0 0 ] T , A =           -µ x 0 0 0 -µ 1 -γ 1 0 0 γ 1 -µ 2 -γ 2 0 0 γ 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -µ 3 -γ 3 0 0 0 γ 3 -µ 4 -γ 4 0 0 0 γ 4 -µ 5 -γ 5 0 0 0 rγ 5 -µ m          
, where µ x > 0 is the natural death rate of healthy cells; µ i > 0 is the natural death rate of i th stage of infected cells, γ i > 0 is the transition rate from i th stage to (i + 1) th stage of infected cells, i = 1, . . . , 5; r > 0 is the number of merozoites released by the late stage of infected cells, µ m > 0 is the natural death rate of merozoites. For different patients the values of the parameters µ x , µ i , γ i , r and µ x are different and they are varying with time for a patient, that is why we assume that A ≤ A(t) ≤ A for some known A, A ∈ R 7×7 and the instant value of A(t) is unavailable. Such a hypothesis is rather natural since on-line identification of these parameters is complicated. Similarly for the healthy RBC recruitment Λ(t), the values Λ, Λ ∈ R + are given such that

Λ ≤ Λ(t) ≤ Λ ∀t ≥ 0.
It is assumed that for β(t) there is no confidence interval (its instant value is highly uncertain).

We suppose that Y = y 1 + y 2 + v(t), i.e. the circulating Plasmodium concentration can be measured with a noise v, while it is required to estimate the sequestered one Z = y 3 + y 4 + y 5 .

IV. INTERVAL OBSERVER DESIGN

We define w(t) = β(t)x(t)m(t) as a new unmeasurable variable, which can be considered as a new uncertain input for [START_REF] Ebihara | L1 gain analysis of linear positive systems and its application[END_REF]. Following [START_REF] Bichara | On the estimation of sequestred infected erythrocytes in plasmodium falciparum malaria patients[END_REF] and using the equation ( 5), we can find:

w = ((CE) T CE) -1 (CE) T ( . Y -CAz -Ce 1 Λ)
where Ẏ is the derivative of the output. Using Lemma 2 and differentiator (4), an estimate Ŷ of Ẏ can be calculated such that for all t ≥ 0:

Ẏ (t) = Ŷ (t) + v (t),
where ||v || ∞ < V for some known V > 0.

Note that CE = 1 in (5), let 0 ≤ z(t) ≤ z(t) ≤ z(t) for all t ≥ 0 and some z, z ∈ R 7 , then using Lemma 1 we obtain the following relations for all t ≥ 0:

w(t) ≤ w(t) ≤ w(t), where w = Ŷ -V -Ce 1 Λ -(CA) + z + (CA) -z and w = Ŷ + V -Ce 1 Λ -(CA) + z + (CA) -z.
Following [START_REF] Chebotarev | Interval observers for continuous-time LPV systems with l 1 /l 2 performance[END_REF], the equations of an interval observer for (5) take the form:

.

ζ(t) = Aζ(t) + e 1 Λ + E + w(t) -E -w(t) + L(Y (t) -Cζ(t)) -|L|V, . ζ(t) = Aζ(t) + e 1 Λ + E + w(t) (6) -E -w(t) + L(Y (t) -Cζ(t)) + |L|V, z(t) = max{0, ζ(t)}, z(t) = max{0, ζ(t)},
where z ∈ R 7 and z ∈ R 7 are respectively the lower and the upper interval estimates for the state z; ζ, ζ ∈ R 7 is the state of [START_REF] Efimov | Interval estimation for LPV systems applying high order sliding mode techniques[END_REF]. The following restrictions are imposed on (6).

Assumption 1. There exist matrices L ∈ R 7×1 , L ∈ R 7×1 such that the matrices (A-LC) and (A-LC) are Metzler.

Assumption 1 fixes the principal conditions to satisfy for positivity of the error dynamics (due to the structure of A this condition is always satisfied for L = L = 0). Theorem 1. Let Assumption 1 be satisfied. Then for all t ∈ R + the estimates z(t) and z(t) given by ( 6) yield the relations:

0 ≤ z(t) ≤ z(t) ≤ z(t) ∀t ≥ 0, (7)
provided that 0 ≤ z(0) ≤ z(0) ≤ z(0). If in addition, there exists a diagonal matrix P ∈ R 7 and γ > 0 such that

A T P + P A + P (γ -2 I 7 + F F T )P + 2I 7 0 for A = (A -LC) 0 0 (A -LC) , F = E + (CA) -+ E -(CA) + -E + (CA) + -E -(CA) - -E + (CA) + -E -(CA) - E + (CA) -+ E -(CA) + , then z, z ∈ L 7 ∞ (ζ, ζ ∈ L 7
∞ and the transfer function

e 1 Λ + LY (t) -|L|V e 1 Λ + LY (t) + |L|V, → ζ ζ has L ∞ gain less than γ).
Proof. Note that z(t) ≥ 0 for all t ≥ 0 and z(t) is also bounded [START_REF] Bichara | On the estimation of sequestred infected erythrocytes in plasmodium falciparum malaria patients[END_REF]. The equation ( 5) can be rewritten as follows:

.

z = (A -LC)z + (A(t) -A )z + Ew + e 1 Λ + LY -Lv
for some A ∈ R 7×7 (A or A) and L ∈ R 7×1 (L or L), then the dynamics of the errors e(t) = z(t) -ζ(t), z(t) = ζ(t) -z(t) obey the equations:

. e(t) = (A -LC)e(t) + g(t),

.

e(t) = (A -LC)e(t) + g(t), (8)
where

g = (A(t) -A)z + Ew -E + w + E -w + Lv + |L|V, g = (A -A(t))z + E + w -E -w -Ew -Lv + |L|V.
Under the introduced conditions, it can be inferred from Lemma 1 that g(t) ≥ 0, g(t) ≥ 0 ∀t ≥ 0. From Assumption 1, we conclude that e(t) ≥ 0 and e(t) ≥ 0 (g, g have the same property and e(0) ≥ 0 and e(0) ≥ 0 by conditions). That implies that the order relation

ζ(t) ≤ z(t) ≤ ζ(t)
is satisfied for all t ≥ 0, then ( 7) is true by construction of z, z and due to nonnegativity of z.

In order to prove boundedness, let us define:

ζ = [ζ T ζ T ] T , = [ T T ] T , = Ce 1 (E -Λ -E + Λ) + e 1 Λ + E Ŷ -|E|V + LY -|L|V, = Ce 1 (E -Λ -E + Λ) + e 1 Λ + E Ŷ + |E|V + LY + |L|V,
then dynamics of interval observer takes the form:

ζ = Aζ + F max{0, ζ} + ,
where the matrices A and F are defined in the theorem formulation and ∈ L 14 ∞ by construction. Consider a Lyapunov function

V (ζ) = ζ T P ζ, then V = ζ T (A T P + P A)ζ + 2ζ T P [F max{0, ζ} + ] ≤ ζ T [A T P + P A + P (γ -2 I + F F T )P + I]ζ + γ 2 T ≤ -ζ T ζ + γ 2 T
and the needed stability conclusion follows.

The obtained interval estimates z, z are nonnegative as the state z is. Note that the presented approach can be easily extended to higher/lower order models of parasitized erythrocytes (when age partition of erythrocytes has more/less than 5 levels as in ( 5)).

Remark 1. Using the results of [START_REF] Briat | Robust stability analysis of uncertain linear positive systems via integral linear constraints: l 1 -and l∞-gain characterizations[END_REF], [START_REF] Ebihara | L1 gain analysis of linear positive systems and its application[END_REF], the stability conditions for (6) can be expressed as a linear programming (LP) problem, then it is necessary to find λ ∈ R 14 + \ {0} such that the following LP problem is feasible:

Aλ + |F |E 14
λ -E 14 < 0.

V. SIMULATION OF THE INTERVAL OBSERVER

In this section the result of numerical experiments are reported for the model [START_REF] Ebihara | L1 gain analysis of linear positive systems and its application[END_REF] and with a real data measured for an infected patient.

A. Model experiments

In [START_REF] Bichara | On the estimation of sequestred infected erythrocytes in plasmodium falciparum malaria patients[END_REF], for a patients without fever, i.e. at 37 • C, the parameters of the matrix A have the following constant values: γ 1 = 1.96, γ 2 = 3.78, γ 3 = 2.85, γ 4 = 1.76, γ 5 = 3.26;

µ 1 = 0, µ 2 = 1.86, µ 3 = 0, µ 4 = 0.1, µ 5 = 0; (9)
µ x = 1 120 , r = 16, µ m = 72.
Assume that admissible deviations of these parameters from the nominal values given above are σ%, then we can calculate the matrices A and A. The nominal value of healthy RBC recruitment is Λ 0 = 5×10 6 120 cellsµl -1 day -1 (as in [START_REF] Bichara | On the estimation of sequestred infected erythrocytes in plasmodium falciparum malaria patients[END_REF], the unit of volume is micro-liter (µl) and the unit of time is day) with admissible deviations ±20%, i.e. For simulations we selected:

0.8Λ

0 = Λ ≤ Λ(t) ≤ Λ = 1.2Λ 0 ∀t ≥ 0.
Λ(t) = Λ 0 (1 + 0.2 sin(2t)),
β(t) = 10 -6 (1 + 0.25 sin(5t))e mod 2 (t,2.5+0.5 sin(0.25t)) , v(t) = V sin(35t), V = 10, A(t) = sin 2 (3t)A + cos 2 (3t)A.

B. Real data experiments

The measured concentration of circulating parasites (peripheral parasitaemia, that is density of parasites in the blood) is taken from the data collected by the US Public Health Service at the National Institutes of Health laboratories in Columbia, South Carolina and Milledgeville, Georgia [START_REF] William | A retrospective examination of sporozoite-and trophozoite-induced infections with plasmodium falciparum in patients previously infected with heterologous species of plasmodium: effect on development of parasitologic and clinical immunity[END_REF], [START_REF] Eichner | Genesis, sequestration and survival of plasmodium falciparum gametocytes: parameter estimates from fitting a model to malariatherapy data[END_REF], where malaria was used for therapy of neurosyphilis, and patients were inoculated through mosquito bite or infected blood (see also [START_REF] Bichara | On the estimation of sequestred infected erythrocytes in plasmodium falciparum malaria patients[END_REF] for more details). The data of patient named S1204 is used in this work. As in [START_REF] Bichara | On the estimation of sequestred infected erythrocytes in plasmodium falciparum malaria patients[END_REF] the parameters of the model have been selected the same as in [START_REF] Farina | Positive Linear Systems: Theory and Applications[END_REF]. For differentiator (4), α = 3 × 10 4 , = 3α and χ = 0.25α, with V = 100 and V = 6000. The interval of uncertainty for Λ(t) has been selected the same as previously. The level of the matrix A uncertainty σ has been adjusted in order to have correct interval estimates for the measured variable Y (t), the obtained value is σ = 31, then the results of interval estimation for both variables, Y (t) and Z(t), are shown in Fig. 2 (the interval estimates for sequestered parasites, Z(t) and Z(t), are plotted in the logarithmic scale). As we can conclude, in order to have adequate interval estimates for Y (t), the value of σ has been selected sufficiently big, which signalizes that the used model [START_REF] Ebihara | L1 gain analysis of linear positive systems and its application[END_REF] with the parameters (9) are rather far from reality, then the interval estimates are rather approximate. From another side, this approach also allows the model accuracy to be evaluated by the measured data.

VI. CONCLUSION

An interval observer is proposed in this work in order to estimate the sequestered parasite population from the measured circulating parasites. It is assumed that almost all parameters and inputs of the model are uncertain (just intervals of admissible values are given) and the measurements are obtained with a noise. Despite of that the proposed observer demonstrates a reasonable accuracy of interval estimation, which is confirmed by numerical experiments.

In a future work the sampled-time kind of measurements can be taken into account.

Figure 1 .

 1 Figure 1. The results of interval estimation for sequestered parasites

  Let z(0) =1 3 z(0) = [500 100 150 50 50 50 50] T . For differentiator (4), α = 2 × 10 3 , = 3α and χ = 0.25α, then V = 80V and L = 1γ 1 0 0 0 0] T have been selected. The results of interval estimation of sequestered Plasmodium Z(t) are shown in Fig.1for σ = 5 and σ = 10. As we can conclude, the size of uncertainty σ influences directly the estimation accuracy.

Figure 2 .

 2 Figure 2. The results of interval estimation on real data

Federation (Grant 074-U01) and the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031).


Khalid Daoudi
email: khalid.daoudi@inria.fr

Jaya Ashwini
email: ashwini.jaya-kumar@inria.fr

Kumar

Pitch-based speech perturbation measures using a novel GCI detection algorithm: Application to pathological voice classification

Keywords: Perturbation measures, Pitch marks, Jitter, Shimmer, Pathological speech classification

Classical pitch-based perturbation measures, such as Jitter and Shimmer, are generally based on detection algorithms of pitch marks which assume the existence of a periodic pitch pattern and/or rely on the linear source-filter speech model. While these assumptions can hold for normal speech, they are generally not valid for pathological speech. The latter can indeed present strong aperiodicity, nonlinearity and turbulence/noise. Recently, we introduced on a novel nonlinear algorithm for Glottal Closure Instants (GCI) detection which has the strong advantage of not making such assumptions. In this paper, we use this new algorithm to compute standard pitch-based perturbation measures and compare its performances to the widely used tool PRAAT. We address the task of classification between normal and pathological speech, and carry out the experiments using the popular MEEI database. The results show that our algorithm leads to significantly higher classification accuracy than PRAAT. Moreover, some important statistical features become significantly discriminative, while they are meaningless when using PRAAT (in the sense that they have almost no discrimination power).

Introduction

Most of the classical approaches in speech processing are based on linear techniques that may not adequately capture the complex dynamics of speech. Indeed, in normal speech, some of the well-known examples of non-linear phenomena include: the existence of turbulent sound source in production process of unvoiced sounds, the existence of a time spread and turbulent component for the excitation source of plosives (which is idealized as an impulse in the linear framework) and the evidences regarding characterization of voiced sounds by highly complex air flows like jets and vortices. In pathological speech, linear methods are definitely not enough to characterize the strong aperiodicity, the turbulence and breathy noise that can be present in such signals.

In particular, dysphonia analysis is generally performed using pitch-based perturbation measures such as Jitter, Shimmer and their variants. Most of the algorithms to compute these measure use linear methods at the frame-basis level and implicitly assume some kind of periodicity constraints. As a canonical example, the widely used tool PRAAT [START_REF] Boersma | Praat: doing phonetics by computer[END_REF] uses a short-term autocorrelation method [START_REF]Praat: a system for doing phonetics by computer[END_REF] for the detection of pitch marks. Such implicit periodicity assumptions can not hold for pathological speech, particularly speech with severe dysphonia.

In [START_REF] Khanagha | Detection of glottal closure instants based on the microcanonical multiscale formalism[END_REF], we recently developed a nonlinear algorithm for Glottal Closure Instants (GCI) detection. This algorithm has the strong advantages of not making any assumption about the signal (such as periodicity of the source), it does not operate a frame-basis and is computationally efficient. The algorithm has been compared to state-of-the-art methods and showed similar performances for clean (normal) speech and significantly better ones on noisy (normal) speech. Precise detection of GCI has many applications in speech technology: closed phase Linear Prediction (LP) analysis [START_REF] Ananthapadmanabha | Epoch extraction from linear prediction residual for identification of closed glottis interval[END_REF][START_REF] Pinson | Pitch synchronous time domain estimation of formant frequencies and bandwidths[END_REF][START_REF] Steiglitz | The use of time-domain selection for improved linear prediction[END_REF][START_REF] Naylor | Estimation of glottal closure instants in voiced speech using the dypsa algorithm[END_REF], pitch synchronous speech processing for converting the pitch and duration of speech [START_REF] Ewender | Accurate pitch marking for prosodic modification of speech segments[END_REF], prosody modification [START_REF] Moulines | Pitch-synchronous waveform processing techniques for text-to-speech synthesis using diphones[END_REF], synthesis [START_REF] Drugman | A deterministic plus stochastic model of the residual signal for improved parametric speech synthesis[END_REF][START_REF] Stylianou | Applying the harmonic plus noise model in concatenative speech synthesis[END_REF], dereverbertion [START_REF] Gaubitch | Spatiotemporal averagingmethod for enhancement of reverberant speech[END_REF], casual-anticasual deconvolution [START_REF] Bozkurt | Mixed-phase speech modeling and formant estimation, using differential phase spectrums[END_REF][START_REF] Drugman | Advances in glottal analysis and its applications[END_REF] and glottal flow estimation [START_REF] Wong | Least squares glottal inverse filtering from the acoustic speech waveform[END_REF]. In this paper, we use our GCI detection algorithm in the framework of pathological voice classification. Indeed, motivated by our good results in noisy (normal) speech, we could fairly expect that our algorithm may also have an interesting behavior in the case of dysphonic speech. We thus proceed to compute classical pitch-based perturbation measures using the GCI locations detected by our algorithm as pitch marks. We then analyze the results we obtain and compare them with one obtained using PRAAT on the task of classification between normal and pathological speech. We carry out the experiments using the popular MEEI-KayPENTAX Voice Disorders database (KPdb) [START_REF]Voice disorders database[END_REF]. The results clearly show that going beyond linear method can indeed lead to much better insight and performances in this task.

The paper is organized as follows. In the next section we give a brief description of our GCI detection algorithm. Then, in section 3, we list the perturbation measures we consider in this paper. The core of the paper is in section 4 where we present the experimental results and their analysis. We then provide a conclusion in section 5.

GCI detection

In this section, we briefly recall the description of our GCI detection methodology and then , in order to make the paper selfcontained, provide the algorithm. The details can be found in [START_REF] Khanagha | Detection of glottal closure instants based on the microcanonical multiscale formalism[END_REF].

Our GCI detection algorithm is based on a novel framework called the Microcanonical Multiscale Formalism (MMF) [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. part 1: singularity analysis[END_REF]. MMF allows the study of local geometrico-statistical properties of complex signals from a multiscale perspective. It is based on precise computation of local parameters called the Singularity Exponents (SE) at every point in signal domain. When correctly defined and estimated, these exponents alone can provide valuable information about local dynamics of complex signals and has recently proven to be promising in many signal processing applications ranging from signal compression to inference and prediction in a quite diverse set of scientific disciplines such as satellite imaging [START_REF] Grazzini | Edge-preserving smoothing of high-resolution images with a partial multifractal reconstruction scheme[END_REF][START_REF] Grazzini | Multifractal Formalism for Remote Sensing: A Contribution to the Description and the Understanding of Meteorological Phenomena in Satellite Images[END_REF][START_REF] Yahia | Motion analysis in oceanographic satellite images using multiscale methods and the energy cascade[END_REF][START_REF] Yahia | High-resolution ocean dynamics from microcanonical formulations in non linear complex signal analysis[END_REF], adaptive optics [START_REF] Maji | Inferring information across scales in acquired complex signals[END_REF][START_REF] Maji | Towards multiscale reconstruction of perturbated phase from hartmann-shack acquisitions[END_REF], computer graphics [START_REF] Badri | Computer graphics effects from the framework of reconstructible systems[END_REF] and natural image processing [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF][START_REF] Turiel | The multi-fractal structure of contrast changes in natural images: from sharp edges to textures[END_REF]. In the field of speech processing, besides GCI detection [START_REF] Khanagha | Detection of glottal closure instants based on the microcanonical multiscale formalism[END_REF] , we have also successfully used MMF in phonetic segmentation [START_REF] Khanagha | A novel textindependent phonetic segmentation algorithm based on the microcanonical multiscale formalism[END_REF][START_REF]Improving text-independent phonetic segmentation based on the microcanonical multiscale formalism[END_REF][START_REF]Phonetic segmentation of speech signal using local singularity analysis[END_REF]. An important subset of SEs is called the Most Singular Manifold (MSM), which is defined as the set of points having the smallest SE values. Indeed for a given point, the smaller the value of SE is, the higher unpredictability is at this point [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. part 1: singularity analysis[END_REF]. It has been established that the critical transitions of the system occur at these points, and this fact has been used in many signal processing applications [START_REF] Yahia | Motion analysis in oceanographic satellite images using multiscale methods and the energy cascade[END_REF][START_REF] Pont | Description, modeling and forecasting of data with optimal wavelets[END_REF]. MSM constitutes the core of our algorithm for GCI detection which consists in the following steps. Let s[n] be a discrete time speech signal with a sampling frequency fs. The multi-scale integral of the following scale ri dependent functional is first defined as:

Γr i (s[n]) = |2s[n] -s[n -ri] -s[n + ri]| (1)
Then the singularity exponent h[n] at the discrete time instance n is computed as:

h[n] = I i=1 hi[n] (2)
where

hi[n] = log(Γi(s[n])) log(ri) (3)
I is the number of scales used for estimation and ri = i/fs.

We then define the regularity-drop functional DL[n] that measures the change in local average of SEs before and after any time instant n, on two adjacent windows of length TL:

DL[n] = n-1 k=n-T L h[k] - n+T L k=n h[k] (4)
The final algorithm1 is given below:

Algorithm 1 : GCI detection

n i |ni -n peak | 6: ngci ← [(n peak + nmsm)/2]
This algorithm has two free parameters: I and TL which can be tuned on development data.

Perturbation measures

The features we consider in this paper are all derived from pitch marks, that we denote T o (i) , which is the duration of the ith (estimated) pitch period. We use the default parameters for PRAAT. That is, the Period Floor and Period Ceiling parameters are set at 50Hz and 10000Hz respectively. For our feature computation algorithm, we do not use any threshold. The features we compute are Jitter, Shimmer and their classical variants.

1. Jitter is computed by:

jitter = 1 N -1 N -1 i=1 T o (i) -T o (i+1) T o ( 5
)
where N is the number of extracted pitch marks and

T o = 1 N N i=1
T o (i) .

RAP, PPQ and sPPQ are computed similarly to Jitter but with 3, 5 and 55 pitch cycles respectively (see PRAAT website www.praat.org) .

2. The T o (i) standard deviation σ is computed by:

σ = 1 N N -1 i=1 (T o -T o (i) ) 2 (6)
where:

3. The T o (i) variance vT o is computed by:

vT o = σ 1 N N i=1 1 N N j=1 T o (j) -T o (i) 2 T o (7)
4. The Skewness is computed by:

skewness = 1 N N i=1 (T o (i) -T o) 3 ( 1 N N i=1 (T o (i) -T o) 2 ) 3 (8)

Shimmer is computed by:

shimmer = 1 N -1 N -1 i=1 A (i) -A (i+1) 1 N N i=1 A (i) (9)
where: A (i) , i = 1, 2, . . . , N are the extracted peakto-peak amplitude and N is number of extracted pitch marks. APQ and sAPQ are computed similarly to Shimmer but with 5 and 55 peak-to-peak cycles.

Experimental results

The MEEI-KayPENTAX Voice Disorders database (KPdb)

The MEEI-KayPENTAX Voice Disorders database [START_REF]Voice disorders database[END_REF] was released in 1994 and has been developed by the MEEI Voice and Speech lab and the Kay Elemetrics (now KayPENTAX) Corp. The recordings consist in sustained phonation of the vowel /ah/ (53 normal and 657 pathological) and utterance of the first sentence of the rainbow passage (53 normal and 662 pathological). All normal vowels and 77 pathological vowels are sampled at 50 kHz, while the remaining 580 pathological vowels are sampled at 25 kHz. 36 of the normal rainbow sentences are sampled at 25 kHz and 17 at 10 kHz. 648 of the pathological sentences are sampled at 25 kHz, 13 at 10 kHz and one at 50 kHz. More details about KPdb can be found in [START_REF]Voice disorders database[END_REF]. In the last years, KPdb has been the most widely used dataset for research in pathological voice classification. In this paper, we use the full dataset of sustained vowels in the experiments.

Results analysis

We use a Random Forest classifier with leave-one-out cross validation method. The performance of the classifier are evaluated in terms of percentage of true positives (TP), i.e, when a disordered subject is correctly assigned to the disordered class and true negatives (TN), i.e, when a normal subject is correctly assigned to the normal class. The overall performance is the percentage of correctly classified subjects in both classes. Classification scores are presented in Table 1. The first remark is that GCI significantly outperforms PRAAT for Jitter and its variants (RAP, PPQ and sPPQ) in both TP and TN (and hence overall scores). For Shimmer, PRAAT outperforms GCI in TN but the GCI's overall score is higher. For Shimmer variants (APQ and sAPQ), GCI significantly outperforms PRAAT. This suggests that GCI is definitely is a better choice than PRAAT for the most widely used pitch perturbation features.

More importantly, the second and third order statistics show a very interesting behavior. It is worth noting that these statistics are generally not used as perturbation measures. The explanation become obvious when looking at PRAAT's scores with these features (vT 0, σ and skewness) . Indeed, they have very low scores on TN (besides for σ which has an acceptable one) which implies that they have no discriminative power. On the other hand, when computed by GCI, these features become discriminative and even significantly outperform Jitter and its variants sometimes. In particular, σ yields a (relatively) very high score and the best one. This shows that GCI captures the expected deviation from the mean of pitch periods in pathological voices. Moreover skewness, which is a measure of asymmetry in the distribution of T o (i) , yields a good score which is significantly higher than Jitter and its variants (in TN). This shows that GCI reveals an asymmetry in the distribution of pitch periods around the mean in pathological voices. To the best of our knowledge, this property has not been observed before in pathological voice analysis. This also suggests that other statistics may be worth analyzing or other new features can be conceived in order to reveal other interesting properties or to achieve higher discriminative power.

Another observation is that PRAAT is unable to process some speech signals sometimes, because of the implicit constraints it imposes in the computation of perturbation measures. On the other hand, thanks to the fact that no constraints or assumptions are imposed by our method, GCI can process any signal. Table 2 presents the number of files missed by PRAAT for each feature. The first 2 columns in "processed files" present the number of files which have been used to compute the scores: the files unprocessed by PRAAT was not considered in score calculation for a fair comparison.

Finally, we emphasize these interesting results were obtained by a "brut" version of our GCI detection algorithm. Indeed, pathological speech was not considered in the conception of the latter. Therefore, we can fairly expect that if we take this option in consideration, we could come up with a more accurate pathological GCI detection algorithm, and hence a better estimation of perturbation measures. Also, a finer analysis would allows us to choose appropriate thresholds to remove outliers (as in PRAAT) and improve classification accuracy.

Conclusion

Using our recently developed GCI detection algorithm, we evaluated some standard pitch-based perturbation measures on the task of classification between normal and pathological speech. The results first showed that Jitter/Shimmer and their variants have significantly better classification performances with our GCI than PRAAT. More importantly, the results showed that second and third order statistics (of pitch periods) are highly discriminative and informative when using GCI, while they have almost no discriminative power and convey no information when using PRAAT (at least on this database and this task). This suggests in particular that other new features can be conceived in order to reveal other interesting properties and/or to achieve higher discriminative power. This is the objective of our ongoing work, as well as the improvement of our GCI detection algorithm to make it more adapted and robust to pathological speech.

Table 1 :

 1 Classification scores on KPdb

		TP [%] TP [%] TN [%] TN [%] Overall [%] Overall [%]
		GCI	PRAAT GCI	PRAAT GCI	PRAAT
	jitter	89.94 87	47.17	45.23	86.7	83.83
	RAP	90.87 88.85	52.83	50.94	87.98	85.98
	PPQ	89.77 88.37	56.6	50.94	87.25	85.53
	sPPQ	87.18 84.03	49.06	41.5	84.1	80.58
	shimmer 92.57 89.16	52.83	58.49	89.56	86.84
	APQ	90.03 90.8	60.38	50.94	87.77	87.77
	sAPQ	89.77 88.11	54.72	47.12	86.95	84.83
	vT o	87.13 84.96	54.72	22.64	84.67	80.23
	σ	94.89 93.03	73.58	58.49	93.27	90.41
	skewness 89.77 83.41	60.38	9.43	87.53	77.79

Table 2 :

 2 Number of files processed/missed by GCI and PRAAT

		processed files		missed files
		total files total files GCI	PRAAT GCI	PRAAT
		normal	pathol	normal normal	pathol pathol
	jitter	53	646	0	0	0	11
	RAP	53	646	0	0	0	11
	PPQ	53	645	0	0	0	12
	sPPQ	53	601	0	0	0	56
	shimmer 53	646	0	0	0	11
	APQ	53	642	0	0	0	15
	sAPQ	53	606	0	0	0	51
	vT o	53	645	0	0	0	12
	σ	53	646	0	0	0	11
	skewness 53	645	0	0	0	12

A Matlab implementation of this algorithm is made publicly available in http://geostat.bordeaux.inria.fr/