Dehui Du

Ping Huang

Kaiqiang Jiang

Frédéric Mallet

Mingrui Yang

MARTE/pCCSL: Modeling and Refining Stochastic Behaviors of CPSs with Probabilistic Logical Clocks

Keywords: cyber-physical systems, MARTE, pCCSL, stochastic hybrid automata, energyaware building, statistical model checking

are networks of heterogeneous embedded systems immersed within a physical environment. Several ad-hoc frameworks and mathematical models have been studied to deal with challenging issues raised by CPSs. In this paper, we explore a more standard-based approach that relies on SysML/MARTE to capture different aspects of CPSs, including structure, behaviors, clock constraints, and non-functional properties. The novelty of our work lies in the use of logical clocks and MARTE/CCSL to drive and coordinate different models. Meanwhile, to capture stochastic behaviors of CPSs, we propose an extension of CCSL, called pCCSL, where logical clocks are adorned with stochastic properties. Possible variants are explored using Statistical Model Checking (SMC) via a transformation from the MARTE/pCCSL models into Stochastic Hybrid Automata. The whole process is illustrated through a case study of energy-aware building, in which the system is modeled by SysML/MARTE/pCCSL and different variants are explored through SMC to help expose the best alternative solutions.

Introduction

Cyber-Physical Systems (CPSs) combine digital computational systems with surrounding physical processes and can be viewed as a network of embedded systems where a (large) number of computational components are deployed within a physical environment [START_REF] Lee | Cyber physical systems: Design challenges[END_REF]. Each component collects information about and offers services to its environment (e.g., environmental monitoring, healthcare monitoring and traffic control). This information is processed either within the component, in the network or at a remote location (e.g., a base station), or in any combination of these. The prominent characteristic of CPSs is that they have to meet a multitude of quantitative constraints, e.g., timing constraints, energy consumption, memory usage, communication bandwidth, QoS, and often with uncertain environment or user behaviors. So how to model and verify CPSs still remains a challenging problem. Now CPSs are spreading in numerous applications with a large societal impact, ranging from automotive systems, manufacturing, medical devices, automated highway systems (AHS) [START_REF] Deshpande | Design and evaluation tools for automated highway systems[END_REF], air traffic control systems [START_REF] Tomlin | Hybrid control models of next generation air traffic management[END_REF], personal and medical devices [START_REF] Jiang | Modeling and verification of a dual chamber implantable pacemaker[END_REF], Smart Grids and Smart Building, etc.

In the literature we find that both a variety of engineering modeling languages (lots of them are UML/SysML/MARTE-based) and a bunch of formal models (e.g., timed automata, hybrid automata, Petri nets, synchronous languages) provide a good support for formal verification. However, the integration of industry standards with verification frameworks is still in its infancy. For instance, a classical flow consists in adorning a UML state machine with some annotations and then transforming it into a timed automata for verification. We intend to go further by combining together several models of various kinds to cover heterogeneous aspects of the systems before transforming them into a language amenable to verification. While several frameworks inspired by Ptolemy [START_REF] Ptolemaeus | System Design, Modeling, and Simulation: Using Ptolemy II. Ptolemy.org[END_REF] address the important issue of heterogeneity, most of them propose an ad-hoc environment and notation, while here we start from UML/SysML/MARTE models. The aim of our work is to facilitate the modeling CPSs with standard-based modeling language.

We consider that the UML, as a general-purpose language, provides a variety of models to cover lots of aspects of CPSs, structural aspects with structured classifiers and components, state-based models with state machines and data-flow models with activities. Because we target embedded systems, we use the MARTE profile, which appears as the best choice [START_REF] Boutekkouk | UML2.0 profiles for embedded systems and systems on a chip (socs)[END_REF] for a UML-based solution.

In particular we focus on its subprofiles covering time, allocation, non-functional properties (NFP) (like power or energy consumption) and Performance Analysis Modeling (PAM). Because CPSs combine discrete and continuous aspects, we follow the lead of other works [START_REF] Selic | Modeling and Analysis of Real-Time and Embedded Systems with UML and MARTE[END_REF] and combine MARTE with SysML. In particular, we use SysML parametrics to capture the equations that link the energy to time and power.

As in any UML-based models, the relationships among models and the consistency is of paramount importance. We claim that logical clocks [START_REF] Lamport | Time, clocks, and the ordering of events in a distributed system[END_REF], just like tagged structures [START_REF] Lee | A framework for comparing models of computation[END_REF], provide a good abstraction to link different models together. Indeed, logical clocks can be used as activation conditions [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] of different models. Clock constraints then define a coordination model to constrain the joint execution of these models. The time subprofile of MARTE extends the UML with logical clocks that can then be used to control the different interactions between the models, e.g., the relationships between a state transition and a part in a structured classifier, or the start of an action and the sampling step to integrate the energy and compute the power consumption, but also relationships between UML and SysML models.

While MARTE extends UML with the notion of logical clock, its companion language CCSL [1] offers a syntax to build clock constraints. CCSL as a declarative language helps build a specification that can be refined when new constraints from the environment, or the platform or the application become available. There may be several implementations that satisfy a CCSL specification, a classical approach [START_REF] Deantoni | Timesquare: Treat your models with logical time[END_REF] consists in defining a policy (e.g., non-deterministic choice or as soon as possible) to decide which solution to retain. Another solution, which is explored here, consists in extending CCSL with stochastic constructs and probability. Such constructs help pick one solution instead of another one by giving the likelihood that a clock ticks.

Once the SysML/MARTE model is built, we propose to analyze the resulting model through a transformation to Stochastic Hybrid Automata and to use Statistical Model Checking (SMC) [START_REF] Bohlender | A review of statistical model checking pitfalls on real-time stochastic models[END_REF][START_REF] Legay | Statistical model checking: An overview[END_REF][START_REF] Sen | Statistical model checking of black-box probabilistic systems[END_REF][START_REF] Younes | Statistical probabilistic model checking with a focus on time-bounded properties[END_REF] to explore different solutions. To illustrate our approach, we take the example of a energy-aware building and show how to explore and compare alternative solutions.

Our contributions are 1) to propose an extension of CCSL, called pCCSL, with stochastic constructs and probability to drive the exploration of alternative solutions when building clock specifications; 2) To show how pCCSL specifications can augment UML/SysML/MARTE models to link together several models; 3) To propose a structural transformation into SHA according to the semantics of pCCSL so as to perform evaluation on the alternative solutions with statistical model checking, which helps designers refine system models. The whole process is illustrated on the example of an energy-aware building.

The remainder of this paper is organized as follows. Section 2 introduces CCSL and our proposed extension, pCCSL. Section 3 introduces our case study and makes a joint use of UML, SysML and MARTE to capture different aspects of this model. Section 4 proposes some transformation rules to transform MARTE/pCCSL into SHA. Finally we position our work with respect to related works before concluding and discussing possible future extensions.

pCCSL: the Probabilistic Extensions of CCSL

Syntax of pCCSL

We first recall the basic constructs of CCSL [START_REF] André | Syntax and semantics of the Clock Constraint Specification Language (CCSL)[END_REF] and then further describe the proposed extensions (see Fig. 1). Core constructs of CCSL. A specification is made of clock relations and declarations. Relations prevent some clocks from ticking depending on configurations. Declarations are meant to declare new clocks, either to capture events of the system or to build intermediate clocks based on other ones. The two basic CCSL relations are subclocking and causality. The former one is inspired by synchronous languages and prevents one clock (the subclock) from ticking when its superclock cannot tick. The superclock is said to be coarser and the subclock is finer than the other one. The latter one is akin to the causality relation in event structures. It prevents a clock (the effect) from ticking when another clock (its cause) has not ticked. The cause is said to be faster than the effect (which is then slower). This typically represents first-in-first-out (FIFO) constructs, e.g., one cannot read a data before it has been written. The other relations (precedence, synchrony, exclusion) are derived.

When declaring a clock, one can specify (if required) how it is related to other clocks. This is done using clock definitions, e.g., filters, or expressions.

Filters allow for precisely defining a subclock based on its superclock. For instance, a Periodic-Filter makes one clock tick every p ticks of its superclock. The Select expression selects an interval (possibly not closed) of ticks to the superclock at which the subclock ticks.

Function-expression builds a new clock based on two (or more clocks). Union builds the coarsest clock that is a super clock of two other clocks, whereas Intersection builds the finest subclock of two clocks. Inf (resp. Sup) builds the slowest (resp. fastest) clock that is faster (resp. slower) than two clocks.

The Sampling takes a triggering and sampling clock and builds the fastest clock slower than the triggering clock and subclock of a sampling clock.

One can also declare a period and optionally a jitter for clock making reference to an ideal physical clock. This information related to physical time is not used in CCSL clock calculus but is a mere annotation to display results with a scale that is meaningful to the user. Proposed extensions for pCCSL. The grammar of pCCSL is shown in Fig. 1. There are three proposed extensions marked with red fonts. The first consists in adding a rate parameter to the subclock relation. It specifies the rate at which the subclock ticks compared to its superclock. It corresponds to a probability for the subclock to tick when its superclock ticks. This is an intermediate solution between saying nothing, which may imply that the subclock never ticks, and completely deciding when the subclock ticks, with a filtering expression for instance.

The second proposed extension consists in assigning a probability to a clock, to replace an expression. Rather than giving a deterministic expression that says when the clock may or cannot tick, we give the probability that the clock ticks. As we show later, the consistent integration of this probabilistic clock with the semantics of the other relations is not trivial and requires some particular caution. However, this is useful to help CCSL clock calculus picking one solution instead of another one when several solutions are possible.

The third and last extension consists in giving a parameter λ to control when a clock ticks or not. The lambda parameter makes a clock tick according to an exponential distribution. The probability that the clock ticks at the time less than x is 1 -e -λx . x here is either relative to the ticks of another clock or, if no reference clock is provided, relative to an absolute ideal physical-time clock. In the latter case, the information is a mere annotation only used to interpret the ticking of the clock according to an ideal physical time reference but is not used in the clock calculus process.

Note that these three extensions cannot add new solutions, they are just meant to reduce the set of solutions in case of several possible solutions. The probability and distribution constructs help us decide on the likelihood for a clock to tick (when allowed by the other constraints).

Semantics of pCCSL

A constraint specification spec = C, Cons , is a tuple where C is a set of clock names and Cons is a set of constraints. The semantics of each individual constraint is given by a special form of transition system called clock-labeled transition systems.

Definition 1 (Labeled Transition System). A labeled transition system (lts) is a structure S, A, -→ consisting of a set S (of elements, s, called states), a set A of labels, and a transition relation

-→ ⊆ S × A × S. s a -→ s is used to denote (s, a, s ) ∈-→.
A clock-labeled transition system is a lts where each label is a set of clocks. From each state, there are maximum 2 n outgoing transitions, where n = |C| is the number of clocks. Each transition corresponds to a particular configuration of ticking clocks. Transition systems may have an infinite number of states.

Definition 2 (Clock-Labeled Transition System). A clock-labeled transition system (clts) is a structure S, C, -→ where S, 2 C , -→ is a labelled transition system.

To capture the semantics of the proposed extension we extend the clts by adding a probability to each transition.

Definition 3 (Probabilistic CLTS).

A probabilistic clock-labeled transition system (pclts) is a clts with an extended transition relation -→ ⊆ S × 2 C × P × S, where P ⊆ Q is a real number between 0 and 1 (i.e., a probability).

For a given transition t = (s, Γ, p, s ) ∈-→, π(t) = p denotes the probability p that the transition t is fired.

For a pclts S, C, -→ , we call s • the set of all transitions whose source is s:

s • = {(s, Γ, p, s ) ∈-→}
Note that s • can never be empty since it is always possible to do nothing in CCSL, i.e., (s, ∅, p, s) is always in -→ for all s ∈ S and for some value p.

Given a clock c ∈ C, let us call s • c the set of all transitions whose source is s and such that the clock c ticks:

s • c = {(s, Γ, p, s ) ∈-→| c ∈ Γ }
For a pclts to be well-formed, it must satisfy the two following conditions:

∀s ∈ S, t∈ s • π(t) = 1 (1) ∀s ∈ S, ∀c ∈ C, t∈ s • c π(t) = p c (2)
In Eq. 2, for each clock c ∈ C, the probability p c is either manually assigned by the user with a declaration 'Clock c probability p', or derived using the rate in a subclocking relation or assigned to the default value 1/|s • | otherwise.

A 'normal' clts can be seen as a probabilistic clts where all the probabilities are assigned with default values 1/|s • | for all the states s ∈ S. Subclocking and synchrony. Let a, b ∈ C be two clocks and r ∈ Q a real number such that 0 ≤ r ≤ 1. The subclocking relation (see Fig. 2(a)), b ⊆ a rate r is defined as a pclts {s0}, {a, b}, -→ ⊆ , such that -→ ⊆ = {(s0, {}, 1 -p a , s0), (s0, {a, b}, p a * r, s0), (s0, {a}, p a * (1 -r), s0)}, where p a ∈ Q is the probability assigned to clock a. Let us note that Eq. 1 is satisfied since

t∈ s0 • π(t) = (1 -p a ) + (p a * r) + (p a * (1 -r)) = 1. Eq. 2 is also satisfied since t∈ s0 • b π(t) = p a * r = p b and t∈ s0 • a π(t) = (p a * r) + (p a * (1 -r)) = p a .
If no probability was assigned then the default is 2/3. If no rate is assigned, then r defaults to 1/2. With default values, each one of the three transitions has a probability of 1/3, i.e., each transition has the same probability to be fired. The transition {b} however has a probability of 0 since it would contradict the subclocking relation.

Note that if both the probability of a is given and the rate of b relative to a is given, then p b = p a * r. In any other cases, the specification is ill-formed.

The synchrony constraint is a special case of subclock defined as follows

a = b ≡ b ⊆ a rate 1, which implies p a = p b . s0 {a} [pa * r]{a, b} [1 -pa]{} [pa * (1 -r)] s0 [pa]{a} [1 -pa -p b ]{} [p b ]{b} (a) b ⊆ a rate r (b) a # b
Fig. 2: pclts for subclocking and exclusion The exclusion constraint (see Fig. 2(b)), a # b is defined similarly using the following transition relation -→ # = {(s 0 , {}, 1-p b -p a , s 0 ), (s 0 , {b}, p b , s 0 ), (s 0 , {a}, p a , s 0 )}. Again the two consistency rules are satisfied. Causality. Let a, c ∈ C be two clocks. The causality relation (see Fig. 3), a ≤ c is defined as a pclts N, {a, c}, -→ ≤ , such that -→ ≤ = i∈N s i , where N is the set of natural numbers, s 0 = {(0, {}, 1 -p a , 0), (0, {a, c}, p c , 0), (0, {a}, p a -p c , 0)}, and for all i ∈ N, i > 0,

0 start 1 i ∈ N [pa * (1 -pc)]{c} [pa -pc]{a} [1 -pa]{} [pc]{a, c} [(1 -pa) * pc]{a} [(1 -pa) * (1 -pc)]{} [pa * pc]{a, c} . . . . . . . . .
s i = {(i, {}, (1 - p a ) * (1 -p c ), i), (i, {a, c}, p a * p c , i), (i, {a}, p a * (1 -p c ), i + 1), (i, {c}, (1 -p a ) * p c , i -1)}. p a , p c ∈ Q
are the probabilities assigned to clocks a and c respectively.

Note that whatever the values of p a and p c , then a will occur more frequently than c since c cannot occur alone in state 0 and to reach state 0, a and c must have occurred exactly the same number of times.

When a size is associated with the causality constraint, then the transition system becomes finite and is such that the number of states equals to size + 1. With size 1, there are two states, with size 2 there are three states and so on.

The precedence is very similar to causality except that in state 0 the simultaneous occurrence of a and c is forbidden. The whole semantics of the operators is available in [START_REF] Mallet | Correctness issues on MARTE/CCSL constraints[END_REF] and we just give here the ones used in this paper and that have been extended with stochastic constructs.

Composition

Each constraint is expressed as a clts and a specification is then captured as the synchronized product [START_REF] Arnold | The altarica formalism for describing concurrent systems[END_REF] of all the clts. The process for standard CCSL constraints is explained in detail in [START_REF] Mallet | Correctness issues on MARTE/CCSL constraints[END_REF].

With the proposed extensions, for each clts L = S, C, -→ L that results from the synchronized product and such that C is the union of all the clocks of all the composed clts, we derive a pclts P = S, C, -→ P such that the two consistency rules (Eq.1-2) are satisfied.

For instance, let us consider the following pCCSL specification:

1. Clock a with probability p a 2. Clock b ⊆ a rate: r 3. Clock a ≤ c size: 1

The semantics of each constraint is captured with a clts which are then composed through a synchronized product. Then the probabilities are added to the transitions resulting in the pclts shown in Fig. 4.

s0 start s1 [pa * (1 -r)]{a} [1/2]{c} [1 -pa]{} [pa * r]{a, b} [1/2]{}
Fig. 4: Example of composition through synchronized product.

Note that Eq. 1 is satisfied since

t∈ s0 • π(t) = (1 -p a ) + (p a * (1 -r)) + (p a * r) = 1 and t∈ s1 • π(t) = 1/2 + 1/2 = 1. Eq. 2 is also satisfied since t∈ s0 • a π(t) = (p a * (1 -r)) + (p a * r) = p a and t∈ s0 • b π(t) = p a * r = p b .
Since the probability of c is not given and c is the only clock that can tick from state s1 then the outgoing transitions are assigned with the default value 1/2 since

|s1 • | = 2.

Modeling and Refining Energy-Aware Building with MARTE/pCCSL

This section describes the application of MARTE/pCCSL to a hybrid system benchmark, which addresses the control problem of the temperature of rooms in a given building with limited heaters. With MARTE/pCCSL, we can flexibly model every part of the system in multi-views (structure, equation, behavior, and clock constraints) with quite loose coupling between them. MARTE facilitates the modeling of Time and Non Functional Properties (NFP, here energy and temperature), and pCCSL provides precise and probabilistic time control through definition and constraints of logical clocks used between abstract specification and its refinements to help us expose the best alternative solutions. Next we briefly introduce the case of energy-aware building and then present its SysML/MARTE/pCCSL models in multi-views.

Energy-Aware Building Setup

Fig. 5: Representation of building layout, outside temperature and dynamic user profile As a starting point, consider a setup proposed by [START_REF] David | An evaluation framework for energy aware buildings using statistical model checking[END_REF] as an extension of the challenging benchmark for hybrid systems model-checkers addressed in [START_REF] Fehnker | Benchmarks for hybrid systems verification[END_REF]. The benchmark consists of a building layout with temperature dynamics, autonomous heaters and a central controller deciding which room gets a heater. The room temperature dynamics is described by a differential equation:

T i = j =i a i,j (T j -T i ) + b i (u -T i ) + c i h i
where T i and T j are the temperatures in rooms i and j respectively, u is the environment temperature and h i is equal to 1 when the heater is turned on in the room i and 0 otherwise. The building layout is encoded by an adjacency matrix a where a i,j is a heat exchange coefficient between rooms i and j. The heat exchange with an environment is encoded in a separate vector b, where b i is the energy loss coefficient for room i. An energy supply from a heater is encoded as a vector c, where c i is a power coefficient for room i. Fig. 5(a) shows a benchmark building configuration instance (HEAT15 in [START_REF] Fehnker | Benchmarks for hybrid systems verification[END_REF]) with rooms and heaters, where the wall thickness corresponds to an isolation defined by a and b. The definition of matrix a, vectors b and c can be found in [START_REF] David | An evaluation framework for energy aware buildings using statistical model checking[END_REF].

Each heater is equipped with a bang-bang controller configured to turn on the heating (h i := 1) when the temperature T i is below threshold on i and turn off (h i := 0) when the temperature T i is greater than off i . Whenever the heating is turned on, the heaters consume an amount of power denoted by vector pow. The central controller can switch-over the heating from one room to another. The room is said to be needing a heater if the temperature drops below its get threshold and it is said to be outside comfort zone if the temperature drops below low. To decide when the heating can be switched over, we consider a control strategy, which is based on heuristics that the temperature difference between rooms should not be too high. To be precise, if room i whose temperature T i ≤ get i has no heater while room j has a heater, the heater can move from room j to room i when the difference T j -T i ≥ dif i . To reduce the non-determinism further, we consider probabilistic choices between the possible room destinations denoted by probabilistic weights imp. The temperature thresholds for each room used above also refer to [START_REF] David | An evaluation framework for energy aware buildings using statistical model checking[END_REF].

Further we propose to augment this setup with specific weather conditions and different user profiles to make a more realistic case for optimizing the energy consumption. The benchmark [START_REF] Fehnker | Benchmarks for hybrid systems verification[END_REF] assumes that the environment temperature is within a range between 0 • C and -2 • C without any specific dynamics. Fig. 5(b) shows one daily cycle and Fig. 5(c) shows the user profile with dynamic temperature thresholds in a day. To this user profile, we have two different refinement versions, one of which is a more complex user profile with probabilistic choice that is explained later to show the usage of MARTE/pCCSL and probabilistic clock and to demonstrate how to refine the models. Structure Modeling As described above, the structure of our energy-aware building is modeled with SysML Internal Block Diagrams shown in Fig. 6. BuildingContext represents the whole context that consists of five rooms and a central controller as well as the environment such as weather. Room represents the common template of five rooms, each of which contains a heater, a TempSensor, and a UserProfile for user-defined temperature conditioning thresholds. We limit the maximum number of heaters turned on at a time to three, though each room has a heater. Each part in the energy-aware building is stereotyped by clock that means their behaviors should be controlled by certain clock constraints (specified in pCCSL). For instance, controller monitors the information of rooms with the constraint monitor = sensor i (i is from 1 to 5), representing that monitor clock in controller and sensor clock of each room is synchronous. Besides such discrete control behavior, continuous behavior such as temperature change of each room as well as outside environment (weather ) is constrained by clock, that is, continuous behavior is discretized by clocks. We use several nfpTypes to define related variables like Weather.T and Room.T as NFP Temperature, Heater.energy as NFP Energy, Heater.power as NFP Power so that NFP properties like unit can be used. Several nfpTypes are imported from MARTE Library (NFP Energy and NFP Power ), and NFP Temperature is defined by us using MARTE constructs (units, dimension and nfpType stereotypes). The definition is simple and similar to most nfpTypes defined in MARTE Library, so we do not discuss the details. Time/Clock Domain Fig. 7 shows the time domain of the whole system. Seven clocks are defined to constrain the time of the system and three of them are defined as probabilistic ones with the keyword rate. We define a new clockType called BuildingClock which serves as a base time of the whole system. BuildingClock is discrete and owns a read-only attribute resolution. sysClk, hour, stepClk are three instances of BuildingClock with their own time units and resolutions for different usage. sysClk with high resolution (0.01s) is used in precise control process; hour with low resolution (1h = 3600s) is used to specify user profile. stepClk is similar to sysClk, but is mainly used for discretizing continuous behavior. Then we import idealClk from MARTE TimeLibrary to constrain our clock instances with pCCSL. Besides, we define three probabilistic subclocks substepClk, refinehour1 and refinehour2. substepClk is the subclock of stepClk with rate 2/3. refinehour1 and refinehour2 is the subclock of hour, where refinehour1 is used to describe user behaviors when the user may need to go out for meeting in the afternoon as well as refinehour2 is used when the user will not go out and just work. 8. We present these activities stereotyped by timedProcessing, of which the first represents monitoring and scheduling behaviors of the system and the rest represent pre-defined user behaviors (user profiles). System behavior references sysClk and the whole control process is triggered by a timedEvent called monitorEvent that occurs every 120 seconds. Each time the system readSen-sorData (e.g. temperature of rooms) and also readProfileData (e.g. whether user stays in the room or not). Then the system starts Schedule process made of three sequential actions: getNeedList (to collect need and heating information), chooseMostImportant (to find the room that needs to be heated most according to the strategy if more than one room is in the list), and changeHeatingConfig (that changes heating parameter of rooms or mode of heaters). Each action is also stereotyped by timedProcessing and may include complex operations implemented in code. The chooseMostImportant contains the code implementation of scheduling strategy described in the last subsection. The changeHeatingConfig contains the code that will set necessary parameters such as Heater.on (for denoting the running mode of heater), HeatingVector (for denoting the heated state of room), and Room.cold (for denoting the discomfort state of user) and so on. The pink part in Fig. 8 is the pCCSL specification including three kinds of constraints:

MARTE/pCCSL Models of Energy-Aware Building

distribution clock (e.g. Clock getNeedList with distribution 0.03 that means getNeedList occurs with an exponential time delay whose parameter λ = 0.03 on this clock) to describe the possibility of the unstable performance of sensors; clock synchrony (e.g. monitorEvent = readSensorData that means monitorEvent coincides with readSensorData); clock causality (e.g. getNeedList ≤ chooseMostImportant that means getNeedList is always followed by chooseMostImportant).

Refinement of Behavior Models

The refinement processes are shown in Fig. 10 and Fig. 11. The refined version shown as Fig. 10 uses subclock refinehour1 with rate 1. User arrives at the building between 8:00 and 9:00. After 4 hours, user may go for lunch between 12:00 and 13:00. After lunch, user can either go on working in the building (with a probability of 0.6) or go to a meeting (with a probability 0.4). If working in the afternoon, user may leave the building after 4 hours (i.e. between 17:00 and 18:00). If not, user will not go back to the building until next day. When user stays in the building, temperature thresholds are relatively higher than that when user does not, so that energy can be saved due to low temperature duration. As shown in the pink part of Fig. 10, pCCSL specifies the probabilistic clocks (e.g. Clock outMeeting with probability 0.4 that means outMeeting may occur with a probability of 0.4). We also use MARTE stereotype GaStep from the Performance Analysis Modeling (PAM) subprofile to report that information on the UML model. The refined version shown as Fig. 11 uses subclock refinehour2 with rate 3/4. The user doesn't need to work 8 hours a day, he/she only needs to work determined time given by subclock in the morning and afternoon, and the user can leave in the rest time and we neglect it. User may also arrive at the building between 8:00 and 9:00, and then, in 4 hours of morning time, according to the rate, user may work 3 hour randomly and the other time is rest time (neglected). In the afternoon, user may work 3 hour randomly, too. Thus, as we can see, the abstract specification and its refinements are connected logically by clock hour and its subclock refinehour1 and refinehour2. Why we use user profiles to model is that user profiles are stochastic behaviors and we can refine these behaviors. The relation between Fig. 9, Fig. 10 and Fig. 11 is that Fig. 9 is the abstract specification of user profiles, while Fig. 10 is refinement1 and Fig. 11 is the refinement2 of the abstract specification. Lay down an abstract specification of some behaviors first, and then we can refine it with probabilistic clocks to model stochastic behaviors and various solutions.

Equation Modeling Continuous behaviors described by differential equations also play an important role. Most continuous variables need to be monitored in effective real-time control, such as temperature and energy consumption in this case. To model continuous behaviors of CPSs, the parametric models are used to present the Ordinary Differential Equations (ODE). As shown in Fig. 12, we consider four equations: RoomTemperatureEquation (that determines the temperature change of each room according to the adjacent room, weather, and heater), HeaterEnergyEquation (for monitoring energy consumption of each heater), DiscomfortMonitorEquation (for monitoring discomfort value of users), WeatherEquation (that describes outside temperature curve). The former three equations use derivatives like d(T)/d(t) where t refers to the clock substepClk defined already as discretization step. Fig. 12: Equations of continuous behavior of energy aware building (modeled as a SysML parametrics)

From this case study, we can find the reason why we use pCCSL rather than CCSL. Since pCCSL has the concept of rate, it can present the stochastic behaviors of the system effectivly and simply. CCSL does has logical time and ticks, but its tick position is certain. For instance, if A is a subclock of B and A ticks every four B, then the tick positions of A must be the fist position of every four consecutive B, that is the first, the fifth, the ninth and so on. But for pCCSL, if A is a subclock of B with rate 1/4, the tick position of A is stochastic. We can only ensure that A will tick only once in four consecutive B, but the certain position is unknown. Thus, it can be seen that pCCSL gives a global rate rather than a precise position as before, which can be used to model stochastic behaviors of CPS.

Transformation from MARTE/pCCSL to SHA for Evaluation

This section discusses the transformation approach from MARTE/pCCSL specification to Stochastic Hybrid Automata (SHA) used in Uppaal-smc [START_REF] Bulychev | Uppaalsmc: Statistical model checking for priced timed automata[END_REF]. The transformation mainly targets stochastic behaviors and continuous behaviors (ODE) in terms of SHA. The core elements of MARTE/pCCSL can be directly mapped to the elements of SHA. At last, we compare the energy consumption of users using Uppaal-smc aiming to help designers refine the design models.

Stochastic Hybrid Automata

Definition 4 (Stochastic Hybrid Automata). A Stochastic Hybrid Automaton (SHA) H is a tuple H = (L, l 0 , X, Σ, E, F, I), where -L is a finite set of locations, l 0 ∈ L is the initial location, -X is a finite set of continuous variables, -Σ is a finite set of actions, and Σ = Σ i + Σ o , where Σ i is the set of input actions, Σ o is the set of output actions, -E is a finite set of transitions, for each transition denoted by (l, g, α, ϕ, l ),l, l ∈ L, g is a predicate defined by R X and action label α ∈ Σ, ϕ is a binary relationship on R X . -F (l) is a time delay function for each location l ∈ L, -I is a finite set of invariants. Uppaal-smc supports the analysis of stochastic hybrid automata (SHA) that are timed automata whose clock rates can be changed to be constants or expressions depending on other clocks, effectively defining ODEs. This generalizes the model used in previous works [START_REF] David | Statistical model checking for networks of priced timed automata[END_REF][START_REF] David | Time for statistical model checking of real-time systems[END_REF] where only linear priced automata were handled.

Definition 5 (Semantics of SHA). The semantics are denoted by a timed Labeled Transition System, which supports the seamlessly transformation from MARTE/pCCSL to SHA. For a SHA H, each location (l, v) ∈ L × R X and v |= I(l), where the transitions may be of two kinds:

-timed transition: (l, v) d -→ (l, v ), where d ∈ R ≥0 and v = F (d, v) -event-triggered transition: (l, v) α -→ (l , v ) , which means the transition is enabled with v |= g and ϕ ∈ (v, v ).
For timed transition, the probability density distribution of delay is a uniform distribution or an exponential distribution depending on the invariant of l and µ s (t) • dt = 1 (µ s is the probabilistic density function). Let E l denotes the disjunction of guards g such that (l, g, o, -, -) ∈ E j for some output o. Let d(l,v) denotes the infimum delay and D(l,v) denotes the supremum delay, i.e. d(l

,v)= inf{d ∈ R ≥0 : v + R j • d |= E l } and D(l,v)= sup{d ∈ R ≥0 : v + R j • d |= I j l }. If D(l, v
)<∞, the probabilistic density function µ s means a uniform distribution on [d(l, v), D(l, v)]. Otherwise, µ s means an exponential distribution with rate P(l), where I j l dosen't put an upper bound on the possible delays out of (l,v) and P : L j → R ≥0 is an additional distribution rate. For action transition, the output probability function γ s over Σ j o is the uniform distribution over the set {o : (l, g, o, -, -) ∈ E j ∧ v |= g} whenever this set is non-empty. The detailed semantics of SHA is referred to [START_REF] David | Statistical model checking for networks of priced timed automata[END_REF].

Transformation from MARTE/pCCSL Models to SHA

Transformation from Mode/State-based MARTE/CCSL behavior to Timed Automata (TA) has been discussed in our previous work [START_REF] Suryadevara | Verifying MARTE/CCSL mode behaviors using UPPAAL[END_REF]. The transformation approach to be presented in this section will reinforce it by encoding activity behavior, stochastic behavior and continuous behavior into SHA. First, we need to unify all chronometric clocks by converting all physical time units into one. For instance, if we choose the time unit of highest resolution s as the base unit, the conversion has to be conducted through multiplying all clocks in userProfile using the unit h by 3600. Then we can encode behaviors (modelled by Activity Diagram or State machine) with pclts according to the semantics of pCCSL. After that, we can map pclts without probabilistic time delay and choice into TA, which addresses two classes of timing constraints: deterministic delay and non-deterministic duration. For instance, the duration of absent is [START_REF] Boutekkouk | UML2.0 profiles for embedded systems and systems on a chip (socs)[END_REF][START_REF] Brosig | Quantitative evaluation of model-driven performance analysis and simulation of component-based architectures[END_REF](h) that is mapped to an invariant (t ≤ 9) and a guard (t ≥ 8). Then, we add probabilistic aspects in transformed TA. For the activity who contains several actions constrained by a clock with exponential distribution, it will be mapped to an location with its exponential parameter λ equalling to the exponential distribution. When we use GaStep with the property prob to describe the probabilistic choice, it is mapped to two branches with their probabilistic weights. The mapping rules are summarized in Fig. 13. Due to the space limitation of the paper, the detailed correctness demonstration is not given here. But, the correctness of mapping is ensured with the semantics of pCCSL in [START_REF] Suryadevara | Verifying MARTE/CCSL mode behaviors using UPPAAL[END_REF][START_REF] Zhang | Timed automata semantics of spatial-temporal consistency language STeC[END_REF].

Based on the above mapping rules, we can get corresponding SHA models. And then, we can compare and evaluate different refinement versions with Uppaal-smc to expose the best alternative one.

Evaluation with UPPAAL-SMC

With transformed SHA, we conduct evaluation experiment in Uppaal-smc. We can get a specific and quantitative analysis using SMC and the evaluation analysis results help us to expose the proper one from different refinement solutions.

To compare and evaluate energy consumption of two refined versions, the following query is used:

P r[energy ≤ 1000000]( time ≥ 5 * 24)
The experiment result is shown in Fig. 14. As we can see, energy consumption of refinement1 is about 50% more than that of refinement2. From the model point of view, the time of staying at the building of refinement1 is 4 + 4 * 0.6 = 6.4hours, which is a little more than that of refinement2. Therefore, the energy consumption of refinement2 should be less than that of refinement1 and our experiment result dose confirm this. So, if we just consider energy-saving, the user profile of refinement version 2 is more suitable.

The details of our experiments can be found in https://github.com/ECNU-MODANA/AL-Modana. git.

Related Work

This section compares our approach to related works for 1) the modeling of Cyber-Physical Systems in general, 2) the use of MARTE to conduct various kinds of analysis, 3) the use of Statistical Model Checking along with UML-based models.

The main challenges [START_REF] Lee | Cyber physical systems: Design challenges[END_REF] of designing CPSs resides in combining physical and computational models to deal with both the digital embedded systems, their environment and their interactions in a close-loop fashion. In [START_REF] Lee | The past, present and future of cyber-physical systems: A focus on models[END_REF], the author acknowledges the need to extend deterministic paradigms with probabilistic constructs to capture the unknown (and unknowable) behaviors in a more precise manner. PTides [START_REF] Zou | Execution strategies for PTIDES, a programming model for distributed embedded systems[END_REF] is mentioned as a successful example of simulation framework that can deal with these challenges. While PTides gives an ad-hoc solution, we strive to rely on industry-standard models as a front-end. We think that UML has a large-enough acceptance in industry to be a good base, and we use some of its extensions to deal with non-functional properties (with MARTE) and systems engineering properties (Parametrics of SysML). PTides, as an analysis tool, could very well be used as a back-end provided that a semantic-preserving transformation is given.

Combining MARTE and SysML to address issues of CPSs was initially suggested in [START_REF] Selic | Modeling and Analysis of Real-Time and Embedded Systems with UML and MARTE[END_REF]. MARTE was also experimented on industrial-size case studies [START_REF] Demathieu | First experiments using the UML profile for MARTE[END_REF] and proved [START_REF] Iqbal | Experiences of applying UML/MARTE on three industrial projects[END_REF] first, to be up to the task, second, to require additional ad-hoc extensions depending on the kind of analysis that is targeted. MARTE has also been used standalone on academic examples to address some aspects of CPSs, notably in [START_REF] Liu | Spatio-temporal UML statechart for cyber-physical systems[END_REF] where an extension of StateCharts is proposed and in [START_REF] Liu | Hybrid MARTE statecharts[END_REF] where Fig. 14: Energy consumption comparison between two refinement versions a hybrid version is explored. MARTE is also often used for its ability to describe non-functional properties [START_REF] Berardinelli | UML profiles for non-functional properties at work: Analyzing reliability, availability and performance[END_REF] or as a base to conduct performance analysis [START_REF] Espinoza | Annotating UML models with non-functional properties for quantitative analysis[END_REF][START_REF] Tribastone | Automatic extraction of PEPA performance models from UML activity diagrams annotated with the MARTE profile[END_REF] or other sort of non-functional analysis like dependability analysis [START_REF] Bernardi | A dependability profile within MARTE[END_REF]. In literature [START_REF] Shorin | Transforming uml state machines into stochastic petri nets for energy consumption estimation of embedded systems[END_REF], they use MARTE state machines to model embedded system and transform it to stochastic petri net for energy consumption estimation. In all these examples, only one model of UML is used (like state machine or sequence diagram) while this paper focuses more on finding a practical solution to combine together several UML/SysML/MARTE models. We believe that logical clocks, in a similar way than tagged systems [START_REF] Lee | A framework for comparing models of computation[END_REF], provides the good abstraction level to do that. CCSL [START_REF] André | Syntax and semantics of the Clock Constraint Specification Language (CCSL)[END_REF] then provides a concrete syntax to handle logical clocks. While CCSL is a deterministic declarative polychronous specification language, its pCCSL extension provides new constructs to model the unknown(able) behaviors. Our proposed transformation into stochastic hybrid automata is largely inspired by our previous work to transform CCSL to timed automata [START_REF] Suryadevara | Verifying MARTE/CCSL mode behaviors using UPPAAL[END_REF][START_REF] Zhang | Timed automata semantics of spatial-temporal consistency language STeC[END_REF] however we deal here with new operators of CCSL and with its stochastic extension.

Other works, either based on MARTE [START_REF] Shorin | Formal description of an approach for power consumption estimation of embedded systems[END_REF] or just on UML [START_REF] Brosig | Quantitative evaluation of model-driven performance analysis and simulation of component-based architectures[END_REF] also use stochastic models to analyze the energy consumption. In this case, the authors rely on Stochastic Petri Nets while in our work we use Stochastic Hybrid Automata. The key difference here again is in the use of probabilistic logical clocks to coordinate the models. However, it is rather a practical matter for us and Stochastic Petri Nets could also be used as a backend instead.

Statistical model checking is a highly scalable simulation-based approach which is useful to bound the probability of making an error by increasing the simulation effort. SMC gets widely accepted in various research areas such as software engineering, in particular for industrial applications [START_REF] Boudjadar | Schedulability and energy efficiency for multi-core hierarchical scheduling systems[END_REF][START_REF] Fang | Formal verification and simulation for platform screen doors and collision avoidance in subway control systems[END_REF], or even for solving problems originating from systems biology [START_REF] David | Statistical model checker for biological systems[END_REF][START_REF] Jiang | Modeling and verification of a dual chamber implantable pacemaker[END_REF]. Inspired by our previous work [START_REF] David | An evaluation framework for energy aware buildings using statistical model checking[END_REF] which particularly focuses on the analysis of energy consumption of energy-aware building case in different cases to find the most significant factors that influence energy consumption, we attempted to apply MARTE/pCCSL to this case and further evaluate NFP performance of CPSs with SMC.

Conclusion

This paper explores an idea of using probabilistic logical clocks in relation to UML models to capture several aspects of CPSs including the stochastic behavior of unknowable events. We believe that using industry-standard is important for a better adoption and we show how to make a consistent use of SysML and MARTE. Clocks through pCCSL are used to coordinate multi-view models. The underlying pCCSL specification is key for the semantic interpretation of models and the transformation into Stochastic Hybrid Automata.

To summarize our contributions 1) MARTE/CCSL has been extended to model stochastic behaviors. The syntax and semantics models of pCCSL are introduced; 2) To show how pCCSL specifications can augment SysML/MARTE models to link together several models and refine the abstract specification; 3) According to the semantics model of pCCSL, we present the transformation process from MARTE/pCCSL to SHA, whose aim is to perform evaluation on the alternative solutions with statistical model checking, which helps designers expose the best solution. The resulting models are analyzed with Uppaal-smc. The process is demonstrated on a simple case study of energy-aware building.

MARTE proposes a subprofile, called Performance Analysis Modeling (PAM), dedicated to performance analysis. We think this contribution is an important step for the integration of MARTE Time model with PAM models. Future work shall consider a more extensive use of PAM, while pCCSL may be kept as much as possible as a lower level intermediate semantic model to ease the transformation from pure UML to other formal models.

Fig. 1 :

 1 Fig. 1: Grammar of pCCSL

Fig. 3 :

 3 Fig. 3: pCCSL causality (infinite state pclts): a ≤ c.

Fig. 6 :

 6 Fig. 6: Structure of an energy-aware building (Internal Block Diagram)

Fig. 7 :

 7 Fig. 7: Time domain of an energy-aware building

Fig. 8 :

 8 Fig. 8: Activity diagram for system behaviors

Fig. 9 :

 9 Fig. 9: Abstract activity diagram for user behaviors

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Refined version 1 of activity diagram for user behaviors (refinehour1⊆hour, rate=1)

Fig. 13 :

 13 Fig. 13: Mapping rule from MARTE/pCCSL to SHA

Acknowledgment

This work is partly supported by NSFC under Grant No. 61472140,61170084, NSF of Shanghai under Grant No. 14ZR1412500, and the Danish National Research Foundation Grant No. 61361136002.


Sébastien Fenet
email: sebastien.fenet@zoho.com

Roland Badeau

Gaël Richard

Reassigned Time-Frequency Representations of Discrete Time Signals and Application to the Constant-Q Transform

Keywords: Reassigned time-frequency representations, reassigned CQT, discrete signals, synchrosqueezing 5 models

In this paper we provide a formal justification of the use of time-frequency reassignment techniques on time-frequency transforms of discrete time signals.

State of the art techniques indeed rely on formulae established in the continuous

case which are applied, in a somehow inaccurate manner, to discrete time signals.

Here, we formally derive a general framework for discrete time reassignment. To illustrate its applicability and generality, this framework is applied to a specific transform: the Constant-Q Transform.

Introduction

Time-frequency reassignment has received great attention over the last decades, especially for the task of sinusoidal parameter estimation in noisy data.

Numerous methods have been developed based on Fourier analysis [START_REF] Auger | Improving the readability of time-frequency and time-scale representation by the reassignment method[END_REF][START_REF] Abe | Design criteria for simple sinusoidal parameter estimation based on quadratic interpolation of FFT magnitude peaks[END_REF][START_REF] Betser | Review and discussion on STFT-based frequency estimation methods[END_REF], on subspace decomposition [START_REF] Roy | Esprit -a subspace rotation approach to estimation of parameters of cisoids in noise[END_REF][START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] or on more general models such as AM/FM some time-frequency point in the signal to a different point in the time-frequency plane that depends on the window used for the spectral computation. Time-frequency reassignment methods emerge from the idea first proposed by Kodera [8]. This original approach uses the phase information of the timefrequency representation and remains difficult to use in practice. Later on, Auger and Flandrin [START_REF] Auger | Improving the readability of time-frequency and time-scale representation by the reassignment method[END_REF] proposed a new closed-form solution to this problem which applies to a wide variety of time-frequency representations and relies on much more straightforward computations of the reassigned indexes. This work has opened the door to the use of time-frequency reassignment in numerous domains such as physics [START_REF] Niethammer | Time-frequency representations of lamb waves[END_REF], radar imaging [START_REF] Ram | Doppler-based detection and tracking of humans in indoor environments[END_REF] or audio [START_REF] Peeters | Template-based estimation of time-varying tempo[END_REF]. It has also led to the development of numerous extensions and adaptations of the original method [START_REF] Holighaus | Reassignment and synchrosqueezing for general time-frequency filter banks, subsampling and processing[END_REF][START_REF] Auger | Making Reassignment Adjustable: the Levenberg-Marquardt Approach[END_REF][START_REF] Oberlin | Second-order synchrosqueezing transform or invertible reassignment? towards ideal time-frequency representations[END_REF].

Another solution to the problem, named synchrosqueezing, has been proposed by Daubechies and Maes in [START_REF] Daubechies | A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models[END_REF]. Notably because it offers the ability to reconstruct the time signals, the technique has drawn a lot of interest and has become the root of multiple applications and enhancements [START_REF] Daubechies | Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool[END_REF][START_REF] Li | A generalized synchrosqueezing transform for enhancing signal timefrequency representation[END_REF][START_REF] Ahrabian | Synchrosqueezing-based time-frequency analysis of multivariate data[END_REF] .

Although synchrosqueezing was initially presented as a distinct technique from Auger and Flandrin's reassignment method, the strong connection that exists between the two has been clarified in [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF].

Traditionally, reassignment calculations are carried out in the context of continuous time signals (see [START_REF] Hainsworth | Time frequency reassignment: A review and analysis[END_REF] for a review) while most applications involve discrete time signals. In practice, all results obtained in the continuous time case are applied, without detailed justification, to the discrete case.

In this work, we propose a formal framework for the computation of the reassigned transforms which fully takes into account the discrete time aspect. Our approach consists of expressing the magnitude of the time-frequency transform of a discrete signal as a mass function in the time-frequency plane and in assigning the energy to the centre of mass of this representation. Interestingly, the obtained mathematical expressions are very similar to the classical expressions proved in the continuous time case. The main advantage of our approach is that we obtain exact expressions of the solution when discrete signals are considered. This opens the door to the implementation of exact solutions and, should the implementation require an approximate solution, we are able to characterise the introduced error. To some extent, this work also gives a formal justification of the common approximation made when applying continuous time formulae to discrete time signals.

The paper is organised as follows. The mathematical model and the derivations of closed-form expressions for reassigned time and frequency indexes are provided in section 2. We then discuss in the subsequent section the merit of the proposed solution. An application of the framework to a specific timefrequency representation, the Constant Q Transform (CQT), is finally proposed in section 4.

Mathematical model

Traditional time-frequency representations

As stated in the introduction, our aim is to derive the mathematical formulation of the reassigned time-frequency representation of a discrete numerical signal. We thus consider a discrete time signal x. Such a signal maps any discrete index n P Z to a complex value x n P C. Its frequency content Xpξq is defined for any normalised frequency ξ P R{Z by the Discrete Time Fourier Tranform (we use the standard notation R{Z to indicate that the normalised frequency is defined modulo 1).

The information that one reads in a time-frequency representation of x is the amount of energy in x at time t P R and normalised frequency ν P r0; 1s.

In order to evaluate this energy, the scalar product between x and a kernel is computed. The kernel consists of a windowing sequence h t,ν , centred on t, multiplied by the harmonic function of frequency ν. Let us note that the pt, νq exponent makes it explicit that the windowing sequence depends on the time of interest and may also depends on the frequency of interest. More precisely, the transforms of x that fall within the scope of this paper can be written in the following form:

T x,h t,ν pt, νq " ÿ nPZ h t,ν n x n e ´j2πνn (1)
The time-frequency representation at time t and frequency ν is finally obtained by considering the squared magnitude of the transform:

spt, νq " | T x,h t,ν pt, νq | 2 . ( 2
)
It is interesting here to recall the Heisenberg-Gabor limit [START_REF] Gabor | Theory of communication[END_REF] that constrains the design of the windowing sequence h t,ν . More precisely, the Gabor limit states that there is a trade-off between the temporal and spectral resolutions when representing a signal in the time-frequency plane. In practice, adjusting the support of the windowing sequence is a direct way to tune this trade-off. A wide support will result in a precise frequency resolution with a poor temporal resolution. Conversely, a narrow support will provide a good temporal resolution at the cost of the frequency resolution. In order to ensure consistency, we consider that the windowing sequences are of finite support and that they are normalised by the size of their supports. For instance, h being a continuous window function of finite temporal support, the windowing sequence is defined by h t,ν n " h pn ´tq in the case of a short-term Fourier Transform (STFT) or by h t,ν n " νh pνpn ´tqq in the case of a constant-Q Transform (CQT) for a set of frequencies ν within " 0; 1 2 ‰ (see Section 4 for more details).

In addition to constraining the design of the window, the choice of a given time-frequency transform also determines the set of time-frequency points pt, νq at which the representation is evaluated. Typically, a Short-Time Fourier Transform with a temporal hop size ∆ t and a spectral hop size ∆ ν is obtained with the following set of points: tpt 0 `k∆ t , ν 0 `k1 ∆ ν q for pk, k 1 q P N 2 qu. In contrast, the Constant-Q Transform, whose spectral geometric progression is often denoted by 2 1 r (r being the number of bins per octave), is obtained with the set tpt 0 `k∆ t , ν 0 2 1 r k 1 q for pk, k 1 q P N 2 u. In these expressions, t 0 naturally denotes the lowest time index of the representation and ν 0 the lowest frequency bin.

Let us make explicit here that in the following derivations we will be using the notation z for the complex conjugate of z and the symbol ˙for the discrete convolution operator.

Time-frequency representations as mass functions

Reassignment techniques rely on the idea that time-frequency representations, at a given point pt, νq, can be written as the sum of a mass function defined on the time-frequency plane pn, ξq. Given our context, which involves a discrete time axis and a periodic frequency axis, we have pn, ξq P Z ˆpR{Zq.

We thus look for an expression of the form:

spt, νq "

ÿ nPZ ż ν`1 2 ξ"ν´1 2 Φ t,ν pn, ξqdξ (3)
where the function Φ t,ν is real-valued.

Let us define W , a discrete version of Rihaczek's ambiguity function [START_REF] Rihaczek | Signal energy distribution in time and frequency[END_REF], for any sequence ϕ P 1 pZq, time index n P Z and frequency ξ P R{Z by:

W ϕ pn, ξq " ÿ τ PZ ϕ n`τ ϕ n e ´j2πξτ . (4)
Proposition 1. The time-frequency representation spt, νq of a discrete time signal x P 1 pZq, as defined in equation [START_REF] Abe | Design criteria for simple sinusoidal parameter estimation based on quadratic interpolation of FFT magnitude peaks[END_REF], can be written as the sum of a mass function, as in equation ( 3), with:

Φ t,ν pn, ξq " ! W h t,ν pn, ν ´ξq W x pn, ξq
) .

Putting things together, this means that the time-frequency representation of x at point pt, νq can be written in the following form:

spt, νq " ÿ nPZ ż ν`1 2 ξ"ν´1 2 ! W h t,ν pn, ν ´ξq W x pn, ξq ) dξ (5)
Proof. Let us evaluate the following expression:

Ept, νq " ÿ nPZ ż ν`1 2 ξ"ν´1 2 W h t,ν pn, ν ´ξq W x pn, ξq dξ .
We have:

Ept, νq " ÿ nPZ ż ν`1 2 ξ"ν´1 2 # ÿ τ1PZ h t,ν n`τ1 h t,ν n e ´j2πpν´ξqτ1 ÿ τ2PZ
x n`τ2 x n e ´j2πξτ2 + dξ .

Knowing that h t,ν is of finite support and that x is in 1 pZq, Fubini's theorem ensures that the summations can be permuted:

Ept, νq " ÿ n,τ1,τ2PZ # h t,ν n`τ1 h t,ν n x n`τ2 x n e ´j2πντ1 ż ν`1 2 ξ"ν´1 2 e ´j2πξpτ2´τ1q dξ + .
Knowing that, for any integer k, we have:

ż ν`1 2 ξ"ν´1 2 e ´j2πξk dξ " $ & % 1 if k " 0 0 otherwise
the above expression can be rewritten with respect to a single shift variable τ :

Ept, νq " ÿ n,τ PZ h t,ν n`τ h t,ν n x n`τ x n e ´j2πντ .
By the substitution τ Þ Ñ m ´n we get:

Ept, νq " ÿ n,mPZ h t,ν m h t,ν n x m x n e ´j2πνpm´nq " # ÿ mPZ h t,ν m x m e ´j2πνm +# ÿ nPZ h t,ν n x n e ´j2πνn + " | T x,h t,ν pt, νq| 2 .
Altogether, we have proved that:

110 | T x,h t,ν pt, νq| 2 " ÿ nPZ ż ν`1 2 ξ"ν´1 2 W h t,ν pn, ν ´ξq W x pn, ξq dξ . (6)
By applying the real part operator on both sides of the above equality, using its linearity and knowing that | T x,h t,ν pt, νq| is real, we get:

| T x,h t,ν pt, νq| 2 " ÿ nPZ ż ν`1 2 ξ"ν´1 2 ! W h t,ν pn, ν ´ξq W x pn, ξq ) dξ.
l

We have thus proved that the time-frequency representation of x at a given time t and frequency ν, which we have defined as the squared magnitude of the transform, can be expressed as the sum of a real-valued function in the time-frequency plane. In accordance with the context of discrete time signals,

we have a discrete time axis (leading to a discrete summation over the time index) and a continuous and periodic frequency axis (leading to a continuous and finite summation over the frequency variable). It is also interesting to note that the contributions of the window h t,ν and the signal x are quite separated and symmetric. They are induced through the same transform W that applies separately on the window and on the signal. We may now consider the function that is summed as a mass function and, as such, look for its centre of mass in the time-frequency plane.

Reassigned time index

Given the formulation (3) of the time-frequency representation at point pt, νq as the two-dimensional sum of a mass function, it is possible to evaluate the centre of mass of this distribution in the time-frequency plane. By definition, the first coordinate of the centre of mass is given by the following formula: tt,ν "

ř nPZ ş ν`1 2 ξ"ν´1 2 nΦ t,ν pn, ξqdξ ř nPZ ş ν`1 2 ξ"ν´1 2 Φ t,ν pn, ξqdξ . (7)
Proposition 2. The temporal coordinate of the centre of mass in the timefrequency plane of the time-frequency representation evaluated at point pt, νq, defined in equation [START_REF] Betser | Estimation of frequency for AM/FM models using the phase vocoder framework[END_REF], can also be expressed as:

tt,ν " t ` # T x, r t h t,ν pt, νq T x,h t,ν pt, νq + (8)
where r t is the unit ramp centred on t: r t n " n ´t.

Proof. Let us evaluate the following quantity:

E 1 pt, νq " ÿ nPZ ż ν`1 2 ξ"ν´1 2 nΦ t,ν pn, ξqdξ . (9)
We have:

E 1 pt, νq " # ÿ nPZ ż ν`1 2 ξ"ν´1 2 ÿ τ1PZ nh t,ν n`τ1 h t,ν n e ´j2πpν´ξqτ1 ÿ τ2PZ
x n`τ2 x n e ´j2πξτ2 dξ + .

By using the same grouping and the same substitution as in the proof of Proposition 1, the expression becomes:

E 1 pt, νq " # ÿ n,mPZ n h t,ν m h t,ν n x m x n e ´j2πνpm´nq + .
Expanding n as t `n ´t, we have:

E 1 pt, νq " t # ÿ nPZ h t,ν n x n e ´j2πνn ÿ mPZ h t,ν m x m e ´j2πνm + ` # ÿ nPZ pn ´tq h t,ν n x n e ´j2πνn ÿ mPZ h t,ν m x m e ´j2πνm

+

This can be more concisely written as:

E 1 pt, νq " t | T x,h t,ν pt, νq| 2 ` ! T x,r t h t,ν pt, νq T x,h t,ν pt, νq ) (10)
Given the definitions ( 7) and ( 9), we know that dividing E 1 pt, νq by

ř nPZ ş ν`1 2
ξ"ν´1 2 Φ t,ν pn, ξqdξ yields tt,ν . Besides, we know from equations ( 2) and

(3) that:

ÿ nPZ ż ν`1 2 ξ"ν´1 2 Φ t,ν pn, ξqdξ " | T x,h t,ν pt, νq| 2 (11)
" T x,h t,ν pt, νq T x,h t,ν pt, νq .

Thus, dividing equation [START_REF] Ram | Doppler-based detection and tracking of humans in indoor environments[END_REF] by

ř nPZ ş ν`1 2 ξ"ν´1 2 Φ t,ν pn, ξqdξ, we finally get: tt,ν " t ` # T x, r t h t,ν pt, νq T x,h t,ν pt, νq + l
Let us recall that we are currently looking for the centre of mass of the mass function Φ t,ν from which the time-frequency representation spt, νq stems. What we have done here is establishing a simple expression of the time coordinate tt,ν of the centre of mass. The original definition of the centre of mass, shown in equation ( 7), indeed leads to quite complex computations while our calculations have showed that the time coordinate of the centre of mass can be expressed mostly as the division of two transforms. The benefits of this expression will be further discussed in section 3 but it is already remarkable that we could express this coordinate in a nice and simple way.

Reassigned frequency

Similarly to what has been done for the time coordinate, the expression of the time-frequency representation at point pt, νq as the sum of a mass function Φ t,ν , proved in Proposition 1, allows to determine the frequency coordinate of the centre of mass of Φ t,ν . The latter is defined by: νt,ν "

ř nPZ ş ν`1 2 ξ"ν´1 2 ξΦ t,ν pn, ξqdξ ř nPZ ş ν`1 2 ξ"ν´1 2 Φ t,ν pn, ξqdξ . ( 12
)
Proposition 3. The frequency coordinate of the centre of mass in the timefrequency plane of the time-frequency representation evaluated at point pt, νq, defined in equation ( 12), can also be expressed as:

νt,ν " ν ´1 2π # T x,h t,ν ˙g pt, νq T x,h t,ν pt, νq + ( 13
)
where g is the differentiator filter, i.e. g n " p´1q n n for n ‰ 0 and g 0 " 0.

Proof. Let us evaluate the following expression:

E 2 pt, νq " ÿ nPZ ż ν`1 2 ξ"ν´1 2 ξΦ t,ν pn, ξqdξ . (14)
E 2 pt, νq " # ÿ nPZ ż ν`1 2 ξ"ν´1 2 ξ ÿ τ1PZ h t,ν n`τ1 h t,ν n e ´j2πpν´ξqτ1 ÿ τ2PZ
x n`τ2 x n e ´j2πξτ2 dξ + .

By the substitution ξ Þ Ñ ξ 1 `ν we get:

E 2 pt, νq " ν | T x,h t,ν pt, νq| 2 ` # ÿ nPZ ż 1 2 ξ 1 "´1 2 ξ 1 ÿ τ1PZ h t,ν n`τ1 h t,ν n e j2πξ 1 τ1 ÿ τ2PZ x n`τ2 x n e ´j2πpξ 1 `νqτ2 dξ 1 + looooooooooooooooooooooooooooooooooooooooooooooooomooooooooooooooooooooooooooooooooooooooooooooooooon E 2 2 pt,νq .
In the above expression, let us call the second operand E 2 2 pt, νq. We now have:

E 2 pt, νq " ν | T x,h t,ν pt, νq| 2 `E2 2 pt, νq . (15)
Adequately reordering the summations in E 2 2 pt, νq, the latter can be written as:

E 2 2 pt, νq " # ÿ n,τ1,τ2PZ h t,ν n`τ1 h t,ν n x n`τ2 x n e ´j2πντ2 ż 1 2 ξ 1 "´1 2 ξ 1 e j2πξ 1 pτ1´τ2q dξ 1 + .
Let us define the sequence g in the following way:

g n " ż 1 2 ξ 1 "´1 2 2jπξ 1 e j2πξ 1 n dξ 1 " $ & % p´1q n n if n ‰ 0 0 if n " 0 . ( 16
)
The expression of E 2 2 pt, νq consequently becomes:

170 E 2 2 pt, νq " # j 2π ÿ n,τ1,τ2PZ h t,ν n`τ1 h t,ν n x n`τ2
x n e ´j2πντ2 g τ2´τ1 + .

We can then use the substitution:

$ & % τ 2 Þ Ñ l ´n τ 1 Þ Ñ l ´n ´m
to get the following equation.

E 2 2 pt, νq " # j 2π ÿ n,m,lPZ h t,ν l´m h t,ν n x l x n e ´j2πνpl´nq g m + .
As far as the terms in m are concerned, we observe that:

ÿ mPZ h t,ν l´m g m " ph t,ν ˙gq l .
The whole expression can thus be rewritten:

E 2 2 pt, νq " # j 2π ÿ nPZ h t,ν n x n e ´j2πνn ÿ lPZ ph t,ν ˙gq l x l e ´j2πνl + .
By identifying the transforms in the above expression of E 2 2 pt, νq and reinjecting in equation ( 15), we get:

E 2 pt, νq " ν| T x,h t,ν pt, νq| 2 ´1 2π ! T x,h t,ν pt, νq T x,h t,ν ˙g pt, νq
) .

Taking into account equation [START_REF] Peeters | Template-based estimation of time-varying tempo[END_REF] and given the definition ( 14) of E 2 pt, νq, it is only needed to divide the latter by | T x,h t,ν pt, νq| 2 to obtain νt,ν , as defined in equation( 12). This finally leads to:

νt,ν " ν ´1 2π # T x,h t,ν ˙g pt, νq T x,h t,ν pt, νq + l
Similarly to what we have done in Section 2.3 for the time coordinate, we have obtained here a simple expression of the frequency coordinate of the centre of mass of Φ t,ν pn, ξq. Although the mathematical steps that led to this expression are rather different from the ones performed in the case of the first coordinate, the final result in equation ( 13) is quite symmetric with the one in equation [START_REF] Kodera | Analysis of time-varying signals with small BT values[END_REF]. In a similar fashion, the frequency coordinate is expressed as a simple combination of two transforms. This expression of νt,ν will be further discussed in the next section.

Discussion

In the previous section, we have obtained exact expressions of the coordinates p tt,ν , νt,ν q of the centre of mass of Φ t,ν . The latter is a mass function whose summation in the time frequency plane gives spt, νq i.e. the value of the time-frequency representation of x at time t and frequency ν. The obtained expressions are valid for any time-frequency representation of a discrete signal with any kind of windowing. Time-frequency reassignment consists of building an enhanced time-frequency representation of x simply by reassigning the value spt, νq to the time-frequency point p tt,ν , νt,ν q.

Practical interest of the expressions

The major advantage of the presented formulae for tt,ν and νt,ν is that they provide a way of computing the reassigned index with the computation of only three transforms, one of which is the 'original' one (T x,h t,ν pt, νq). The presented reassignment tool is thus very powerful since it only involves an increase of the computation time by a factor 3 compared to the original representation, which does not seem a very high price to pay considering the tremendous gain of resolution.

In detail, in order to compute the reassigned time index as in equation ( 8), one has to compute an additional transform of x with a window function which consists of the original window function multiplied by the unit ramp function. The real part of the division of this ramp-windowed transform by the original transform provides the time shift. In order to compute the reassigned frequency index as in equation ( 13), one has to compute the transform of x with a window sequence which is the original window convolved with some discrete filter g. Let us note that this modified window can be pre-computed once so that calculating the transform of x with this window has the same complexity as with the original window. The frequency shift is proportional to the imaginary part of the transform with the convolved window divided by the original transform.

Comparison with previously established expressions

The transition from the mathematical proof of the time-frequency reassignment formulae in the context of continuous time signals [START_REF] Auger | Improving the readability of time-frequency and time-scale representation by the reassignment method[END_REF] to the one for discrete time signals is not straightforward. The continuous proof indeed relies on the rewriting of the spectrogram as a summed product of Wigner-Ville transforms [START_REF] Wigner | On the quantum correction for thermodynamic equilibrium[END_REF][START_REF] Ville | Théorie et Applications de la Notion de Signal Analytique[END_REF][START_REF] Classen | The Wigner Distribution: A Tool for Time-Frequency Signal Analysis -Part III: Relations with other time-frequency transformations[END_REF]. However, since the Wigner-Ville transform involves non-integer shifts, it cannot be directly applied to discrete-time signals. We have overcome the issue by introducing in equation ( 4) the transform W , which can be seen as a discrete version of Rihaczek's ambiguity function [START_REF] Rihaczek | Signal energy distribution in time and frequency[END_REF]. This transform is fully applicable to discrete time signals but it lacks the symmetric design of the Wigner-Ville transform. As a result, the transform is not bounded to real values. Fortunately, we could still step back to the expression of a real-valued mass function by applying the real operator on equation [START_REF] Abe | AM/FM rate estimation for time-varying sinusoidal modeling[END_REF]. Another pitfall when switching to discrete time signals is the choice of ξ's integration interval in equation ( 3). Given the spectral periodicity of discrete signals, it is clear enough that the interval should be of unit length. However, the localisation of the interval is not obvious. It does not intervene in the calculation of the reassigned time index but it has a meaningful impact when deriving the reassigned frequency index. Our calculations have shown that, in a rather expectable and logical way, the interval of integration should be centred around ν when studying the time-frequency representation at point pt, νq.

In spite of the meaningfully dissimilar calculations steps that led to them, the formulae that we finally obtain for tt,ν and νt,ν are very similar to the ones obtained in the continuous case by Auger and Flandrin [START_REF] Auger | Improving the readability of time-frequency and time-scale representation by the reassignment method[END_REF]. The expression of the reassigned time index is indeed the same, with the same appearance in the calculation of the window multiplied by the ramp function. One meaningful contribution of our work, though, is that our expressions have been derived in the context of discrete time signals. There is thus no approximation when applying it to real signals and we have formally proved that the expression is valid. As far as the reassigned frequency index is concerned, the continuous proof, which relies on an integration by parts, leads to an expression that brings the derivative of the original window into play. On the other hand, our calculations lead to a new window which is the original one convolved with a filter g. For that matter, it is very interesting to note that the filter g that we obtain is the ideal differentiator filter. It is indeed clear in the definition of g, given in equation [START_REF] Daubechies | Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool[END_REF], that it is the inverse Fourier Transform of 2jπξ 1 which is the frequency response of the perfect differentiator. The similarity between our expression and the one from Auger and Flandrin then becomes clear. To some extent, we may indeed consider that the "derivative" of h approximates the convolution of h by g. We however take full advantage here of our discrete time approach. We indeed have an exact formulation for the reassigned frequency index whereas Auger and Flandrin propose to compute the "derivative" of h. The derivative is quite an ill-defined concept for discrete time signals and applying their formula raises the question of the best way to estimate it. Our work brings an answer to this, until now, open question. The technically correct way of using Auger and Flandrin formula is to estimate the derivative with an ideal differentiator filtering.

From reassignment to synchrosqueezing

As stated in the introduction the pioneering approach of Kodera [START_REF] Kodera | Analysis of time-varying signals with small BT values[END_REF] suggests to compute the reassigned time and frequency indexes by using the partial derivatives of the phase of the time-frequency representation. In a similar spirit, the synchrosqueezing method [START_REF] Daubechies | A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models[END_REF] makes use of the partial derivatives of the phase to compute the reassigned frequency index. A notable difference of synchrosqueezing, however, is that only the frequency index is reassigned.

Our method, on the other hand, looks for the centre of mass of the distribution Φ t,ν related to the time-frequency representation spt, νq. The equivalence between the computation of the reassigned indexes thanks to the partial derivatives of the phase and the computation of the same indexes through the centre of mass has however been proved in [START_REF] Auger | Improving the readability of time-frequency and time-scale representation by the reassignment method[END_REF].

As a result, as explained in [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF], the synchrosqueezing method can be computed within our framework in a straightforward way. First, the reassigned frequency νt,ν is determined according to equation [START_REF] Auger | Making Reassignment Adjustable: the Levenberg-Marquardt Approach[END_REF]. Finally, the complex value T x,h t,ν pt, νq is additively reassigned to the point pt, νt,ν q. The differences with the reassignment method that we detailed in the first place are thus clear and limited: only the frequency index is reassigned and the reassignment is operated on the complex value of the transform instead of its squared magnitude.

Application: A reassigned Constant-Q Transform

As explained, the formal framework derived above is valid for a wide variety of time-frequency transforms. To illustrate its applicability, we propose to demonstrate the merit of the reassignment strategy on the Constant-Q Transform. The interested reader will find a Python implementation of this reassigned Constant-Q Transform on the authors' homepages1 .

Presentation of the transform

The Constant-Q Transform is a common tool in the field of audio processing.

Its numerous fields of application include main melody extraction [START_REF] Fuentes | Probabilistic model for main melody extraction using constant-Q transform[END_REF], audiofingerprinting [START_REF] Fenet | A scalable audio fingerprint method with robustness to pitch-shifting[END_REF], chords detection [START_REF] Nawab | Identification of musical chords using constant-Q spectra[END_REF][START_REF] Harte | Detecting Harmonic Change in Musical Audio[END_REF]. It was first proposed in [START_REF] Brown | Calculation of a constant Q spectral transform[END_REF] by Brown, who aimed at designing a time-frequency representation that mimics our perception of sound. To this end, the Constant-Q Transform has frequency bins with geometrically spaced centre frequencies. Moreover, their frequency resolution is inversely proportional to the centre frequency. Its direct calculation, however, is computationally heavy compared to the Short Time Fourier

Transform which can rely on the Fast Fourier Transform. Hence the proposition by Brown and Puckette of a faster implementation in [START_REF] Brown | An efficient algorithm for the calculation of a constant Q transform[END_REF] and a follow-up in [START_REF] Schörkhuber | Constant-Q transform toolbox for music processing[END_REF] by Klapuri. Another pitfall of the initial transform is its non-invertibility. As a consequence, several authors have focused on designing an invertible version of the Constant-Q Transform [START_REF] Velasco | Constructing an invertible constant-Q transform with non-stationary Gabor frames[END_REF][START_REF] Fitzgerald | Towards an inverse constant Q transform[END_REF][START_REF] Holighaus | A Framework for Invertible, Real-Time Constant-Q Transforms[END_REF][START_REF] Nagathil | Optimal Signal Reconstruction from a Constant-Q Spectrum[END_REF]. Very recently, the concept of reassigned Constant-Q Transform has appeared in [START_REF] Pruša | Real-Time Audio Visualization With Reassigned Non-uniform Filter Bank[END_REF] and [START_REF] Holighaus | Reassignment and synchrosqueezing for general time-frequency filter banks, subsampling and processing[END_REF]. The point of view is however slightly different from the one adopted here since both works present the Constant-Q Transform as an application example of the theory derived in [START_REF] Holighaus | Reassignment and synchrosqueezing for general time-frequency filter banks, subsampling and processing[END_REF]. In the latter, reassignment formulae are specifically derived for filter banks.

Since they rely on the work by Auger and Flandrin [START_REF] Auger | Improving the readability of time-frequency and time-scale representation by the reassignment method[END_REF], the derivations are only carried out for continuous-time signals.

More precisely, in the context of the Constant-Q spectrogram, one evaluates the amount of energy in x at any time t for a set of frequencies tν k u. The ν k 's are distributed according to the following law:

ν k " 2 k r ν 0
The parameter r specifies the resolution of the transform (the bigger r is, the more numerous the frequency bins are) and ν 0 is the centre of the lowest frequency bin of the transform. A meaningful interest of the Constant-Q Transform is that the ν k 's can be aligned on the note frequencies of the Western scale. The parameter r then corresponds to the number of bins per octave, or equivalently, r{12 to the number of bins per note.

As far as the windowing functions are concerned, the Constant-Q Transform scales them with respect to frequency:

h t,ν n " νhpνpn ´tqq
where h is a window function. The above relationship implies that the temporal support of the window decreases as ν k increases. This equivalently means that the temporal resolution of the transform increases with frequency, and, according to the Heisenberg-Gabor principle [START_REF] Gabor | Theory of communication[END_REF], that its frequency resolution decreases.

The results presented thereafter have been obtained with the following set of parameters:

$ & % r " 36 ν 0 " 2.

45ˆ10

´3 and h t,ν0 is a Hann window whose length is about 22500 samples. The analysed sound has a sampling rate of 44100 Hz, which means that the lowest frequency bin of the transform corresponds to the continuous frequency 107.9 Hz and that the support of h t,ν0 is half a second long. With this setup the frequency bins of the transform are aligned on the note frequencies of the Western scale, with a resolution of 3 bins per note. More specifically, the three lowest bins of the transform are centre on 110 Hz, which corresponds to the frequency of the musical note A 2 .

Considerations on the resolution

Figure 1 shows the spectrogram of a piano sound obtained by means of the traditional Constant-Q Transform. We can see that the frequency resolution of the transform is well adapted to the musical sound since every frequency of the sound is well resolved. As mentioned, this transform has a better frequency resolution in low frequencies, where the musical frequencies are closer to each other. Of course, the counterpart of this better frequency resolution is a loss in temporal resolution. As a result, in spite of the fact that all frequencies of the piano sound have been played simultaneously, one can observe in the spectrogram that the energy in low frequencies starts earlier than in the higher frequencies.

Similarly, the energy in low frequencies vanishes later. This is a typical artifact of the Constant-Q Transform due to these resolution considerations.

Application of the reassignment tool

It is quite straightforward to apply the reassignment formulae to the Constant-Q Transform. When evaluating the energy in x at time t and frequency ν k , the process is the following. We first compute the magnitude of the Constant-Q Transform at pt, ν k q. We then compute the Constant-Q Transform with the ramp-multiplied windows. We finally compute the Constant-Q Transform with the g-convolved windows. Thanks to equations ( 8) and ( 13), we obtain the indexes tt,ν k and νt,ν k . In the end, the energy computed at point pt, ν k q is additively reassigned to the point p tt,ν k , νt,ν k q in the time-frequency plane.

Figure 2 shows the reassigned Constant-Q Transform of the same piano sound as in Figure 1. One can see that the temporal resolution has much improved in the second figure. The typical artifacts of the widening energy bursts in low frequencies seem to be accurately handled in the reassigned spectrogram.

Besides, one can see that the frequency resolution has also improved. Even if the original tool's design ensures a suitable frequency resolution for sound applications, we can observe that the reassigned transform has much more localised (in terms of frequency) energy bursts.

Conclusion

In To illustrate the applicability of the derived framework as well as the benefits of time-frequency reassignment, the obtained results have then been applied to the Constant-Q Transform. This transform, commonly used in audio processing, has geometrically-spaced frequency bins with a variable resolution. This makes it a well adapted tool for the analysis of sound. This frequency setup, however, comes at the cost of time resolution, which has the serious drawback of being variable with frequency. Once reassigned though, the representation provides a much better resolution, in time and frequency. A Python implementation of this enhanced tool is available on the authors' homepages2 . It is our belief that the proposed reassigned Constant-Q Transform has a strong potential for numerous audio processing applications, including pitch/multipich estimation, onset detection and chords transcription.

Figure 1 :

 1 Figure 1: Spectrogram of a piano sound, sampled at 44100 Hz, obtained by means of the traditional Constant-Q Transform. The X axis is the time, it is labelled with the integer indexes of the sampled signal. The Y axis represents the frequency. It is noticeable that the latter follows a logarithmic scale.

Figure 2 :

 2 Figure 2: Spectrogram of the same piano sound as in Figure 1 obtained by means of the reassigned Constant-Q Transform.

  this work we have derived a formal framework for the time-frequency reassignment of discrete time signals. The principle of time-frequency reassignment is to express the time-frequency transform of interest as the summation of a real-valued mass function in the time-frequency plane. The energy is then reassigned to the centre of mass of this function. The first derivations carried out in this paper result in an explicit expression of the mass function which is valid for a wide variety of time-frequency transforms. Based upon the latter, an expression of the coordinates of the centre of mass could be obtained. Interestingly, it consists of very simple combinations of three time-frequency transforms, each of them using a specific window. This result makes time-frequency reassignment a rather generic, simple and easily computable technique.

googledrive.com/host/0B1DUW3X2T63neXhKWmxrTW5nUlE/ressources.html

googledrive.com/host/0B1DUW3X2T63neXhKWmxrTW5nUlE/ressources.html

We have:


Toward variational data assimilation for coupled models: first experiments on a diffusion problem R émi Pellerej, Arthur Vidard and Florian Lemari é remi.pellerej@imag.fr Inria, Univ. Grenoble-Alpes, CNRS, Laboratoire Jean Kuntzmann, F-38000 Grenoble

Context

Ocean-atmosphere coupled models have a key role in weather forecast nowadays The coupling methods may severely impact the model solution. An exact solution of the coupling problem can be obtained using a Global-in-Time Schwarz method (Lemarié et al.

[2014])

The initialisation of coupled models also has a major impact on the forecast solution [START_REF] Mulholland | Origin and impact of initialization shocks in coupled atmosphere-ocean forecasts[END_REF])

Few coupled DA methods started to be developed [START_REF] Smith | Exploring strategies for coupled 4d-var data assimilation using an idealised atmosphere-ocean model[END_REF], [START_REF] Laloyaux | A coupled data assimilation system for climate reanalysis[END_REF]...) for coupled systems, and showed promising results

Our approach

The dynamical equations of our system are coupled using an iterative Schwarz domain decomposition method [START_REF] Gander | Schwarz methods over the course of time[END_REF])

We are using variational DA techniques, which require minimization iterations and we are looking to take benefit of the minimization iterations to converge toward the exact solution of the coupling problem: the minimisations iterations substitute the Schwarz iterations Three general variational DA algorithms, are presented here and applied to a simple coupled system [START_REF] Pellerej | Toward variational data assimilation for coupled models: first experiments on a diffusion problem[END_REF])

Model problem and coupling strategy

Let us define two models on each space-time domain Ω d × [0, T ] (d = 1, 2), with a common interface Γ = {z = 0}.

Problem: How to strongly couple the two models at their interface Γ ? ⇒ We propose to use a global-in-time Schwarz algorithm [START_REF] Gander | Schwarz methods over the course of time[END_REF])

Ω 2 Ω 1 Γ z = 0 z = L 2 z = -L 1 L 2 u 2 = f 2 L 1 u 1 = f 1
For a given initial condition u 0 ∈ H 1 (Ω 1 ∪ Ω 2 ) and first-guess u 0 1 (0, t), the coupling algorithm reads

     L 2 u k 2 = f 2 on Ω 2 × T W u k 2 (z, 0) = u 0 (z) z ∈ Ω 2 G 2 u k 2 = G 1 u k-1 1 on Γ × T W      L 1 u k 1 = f 1 on Ω 1 × T W u k 1 (z, 0) = u 0 (z) z ∈ Ω 1 F 1 u k 1 = F 2 u k 2 on Γ × T W (1)
F d and G d are the interface operators, k is the iteration number, T W = [0, T ], and

f d ∈ L 2 (0, T ; L 2 (Ω d )) is a given right-hand side
At convergence, this algorithm provides a mathematically strongly coupled solution which satisfies

F 1 u 1 = F 2 u 2 and G 2 u 2 = G 1 u 1 on Γ × T W
The convergence speed of the method greatly depends on the choice for F d and G d operators, and the choice of the first-guess

Classic data assimilation

Let us introduce the classic cost function for variational data assimilation in the uncoupled case, for a domain

Ω d x 0,d = u 0,d (z) = u 0 (z), z ∈ Ω d (d = 1, 2) is the controlled state vector J U cpl (x 0,d ) = J b (x 0,d ) x 0,d -x b d , B -1 (x 0,d -x b d ) Ω d + J o (x 0,d ) T 0 y -H (x d ) , R -1 (y -H (x d )) Ω d dt (2)
where • Σ is the usual Euclidian inner product on a spatial domain Σ.

Toward a coupled variational data assimilation

If the DA process is done separately on each subdomain, the initial condition u 0 = (x a 0,1 , x a 0,2 ) T obtained on Ω does not satisfy the interface conditions. The interface imbalance in the initial condition can severely damage the forecast skills of coupled models [START_REF] Mulholland | Origin and impact of initialization shocks in coupled atmosphere-ocean forecasts[END_REF])

Objective: properly take into account the coupling in the assimilation process Full Iterative Method (FIM) x 0 = u 0 (z), z ∈ Ω We iterate the models till convergence of the Schwarz algorithm (k cvg iterations) The first-guess u 0 1 in ( 1) is updated after each minimization iteration

J F IM (x 0 ) = J b (x 0 ) + T 0 y -H (x cvg ) , R -1 (y -H (x cvg )) Ω dt (3)
where

x cvg = (u k cvg 1 , u k cvg 2 ) T

Truncated Iterative Method (TIM)

x 0 = (u 0 (z), u 0 1 (0, t)) T The Schwarz iterations are truncated at k max iterations Extended cost function (misfit in the interface conditions) [START_REF] Gejadze | On a 2d 'zoom' for the 1d shallow water model: Coupling and data assimilation[END_REF])

J T IM (x 0 ) = J b (x 0 ) + T 0 y -H x trunc , R -1 (y -H x trunc ) Ω dt + J s (4)
where

J s = α F F 1 u 1 (0, t) -F 2 u 2 (0, t) 2 [0,T ] + α G G 1 u 1 (0, t) -G 2 u 2 (0, t) 2 [0,T ] with a 2 Σ = a, a Σ and x trunc = (u k max 1 , u k max 2 ) T

Coupled Assimilation Method with Uncoupled models (CAMU)

x 0 = (x 0,1 , x 0,2 ) T with x 0,d = (u 0 | z∈Ω d , u 0 d (0, t)) We suppress the coupling between both models The cost function for the CAMU is

J CAM U (x 0 ) =    2 d=1 (J b (x 0,d ) + J o (x 0,d ))    + J s
(5)

Algo Control vector # of coupling iterations extended cost function

Adjoint of the coupling

Coupling The originality of these algorithms is the use of a Schwarz algorithm to couple our models jointly to the DA process with an extended cost function.

FIM (u 0 (z)) k cvg no yes strong TIM (u 0 (z), u 0 1 ) T k max yes yes ∼strong CAMU (u 0 (z), u 0 1 , u 0 2 ) T 0 yes no weak

Application to a 1D diffusion problem

Previous algorithms are applied on a 1D linear diffusion problem. We consider: Computational Cost (Model Itegrations)

L d = ∂ t + ν d ∂ 2 z ν 1 = ν 2 the diffusion

Single column observation experiment results

Algo

α G α F k max # of
10 -1 10 0 10 1 J FIM TIM α G =0, k max =5 TIM α G =0.01, k max =1 CAMU α G =0.001, α F =4 0
J o FIM TIM α G =0, k max =5 TIM α G =0.01, k max =1 CAMU α G =0.001, α F =4 0

Interface imbalance indicator

FIM TIM α G =0, k max =5 TIM α G =0.01, k max =1 CAMU α G =0.001, α F =4

Conclusions and perspectives

In the framework of an iterative coupling, we set up few data assimilation algorithms. Adding a physical constraint on the interface conditions in the cost function can have a beneficial effect on the performance of the method and allow to save coupling iterations An approach which only requires the adjoint of each individual model but not the adjoint of the coupling showed promising results

The methods are very sensitive to the parameters choices We only test the algorithms on a simple linear problem

Perspectives

Algorithm convergence and conditioning problem when J S is part of the cost function will be studied Since the objective is to apply such methods to ocean-atmosphere coupled models, increasingly complex models including physical parameterisations for subgrid scales, and non-linearities will be considered

  coefficients in each subdomain F d = ν d ∂ z and G d = Id the interface operators on Γ (Dirichlet-Neumann) u d (z, t) = U 0 4 e -|z| α d 3 + cos 2 3πt τ on Ω d × T W the analytical solution Single column observation experiment: Observations are available in Ω \ {Γ} at the end of the time-window (i.e. at t = T ) We define the interface imbalance indicator, equal to J s with α G = 0.01 and α F = 40

Figure 1 :

 1 Figure 1 : Evolution of different terms with respect to number of model iterations for few configurations

Table 1 :

 1 Overview of the properties of the coupled variational DA methods described

Table 2 :

 2 Results obtained for the three coupled variational DA methods

	0	50	100 150 200 250 300 350 400


Minh Khue

Phan Tran
email: m.phantran@genious.com

Philippe Robert

Francois Bremond

An assistive system to improve game

Keywords: assistance, serious game, older adults, support system, engagement, performance

An assistive system to improve game usability for patients with cognitive disorders.

M.K. PHAN TRAN, P. ROBERT, F. BREMOND. Gerontechnology Nice-2016;doi:10.4017/ gt.2016.15.s.758.00

Purpose

With the ageing population, Alzheimer's disease and related pathologies are a major public health challenge. In parallel, new technologies of information and communication take an increasingly important role in our daily life and can be a support for both the evaluation and the direct support of users. Serious games 1 are computer applications, whose intention is to combine with consistency, both serious aspects (Serious) as teaching, learning, communication, rehabilitation, with playful aspects from the video game (game). Under the AZ@GAME project 2 winner AAP e-health #1 of Future Investments, games have been developed with the aim of stimulating cognitive and physical abilities of the patient. One of perceived problems concerns the engagement of patients to understand and to practice these games. Beside of lack of game culture and limited acceptance towards new technologies, memory and cognitive impairments provoke many difficulties of use and influence negatively on their game performance. Thus, engagement of patients decreases progressively and patients cannot benefit all the positives effects of games.

Method

We present an automated support system (Figure 1) containing a virtual avatar that can interact and provide help in different appropriate moments thanks to engaging strategies based on patients activities and game states. The system has been tested in three experiments with patients (94 participants including 51% man and 49% woman with an average age of 76 years old) distributed in three cognitive states (Mnesic Complaints, Mild Cognitive Impairment and Alzheimer's disease) who participated in a training session using the game to boost concentration attention. The games were offered with or without the help of the system.

Results & Discussion

The first results [START_REF] Tran | Comment intéresser les personnes âgées aux Serious Game? [How to interest older adults in the serious game ?[END_REF] have remarked that most patients can interact with the avatar and accept its invitation to play games. The second results [START_REF] Tran | Assistance for Older Adults in Serious Game using an Interactive System. The Games and Learning Figure 1 -Support system structure Alliance Conference[END_REF] have emphasized usability of the system. Patients have better performance when playing with support of the system than when playing alone, especially among Alzheimer's patients. The last results [START_REF] Tran | A virtual agent for enhancing performance of older adults with dementia in Serious Game[END_REF] have confirmed efficiency of the system compared with a therapist. Alzheimer's patient can have similar performance thanks to the system, as compared to playing with a therapist. All of three experiments have been realized in a hospital in a short time. A longitudinal study in other real contexts (e.g. living home) is necessary for analyzing the use of the system with patients in long-term.