Jae-Gon Lee
email: jglee@vslab.kaist.ac.kr

Moo-Kyoung Chung

Ki-Yong Ahn

Sang-Heon Lee
email: shlee@vslab.kaist.ac.kr

Chong-Min Kyung
email: kyung@vslab.kaist.ac.kr

A Prediction Packetizing Scheme for Reducing Channel Traffic in Transaction-Level Hardware/Software Co-emulation

This paper presents a scheme for efficient channel usage between simulator and accelerator where the accelerator models some RTL sub-blocks in the accelerator-based hardware/software co-simulation while the simulator runs transaction-level model of the remaining part of the whole chip being verified. With conventional simulation accelerator, evaluations of simulator and accelerator alternate at every valid simulation time, which results in poor simulation performance due to startup overhead of simulatoraccelerator channel access. The startup overhead can be reduced by merging multiple transactions on the channel into a single burst traffic. We propose a predictive packetizing scheme for reducing channel traffic by merging as many transactions into a burst traffic as possible based on 'prediction and rollback.' Under ideal condition with 100% prediction accuracy, the proposed method shows a performance gain of 1500% compared to the conventional one.

1 Experimental results obtained with iPROVE TM tested with Pentium-4 2.8 GHz with 512 Mbytes of RAM and 32-bit PCI bus running at 33 MHz.

Introduction 1.Transaction-Level Modeling

Transaction-level modeling (TLM), usually described in SystemC, is a modeling style for SoC design with its focus on the external functional behavior of each block and interblock communications without imposing excessive implementation details [START_REF] Grotker | System Design with SystemC TM[END_REF][START_REF]Functional Specification for SystemC TM 2[END_REF]. According to the modeling of time, TLM is divided into two categories: architectural TLM and micro-architectural TLM [START_REF] Muller | SystemC TM : Methodology and Applications[END_REF][START_REF] Pasricha | Transaction level modeling of SoC with Sys-temC 2.0[END_REF]. In architectural TLM, simulation time is only roughly modeled. This method is suitable for early stage prototyping of SoC. In microarchitectural TLM, the simulation time is modeled in a fully cycle-accurate manner. With micro-architectural TLM, we can verify SoC design in the early design stage with 100% cycle accuracy [START_REF] Arm | AMBA AHB Cycle Level Inteface (AHB CLI) Specification[END_REF][START_REF] Caldari | Transaction-Level Models for AMBA Bus Architecture Using SystemC 2.0[END_REF]. This paper deals with simulation acceleration of micro-architectural TLM.

The simulation speed of TLM is much faster than that of RTL simulation. It is reported that the micro-architectural transaction-level models run at least two orders of magnitude faster than RTL models; simulation speeds of at least 100 kHz for a complete system simulation are readily achievable [START_REF] Arm | ARM System-Level Modeling[END_REF][START_REF] Clouard | Experiences and Challenges of Transaction-Level Modeling with SystemC 2[END_REF][START_REF] Muller | SystemC TM : Methodology and Applications[END_REF]. It is also possible to mix transaction level models with RTL models for gradual refinements but low simulation speed of RTL blocks limits the total simulation performance. In simulation accelerators, introduced to increase the RTL simulation speed, the limited throughput of the channel between simulator and simulation accelerator often restricts the overall performance gain.

Characteristics of the Simulator-Accelerator Channel

The channel between the simulator and the accelerator is composed of layers of API (Application Program Interfaces), device driver, and physical media each with static startup overhead. When PCI-based built-in simulation accelerator is used, experimental results show that startup overhead time is as big as 12.2 usec for each channel access whereas the payload times for simulator-to-accelerator and accelerator-to-simulator transfers are 49.95 nsec/word and 75.73 nsec/word each. 1 To utilize the channel more efficiently, we need to send lots of data at a single time. But that's not possible with conventional simulation accelerators where the progress of simulator and accelerator are synchronized every simulation time. The amount of data to be sent between two simulation domains for a single simulation cycle is often too small to justify the startup overhead time of the channel; for an SoC design where each building block is interconnected via system bus, the amount of data does not exceed five words at a time. As a result, the transfer on the channel is composed of a series of short bidirectional transfers.

If we can remove one of the two types of transfers, we can merge remaining data transfers as a single transfer to minimize the number of channel accesses. This requires prediction of the states of the other simulation domain and recoveries from bad predictions. This paper deals with an optimistic simulator-accelerator channel usage based on 'prediction and rollback' to maximize the simulation speed when transaction-level models are executed by simulator while the behavior of RTL blocks is represented by accelerator.

Related Works

Parallel discrete event simulation (PDES), sometimes called distributed simulation, refers to the execution of a single discrete event simulation program on a parallel computer [START_REF] Fujimoto | Parallel discrete event simulation[END_REF]. The concept of 'prediction and rollback' was first developed in PDES to extract maximum parallelism from a given problem. In the conservative method the progress of each process is synchronized at every simulation time. On the other hand, in optimistic method, each process proceeds assuming that there are no incoming messages. When a process receives a message with a time stamp smaller than the current simulation time, the process rollbacks to a previous state and sends negative messages to negate incorrect messages sent by the process.

The 'prediction and rollback' concept was first applied to SoC design by Yoo [START_REF] Yoo | Fast Hardware-Software Coverification by Optimistic Execution of Real Processor[END_REF], where ARM prototype board models ARM processor and high level simulator models the behavior of hardware IPs. Yoo tried to minimize the channel accesses between the prototype board and simulator. But, his approach is based on high-level simulation without cycle-accurate behavior, which limits its practical use.

In [START_REF] Lee | Simulation Acceleration of Transaction-Level Models for SoC with RTL sub-blocks[END_REF], we proposed an "optimistic" simualtoraccelerator channel usage scheme based on 'prediction and rollback'. With the proposed method, one of the two verification domains, i.e, simulation domain and acceleration domain, leads the other. The leader predicts the responses of the lagger so as to remove lagger-to-leader transfers on the channel. This allows to merge multiple leader-to-lagger transfers on the channel to minimize startup overhead of channel access. The proposed method has two operating modes, i.e., Simulator Leading Accelerator (SLA) and Accelerator Leading Simulator (ALS) depending on which of the two verification do-mains leads the other. The duration of each SLA or ALS phase is named as a transition and a single transition is composed of four steps.

• Run-Ahead step (RA step) where leader proceeds predicting responses of lagger. Output of leader are not sent to lagger until RA step is over. Instead they are stored in Leader Output Buffer (LOB).

• Follow-Up step (FU step) where lagger follows up with the leader.

• Optional RollBack step (RB step) where leader rolls back to a previous state in case of prediction error detected.

• Optional Roll-Forth step (RF step) where leader runs again from the previous state to reach the progress of lagger.

Even though the proposed method can have great performance improvement compared to the conventional method, the paper focused on cases where two blocks residing in two different verification domains communicate over static interconnections between them and failed to handle cases where multiple blocks are dynamically interconnected with a bus.

Problem Definition

We applied the "optimistic" channel usage pattern to SoC model where building blocks with different abstraction levels, i.e., TL and RTL, are interconnected with a system bus. It is assumed that the system bus follows the Advanced High-performance Bus (AHB) specification. 2 The bus dynamically utilizes interconnections among bus components, i.e., bus masters and bus slaves, to use the common resource for multidirectional data transfers. The dynamic utilization includes dynamic decision of data flow direction, dynamic decision of active bus components, dynamic decision of active interconnections, etc., which complicate the application of the "optimistic" channel usage pattern to SoC verification. The complications are summarized as follows.

1. How to split a single bus model into two sub-bus models without combinatorial half loops between them?

2. How to limit number and types of signals between the two sub-bus models so that we can predict contents of at least one of the bidirectional data transfers between them?

3. Dynamic decisions on how to packetize bus signal values between the two sub-bus models.

4. Dynamic decisions among SLA, ALS and "conservative" operating modes.

We can meet the first constraint by letting each bus component to be present only in one of the two verification domains, i.e., simulation domain and acceleration domain. As most bus specifications limit communication between bus components to take place only at edges of a clock signal, there can be no combinatorial half loop if each component is present only in a single verification domain.

We can meet the second constraint by limiting the subject of data transfers between simulator and accelerator only to elements of minimal set of active bus signals. Set 4 Except for the arbitration request signals, all the elements of MSABS are related to data transfer on the bus and we call it transaction bus signals. In short, MSABS can be divided into two subsets of set of transaction bus signals and set of arbitration request signals. Among elements of set of transaction bus signals, values of address and control signals of active bus master can be deduced from their values at the start of a burst transfer on the target bus as their values either increase linearly over time or remain constant throughout a single burst transaction on the target bus. In other words, they are "predictable." Responses of active bus slave are also "predictable" as they just represent whether the active bus slave can handle bus transaction at a particular target time, which can be modeled with a simple producer-consumer model. This leaves read data and 3 Active bus master refers to a bus master that is granted for bus access and active bus slave refers to a bus slave that is accessed by the active bus master. 4 Signal name in parentheses denotes the name of the corresponding signal under AHB specification [START_REF] Arm | AMBA Specification[END_REF]. It is assumed that arbitration priority and address maps of bus slaves are statically defined. This removes output signals of arbiter and decoder from the minimal set of active bus signals whose states (values) can be deduced from arbitration request signals and address signals.

Set of transaction bus signals

Set of bus signals

Set write data whose values cannot be effectively predicted, i.e., "non-predictable." As only one of the two data signals is active at a time, we set the source of the data flow as leader and the sink of it as lagger so that we do not need to predict values of data signals.

Elements of set of arbitration request signals are "nonpredictable" as we cannot predict whether a bus master will request for bus access at a particular target time. We cannot apply the same solution we used for data signals here as there can be sources of arbitration request signals, i.e., bus masters, on either side of simulator-accelerator channel. In other words, we should be able to predict the values of bus request signals driven by bus masters residing in lagger. But as the arbitration request signals contribute only to the generation of arbitration result signal, it suffices to be able to predict arbitration result signal value. In SoC designs where large amount of data flow in bursts between building blocks, the arbitration result tends to change only occasionally and we can effectively predict its value from its previous one. Figure 1 summarizes the above explanation. When there is any signal other than bus signals interconnecting two building blocks residing in different verification domains, interrupt signal to be one of the most common examples, it should be treated the same as elements of MSABS and should be a subject of prediction, too.

The remaining problems related to dynamic use of simulator-accelerator channel are handled in Sect. 5.

Bus modeling

Figure 2 shows how a single SoC model is split into two verification domains according to abstraction level of each component. For this purpose, a single bus model should be

Operations of Channel Wrapper

The operation of CW is realized with a state diagram shown in Fig. 3. Each of the two CWs holds its own state and their combination represents each step of a transition. Operations of CW can be grouped into six paths denoted as F (roll-Forth), P (Prediction), S (Synchronization), L (Lagger), R (Report), and C (Conservative) paths in Fig. 3. Each path represents operations of a CW for a single simulation cycle. Except for the case when operation of a CW is blocked at some blocking read operations, state of the CW flows from START to END at every positive edge of a clock signal. This is called a unit cycle operation of CW. Let's assume that two CWs start running in conventional operating modes: "conservative" cycle-by-cycle synchronization. During this time, two CWs take C-path together leapfrogging each other. At some moment, one of the two CW realizes that it can predict the response of the other and takes P-path instead of C-path. This is the start of a new tran-sition. Even though leader takes P-path, it does not run in "optimistic" way (no predictions are made) when this is the first time to take P-path for a transition. This is to store the state of leader before taking "optimistic" operations for possible rollbacks in the future. Leading CW takes rb store state (denoted as P-5 in Fig. 3) to register a state store after current unit cycle operation is over. 5 After that leader takes C-path for "conservative" operation (P-6). The "optimistic" channel usage starts when leader comes back to P-path again at the very next unit cycle operation. 6 Now leader takes different path in the P-path to write output values of leader to LOB and to predicts the response of lagger. The prediction result is stored in the LOB along with output data of leader. The prediction results is reflected on the current verification domain before current unit cycle operation is over; predicted signal value is sent to the current verification domain as if it was read from the other verification domain. This continues until leader cannot predict the response of lagger. When leader cannot proceed without synchronization, leader takes S-path instead of P-path. In S-path, the contents of LOB are flushed to lagger (S-2) and leader waits for the reply from lagger in Get response state (S-3).

Until then, lagger waits for leader to write data in Read input data state (C-3). Now that LOB is flushed, lagger can get out of the blocking read operation to finish the simulation cycle and to take L-path at the next simulation cycle. Lagger checks a single prediction every time it reaches Prediction check state (L-1). If the prediction coincides with the actual response of lagger, lagger reads another leaderto-lagger data and finishes the unit cycle operation. When all leader-to-lagger data are consumed, i.e., all the predictions are correct, lagger takes R-path, and reports this to leader, which has been waiting for this response in S-3. 7 In R-path, the output of lagger is sent to lagger immediately (R-2) and lagger once again waits for the leader-to-lagger at Read input data state (R-3). And this is the end of a successful transition.

When a prediction failure is detected in Prediction check state in L-path, it is reported to leader immediately (L-5). After reporting the prediction failure, lagger waits for the 5 For simulator, the state is stored after all operations for the current simulation time is over and all the variables are stabilized. This is to save memory requirements for the state storage. For accelerator, the state is stored as soon as ACW reaches P-5 as the signal values of accelerator is stabilized as soon as clock signal toggles [START_REF] Lee | Simulation Acceleration of Transaction-Level Models for SoC with RTL sub-blocks[END_REF]. 6 If leader cannot predict the response of lagger at the next simulation cycle, the transition is over and there is no optimistic channel usage and the proposed method works just the same as conventional method with unnecessary state store overhead spent at the previous unit cycle operation. 7 Lagger can figure out that the last leader-to-lagger data is reached as the last leader-to-lagger data does not contain prediction. The last unit cycle operation of leading CW does not predict the state of lagger as it tries to read it from lagger as conventional method does. leader at L-6. Upon receiving this message, leader takes prediction failure path in S-path. First, leader stores the last response of lagger (S-5), which leader failed to predict, and requests state restore (S-6). After the unit cycle operation of leading CW is over, the state of leader is rolled back to a state P-path stored in the past. 8 Now leader takes F-path for roll-forth operations. The operations of F-path resemble those of P-path except that F-path does not write output signal values to LOB. After iterating taking F-path for the number of successful predictions, the transition is over. Operations of each path and their relations to transition steps are summarized in Tbl. 1.

Experimental Results

The performance of the proposed idea is highly sensitive to the prediction accuracy. Low prediction accuracy degrades performance by increasing not only the number of rollbacks but also the number of state restores and state stores. But the biggest degradation comes from the increased number of clock cycles to be processed by leader and channel accesses as shown in Tbl. 2. Table 2 shows the effect of prediction accuracy to the duration of time spent by each component operation under ALS operating mode. We assumed simulator speed of 1,000 kcycles/sec, accelerator speed of 10 Mcycles/sec, LOB depth of 64 and 1,000 rollback variables. T sim. and T acc. stand for the average time spent by simulator and accelerator to model the behavior of a target SoC model for a single target clock cycle, respectively. T store and T restore stands for the time spent in stor-8 The state of CW is not rolled back.

ing and restoring the state of leader. T store stands for the time spent in accessing simulator-accelerator channel. Conventional method has a simulation speed of 38.9 kcycles/sec under the same environment. The proposed method has performance gain of 16.75 when all the predictions are correct. The performance of ALS drops as the prediction accuracy drops. When it equals to 10%, the performance of ALS is about the same as that of conventional method.

Figure 4 shows the performance estimation of ALS under four different configurations with two different simulator speed and two LOB depths. As the performance of hardware-based simulation accelerator is independent of design sizes, we kept accelerator speed to be constant. The bigger the simulator performance gets, we get the more performance gain from the proposed method. LOB depth decides the maximum number of predictions and tends to accelerate the performance gain of the proposed idea when the prediction accuracy is high. On the other hand, it degrades the performance gain when the prediction accuracy is low.

The performance of SLA has similar tendencies: maximum performance gain of 3.25 and 15.34 for simulation performance of 100 kcycles/sec and 1,000 kcycles/sec, each. But it was found that SLA suffers more from low prediction accuracies. This is because relatively low operating speed of simulation domain compared to that of acceleration domain enlarges the effect of the most dominant factor of performance degradation, i.e., time spent by leader. SLA has the same simulation performance as the conventional method when the prediction accuracy is 98%[70%] assuming that the simulation performance is 100 kcycles/sec[1,000 kcycles/sec].

Prob.

1.000 0.990 0.960 0.900 0.800 0.600 0.300 0.100 T sim. 1.0e-6 1.0e-6 1.0e-6 1.0e-6 1.0e-6 1.0e-6 1.0e-6 1.0e-6

T acc.

Conclusion

Micro-architectural transaction-level modeling enabled early stage SoC verification with full cycle accuracy and fast simulation speed. Usually, transaction-level models are gradually refined to RTL models. But as the simulation

Figure 1 .

 1 Figure 1. Grouping of bus signals.

  of bus signals is a set of signals holding all the signals present in a bus specification. Set of active bus signals is a subset of set of bus signals, all of whose elements influences the operation of the bus. Signals driven by either by active bus master or active bus slave are elements of set of active bus signals. 3 Signals driven by arbiter and decoder are also elements of set of active bus signals as well as arbitration request signals driven by any of bus masters. Minimal set of active bus signals (MSABS) is a subset of set of active bus signals, values of whose elements can exclusively define the operations and states of the bus without redundant elements whose states can be deduced by combinations of states of the other elements of MSABS. Specifically, MSABS includes address (HADDR), control signals (HTRANS, HWRITE, HSIZE, HBURST, and HPROT), write data (HWDATA) of active bus master, read data (HRDATA), responses (HRESP, HREADY, and HSPLITx) of active bus slave, and arbitration request signals (HBUSREQx) of all bus masters.

Table 1 . Roll steps and CW states.

 1

	Roll Leader Lagger Description
	step	state	state
	RA	P	L,R,C Leader predicts responses of
	step	path	path	lagger to remove read trans-
				action except for the first
				time in P-path.
	FU	S	L	Lagger follows up leader
	step	path	path	until either they are synch-
				ronized or prediction error
				is found.
			R	Lagger reports that all the
				predictions were correct.
	RB	S	L	States of leader gets
	step	path	path	rolled back to a prev. state
	RF	F	L	Leader follows up the prog-
	step	path	path	ress of lagger, which is
				waiting for leader.

Table 2 . Performance of ALS.

 2

		1.0e-7 1.6e-7 2.9e-7 4.9e-7 8.1e-7 1.5e-6 2.4e-6 3.0e-6
	T store	4.69-10 7.6e-10 1.6e-9 3.3e-9 6.2e-9 1.2e-8 2.1e-8 2.7e-8
	T rest.	0 2.9e-10 1.2e-9 2.9e-9 5.7e-9 1.2e-8 2.0e-8 2.6e-8
	T ch.	4.3e-7 6.8e-7 1.5e-6 2.9e-6 5.4e-6 1.1e-5 1.8e-5 2.3e-5
	Perform.	652k	543k 363k 226k 138k 76.7k 46.1k 36.7k
	Ratio	16.75 13.97 9.33 5.80 3.56 1.91 1.19 0.94

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

AHB specification is one of the most popular and widely-used bus standards for embedded systems proposed by ARM [1]. Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

F-path P-path S-path L-path R-path C-path

speed of RTL models are too slow, the verification speed degrades as the proportion of RTL blocks increases. We can alleviate this problem with simulation accelerator, but now the throughput of simulator-accelerator channel limits the simulation performance. This gets worse as transactions on the channel are composed of series of short transfers, which suffers from static startup overhead of the channel. To minimize the effect of startup overhead and to get maximum simulation speed, we introduced the concept of 'prediction and rollback' to the synchronization between simulator and accelerator. With this method, the progress of simulator and accelerator are not synchronized at every simulaton time, but they are synchronized only when it is inevitable for cycle accurate behavior. We adopted the concept to a system bus model to get a profound performance gain when prediction accuracy is high.


Ambar A Gadkari

S Ramesh
email: ramesh@cse.iitb.ac.in

Automated Synthesis of Assertion Monitors using Visual Specifications

Automated synthesis of monitors from high-level

properties plays a significant role in assertion-based verification. We present here a methodology to synthesize assertion monitors from visual specifications given in CESC (Clocked Event Sequence Chart). CESC is a visual language designed for specifying system level interactions involving single and multiple clock domains. It has well-defined graphical and textual syntax and formal semantics based on synchronous language paradigm enabling formal analysis of specifications. In this paper we provide an overview of CESC language with few illustrative examples. The algorithm for automated synthesis of assertion monitors from CESC specifications is described. A few examples from standard bus protocols (OCP-IP and AMBA) are presented to demonstrate the application of monitor synthesis algorithm.

Introduction

Assertion-based verification [START_REF] Foster | Assertion-Based Design[END_REF][START_REF]Assertion-Based Verification[END_REF][3][START_REF]Verisity Design[END_REF] is gaining popularity in system level design due to its capability to combine formal methods with simulation based validation techniques. Precise and unambiguous specification is central to the success of system design and verification. Capturing high-level assertions using specification languages such as PSL/Sugar [START_REF]Property Specification Language Reference Manual[END_REF] or temporal logic becomes complex for interactions involving long event sequences [START_REF] Fisler | Towards Diagrammability and Efficiency in Event Sequence Languages[END_REF]. On the other hand, manual construction of assertion monitors using native languages is error-prone and does not scale well. This paper provides a novel methodology for automated synthesis of assertion monitors for complex SoC Designs. The methodology makes use of a visual specification language called CESC (Clocked Event Sequence Chart).

CESC is a visual language designed for specifying interactions between different modules in the system with single as well as multiple clock domains typically found in SoC designs. Properties or assertions related to interactions can be captured intuitively using CESC. CESC has a precisely defined abstract textual syntax. In [START_REF] Gadkari | CESC: A Visual Formalism for Specification and Verification of SoCs[END_REF], we defined CESC and illustrated its use in the context of several real-life examples including single and multiple clock domains.

In this paper, we present an automated procedure to synthesize assertion monitors from CESC. An important feature of the procedure is that the monitor synthesized consists of a number of local monitors one for each clock domain in the given input CESC specification; the monitors communicate and synchronize with each other exchanging the information about the local states using a scoreboardlike data structure. The scoreboard dynamically maintains the information about event occurrences, which is used in implementing the causality checks within the same and across different clock domains. A formal semantics has been defined for CESC based on the synchronous language paradigm [START_REF] Berry | The Esterel V5 Language Primer[END_REF][START_REF] Halbwachs | Synchronous Programming of Reactive Systems[END_REF], which forms the basis for the correctness of the synthesis procedure. The methodology to automatically synthesize monitors from CESC specifications helps in smooth integration of formal specification with simulation-based verification, thus reducing the cycle time and errors involved in the manual development of monitors.

This paper is organized as follows: Section 2 provides the background and related work. Section 3 gives a quick overview of CESC language. Section 4 describes the role of monitors in assertion-based verification and formally defines the monitor. Section 5 describes the algorithm for synthesis of monitors from CESC specifications. Section 6 includes a few examples of application of the monitor synthesis algorithm on standard bus protocols. Section 7 concludes the paper.

Background

With the increasing adoption of formal verification, many property specification languages having notations similar to temporal logics or regular expressions have emerged for hardware verification. These include Acellera's Sugar [START_REF]Property Specification Language Reference Manual[END_REF], Synopsys' OVA [START_REF]Assertion-Based Verification[END_REF] and Intel's FTL [START_REF] Armoni | The ForSpec temporal logic: A new temporal property specification language[END_REF]. These languages describe various event sequences and their temporal relationships in abstract symbolic notations unfamiliar to practicing engineers. The industrial specifications often include intuitive and visual constructs such as timing diagrams, sequence diagrams, message sequence charts, which are easier to use. For example, the specifications for protocol standards such as AMBA CLI [START_REF] Arm Ltd | AMBA AHB Cycle Level Interface Specification[END_REF] and OCP-IP [START_REF]OCP Specification version 1.0[END_REF] use notations like sequence diagrams and timing diagrams to explain the cycle accurate behavior of transactions. But these notations are often informal and ambiguous. Recently, there have been efforts to formalize MSC-like notations and use them in system design [START_REF] Bunker | Using Live Sequence Charts for Hardware Protocol Specification and Compliance Verification[END_REF][START_REF] Krüger | Distributed System Design with Message Sequence Charts[END_REF]. CESC language differs from these formalisms in many ways: It offers features such as synchronizing clocks (grid lines) suited for formally specifying the interactions in hardware systems; besides, it has support for structural constructs and multiple clocks thereby making it a good choice for the system level specification of SoCs which are Globally Asynchronous Locally Synchronous (GALS) systems.

There have been efforts on formalizing other visual notations such as timing diagrams [START_REF] Fisler | Towards Diagrammability and Efficiency in Event Sequence Languages[END_REF][START_REF] Amla | Model Checking Timing Diagrams[END_REF]. But compared to CESC, these are limited in capability for specifying the complex event behaviors such as causality relationships, repetitive event sequences and multiple clocked synchronizations that are commonly found in SoC design context. Also, the structural constructs such as sequential and parallel compositions and alternative constructs, which are useful for systematic building of specifications for large systemlevel designs, are missing in those notations.

With assertion-based verification approaches gaining wider acceptance in industry [16, 3-5], few attempts have been made for automatic construction of monitors and protocol checkers from temporal language specifications [START_REF] Geilen | On Construction of Monitors for Temporal Logic Properties[END_REF][START_REF] Abarbanel | FoCs: Automatic Generation of Simulation Checkers from Formal Specifications[END_REF]. They involve generation of automata equivalent to the specified property in temporal logic. The work described in this paper attempts to take this one step further by providing a methodology to synthesize assertion monitors directly from the visual specifications given in CESC.

Overview of CESC

Clocked Event Sequence Chart (CESC) language is designed for specification of the interaction behaviors in clocked systems. The simplest CESC is known as SCESC (Single Clocked Event Sequence Chart) and represents a finite duration event sequence or interaction scenario. The basic constructs provide representation for agents (referred to as instances pictured as vertical lines), synchronizing clocks (shown as horizontal grid lines) and the presence or absence of events. The events can be guarded by conditions given as propositions formed over system variables or events. Connecting arrows show the causality relationship between the events. Events occurring on external agents are shown on the frame of the chart and are referred to as environment events. Figure 1 shows an SCESC corresponding to the event sequence related to a typical read protocol within single clock domain. Various structural constructs are provided to enable hierarchical specification of complex interaction scenarios. Such constructs include sequential and parallel composition, loop, alternative, and implication. CESC constructs also include a special construct for asynchronous parallel composition to allow specification of interactions involving multiple clocks.

Figure 2 shows the CESC corresponding to the read protocol involving multiple clocks.

Semantics

The semantics of CESC is based on clocked traces or runs. Each run describes the valuations of conditions and events along sequence of clock ticks. For formally defining the notion of a run we provide the following two definitions.

Definition: State (s)

A state s STATES is defined as the assignment of truth-values to the propositions and the events. It is given as, s = {(f1, f2) | f1: PROP Boolean; f2: EVENTS Boolean}. The notation i is used to denote the projection onto i-th element in a tuple. Thus, 1 (s) = f1 and 2 (s) = f2 for a state s.

Definition: Run (r)

A run r RUNS is defined as r: N STATES, a mapping that associates a sequence of clock ticks represented as integers with a sequence of states. r(n) denotes a state at n-th clock tick. Thus, 1 (r(n)) denotes f1 and 2 (r(n)) = f2 at clock tick n. We further denote the following: p PROP (the set of all propositions),

r(n) p if 1 (r(n))(p) = True; r(n) p if 1 (r(n))(p) = False. Also, e
EVENTS (the set of all events),

r(n) e if 2 (r(n))(e) = True; r(n) e if 2 (r(n))(e) = False;
The semantic domain consists of a set of all possible runs. Intuitively, such a set represents an unconstrained set of behaviors that any system can exhibit. The specifications given in form of CESC constrain the permitted behaviors. Each SCESC chart specifies the constraints on the event occurrences within a finite interval on any run. The length of this interval is determined by number of grid lines denoted within a given chart. Thus, for CLK = {1,..., n}, given SCESC specifies the event occurrence pattern over 'n' clock ticks on any run. Intuitively, it can be seen (Figure 3) that for every run associated with an SCESC there is a finite interval in which the events occur according to the ordering specified by the SCESC. It may be noted that since SCESC does not provide any absolute notion of clock, the starting point of this interval is arbitrary.

Assertion Monitors

A typical system-on-chip (SoC) verification flow is shown in Figure 4 (exclude the grey boxes). The verification plan consisting of different scenarios specified as CESCs is easier to develop and review and can be formally analyzed for specification inconsistencies. The automated synthesis of monitors from CESC eliminates the step of manually developing checkers and monitors, thus reducing the cycle time and making the process less prone to human errors.

Informal specification documents

Stimuli

The monitor corresponding to a CESC is a finite state machine that can detect the traces (runs), which exhibit the event sequence behavior specified by the given CESC. The monitor automaton is an extension of the string-matching automaton described in [START_REF] Cormen | Introduction to Algorithms[END_REF]. The monitor automaton operates on clocked event traces. Each element of the input trace is a set valuations of events and propositions given as {(f1, f2) | f1: PROP Boolean; f2: EVENTS Boolean}. The monitor automaton uses a dynamic 'scoreboard' for storing the information regarding the event occurrences, which is helpful in implementing the checks related to causality relationships between events during a run.

Definition: Monitor

A monitor is a finite automata defined as a 5-tuple Q, , , s 0 , s f where, Q is a set of states. s 0 and s f in Q represent the initial and final state respectively. is the finite input alphabet consisting of events (EVENTS) and propositional symbols (PROP). The transition function is a mapping from Q EXP ACT to Q, where EXP is the set of logical expressions formed over EVENTS and PROP using logical connectives ' ', ' ' and ' ' with their standard meaning. ACT is a set of actions given as {Add_evt(), Del_evt(), Null}, which can be performed on the 'scoreboard'.

Following the synchronous model of systems [START_REF] Halbwachs | Synchronous Programming of Reactive Systems[END_REF], the transitions in a monitor are instantaneous and a single clock tick separates two successive transitions. The monitor automaton begins in initial state s 0 and reads one element of input trace in a clock step. If the monitor is in state 's' and reads element 'e' from the input trace, it takes the transition labeled 'exp/act', provided the logical expression 'exp' evaluates to True for the valuations of EVENTS and PROP given by element 'e'. The action 'act' is performed on scoreboard while taking this transition. A sequence of transitions that takes the monitor from initial state to final state is considered as the 'accepting' run and the input trace corresponding to that run represents the finite word in the language of the monitor.

An example of a monitor corresponding to an SCESC is shown in Figure 5.

Synthesis of Monitors

We have developed the algorithm for automated synthesis of monitors from CESC specifications. The algorithm constructs localized monitors for every SCESC, which are then combined using various composition operations. The monitors in each clock domain communicate and synchronize with the monitors in different domains exchanging the information about the local states using the dynamic scoreboard.

Translation Algorithm -Tr

The algorithm for construction of monitor from SCESC is described below. The major subroutines used in main algorithm are described below.

Pattern P of events is an array where each element of the array is a logical expression formed over EVENTS and PROP. The pattern corresponding to given SCESC is extracted using following subroutine.

extract_pattern input: SCESC 'C' output: Array 'P' begin

For clock tick i= 0 to n in C do begin Expression 'exp' corresponding to set of events on each grid line in C is obtained as follows:

Event 'e' translates to exp = (e) Event 'p:e' translates to exp = ( p e) Multiple events 'e1 ... ek' translate to exp = (e1 ..

. ek) P[i] = exp end return P end

The subroutine to compute the transition function for a given pattern P is described below. For describing the relation 'suffix_of' used in the above routine we first define a notion of 'matching' for the pattern. P k represents the sub-pattern of P consisting of first n elements, also referred to as 'prefix' of P of length k. Similarly, T s represents the sub-trace of input trace, after reading which the monitor is in state s. An element of pattern e p is said to be 'matched' by an element e T of input trace if the logical expression e p evaluates to True for the valuation of EVENTS and PROP given by e T . A prefix P k of pattern P 'matches' with a suffix of input trace T s e (written as P k suffix_of T s e), provided they have same length and there exists an element-by-element 'matching' between their respective elements.

The subroutine add_causality_check(e x , e y ) is implemented as follows: For every transition that depends on the occurrence of event e x , an action 'Add_evt(e x )' is associated with the transition. For every transition that depends on occurrence of event e y , an additional guard 'Chk_evt(e x )' is included along with the matching of corresponding elements from the pattern. For all the backward transitions (i.e. those taking monitor from higher numbered state to lower numbered state) all the Add_evt actions appearing on the forward path between these two states are reversed by including Del_evt actions for all the events. The example given in Figure 5 illustrates monitor automata corresponding to the given SCESC with causality arrow.

It can be shown that the finite words accepted by the monitor automaton correspond to the finite prefixes of the valid runs of the system meeting the specification given by corresponding CESC. Result: Let 'M' be the monitor obtained by applying the above-described algorithm 'Tr' on SCESC 'C'. The set of runs (sequence of states) corresponding to C is given as [[C]]. Also, the language of monitor automata M is given as L(M). Then we can show that:

[[C]] = * L(M)
.

Case Study Examples

The effectiveness of the CESC based specification and monitor synthesis approach has been studied by applying it on two standard bus protocols used widely in industry for SoC design, namely, the Open Core Protocol (OCP) [START_REF]OCP Specification version 1.0[END_REF] and Advanced Microcontroller Bus Architecture (AMBA) [START_REF] Arm Ltd | AMBA AHB Cycle Level Interface Specification[END_REF] bus protocols. We present here a few illustrative examples. Example 1. Figure 6 shows the SCESC corresponding to simple read transaction (refer to page 44 of [START_REF]OCP Specification version 1.0[END_REF]) in the OCP; the corresponding monitor, which can detect the depicted scenario is also included. Example 3. We include here an example of application of CESC for specification of interactions in AMBA bus protocol. A master and bus transaction sequence (refer to page 23 [START_REF] Arm Ltd | AMBA AHB Cycle Level Interface Specification[END_REF]) is shown using SCESC with corresponding monitor in Figure 8.

Conclusion

We have presented here a methodology to synthesize assertion monitors from high-level visual specifications. This effectiveness of the methodology has been tested on different parts of specifications of industrial protocols like OCP and AMBA and the initial results are very encouraging. The high level constructs of CESC were found to be very useful in structuring complex specifications and arriving at monitors in a compositional way. As part of the future work, we plan to use the synthesized monitors for checking the implementations of these protocols.

Figure 1 .

 1 Figure 1. Typical read protocol (single clocked)

Figure 2 .

 2 Figure 2. Typical read protocol (multi-clocked)

Figure 3 .

 3 Figure 3. Semantic mapping for SCESCFor defining the semantics of multi-clocked CESCs a global run is defined over a global clock, which is obtained as a union of clock ticks contributed by all the component clocks in the system.

Figure 4 .

 4 Figure 4. SoC verification flow Such verification flows often involve multiple reviews and iterations during the development of verification plan and implementation of monitors / protocol checkers making the entire process timeconsuming. Also the manual development of monitors introduces is prone to large number of errors. The CESC based flow suggested in [7], (refer to Figure 4 by including the grey boxes), helps in multiple ways. The verification plan consisting of different scenarios specified as CESCs is easier to develop and review and can be formally analyzed for specification

Figure 5 .

 5 Figure 5. SCESC and corresponding monitor

  the number of clock ticks in C. Q = {0, ..., n} /* set of states in M */ = EVENTS PROP /* input alphabet */ s 0 = itial and final states */ 0, s f = n /* in P = extract_pattern(C) /* events pattern */ = compute_transition_func(P, ) /* transition function */ For every pair (e x , e y ) in C connected by a causality arrow do call add_causality_check(e x , e y ) return M end /* end of main routine */

Figure 6 . OCP simple read operation Example 2 .

 62 Figure 6. OCP simple read operation Example 2.Figure 7 shows a monitor corresponding to a pipelined burst read operation (refer page 49 [12]) in OCP. Various events and related scoreboard actions are shown. (In Figure 7 NOT of action Add_evt indicates Del_evt).

Figure 7 .

 7 Figure 7. OCP burst operation

Figure 8 .

 8 Figure 8. AMBA AHB CLI transaction

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Acknowledgements

We would like to thank Dr. Rubin A. Parekhji for discussions and valuable suggestions. We thank the Heads of Computer Science departments in IIT Bombay and IISc Bangalore for their support.


Martin Brossard

Jean-Philippe Condomines

Silvère Bonnabel

Tightly Coupled Navigation and Wind Estimation for Mini UAVs

A novel nonlinear filter for the state-space model based on aided Inertial Navigation System, 1D Pitot static tube, angle-of-attack and angle-of-sideslip sensors is proposed. The solution, based on Multiplicative Kalman filtering, estimates in real-time the orientation, the velocity and the position of an Unmanned Aerial Vehicle along with sensor bias and 3D wind components. This paper describes and justifies the designed tightly coupled estimation scheme with both theoretical and experimental considerations. We validate then the whole approach for a mini UAV controlled through the well-known Paparazzi autopilot system which we equip with a set of low-cost sensors (accelerometers, gyros, GPS, magnetometer, barometer, 1D Pitot static tube and angular sensors), by successfully comparing the estimates obtained from real flight data with the 3D wind ground truth provided from a 60-m weather measurement tower.

Nomenclature z

state vector q unit quaternion R q direction cosine matrix parametrized by quaternion q v UAV velocity in inertial frame x UAV position in inertial frame

I. Introduction

The knowledge of the instantaneous wind velocity constitutes a precious information for Unmanned Aircraft Vehicles (UAVs), e.g., for atmospheric energy harvesting, wind disturbance rejection, target tracking, [START_REF] Al-Sabban | Wind-Energy Based Path Planning for Unmanned Aerial Vehicles Using Markov Decision Processes[END_REF][START_REF] Ragi | Dynamic UAV Path Planning for Multitarget Tracking[END_REF][START_REF] Rysdyk | Unmanned Aerial Vehicle Path Following for Target Observation in Wind[END_REF][START_REF] Ware | An Analysis of Wind Field Rstimation and Exploitation for Quadrotor Flight in the Urban Canopy Layer[END_REF][START_REF] Lawrance | Long Endurance Autonomous Flight for Unmanned Aerial Vehicles[END_REF] since wind significantly impacts flight performances, particularly for small and light vehicles. RIt may be useful to utilize fleets of UAVs to map the wind, e.g., for meteorological purposes, gas mapping, [START_REF] Reymann | Adaptive Sampling of Cumulus Clouds with UAVs[END_REF][START_REF] Elston | Overview of Small Fixed-Wing Unmanned Aircraft for Meteorological Sampling[END_REF][START_REF] Mayer | A 'No-Flow-Sensor' Wind Estimation Algorithm for Unmanned Aerial Systems[END_REF][START_REF] Neumann | Autonomous Gas-Sensitive Microdrone: Wind Vector Estimation and Gas Distribution Mapping[END_REF] with several advantages over manned flight as low cost or simplicity of deployment. However, such mapping procedures usually assume the local wind estimation, i.e., the real-time estimation at the current location of the UAV, as an output of a yet established black-box, whereas it strongly a↵ects the mapping quality.

To estimate the local wind onboard an aircraft, one method consists in comparing the aircraft dynamical model and the wind velocity. [START_REF] Lee | Estimation of Maneuvering Aircraft States and Time-Varying Wind with Turbulence[END_REF] Such techniques require accurately knowing the aircraft aerodynamic which is a clear limitation for mini UAVs, especially when the assumed steady state conditions are violated. For these reasons, we focus on the other main approach based on the wind triangle relationship which relates the airspeed, the ground speed and the wind speed. [START_REF] Beard | Small Unmanned Aircraft: Theory and Practice[END_REF][START_REF] Klein | Aircraft System Identification: Theory and Practice[END_REF][START_REF] Langelaan | Wind Field Estimation for Small Unmanned Aerial Vehicles[END_REF] With this approach, Langelaan et al. investigate direct computations of the wind through simulation. [START_REF] Langelaan | Wind Field Estimation for Small Unmanned Aerial Vehicles[END_REF] Van den Kroonenberg et al. [START_REF] Van Den Kroonenberg | Measuring the Wind Vector Using the Autonomous Mini Aerial Vehicle MEAV[END_REF] assume the horizontal wind speed as constant and no vertical wind during one flight loop and proposes then an experimental estimator through Kalman filtering. Similar estimators with no "flow sensor" have been tested, [START_REF] Mayer | A 'No-Flow-Sensor' Wind Estimation Algorithm for Unmanned Aerial Systems[END_REF][START_REF] Neumann | Autonomous Gas-Sensitive Microdrone: Wind Vector Estimation and Gas Distribution Mapping[END_REF] but an accurate 3D wind estimation need additional sensors, 15 e.g., Pitot tube, pressure strip or vanes. [START_REF] Elston | Overview of Small Fixed-Wing Unmanned Aircraft for Meteorological Sampling[END_REF] With an added Pitot tube, Cho et al. [START_REF] Cho | Wind Estimation and Airspeed Calibration Using a UAV with a Single-Antenna GPS Receiver and Pitot Tube[END_REF] design a relatively simple nonlinear Kalman filter for horizontal wind estimation. Recent works [START_REF] Rhudy | UAV Attitude, Heading, and Wind Estimation Using GPS/INS and an Air Data System[END_REF][START_REF] Condomines | Experimental Wind Field Estimation and Aircraft Identification[END_REF][START_REF] Johansen | On Estimation of Wind Velocity, Angleof-Attack and Sideslip Angle of Small UAVs Using Standard Sensors[END_REF] expand upon this concept to include vertical wind, incorporate Angle-of-Attack (AoA) and Angle-of-Sideslip (AoS) which both relate the aircraft body frame to the wind frame, and then uses Kalman filters to estimate attitude, ground speed and wind components. Such techniques require the capability to accurately determining attitude and speeds from various cheap sensors such as Global Positioning System (GPS), Inertial Measurement Unit (IMU) and magnetometer which is related to the UAV navigation problem, a well developed field [START_REF] Martin | Generalized Multiplicative Extended Kalman Filter for Aided Attitude and Heading Reference System[END_REF][START_REF] Trawny | Indirect Kalman Filter for 3D Attitude Estimation[END_REF][START_REF]Invariant Extended Kalman Filter: Theory and Application to a Velocity-Aided Attitude Estimation Problem[END_REF][START_REF] Markley | Attitude Error Representations for Kalman Filtering[END_REF][START_REF] Crassidis | Unscented filtering for spacecraft attitude estimation[END_REF] in which, however, most existing works implicitly assume that the aircraft is always pointing in the direction of its total velocity, i.e., its AoA and AoS are negligible. Finally, validating such wind estimators in experimental flight is challenging due to the di culty in accurately measuring the wind. One can yet use ground weather station data, radiosondes and piloted balloons, or SOnic Detection and Ranging, [START_REF] Mayer | A 'No-Flow-Sensor' Wind Estimation Algorithm for Unmanned Aerial Systems[END_REF][START_REF] Van Den Kroonenberg | Measuring the Wind Vector Using the Autonomous Mini Aerial Vehicle MEAV[END_REF][START_REF] Rhudy | Onboard Wind Velocity Estimation Comparison for Unmanned Aircraft Systems[END_REF][START_REF] Cho | Wind Estimation and Airspeed Calibration Using a UAV with a Single-Antenna GPS Receiver and Pitot Tube[END_REF][START_REF] Rhudy | UAV Attitude, Heading, and Wind Estimation Using GPS/INS and an Air Data System[END_REF] albeit only an approximation to the local wind at the aircraft.

Our approach assumes the low-cost sensor suite for navigation, i.e., GPS, IMU, magnometer and barometer, with additional 1D Pitot static tube and angle vanes. we mounted on our experimental drone. We opt for a sensor fusion scheme where navigation and wind estimation are tightly coupled. The contributions of the present paper are as follows • As explained in the sequel, there are possible choices for the state variables. The one proposed here has not been much used yet for tightly coupled navigation and wind estimation.

• The fusion is performed by an (quaternion based Multiplicative) extended Kalman filter (M-EKF). Part of the initial nonlinear problem is linearized resorting to the virtual measurements technique, [START_REF] Guillaume | Observable Systems Transformable into Implicit A ne Forms[END_REF][START_REF] Tan | Analytical SLAM Without Linearization[END_REF] to handle our specific AoA and AoS measurements.

• The wind model is based on experimentally collected measurements of the wind.

• Monte-Carlo simulations using the software Paparazzi illustrate the good performances of the filter, and shows it competes with wind estimators of previous literature.

• We present validation flight experimental results, where a 60-m measurement tower provides 3D wind "ground truth".

The rest of the paper is organized as follows. Section II defines and describes the considered model followed by the di↵erent approaches to perform filtering. In Section III, we design a MEKF-Wind filter for tightly coupled navigation and wind estimation. In Section IV parameter tuning is discussed, and experimentally obtained wind statistics support the chosen wind stochastic model and allow determining its noise amplitude. Section V illustrates the filter performances on realistic simulations. Section VI presents then experimental results, compared to data provided from a 60-m measurement tower, and we finally draw our conclusions in Section VII. Lower case non bold, lower case bold and upper case bold letters refer to, respectively, scalars, vectors or quaternions, and matrices.

II. Quantities of Interest, Measurements, and Estimation Schemes

This section presents the considered estimation problem, describes the various measurements available on our UAV, and then discusses various estimation schemes.

A. Dynamical Model

Taking into account the wind velocity along with UAV navigation states, we define the dynamic model as [START_REF] Klein | Aircraft System Identification: Theory and Practice[END_REF][START_REF] Rhudy | UAV Attitude, Heading, and Wind Estimation Using GPS/INS and an Air Data System[END_REF]

q = 1 2 q ⇤ (! ! b + n ! ) , (1) v
= q ⇤ (a a b + n a ) ⇤ q 1 + g, (2) ẋ = v, (3)
!b = ! b /⌧ ! + n ! b , (4) ȧb = a b /⌧ a + n a b , (5) ẇ = n w , ( 6
)
where q is the unit quaternion that encodes the rotation of the body w.r.t. the inertial frame, v and x are, respectively, the UAV ground speed and position relative to the inertial frame, ! b and a b are IMU bias, w is the wind speed, inputs ! and a correspond, respectively, to biased noisy angular rate and body axis acceleration which are provided from an IMU, and g denotes the gravity field vector. We benefit here from the common quaternional form, whose basis properties, relations and advantages w.r.t. other representations are given in Appendix A. ( 1)-( 5) are the standard IMU-based navigation dynamic where we suppose the biases follow a first order Gauss-Markov stochastic process in ( 4)-( 5), with time constant ⌧ ! and ⌧ a , see. [START_REF] El-Sheimy | Analysis and Modeling of Inertial Sensors Using Allan Variance[END_REF][START_REF] Allan | Statistics of Time and Frequency Data Analysis[END_REF] our constant wind model model ( 6) is justified in Section IV. We finally define

n = h n T ! n T a n T ! b n T a b n T w i T , (7)
that reflects the thermo-mechanical white and flicker noises on low-cost IMU direct measurements, the noise on low-cost IMU bias dynamics due to temperature e↵ects and calibration errors, and the noise on wind dynamic.

B. Measurements

The UAV disposes of a large panel of low-cost sensors, giving our noisy measurements that we stack in

y = h y T v y T x y T m y b y p y ↵ y i T , (8)
where y v and y x correspond, respectively, to the GPS velocity and position measurements, y m comes from a magnetometer and y b from a barometer. An 1D Pitot static tube obtains the measurement y p of the airspeed in the forward u-direction since the tube is mounted on the UAV along the longitudinal axis, see Figure 1. The vanes (or pressure strips) provide measurements of AoA ↵ and AoS which are represented in Figure 2 and calculated, respectively, as 12

↵ = tan 1 ✓ w a u a ◆ , (9)
= sin 1 v a p u 2 a + v 2 a ! , (10)
where

v a = h u a v a w a i T
denotes the airspeed and with measured parallel to the aircraft w-axis and thus does not contain a vertical airspeed component. The AoA represents the positive angle between the UAV wings and the airspeed, which is needed to generate lift, when the AoS relates to the rotation of the aircraft centerline from the airspeed, which decreases as the lateral acceleration increases during cornering, whereas in steady and level flight, airspeed v a is commonly aligned with the u-axis, meaning is null. The equations

AoA vane

AoS vane Pitot tube The wind triangle relationship between ground speed v, airspeed va and wind w, projected, respectively, onto the vertical (at left) and horizontal (at right) planes, making in evidence the Angle-of-Attack ↵ and the Angle-of-Sideslip . The direction of the ground speed is specified using the course , which represents the angle between the projection of v onto the xy-plane and north (i.e., the direction of travel relative to the Earth's surface) and the flight path which is the angle between the xy-plane and v. The pitch ✓ moves the UAV nose up or down and heading (yaw) is the direction that the UAV is pointed. In the absence of wind, the di↵erence between course and heading, the so-called crab-angle, is the Angle-of-Sideslip.

then read

y v = v + n v , (11)
y x = x + n x , (12)
y m = q 1 ⇤ b ⇤ q + n m (13)
y b = I z x + n b , (14)
y p = I u v a + n p , (15)
y ↵ = ↵ + n ↵ , (16)
y = + n , ( 17
)
where b denotes the magnetic field vector, the line vectors

I z = h 0 0 1 i and I u = h 1 0 0 i select, respectively,
the third and first component of a vector and therefore extract the altitude in (14) and airspeed in the forward u-direction in (15). We finally stack all the sensor noises in

n y = h n T v n T x n T m n b n p n ↵ n i T . ( 18
)
These measurements ( 11)-( 17) arrive in discrete time at di↵erent sampling rates, which are usually much slower than the IMU measurement rate. The proposed filter in Section III matches particularly these asynchronous measurements since each of them can be treated independently or combined, and a loss of on of them (e.g., GPS) is straightforwardly handled by the filter. [START_REF] Beard | Small Unmanned Aircraft: Theory and Practice[END_REF]

C. Estimation Schemes

There are various ways to tackle dynamical model ( 1)-( 6) with measurement (8). Indeed, beyond the choice of a nonlinear filter (many possibilities exist such as [START_REF] Martin | Generalized Multiplicative Extended Kalman Filter for Aided Attitude and Heading Reference System[END_REF][START_REF] Trawny | Indirect Kalman Filter for 3D Attitude Estimation[END_REF][START_REF]Invariant Extended Kalman Filter: Theory and Application to a Velocity-Aided Attitude Estimation Problem[END_REF][START_REF] Markley | Attitude Error Representations for Kalman Filtering[END_REF][START_REF] Crassidis | Unscented filtering for spacecraft attitude estimation[END_REF] ), there are important choices to be made. This is mainly due to the existence of algebraic equations that relate some variables, and particularly the so-called wind triangle relationship [START_REF] Beard | Small Unmanned Aircraft: Theory and Practice[END_REF][START_REF] Klein | Aircraft System Identification: Theory and Practice[END_REF]

v = q ⇤ v a ⇤ q 1 + w, (19)
such that we should explicitly select two variables among v, v a and w, while the remaining variable are implicitly defined. Let us consider in the sequel the two natural couples (v, w) and (v a , w) even if considering (v, v a ) is also theoretically valid. With one couple selected, we can perform the wind estimation in two manners: treating the navigation problem first, which leads to a loosely coupled scheme; or in a tightly coupled approach, i.e., opting for a scheme that incorporates both the navigation and wind estimations. Thus, three main alternative approaches emerge, and we detail them in the rest of this section.

Loosely Coupled Navigation and Wind Estimation

Navigation estimator

Wind estimator

!, a y v , y x , y m , y b y p , y ↵ , y x, !b , âb q, v ŵ Figure 3
. Loosely coupled approach. A first navigation estimator, typically an aided Inertial Navigation System, obtains estimates of navigation states. The estimated quaternion q and ground speed v are then used in a wind estimator building upon equation (19) to obtain the estimated wind ŵ.

This cascaded scheme consists of starting from the navigation estimator with dynamic (1)-( 5) and measurements ( 11)-( 14), to then launch a second estimator to obtain wind estimates by using the wind dynamic (6), the wind triangle (19), and the measurements ( 15)-( 17), as depicted in Figure 3. The second estimator takes the required navigation estimates as input random variables with known mean and covariance to obtain ŵ through, e.g., Kalman filtering. [START_REF] Cho | Wind Estimation and Airspeed Calibration Using a UAV with a Single-Antenna GPS Receiver and Pitot Tube[END_REF][START_REF] Johansen | On Estimation of Wind Velocity, Angleof-Attack and Sideslip Angle of Small UAVs Using Standard Sensors[END_REF][START_REF] Tian | UAV Flight Test Evaluation of Fusion Algorithms for Estimation of Angle of Attack and Sideslip Angle[END_REF] The wind estimation depends strongly of the accuracy of the navigation estimates but have no impact on the navigation itself.

Tightly Coupled Navigation and Wind Estimation with v and w as State Variables

Navigation and wind estimator

!, a y v , y x , y m , y b y p , y ↵ , y x, !b , âb q, v ŵ Figure 4
. Tightly coupled approach. The only estimator obtains both navigation state and wind estimates, from all the available measurements, introducing thus correlation between estimated variables.

It is also possible to couple the whole model ( 1)-( 6) when considered the whole measurement (8). The main interest of such approach compared to a cascaded scheme is the introduction of correlation between the wind estimates and navigation estimates, such that the airspeed related measurements ( 15)-( 17) serve to improve the estimation of the navigation variables. We follow this approach in Section III, that are not been yet considered to the best of our knowledge.

Tightly Coupled Navigation and Wind Estimation with v

a and w as State Variables The ideas are similar to the previous approach by substituting the airspeed v a in place of the ground speed v in the state, an approach yet treated via Unscented Kalman Filter (UKF). [START_REF] Rhudy | Onboard Wind Velocity Estimation Comparison for Unmanned Aircraft Systems[END_REF][START_REF] Rhudy | UAV Attitude, Heading, and Wind Estimation Using GPS/INS and an Air Data System[END_REF][START_REF] Condomines | Experimental Wind Field Estimation and Aircraft Identification[END_REF] We stress that for a same type of filter, the results between this approach and the above di↵er since: i) linearity is introduced in the 1D Pitot measurement (15) in place of GPS speed (11) which is rewritten as

y v = q ⇤ v a ⇤ q 1 + w + n v ; (20)
and ii) error definitions di↵er strongly with the previous approach. Indeed, if we define v a and w errors as Gaussian, the error on v is not Gaussian and conversely.

III. Proposed MEKF-Wind Filter

This section describe the proposed filter, which introduces correlation between the UAV state variables and wind estimates, via a tightly coupled scheme, see Section II.C.

A. Retained Final Model

We include both the ground speed and the wind in our state along with attitude, position a , bias and wind, leading to the definitions of the state vector as

z = h q T v T x T ! T b a T b w T i T , ( 21
)
and state dynamic given in ( 1)- (6). We assume all the measurements (11)- (17) are available in discrete time and rewrite the 1D Pitot measurement (15) in function of the state variables as

y p = I u q 1 ⇤ (v w) ⇤ q + n p . ( 22
)
We choose to consider ground speed rather than airspeed since: i) it makes the retained model as an augmented model of the initial navigation problem, from which many works have been yet consecrated and a The position itself is unneeded for the local wind estimation but requiring if we search to use the local estimation in a mapping scheme. [START_REF] Reymann | Adaptive Sampling of Cumulus Clouds with UAVs[END_REF] thus we can inspire of them and easily observe the consequences of our wind estimation on the navigation (not considering airspeed related measurements ( 15)-( 17) lead directly to the original navigation problem); and ii) it let the GPS speed measurement linear and not attitude dependent in contrast with (20), which preserves a correct correction when attitude estimates are inaccurate. One could argue that the Pitot measurement (22) becomes nonlinear but: ( 22) is a scalar measurement whereas GPS gives three speed measurements (11); and we preserve linearity on the most confident measurement since in practice the noise on Pitot tube is much higher than on the GPS.

B. Transformation of the Angle-of-Attack and Angle-of-Sideslip Measurements

This section represents a major originality of the proposed filter. Indeed, the raw angle measurements ( 16)-( 17) are strongly non linear w.r.t. the state, that can lead to degraded performances. A naive approach to treat these measurements consists of considering the tangent of y ↵ and sinus of y , but it preserves a strong non-linearity. Inspired from works which convert range and angles in/from Cartesian coordinates for chemical processes estimation, [START_REF] Guillaume | Observable Systems Transformable into Implicit A ne Forms[END_REF] robotics, ? or radar, [START_REF] Longbin | Unbiased Converted Measurements for Tracking[END_REF] we choose to transform these measurements in the following way, thus improving the filter performances, see Section V for an illustration.

Let us consider the two modified measurements

ỹ↵ = Cv a , = C q 1 ⇤ (v w) ⇤ q , ( 23
) ỹ = v T a Bv a , = q 1 ⇤ (v w) ⇤ q T B q 1 ⇤ (v w) ⇤ q , ( 24
)
where C and B are defined, respectively, as

C = h sin (y ↵ ) 0 cos (y ↵ ) i and (25)
B = 2 6 4 sin 2 (y ) 0 0 0 cos 2 (y ) 0 0 0 0 3 7 5 . (26)
In the noise-free case, we have

ỹ↵ = sin (↵) u a cos (↵) w a , = 0, ( 27
) ỹ = sin 2 ( ) u 2 a cos 2 ( ) v 2 a , = sin 2 ( ) u 2 a 1 sin 2 ( ) v 2 a , = 0, (28)
and we treat these modified measurements in our following proposed filter.

C. Proposed Filter

We choose to derive a Multiplicative EKF (M-EKF), [START_REF] Trawny | Indirect Kalman Filter for 3D Attitude Estimation[END_REF][START_REF] Markley | Attitude Error Representations for Kalman Filtering[END_REF] called MEKF-Wind, for the considered estimation problem. The M-EKF is the standard algorithm in the inertial navigation systems industry. Albeit increasingly replaced with the UKF, it remains computationally more e cient, [START_REF] Yang | Comparison of Unscented and Extended Kalman Filters with Application in Vehicle Navigation[END_REF] and the name Multiplicative is coined as the quaternion-based error is computed using a quaternion product ⇤ instead of subtraction. Specifically, the three-dimensional error vector ✏ q obtained from the quaternion error after using a small angle approximation is maintained in the error state, resulting in a minimal representation that is safe to use in calculating the covariance. In the following, we make use a continuous-discrete representation [START_REF] Beard | Small Unmanned Aircraft: Theory and Practice[END_REF] for our implementation of the MEKF-Wind.

The filter consists of two steps: propagation and update, which are summarized in Algorithm 1 b . During prediction (lines 1 and 2), the filter produces estimates of the current state variables, along with their Algorithm 1: continuous-discrete MEKF-Wind Input: initial state ẑ, initial covariance matrix P, IMU sampling time t; uncertainties, and we move the mathematical machinery to Appendix B. In the update step (lines 3 to 6), the gain matrix K is divided into K ⇤ for the quaternion error and K + for the remaining state errors. The previous treatment of the anglular measurements results in the definition of the transformed measurement ỹ and transformed measurement function h (.), respectively, as

1 ẑ = " q z# = ẑ + ż t,
K = " K ⇤ K + # = PH T S 1 ; // gain matrix 5 ẑ+ = " (K ⇤ (ỹ h (ẑ))) ⇤ q z + K + (ỹ h (ẑ)) # ; //
ỹ = h y T v y T x y T m y b y p E {ỹ ↵ } E{ỹ } i T , (29)
h (z) = 2
where C, B are given in ( 25)-( 26) and let the observation matrix H be defined as

H = 2
. . . . . . 0 I u RT q (v ŵ) ⇥ I u RT q 0 . . . . . . I u RT q C RT q (v ŵ) ⇥ C RT q . . . . . . . . . C RT q 2v a B RT q (v ŵ) ⇥ 2v a B RT q 0 0 0 2v a B RT
with Rq the (estimated) direction cosine matrix corresponding to the rotation operator q ⇤ (.) ⇤ q 1 , (.) ⇥ the skew symmetric operator representing cross products as matrix multiplications, and where the measurement transformation introduces the bias E {ỹ ↵ } and E {ỹ }, which are null in first order approximation. Finally, Q and N encodes, respectively, the covariances of the dynamic process noise n and measurement noise n y , typically diagonal matrices, and H n is the Jacobian of ỹ w.r.t. the measurement noise defined as the diagonal matrix

H n = diag ✓ 1, 1, 1, 1, 1, @ ỹ↵ @n ↵ , @ ỹ @n ◆ , with (32)
@ ỹ↵ @n ↵ = h cos (↵) 0 sin(↵) i va , (33)
@ ỹ @n = h sin (2 ) 0 sin(2 ) i diag (v a ) va , (34)
and where diag (.) converts its vector arguments to a diagonal matrix with the vectors aligned on the main diagonal. We let the getting of bias and Jacobians of ỹ↵ and ỹ , leading to both H and H n , in Appendix C.

IV. Noise Tuning and Wind Statistics

We must provide our filter with some values about the noise covariances both for the dynamic and measurement models, a↵ecting both filter convergence and performances. Let us divide the various noises into three classes: i) navigation noise incorporating the IMU random process, IMU and navigation measurements noises, i.e, n ! b , n a b , n ! , n a , n v , n x , n m and n b ; ii) airspeed sensor noise containing the noises n p , n ↵ and n of the airspeed related measurements ( 15)-( 17); and iii) wind process noise n w . The two first classes are discussed in the following subsection, the third one is based on real data and discussed in the following subsection.

A. Sensor Noise Identification

We consider in this subsection the parameter tuning of the navigation noise and airspeed sensor noise. Classical techniques and sensor characteristics yield the noise covariances for the IMU and measurements ( 11)-( 14), which can be done simultaneously to bias time constant identification. [START_REF] El-Sheimy | Analysis and Modeling of Inertial Sensors Using Allan Variance[END_REF][START_REF] Allan | Statistics of Time and Frequency Data Analysis[END_REF] Concerning then the airspeed related measurements ( 15)-( 17), we obtain a first covariance of the static noise when the UAV is motionless. We then have the noise due to non-perfect calibration in wind tunnel and inherent noise on the instruments, that are inaccurate for small airspeed v a < 10 m/s, i.e., before and in takeo↵. More precisely, we consider only the airspeed related measurements ( 15)-( 17) when the flight is yet established, which is equivalent to put an infinite noise covariance before flight establishment, that additionally avoids trouble due to wrong initialization.

B. Wind Statistics

We concentrate now on the justification of our wind dynamical model ( 6) and then propose a simple wind noise tuning, by analyzing data provided by the Midi-Pyrénées Observatory c . The Midi-Pyrénées Observatory land in Campistrous d disposes of a 60-m measurement tower which provides wind measurements in a large flat area, see Figure 5. The wind sensors consist of sonic anemometers CSAT3 and GILL, which obtain accurately 3D wind estimates equivalent to a ground truth with high sampling rate of 10 Hz at three di↵erent altitudes: 30 m, 45 m and 60 m. We plot data acquired on 27 April 2017, give a half-hour zoom representing the wind variations for the time of a UAV flight in Figure and wind statistics about these data in, respectively, Figure 6, Figure 7 and Table 1. From these data, we observe clearly weak correlation between each component of the wind in Figure 6 and in the right part of Table 1, i.e., each wind component follows a relatively independent dynamic, leading to choosing a diagonal covariance of the wind noise process. We also notice a strong correlation of the wind speed w.r.t. the altitude for w

measurement tower UAV

x and particularly for w y . However, defining a static wind profile or a Dryden model [START_REF] Langelaan | Wind Field Estimation for Small Unmanned Aerial Vehicles[END_REF] appears challenging as the factors mapping the wind at di↵erent altitude vary, and this becomes unfeasible for the very short period time of a flight 33 as in Figure 7, therefore justifying our random process model.

Tuning of Wind Process Noise Covariance

The wind model being (at least partially) justified, the parameters are still to be tuned, i.e. the covariance of n w . Compared to the measurement tower, the UAV is moving. When the latter is slowly flying or unmoving, it has di culty to precisely estimate the airspeed and thus the wind, consequently we freeze the wind dynamic with zero noise covariance, i.e., n w = 0, such that the wind converge to a static value (which di↵ers from imposing w = 0). During flight, with n w 6 = 0, an instinctive tunning would to increase the noise as the UAV speed increases since wind is also spatially distributed. [START_REF] Tascikaraoglu | A review of Combined Approaches for Prediction of Short-Term Wind Speed and Power[END_REF][START_REF] Bianchi | The Wind and Wind Turbines[END_REF] However, this appears complicated to handle, and we prefer to have two noise covariance values: one null when the airspeed related measurements are inactive (i.e., before and during takeo↵, see Section IV.A); one non-null value otherwise, and let a numerically and more sophisticated tuning for future research. A practical consequence of such tunning is to have a wind fixed to its initial values until flight establishment (where there are no reason to insert correlation between wind and another variables in the initial P).

V. Simulation Results

A. Simulation Setting

A simulation study is conducted to test the proposed filter when an explicit truth is available. We consider the flight simulation environment proposed by the Paparazzi e autopilot system [START_REF] Brisset | The Paparazzi Solution[END_REF] to compare the three following approaches that assumes the same dynamic ( 1)-( 6): i) v a -UKF, an UKF which incorporates Euler angles and the airspeed (19) in place of the quaternion and the ground speed, and considers the direct measurement (8), which is similar to the filter of Rhudy et. al; [START_REF] Rhudy | UAV Attitude, Heading, and Wind Estimation Using GPS/INS and an Air Data System[END_REF] ii), v-UKF, a novel UKF which incorporates Euler angles, ground speed, and considers the proposed measurement (29); and iii) MEKF-Wind, which incorporates the state vector (21) and considers the transformed measurement (29) where we set bias to zero, improving thus time-e ciency without decreasing performances in practical scenarios f . f Indeed, it appears in simulation that the unbiased measurements increase performances only with unrealistic weak noise.

Parameter

! (rad/s 3 2 ) a (m/s 5 2 ) ! b (rad/s 3 2 ) a b (m/s 5 2 ) w (m/s 3 2 ) v (m/s) x (m) Value 3 • 10 -3 3 • 10 -2 3 • 10 -4 10 -2 10 -2 10 -1 1 Parameter m (T) b (m) p (m/s) ↵ (rad) (rad) ⌧ ! (rad/s) ⌧ a (m/s 2 ) Value 10 -5 1 3 • 10 -1 10 -1 10 -1 8 • 10 2 10 3
Table 2. Parameter values for the simulations, with high-moderate noise.

of 30 in xy-plan. The noises are assumed independent, zero mean and we define their standard deviation to moderate-high realistic values given in Table 2 with the time constants. We then recover data from the simulation system, whose resulting trajectory is displayed in Figure 8. We initialize the three filters with the raw estimates corresponding of the first measurements, specifically for the wind we use a direct method. [START_REF] Langelaan | Wind Field Estimation for Small Unmanned Aerial Vehicles[END_REF] We set a moderate uncertainty for the wind covariance in P and a small uncertainty for the remaining state to be coherent with the conditions before real takeo↵, i.e., being more confident for the UAV states than for the wind velocity.

B. Simulation Results

We give these noisy measurements to the three filters and plot the root square error w.r.t. the wind speed in Figure 9. We clearly observe the accuracy of MEKF-Wind and v-UKF compared to v a -UKF once flight establishment, see Section II.C for justifications. The three filters keep a non null error since they suppose wind is quick varying as in practice (i.e., a high value for n w ) and noise is set moderate-high. We then run 200 Monte-Carlo simulations for this setting to compare the filter statistic performances, which are indicated in Table 3. This confirms the superiority of the MEKF-Wind, which better estimates attitude, position and wind speed, and outperforms v-UKF particularly for the attitude, that is certainly due to its quaternional error representation. The bias errors are similar for the three filters and the Estimated-STandard Deviation of MEKF-Wind are coherent with the Root Mean Square Error, which indicates that the filter is not overconfident. An unscented version of the proposed filter, which makes use of an identical quaternional error representation, [START_REF] Crassidis | Unscented filtering for spacecraft attitude estimation[END_REF] would probably obtains similar results. However, one major advantage of the MEKF-Wind over such unscented filters is its computational time-e ciency, particularly important with the relatively large size of the space-state. [START_REF] Yang | Comparison of Unscented and Extended Kalman Filters with Application in Vehicle Navigation[END_REF] We finally stress that the initial uncertainty and noise covariance values a↵ects the convergence of the estimates for all methods, which indicates the need of tuning numerically the filter parameters. However, with di↵erent initializations and realistic covariances, the tendency and observations are identical as for the indicated results.

State component RMSE Estimated-STD MEKF v-UKF v a -UKF MEKF v-UKF v a -UKF Euler angles ( ) 2.0 5.5 6.0 2.5 1.7 1.7 v (m/s) 3.1 • 10 -2 5.1 • 10 -2 - 9.0 • 10 -2 8.8 • 10 -2 - x (m) 6.3 • 10 -2 1.0 • 10 -1 1.2 • 10 -1 5.1 • 10 -1 5.1 • 10 -1 5.1 • 10 -1 ! b (rad/s) 5.1 • 10 -3 4.9 • 10 -3 5.1 • 10 -3 1.2 • 10 -2 5.2 • 10 -2 5.1 • 10 -2 a b (m/s 2 ) 8.0 • 10 -2 8.0 • 10 -2 8.1 • 10 -2 2.0 • 10 -2 2.1 • 10 -2 2.0 • 10 -2 w (m/s) 1.8 • 10 -1 2.2 • 10 -1 2.7 • 10 -1 2.4 • 10 -1 2.4 • 10 -1 2.0 • 10 -1

VI. Application on Flight Experiment

A. Experimental Setting

An experimental study is necessary for validating such wind estimation methods since the discrepancy between real and simulation flights is particularly important for the wind behavior, which can includes, e.g., gusts or thermal ascents. [START_REF] Bianchi | The Wind and Wind Turbines[END_REF] Specific interest is to obtain real-time wind data as close as the wind experimented by the aircraft. Toward this aim, we conduct the flight experiment in the Midi-Pyrénées Observatory, see Figure 10, where we yet obtain data from the 60-m measurement tower described at Section IV. The validation data at altitudes 30 m g are plotted in Figure 11 for the flight time and as remarked in Section IV, since fitting a wind profile (e.g., a power law w.r.t. the altitude) for such very short-period time is irrelevant, data acquired here are clearly accurate compared to, e.g., a simple portable ground weather station. ENAC UAV laboratory h has developed a fixed-wind UAV so-called "Zagi", which is a mini UAV of wingspan 1.1 m for 0.5 kg visible on Figure 1. Its onboard sensor consists of an Apogee v1.00 PCB flight controller with the Paparazzi autopilot, GPS RTK and HMC5983 magnetometers. We mount a self-made 1D Pitot static tube, connected to HCLA12X5EB pressure sensor, on the nose of the aircraft along the longitudinal axis, which comports both AoA and AoS vanes, each associated with an analog-to-digital converter g Anemometers at 45 m and 60 m were breakdown during flight. We schedule new experiments once sensors are fixed and consider these new experiments flights in the final version of the paper. and the whole sensor suite is yet calibrated in wind tunnel. Data are finally logged on a micro SD card with a light Real-Time Operating System called ChibiOS.

Flight experiments are realized on 6 June 2017, starting at 36 past 11 (am), the weather is cloudy and wind is moderate-high, as confirmed by the measurement tower in Figure 11. We define the flight trajectory as ovals around the measurement tower at stable altitude and search to keep close to the measurement power, to then justifiably compare the wind estimates at the measurement tower and at the UAV locations. Giving the raw acquired measurements to the MEKF-wind, the filter replays the flight and obtains wind estimates plotted in Figure 11. As we activate the airspeed measurements once takeo↵ is e↵ected (see Section IV.B), we have null estimates at the two first minutes. We see that the wind estimates matches the wind as measured by the measurement tower in mean, while none of these estimates closely match the "ground truth". This comes as no surprise since the measurement tower is measuring wind at a di↵erent locations, and local turbulence could cause di↵erences between the wind speed measured at di↵erent points and the values of the tuning parameters impact the dynamics of the wind estimates. Another particular interesting aspect for the navigation is to compare the raw measurements about the airspeed with their estimated values, as depicted in Figure 12. We clearly observe similar global shapes between estimated and measured values, confirming the correct functioning of the filter.

B. Experimental Results

VII. Conclusion & Perspectives

A novel tightly coupled estimator has been presented for the nonlinear state-space formulation of the local wind estimation and UAV navigation problem. Based on the Mulpliticative Kalman filtering, the proposed filter make use of the standard low-cost autopilot sensor suite (GPS, magnetometer, barometer, gyros and accelerometers) with additional 1D Pitot static tube, angle-of-attack and angle-of-sideslip sensors to estimate attitude, ground speed, position, sensors bias and wind. Particularly, the filter transforms the angular measurements, improving thus its accuracy and robustness, and simulations results illustrates the superiority of the proposed filter, which is also computationally e cient for future real onboard implementation. We compare then results from experimental flight data of the proposed filter with ground truth of a 60-m measurement tower, which validate the performances of the proposed approach.

Future works concern the definition of an adaptive wind noise covariance depending on the UAV speed with specific attention to the numerical values and the direct inclusion of the proposed local wind estimation in a mapping framework.

A. Quaternions for Representing Rotations

A unit quaternion provides a convenient mathematical notation for representing orientations and rotations of objects in three dimensions. [START_REF] Kuipers | Quaternions and Rotation Sequences : A Primer with Applications to Orbits, Aerospace, and Virtual Reality[END_REF] Compared to Euler angles they are simpler to compose and avoid the problem of gimbal lock. Compared to rotation matrices, they are more compact, more numerically stable and more e cient. Indeed, for any unit quaternion

q = q r + q i = cos ✓ ✓ 2 ◆ + u sin (✓) , (35)
and for any vector v 2 R 3 the action of the operator q ⇤ (.) ⇤ q 1 on v rotates the vector argument about axis u = q i /kq i k through an angle ✓, from which we can deduce a direction cosine matrix R q and Euler angles such that

q ⇤ v ⇤ q 1 = R q v. (36)
To let a quaternion, that is parameterized by a scalar and a three dimensional vector, operates on a vector v 2 R 3 , note that v can be viewed as a pure quaternion whose real part is zero, i.e., " 0 v # . Thus, we can compactly define the following (non commutative) quaternion product ⇤ as

q = " q r q i # , p = " p r p i # , q ⇤ p = " q r p r q T i p i q r p i + p r q i + (p i ) ⇥ q i # , (37)
where the skew symmetric operator (.) ⇥ represents cross products as matrix multiplications.

B. MEKF-Wind Error and Transition Matrix Definitions

We develop here the machinery to obtain the state transition matrix A, and the Jacobian of the state transition w.r.t. the process noise A n , which are needed to the propagation step of MEKF-Wind, see Algorithm 1. According to the multiplicative error methodology [START_REF] Martin | Generalized Multiplicative Extended Kalman Filter for Aided Attitude and Heading Reference System[END_REF][START_REF] Trawny | Indirect Kalman Filter for 3D Attitude Estimation[END_REF] we choose the following attitude error ✏ q defined as q

= " 0 ✏ q # ⇤ q, (38)
where (.) symbolizes estimated variables and which is after a small angular approximation equivalent to

R q = ⇣ I + (✏ q ) ⇥ ⌘ Rq , (39)
and we then define the remaining errors as

✏ v = v v, ✏ x = x x, ✏ ! b = ! b !b , ✏ a b = a b âb and ✏ w = w ŵ. ( 40
)
We thus derive the error dynamics using the nonlinear dynamics in ( 1)-( 6), which are used to propagate then the error covariance matrix P. The full error dynamic is given as

✏q = q ⇤ ( ✏ ! b + n ! ) ⇤ q 1 , (41) ✏v
= q ⇤ (a a b ) ⇤ q 1 ⇥ ✏ q + q ⇤ (✏ a b n a ) ⇤ q 1 , (42) ✏x = ✏ v , (43) ✏! b = ✏ ! b /⌧ ! + n ! b , ( 44
) ✏a b = ✏ a b /⌧ a + n a b , (
We thus are able to obtain the continuous dynamic of the covariance matrix P as

Ṗ = APA T + A n QA T n , (49)
that in discretized in Algorithm 1 and where Q is the covariance matrix of the process noise (7).

C. Expectations & Jacobians of the Transformed Measurements

We develop here the calculus giving the expectations and Jacobians of the transformed measurements ỹ↵ = Cv a , = C q 1 ⇤ (v w) ⇤ q and (50)

ỹ = v T a Bv a , = q 1 ⇤ (v w) ⇤ q T B q 1 ⇤ (v w) ⇤ q , (51)
with C and B defined, respectively, in (25) and (26). Specifically, when the noise encoded in y ↵ and y is present, these transformed measurements can be biased and we should theoretically evaluate the expected bias E {y ↵ } and E {y } to remove it, that gives us the full alternative measurement (29). Consider firstly the measurement ỹ↵ , we start from the Kalman filter assumptions, i.e., the state and measurement noise

for z ⇠ N (ẑ, ), that we use in the sequel. For the expectation, since state error and noise are independent, we separate the expectation in noise expectation and state error expectation as

E {ỹ ↵ } = E z E n {C} q 1 ⇤ (v w) ⇤ q , = E n {C} E z q 1 ⇤ (v w) ⇤ q , ( 56
)
with after using ( 52)-( 53)

E n {C} = e 2 ↵ 2 h sin ( ↵) 0 cos (↵) i , = e 2 ↵ 2 Ĉ, (57)
E z q 1 ⇤ (v w) ⇤ q = E z n⇣ I (✏ q ) ⇥ ⌘ q 1 ⇤ (v ŵ + ✏ v ✏ w ) ⇤ q o , = q 1 ⇤ (v ŵ) ⇤ q E z n (✏ q ) ⇥ q 1 ⇤ (✏ v ✏ w ) ⇤ q o , (58)
= va p, with (59) p = q 1 ⇤ (P

qv P qw ) ⇥ ⇤ q, (60)
and where we remove the first order term in (58) since random variables have zero mean, use then the cross-covariances P qv and P qw , between ✏ q and, respectively, ✏ v and ✏ w , and the counterpart of the cross product for matrix (.) ⇥ , defined as As a consequence, the expectation is the noise-free measurement multiplied by the factor e 2 ↵ 2 plus an additional term, i.e., the expectation is given as

P
E {ỹ ↵ } = e 2 ↵ 2 Ĉ (v a + p) , = e 2 ↵ 2 ⇣ 0 + Ĉp ⌘ , = e 2 ↵ 2 Ĉp. ( 62
)
We compute then the Jacobians w.r.t. the state and the noise, respectively, as

@ ỹ↵ @z = h C RT q (v ŵ) ⇥ C RT q 0 • • • 0 C RT q i , (63)
@ ỹ↵ @n y = h 0 • • • 0 @ ỹ↵ @n↵ 0 i . ( 65
)
Similarly for ỹ but more tricky, we start from the Kalman filter assumptions, with n ⇠ N (0, ). For the expectation, since state and noise are uncorrelated, we have

E {ỹ } = E z n q 1 ⇤ (v w) ⇤ q T E n {B} q 1 ⇤ (v w) ⇤ q o . ( 66
)
Before developing this expectation with the error variables, we define ṽ = v ŵ, (67)

✏ ṽ = ✏ v ✏ w and (68)
D = RE n {B} RT , (69)
leading the expectation about noise after using of (54) and ( 55) as We can then develop

E n {B} =
E {ỹ } = E z n (ṽ + ✏ ṽ) T ⇣ I + (✏ q ) ⇥ ⌘ D ⇣ I (✏ q ) ⇥ ⌘ (ṽ + ✏ ṽ) o , ( 71
) = E z n⇣ ṽT ṽT (✏ q ) ⇥ + ✏ T ṽ ✏ T ṽ (✏ q ) ⇥ ⌘ D ⇣ ṽ + (✏ q ) ⇥ ṽ + ✏ ṽ + (✏ q ) ⇥ ✏ ṽ⌘o , (72)
and preserve only zero, second and fourth order terms since random variables have zero mean, leading to

E {ỹ } = E z ⇢
ṽT Dṽ + ṽT D (✏ q ) ⇥ ✏ ṽ ṽT (✏ q ) ⇥ D (✏ q ) ⇥ ṽ ṽT (✏ q ) ⇥ D✏ ṽ + ✏ T ṽ D (✏ q ) ⇥ ṽ ✏ T ṽ (✏ q ) ⇥ Dṽ ✏ T ṽ (✏ q ) ⇥ D (✏ q ) ⇥ ✏ ṽ .

(73)

We thus recognize the free-noise measurement ṽT Dṽ = vT a Bv a = 0,

and cancel some parts of (73) with their minus transpose, to obtain

E {ỹ } = E z n ṽT (✏ q ) ⇥ D (✏ q ) ⇥ ṽ ✏ T ṽ (✏ q ) ⇥ D (✏ q ) ⇥ ✏ ṽo , (75)
= E z n ✏ T q (ṽ) ⇥ D (ṽ) ⇥ ✏ q ✏ T ṽ (✏ q ) ⇥ D (✏ q ) ⇥ ✏ ṽo , (

= Tr (ṽ) ⇥ D (ṽ

) ⇥ P qq E z n ✏ T ṽ (✏ q ) ⇥ D (✏ q ) ⇥ ✏ ṽo , (77)
' Tr (ṽ) ⇥ D (ṽ) ⇥ P qq pT Dp,

with Tr (.) the trace operator and by considering the variance of the random variable (✏ q ) ⇥ ✏ ṽ as null (fourth order approximation).

Finally, let us compute the Jacobians w.r.t. the state and the noise, respectively, as

@ ỹ @z = h 2v T a B RT q (v ŵ) ⇥ 2v T a B RT q 0 • • • 0 2v T a B RT q i , (79)
@ ỹ @n = h sin (2y ) 0 sin(2y ) i diag (v a ) va and (80)

@ ỹ @n y = h 0 • • • 0 @ ỹ @n i . ( 81)

b

  angular rate bias a b acceleration bias in body frame w wind velocity in inertial frame v a airspeed in body frame ! biased angular velocity from gyro a biased acceleration in body frame from accelerometer ↵ Angle-of-Attack (AoA) Angle-of-Sideslip (AoS) g gravity field vector b magnetic field vector Q process noise covariance matrix N measurement noise covariance matrix

Figure 1 .Figure 2 .

 12 Figure 1. At left, the "Zagi", a fixed-wind UAV with an 1D Pitot static tube mounted along the longitudinal axis. At right, a zoom of the 1D Pitot static tube with two vanes that measure, respectively, Angle-of-Attack and Angle-of-Sideslip. The whole UAV and sensors are designed, manufactured and instrumented at the ENAC UAV laboratory.

Figure 5 .

 5 Figure 5. 60-m measurement tower, equipped with CSAT3 and GILL sonic anemometers at 30 m, 45 m and 60 m, which provide accurate 3D wind measurements at 10 Hz sampling rate.

c

  http://www.obs-mip.fr/ 1. Empirical Justification of the Wind Dynamic(6)

Figure 6 .

 6 Figure 6. Wind speed obtained by the measurement tower at 30 m, 45 m and 60 m for each wind components, as function of time, between 5 am (t = 5 h) and 10 pm, on 27 April 2017, a relatively cloudy day. The correlation of wind w.r.t. altitude appears evident, and we observe weak correlation between wind components.

Figure 7 .

 7 Figure 7. Wind speed obtained by the measurement tower at 30 m, 45 m and 60 m for each wind components, as function of time, between 1 pm and one half past 1 pm, on 27 April 2017. We clearly see that one can not correlate the measurements at di↵erent altitudes through linear regression.

Figure 8 .

 8 Figure 8. Paparazzi simulation trajectory, starting close to the origin. Flight parameters have standard values, constant wind speed is wx= 4 cos(30) m/s, wy = 4 sin(30) m/s and wz = 0.5 m/s. The trajectory represents a starting 300 m ⇥ 200 m-box survey of 150 s, with soaring, curves and straight lines. The simulation flight consists of 150 s of a starting 300 m ⇥ 200 m-box survey, where we combine standard flight parameters with an artificial constant wind of 4 m/s in xy-plan, 0.5 m/s in the z-axis and with direction e http://wiki.paparazziuav.org/.f Indeed, it appears in simulation that the unbiased measurements increase performances only with unrealistic weak noise.

Figure 9 .

 9 Figure 9. Root square error on the wind speed, as a function of time, with realistic noisy measurements. MEKF-Wind obtain best results than v-UKF and va-UKF once flight establishment.

Figure 10 .

 10 Figure 10. Flight trajectory in the Midi-Pyrénées Observatory land (France), consisting of ovals around the 60-m measurement tower, during 13 min. Particular intentions are made to obtain a trajectory close to the measurement tower.

h

  http://www.enac.fr/en/uavs-systems.

Figure 11 .

 11 Figure 11. Wind speed from the measurement tower at 30 m and wind estimates of the MEKF-Wind, as function of time. Null estimates at the two first minutes are a consequence of the tuning of noise parameters.

Figure 12 .

 12 Figure 12. At left, 1D Pitot tube measurement yp and estimated longitudinal airspeed ûa, as function of time. At right, estimated and measured Angle-of-Attack ↵ and Angle-of-Sideslip . The estimates are coherent with the measurements.

2 2

 2 n↵ ⇠ N (0, ↵ ) follow Gaussian independent distributions. Let us then derive the following formula of scalar expectations E {cos (z)} = e

Table 1 .

 1 Correlation factors between the wind speed at di↵erent altitudes (at left) and di↵erent wind components (at right), with high values w.r.t. the di↵erent altitudes and weak values w.r.t. the di↵erent components.

Table 3 .

 3 Root Mean Square Error (RMSE

b A more sophisticated integration (line 1) with exact quaternion propagation or Runge-Kunta method can also be e↵ected.[START_REF] Sosa | Quaternion Kinematics for the Error-State KF[END_REF] At the end of this state propagation, we additionally enforce to the quaternion estimate q to have a unit norm.

d The site is located in Campistrous, France, at N 43 7 0 464 00 E 0 22 0 98 00 .

Acknowledgment

The authors would like to acknowledge Gauthier Hattenberger and Xavier Paris for their expertise about mini UAVs.