Ling Yin

Julien Deantoni

Frédéric Mallet

Robert De Simone

Schedulability analysis by exhaustive state space construction: translating CCSL to transition-based Generalized Büchi Automata

Keywords: CCSL, Scheduling, Automata

In previous work we defined a language (CCSL) made to express real-time temporal scheduling constraints. It uses the notion of partially independent logical clocks (or time threads), of which seemingly physical discrete time is a special case, hence the name Clock Constraint Specification Language. Constraints can represent (asynchronous) causality and precedence relations, or (synchronous) simultaneity ones. A solution to a set of such constraints is called a schedule, as it brings back all the clocks down to a single, totally ordered discrete time (and thus allows to timing events a schedule slot). In the current paper we study the formal semantics of CCSL, by translation into specific automata recognizing infinite sequences of clock steps. We consider the appropriate acceptance criteria for omega-languages (transition-based generalized Büchi), and motivate this choice. We also study how the automata can be minimized, by removing useless states (which is more than the usual check for emptiness of the whole language in classical modelchecking approaches). We feel our work builds a definite link between the scheduling and the formal semantic model-checking communities.

Schedulability analysis by exhaustive state space construction: translating CCSL to transition-based Generalized Büchi Automata

Résumé : In previous work we defined a language (CCSL) made to express real-time temporal scheduling constraints. It uses the notion of partially independent logical clocks (or time threads), of which seemingly physical discrete time is a special case, hence the name Clock Constraint Specification Language. Constraints can represent (asynchronous) causality and precedence relations, or (synchronous) simultaneity ones. A solution to a set of such constraints is called a schedule, as it brings back all the clocks down to a single, totally ordered discrete time (and thus allows to timing events a schedule slot). In the current paper we study the formal semantics of CCSL, by translation into specific automata recognizing infinite sequences of clock steps. We consider the appropriate acceptance criteria for omega-languages (transition-based generalized Büchi), and motivate this choice. We also study how the automata can be minimized, by removing useless states (which is more than the usual check for emptiness of the whole language in classical model-checking approaches). We feel our work builds a definite link between the scheduling and the formal semantic model-checking communities.

Mots-clés : CCSL; Scheduling; Automata

Introduction

Previously we defined a formalism to express time relations between various events and operations of a given (discrete) system, named the Clock Constraint Specification Language [START_REF] André | Syntax and semantics of the clock constraint specification language (ccsl)[END_REF]. ccsl is meant to capture some essential time and timing notions as found in various theories of timed concurrent processes, both from the synchronous and asynchronous sides. In this paper, we explore the use of ccsl for scheduling requirements modeling and schedulability analysis. Historically, ccsl is based on a notion of "discrete clocks", which deserves (at least) a bit of explanation to avoid confusion. In our case a clock is not a physical device that provides regular physical time intervals (as a watch). It is instead an abstract measure unit, possibly user defined, by which durations and dates may be set or measured, as long as a sequence of occurrences ("ticks") on that clock can be observed or produced. In that sense a logical clock could also be named "logical time base", and largely corresponds in fact to just an abstract event (considered as a sequence of event occurrences), simply stressing more its temporal relevance and role in the design. Besides, logical clocks need not in any case be related to regular physical realities, but they can, and often do. To provide a simple example, a modern processor execution cycle, in these days of low-power and thermal budgets, is often quite irregular (it can even be paused on hold at times) with respect to physical time; nevertheless it is still a useful logical time unit when accounting for software complexity. Numerous other examples can be made up in the case of embedded systems, where contextual environment may produce meaningful events that punctuates the system behavior.

Based on clocks, ccsl introduces formal constructs, called clock constraints, both to define new clocks as functions of existing ones and also to specify relations between existing clocks as well. A significant list of primitive ccsl clock constraints will be presented in the paper, accounting for the language expressiveness and scheduling model description.

The logical clocks and clock constraints give ccsl the power to express most of the tasks/events models in (both real-time and non real-time) scheduling theories. A set of clocks and constraints among them can model the scheduling requirements (including functional requirements, platform features, resource limitations, extra timing requirements and etc), and will form a ccsl specification of the system. Solving the specification constraints amounts to computing a schedule of the system.

Under ccsl context, a schedule is a possibly infinite sequence of steps; for each step (index), several clocks (or maybe one) are assigned to tick simultaneously at that step. Then a schedule can be regarded as an infinite words on powerset of clocks. The general scheduling issues, like finding valid schedules or finding if a system under specific scheduling requirements is schedulable, can be related to formal language and automata theory on infinite words. Modeling scheduling concerns in ccsl allows the use of model checking technique and formal analysis to solve issues taken from the scheduling theory.

We do schedulability analysis by translating a ccsl specification into a transition-based Generalized Büchi automata. Proper acceptance criteria is defined to recognize valid schedules. Instead of finding an individual optimal schedule under a specific consideration, the automata representation represents all schedules that satisfy the requirements, and can be used as a basis for further analysis. The schedulability of a specification is equivalent to checking the emptiness of recognizable language of a corresponding automata. We also provide an algorithm to construct a trim canonical version for a given transition-based generalized automata. The trim automata can be used to guide later on scheduling.

The paper is constructed as follows: We introduce ccsl informally in section 2. Then we discuss its relation with existing scheduling practices and theories and the scheduling issues we want to deal with in section 3. In section 4 we show the translation from ccsl to transition-based Generalized Büchi automata. Automata based solutions to the scheduling issues are provided in section 5. Then we conclude this paper, discuss related work and our future plans in section 6.

An introduction to ccsl

Tick and clock ordering relations

As already mentioned, clock ticks amount to event occurrences. Two event occurrences e and e ′ may temporally be: simultaneous/coincident, denoted as e ≡ e ′ , meaning that the two occurrences will occur in the same time step; or one may precede the other, denoted as e ≺ e ′ , introducing a temporal ordering between them. We also note e#e ′ for not (e ≡ e ′ ), and e e ′ for (e ≺ e ′ or e ≡ e ′ ). Note that causality intentions in design will in general allow relation , when instantaneous/combinatorial causality is also legible.

At the clock level (again, a clock being a sequence of ticks/instants), we find that two natural ordering relations extend the former. First, based on synchronous notions, a subClock relation can be defined: we note c ⊏ c ′ iff ∀i, j :

c i c j ∃k, l : (c i ≡ c ′ k ) ∧ (c j ≡ c ′ l ) ∧ (c ′ k c ′ l
), c i being the i th tick of clock c; this definition means that every tick of c has to be made on a tick of c ′ , but not always the other way around, so there are potentially more ticking instants in c ′ . Moreover, this relation is order preserving. Second, based on asynchronous notions, faster than relation can be defined:

we note c ≺ c ′ iff ∀i, c i ≺ c ′ i ; similarly, c c ′ iff ∀i, c i c ′ i (

strictly vs weakly faster).

The subClock relation will in practice be used to define new clocks from existing ones, with a subsampling mechanism. Often, the latter will be a regular pattern, so that periodic subsystems, or even periodically patterned one, will fall into this design style. Also, when original specifications are themselves not synchronously patterned, the results of scheduling methods are often computed of that nature, so that subclock relations with regular patterns will also be a mean to describe (full or partial solutions to) the system ordering; having a way to deal with partial solutions inside the syntax of the specification itself is useful, as it allows solving approaches by incremental time refinements.

The faster than relations can a priori lead to unbounded drifts between clocks (to ensure that the second clock does not tick too soon, we need to monitor how many ticks the first one has performed "in advance", just so it does not drop below 0, and this requires an integer in simulation). In frequent cases the specification will provide an interval between the clocks to bound drifts.

So far we have implicitly assumed that clocks were infinite sequences of ticks. Often, for technical reasons, one has to assume they can also be finite. So we will let |c| denote the length of c, with the convention that |c| = ∞ when the clock is infinite. The previous definitions of clock ordering have to be adjusted so that they only apply where ticks are actually defined. This introduces no real change for the subClock relation, but it should be noted that, in the case of the faster than relation, c ′ i may occur without a prior c i , if i > |c|. This defeats the notion of causality between c and c ′ being encoded as a time ordering relation. To avoid such problem, when later constructs represent causality as such a time ordering, they will always rely on the fact that clock lengths will carefully match. Later, in dynamic operational semantic definitions, we shall use the predicate T erm(c) to mean that the proper termination of clock c, that is it has ticked up to its limit (whose value is generally dynamically reached).

Primitive ccsl clock functions

We now provide a set of primitive functions to build new clocks out of others. They are mostly of synchronous nature (subClock constructed).

Inria

The main function, filteredBy, will use an explicit mask to decide on each new tick (of the superclock) whether it should be preserved or discarded in the subclock definition. Such a mask consists of a binary word (with "1" meaning keep and "0" discard) of length at least equal to the length of the superclock. In practice, for syntactic reasons, the mask will have to be an ultimately periodic pattern, consisting of some initialization prefix part, and a steady/repeated part. Simpler forms can be extracted using suitable masks, like delay_n, with n a constant number of ticks from the superclock; it simply discards the n first ticks of its input clock. Similarly await_n only ticks once, on the original n th tick of its input clock, or until_n, which selects only the first n ticks of its input clock. The upto operator takes two input clock arguments, and ticks with the first clock as long as the second one has not started. The concat sequential composition also takes two input clocks. It will tick with the first until it has reached its length limit, then starts ticking as the second (when the first is terminated). Of course in practice this will require that checking for the first clock lifespan (or its length) is feasible, at least at run time. In simulation it can also mean that the property can be enforced by the simulation policy, but then it has to be checked that it is valid in the future. Union (and intersect) functions work directly at the level of steps, and create new clocks that tick when at least one (or both) input clocks do. The sampledOn function is a mix of synchronous and asynchronous spirit. It takes two input clocks and ticks synchronously with the second, but only in case the first input clock ticked since the previous tick of the second. So, its role is essentially to record whether there was a tick of the first clock since "last time". Inf (and sup) functions also take two clocks as input, and create new clocks that is faster (or slower) than both of the input clocks. The i th tick of the new defined clock is coincident with the earlier (or later) of the i th tick of input clocks.

Primitive ccsl clock relations

Clock relations are provided to specify relations between two existing clocks. Both synchronous and asynchronous aspects are considered. The fundamental synchronous relation is the subClock discussed above, restricting that the subclock cannot tick if its superclock does not. It is one direction specified, and can be extended into two directions, that is the equality relation stating the equality between the two clocks, one ticks iff the other one ticks. The exclusion relation is the opposite of equality. Two clocks are exclusive means that they never tick at the same step. The basic asynchronous relations are the two versions of faster than discussed in last subsection, the strictPre ≺ and causes . For the strict one, the second clock cannot tick when the index difference is 0. While for the weak one, when their index different is 0, the second clock can tick if the first one ticks. They both limit the drift between the two clocks in the range of [0,∞]. BoundedDiff_i_j bounds an integer interval [i,j] to the drift of the two involved clocks, denoted as i ≤ lef tright ≤ j. It can be viewed as an extension of faster than, not restricting which clock is faster, but stating their index difference cannot exceed the range. The relation alternate indicates the alternate occurrences of instants of the two clocks, starting from the left one.

Considering clock termination: the filteredBy operator preserves termination, and may create a finite clock from an infinite one (if the mask is of finite length); The delay_n, await_n and upto operators create finite clocks (unless the second clock of upto never starts); concat will be of real use only if its first clock argument terminates; finite union clock requires the two clock arguments to be finite, while intersect may produce a finite clock even both arguments are infinite; sampledOn is finite if one of its argument is; by convention inf provides infinite clocks as soon as one of their arguments is infinite, while sup provides finite one as soon as one argument is finite. A subclock of a finite clock is finite (while subclock of an infinite clock may be finite too); a clock in bounded drift of a finite one is finite (a clock with a finite faster than clock is finite as well).

3 The issue: schedulability analysis through formal translation to a proper operational semantic domain

We want to play the connection between the scheduling theory and the automata theory through ccsl, modeling the issues taken from scheduling theory by ccsl, and solving them by translating a ccsl specification into a proper automata with acceptance criteria recognizing specific infinite words. First we discuss the relation between ccsl and existing scheduling theories. Then we present the scheduling issues in ccsl context and provide the overview of our automata based solutions to the issues.

Relations between ccsl and existing scheduling theories

Real-time scheduling theories rely on task models that supposedly abstract real applications. Originally they were rather simple (e.g., independent periodic tasks only for Rate Monotonic Analysis). Always more sophisticated models now appear in the literature. They are all based on numerous distinct parameters, providing numerical constraint values for timing aspects (dispatch time, period, deadline, jitter drift,...). Tasks are considered as iterations of jobs (or jobs as instances of tasks). In our view, the successive timing values for characteristic feature of successive jobs can be seen each as a logical clock, and the time constraint relations between such clocks are usually expressed as simple equalities and bounded inequalities that fall well into the range of ccsl constructs descriptive power. As an example, the notion of sporadicity which states that tasks are invoked at random but with a minimum inter-invocation interval, can be handled by ccsl, as for instance, c = (p filteredBy (001

) ω ) sampledOn d, 0 ≤ d -c ≤ 1 and d ≺ (p wait_1
), states at least three ticks of p between two ticks of d; the fact that p ticks are regular in physical time is ignored in ccsl modeling. Note that logical clock in ccsl can also be used to represent realities of different nature, more ad-hoc to the application, with added intuitive expressiveness. For instance, the exclusive access to a critical resource is modeled by exclusion between corresponding clocks, the size limit of a buffer (3 as an example) can be expressed by a special case of boundedDiff relation, making 0 as the left bound and buffer size as the right bound between write and read clocks, 0 ≤ writeread ≤ 3. With such formal specification, ccsl is not able to say anything about the way to find efficient resulting schedule(s) for specific class of problems modeled that way (for instance, Earliest Deadline First plans to schedule first the job with least time remaining, this remaining time being dynamically computed due to past activities). Further comparisons of ccsl and usual real-time scheduling models would be out of the scope of the current paper. But it should still be mentioned that models closer to ccsl, as for instance AADL and AutoSar-related EAST-ADL models, are also now making their way to mainstream scheduling concerns. Classical (non real-time) scheduling, on its side, provides generally models where the initial constraints are less on timing and more on dependencies (as for program instruction scheduling

Automata based solution to ccsl scheduling issue

A ccsl specification formally describes expected requirements and represents a set of valid schedules, i.e. a set of valid sequences of steps; where for each steps (index), clocks are said to tick or not. Solving the specification constraints consists in assigning clock ticks to steps, and consequently to provide one schedule of the system. The purpose of scheduling is to assign resources and time slots to tasks so that all tasks are accomplished under some constraints (e.g., fair share of resources or time, using as little time slots as possible, etc). In this sense, a valid schedule should be a sequence of steps where every clock ticks as far as it meant to (infinitely often or up to its proper termination). Specifications may become unschedulable because of contradictory constraints, or more subtly because it becomes unfeasible to pack a number of ticks in a given time interval to get a valid schedule. We provide a typical example below, which will reveal a number of concerns involved in ccsl schedulability.

Example 1. For three clocks a, b and c, where b is a subclock of a, c is a wait 3, b is strictly faster than c and their index difference is bounded in [0,2]:

b subClock a, c = a wait 3, b ≺ c, 0 ≤ b -c ≤ 2
The automata like representation in Figure 1 is informally used here to give readers an intuitive view of the issue. Formal semantics is given later. The global behaviors of a specification (on the right in Figure 1) have to respect to all its constraints (on the left in Figure 1). One major problem is that during the construction of the schedules, some choices can lead to invalid schedules. For instance, what is expected from our example is: At least one tick of b should occur before the third tick of a, then c ticks together with the third tick of a and then terminates. Clock b terminates when the index difference between b and c reaches 2, and a ticks infinitely. However, none of the constraints forbids a to tick two times without one b. In state < s, t 2 , r 0 > it becomes unfeasible to fire any clock with respect to the conjunction of constraints: every clocks get halted (cannot tick anymore though it is not proper terminated). The run < s, t 0 , r 0 > {a} --→< s, t 1 , r 0 > {a} --→< s, t 2 , r 0 > consequently cannot produce a valid schedule. To decide the schedulability of a specification and to recognize valid schedules from specification runs, we provide a state-based operational semantics for ccsl, ccAutomata. Such ccAutomata will represent all possible runs of a specification, and thus all valid schedules of its clock ticks. Recognizing valid schedules from runs can be achieved by defining acceptance criteria for infinite words on ccAutomata. The schedulability issue, which demands that one such valid schedule exists, will thus amounts to the emptiness checking problem.

4 Formal translation of CCSL into (Generalized Transition Büchi) ccAutomata

We will provide the formal translation of CCSL clock constraints into appropriate state-based models using operational semantics. We will first describe how the "stateful" features of operators expand into states, providing the general skeletons for the expanded models; then we will consider the acceptance criteria issue, so that exactly valid schedules are recognized. In the next section we will show how the resulting models can be used to solve the schedulability issues.

Acceptance criteria for languages of infinite words have been amply studied in the past. This gave rise to recognition patterns such as Büchi (repeated states), certainly the most famous, but also Muller, Rabin, and Street types of acceptance criteria. The Büchi acceptance criterion also supports variants, such as transition Büchi (rather than state Büchi), or Generalized Büchi criteria (in which several sets of repeated states are provided and must be satisfied collectively).

We shall show how CCSL naturally leads to transition Büchi automata translation for individual constraints, and why state Büchi is less natural; we shall also see how the generalized version is needed when specification are composed of several distinct clocks (and there will be one recognition set per clock). It is known from the literature that both variants can be expanded at some cost into plain (state) Büchi automata, but we will consider how direct analysis may prove to be more efficient.

For homogeneity reasons we shall also see how termination states and finite clocks can be unified with infinite ones (by noticing infinitely the termination). We will also have to rely on stuttering transitions so that the construction of a global automaton from the various constraints (and their individual automaton form) can be considered as a synchronous product, even though some may idle in it. These can be seen as "classical tricks" in the model-checking community, but they are to be applied on models and theories meant to the (distinct) scheduling theory community.

State-based operational semantics of CCSL-ccAutomata

Definition 3. A ccAutomata is a tuple A = (S, Clocks, T erms, I, T ) where • S is a nonempty set of states, and I ⊆ S is the set of initial states.

• Clocks is a finite nonempty set of clock names, which can be regarded as sorts of the automata.

• T erms is a set of clock termination predicates T erms = {T erm(c)|c ∈ Clocks}.

• T ⊆ S × 2 Clocks × 2 T erms × S defines the transition relation, each transition is labeled by a subset of Clocks that tick simultaneously in that step, and a subset of Terms indicating the proper terminations of corresponding clocks.

A transition (s, C, T m, s ′ ) ∈ T is also denoted as s

C,T m ----→ T s ′ or simply s C,T m ----→ s ′ if clear from context. A ccAutomata is deterministic iff ∀s ∈ S, (s C,T m ----→ s 1 ∧ s C,T m ----→ s 2 ) ⇒ (s 1 = s 2 ).
A path of A is a sequence of transitions, either finite or infinite, ρ = s 0

C0,T m0 -----→ s 1 C1,T m1 -----→ • • • C k-1 ,T m k-1 ---------→ s k • • • .
A run is a path from an initial state, and its ticking trace is the series of ticking labels: ----→ s since no constraint is violated or changed if no clock ticks. Function : Clocks → {F in, Inf, F ree} is provided to user to specify that a clock is supposed to be finite ( (c) = F in), infinite ( (c) = Inf ), or either ( (c) = F ree). It needs to be in consistent with the constraint semantics, for instance, for c = c 1 await_3, (c) must be F in. The function can be used to simplify the ccAutomata encoding (for scheduling concerns). If we know a clock is specified to be infinite, the the transitions labeled by its termination can be removed since termination of that clock is not expected (transitions labeled by that will not lead to valid schedules). We will show how the is used to simplify the encoding: for each primitive clock constraint, we first show its ccAutomata encoding for the most natural case (usually the case that (c 1 ) = (c 2 ) = Inf ), then add corresponding termination transitions and states for finite specified clocks. These transitions and states are colored for clear view, where "red" means "un-conclusive", some clocks have terminated while others haven't, and "green" means "terminal" that every clock has terminated.

C 0 ; C 1 ; • • • C k ; • • • .
For simplicity, termination predicts are denoted by ↓ clockname, and stutter transitions are ignored in the figures.

• filteredBy, delay_n,await_n,until_n

Constraint c = c 1 filteredBy w builds a sub clock c of clock c 1 according to the binary word w. Usually w is written as u(v) ω , where u is the transient part and v the periodic part with length m, v is the periodic part with n bits, the power ω means v is repeated an unbounded number of times. When the bit of w is 1, the instant of c 1 is selected, otherwise discarded.

(c 1 ) = F in ⇒ (c) = F in.
The most natural case is (c) = (c 1 ) = Inf shown in Figure 2(a), where Constraints delay_n,await_n,until_n can be regarded as special cases of filteredBy, using 0 n (1) ω , 0 n 1(0) ω and 1 n (0) ω as masks respectively.

δ i c = {c} if u i = 1, otherwise δ i c = ∅, similarly µ i c = {c} when v i = 1, µ i c = ∅ when v i =
• upto c = c 1 upto c 2 creates a finite clock unless c 2 never starts c 1 is infinite. Clock c terminates properly when c 2 starts.

• concat Constraint c = c 1 concat c 2 relies on dynamic terminations, c behaves like c 1 until its termination and then behaves c 2 . Usually c 1 should be finite, otherwise c is just equal to c 1 .

• union,intersection

(c) (c 1 ) = Inf, (c 2 ) = F in 1↓234 ↓254↓26 123426 1 ↓ 2 3 6 12425641↓236 1↓256 123426471↓256 123425426 1 2 5 6 1 ↓ 2 3 4 ↓ 2 5 4 ↓ 2 6 1 ↓ 2 5 3 ↓ 2 5 4↓ 2 6 (d) (c 1 ) = (c 2 ) = F in
(c 1 ) = F in, (c 2 ) = Inf ) strict: ∀i, ∃j, k, c2[j -1] c1[i] ≺ c2[j], c[k] ≡ c2[j] weak: ∀i, ∃j, k, c2[j -1] ≺ c1[i] c2[j], c[k] ≡ c2[j]
Please check the difference reflected in ccAutomata encodings: c 1 exclusion c 2 specifies that they cannot tick together. Their terminations do not affect each other.

• BoundedDiff, strictPre, causes

If m ≤ c 1 -c 2 ≤ n (i, j ∈ N ), then (c 1 ) = (c 2 )
. Similar to sup and inf, the index difference δ = i at state s i . The ccAutomata encoding for Inf case is shown in Figure 14(a). For F in case, at s i , if c 2 (resp. c 1 ) terminates, c 1 (resp. c 2 ) can at most ni times before termination.

Constraint causes is a special case of BoundedDiff that m = 0, n = ∞. strictPre is the same as causes except that at s 0 the slower clock cannot tick.

ccAutomata semantics of CCSL specifications-composition of individual ccAutomatas

A specification is a conjunction of constraints so its resulting behavior have to respect all its constraints. Transitions from two ccAutomatas can be composed only if there is no conflict; conflict here means that a common involved clock ticks in one transition while does not in the other, or it terminates in one transition while stays alive in the other. Definition 2. Given a set of ccAutomatas, A i = {S i , Clocks i , T erms i , I i , T i }, where i ∈ ζ = {1, ..., n}, their composition is a proper synchronized product, ||{A i } i∈ζ = {S, Clocks, T erms, I, T }, where

S = S 1 × • • • × S n ; I = I 1 × • • • × I n ; Clocks = Clocks 1 ∪ • • • ∪ Clocks n ; T erms = T erms 1 ∪ • • • ∪ T erms n ; ∀i, j ∈ ζ(i = j), s i Ci,T mi -----→ Ti s ′ i , s j Cj ,T mj -----→ Tj s ′ j , ∀c ∈ Clocks i ∩ Clocks j , c ∈ C i ⇔ c ∈ C j , T erm(c) ∈ T m i ⇔ T erm(c) ∈ T m j (s 1 , ..., s n ) C1∪•••∪Cn,T mi∪•••∪T mn -----------------→ T (s ′ 1 , ..., s ′ n )

Acceptance criteria for valid schedules

A valid schedule is the ticking trace of a run, in which every clock ticks as far as it meant to. That is to say, every infinite specified clock ticks infinitely often and every finite specified clock ticks up to its proper termination.Thanks to the stutter transitions and the preservation of clock terminations, a clock ticks up to its proper termination can be noticed by the infinite visit of its termination predict. A run ends in a state s f where every clock is terminated properly can always be extended to infinite by s f ∅,T erms ------→ s f . Then acceptance criteria for infinite words can be used to recognize valid schedules.

We say path ρ passes infinitely often through a state (or transition in transition-based style) if there are infinitely many occurrences of it in ρ. The set of states (or transitions) which occur infinitely often in ρ is denoted by Inf (ρ). Acceptance criteria in Büchi recognizing pattern is a set of repeated states, i.e., a state-based Büchi ccAutomata (denoted as state BA for short) is a tuple (A, F ), where A is a ccAutomata, F ⊆ S is a finite set of accepting states. A run of A is accepted if Inf (ρ) ∪ F = ∅. For instance, for the ccAutomata of c 1 subClock c 2 with (c 1 ) = (c 2 ) = Inf (Figure 13(a)), the acceptance criteria would be F = {s1}. The fairness between the two clocks cannot be guaranteed. A run with infinite ticks of c 2 but no tick of c 1 would be accepted, while it can not produce a valid schedule. To solve this problem, we can either expand s 1 into two sates (s 11 , s 12 in Figure 16) and make only s 12 accepting, or use transition-based Büchi automata (tBA) instead,

F = {s 1 {c1,c2},∅
------→ s 1 }. Forcing state-based accepting sets would in fact amount to reiterating locally the general expansion from tBA to BA. But counting on transitions comes more natural.

When composing local tBAs obtained from individual constraints into a global one to represent schedules from the whole specification, we want to avoid infinite runs where some clocks are forbidden forever while non terminated. The problem is that such runs are actually accepted by the naive acceptance criterion. For this type of problems, literature has come up in the past with the notion of Generalized BA, in which a collection of acceptance conditions would be provided, each needing to be checked independently. In our case we shall associate such a condition to every specified finite or infinite clock.

C,T m ----→ s ′ ∈ Inf (ρ), T erm(c) ∈ T m}, ρ is accepted if ∀f ∈ F, f ∈ Inf tick)(ρ) ∨ f ∈ Inf term(ρ)
. This means that in an accepted run, every infinite clock ticks infinitely often, and every finite clock terminates properly (guaranteed by the infinitely often visited termination). Under these conditions we avoid runs where clocks are forbidden to tick while they should (we avoid local deadlock of clocks) and consequently the ticking trace of an accepted run would be a valid schedule.

For Example 1, its accepting criteria would be F = {a, T erm(b), T erm(c)}, since both b and c are supposed to be finite while a should be infinite. Runs where at least one tick of b occurs before the third tick of a would be recognized as accepted runs, while runs that lead to state (s, t 2 , r 0 ) would be ruled out.

Note that in practice, the number of clocks considered in defining accepting conditions could be minimized, as some of them can imply others. For instance, if c = c 1 sampledOn c 2 , one only needs to find schedules in which c ticks infinitely often, as that implies c 1 and c 2 do too. Please check the detail rules for minimizing accepting conditions in Table 1.

Inria

Table 1: Acceptance condition minimizing rules

Constraint Implication c = c1 filteredBy w, c = c1 delay_n (c1) = F in ⇒ (c) = F in, (c) = Inf ⇒ (c1) = Inf c = c1 await_n / until_n (c1) = F in ⇒ (c) = F in c = c1 upto c2 (c) = Inf ⇒ (c1) = Inf c = c1 concat/union/inf c2 (c1) = Inf ∨ (c2) = Inf ⇔ (c) = Inf , (c1) = (c2) = F in ⇔ (c) = F in c = c1 intersect/sampledOn c2 (c1) = F in ∨ (c2) = F in ⇒ (c) = F in, (c) = Inf ⇒ (c1) = (c2) = Inf c = c1 sup c2 (c1) = F in ∨ (c2) = F in ⇔ (c) = F in, (c) = Inf ⇔ (c1) = (c2) = Inf c subClock c1 (c) = Inf ⇒ (c1) = Inf c equality/alternate/boundedDiff_i_j c1 (c) = (c1) = Inf c1 strictPre/causes c2 (c2) = Inf ⇒ (c1) = Inf , (c1) = F in ⇒ (c2) = F in

Schedulability checked on automata translation

Once the tGBA built, schedulability analysis is equivalent to finding the existence of accepted infinite words. This (non)emptiness check has also been identified in model-checking context (where automata are generated from temporal logical formulaes). It is straightforward to check the emptiness on plain (state) Büchi automata (existence of a loop or SCC (Strongly Connected Component) including an accepting state, of course amongst the reachable states). But it is a bit more subtle to be checked on extended tGBAs directly, and we will consider this in the current section.

In addition to checking for emptiness, we will even attempt to compute the actual subset of states from which a valid infinite word may be continued. This will then allow to trim the structure and make the automata canonical so that it accepts the same language and no useless (undesired) state remains. This is helpful if one considers an approach where the resulting automata will be used later to guide an actual run-time simulation (similar as in controller synthesis aims). Hopefully, most related model-checking techniques actually do compute such sets of states and decide later language emptiness (or not) by checking whether the initial state was preserved as one of the useful states. The previous definition of transition-based Generalized Büchi Automata can be read as: a state is useful if it can be the root of a path which goes infinitely often through an accepting transition for every accepting state (in our case it is a path where each infinite specified clock ticks infinitely often, and proper termination of every finite specified clock is visited infinitely often), otherwise it is useless. We derive an algorithm from this. It is important to note here that a state may already be useless (in the sense that it leads to no valid infinite behavior) while it takes several steps in any way before bad states where clocks get halted.

For technical reasons, we require the finiteness of state space so that the fixpoints converge (and if not finite-state, the intuitive explanations have to be adapted because an infinite path does not mean a loop or SCC anymore). If a specification contains a constraint that may cause infinite state space (strictPre, causes, sup and inf), there must a be boundedDiff restricting finite drifts 1 . And we want the algorithm to be easily applicable on symbolic state set representations, such as BDDs. The operations used in the algorithm (union, intersection, complementation, Enable, P re), can easily be manipulated on symbolic representation of state sets. For a given specified (infinite or finite) clock c, Enable V (c) computes the set of states that immediately perform a transition where c ticks (if (c) = Inf ) or terminates (if (c) = F in). P re V (Z) computes the set of predecessors of states in Z. Subscript V is the scope where computations are carried.

EnableV (c) =    {s ∈ V |∃s ′ ∈ V, s C,T m ----→ s ′ , c ∈ C}, (c) = Inf , {s ∈ V |∃s ′ ∈ V, s C,T m ----→ s ′ , T erm(c) ∈ T m}, (c) = F in P reV (Z) = {s ′ ∈ V |∃s ∈ Z, s ′ C,T m
----→ s}. Our algorithm (Algorithm 1) consists of three steps: The first step computes the set of potentially useful states, P U . A state is potentially useful if in its future, every infinite clock can tick and every finite clock can properly terminate (may several times, infinitely, or immediately). For each specified clock, we compute independently the set of states that may reach a transition where it ticks (infinite specified) or terminates (finite specified) in a finite but unbounded number of steps, which is a straightforward backward (smallest) fixpoint of functional f (X) = EnableS(c) ∪ P reS(X). P U is then the intersection of all those result sets. P U is closed by co-reachability (if a state is in it, all states that may reach it also are). So if P U does not contain an initial state then we are done, the ccsl specification is not schedulable. Obviously a state that is not potentially useful is useless (from it there are clocks that get halted). But a P U state may not turn out to be actually useful: consider the case of a loop (or more generally a SCC) inside which two distinct clocks never tick, but from which there may be two actual transitions (leaving the SCC at different points), each of them allowing one of the two clocks to tick (and never the other Inria anymore). Then seemingly all states inside the SCC were labeled by both (here, all) clocks, while in fact the tickings cannot occur in the same path. (Case for finite clocks on proper terminations is similar.) Our solution to this is to apply the same algorithm as before, but this time restricting to P U as the set of reachable states (as thus regarding accepting transitions only if both their target and source states are in it). The set OU of once-useful states contains states that remain potentially useful. OU is also co-reachability closed.

The last step is to check the emptiness and get the set of useful states, U . We check if there is at least one loop (or SCC) inside U , by computing a greatest fixpoint on function f (X) = OU ∩ P reS(X). It computes the set of states from which there is an unbounded path inside U , containing SCCs plus extra states that may lead to them (the emptiness problem is therefore preserved between SCCs and this slight extension). U is the set of states preserved by this step. Since the state space is finite, the algorithm eventually ends. Proposition 5.1 A state is in U iff there is an infinite path from it in which each infinite clock ticks infinitely often and the termination of each finite clock is visited infinitely often (or in general terms where an accepting transition is performed infinitely often for each acceptance condition of the tGBA).

Proof (sketch): The ⇐ direction is easy: if there is such an infinite path, then the states found in it are all in OU by definition, and fold itself into a SCC in the end, so that all states remain in U .

For the ⇒ part: if s is in U , then it is in OU . From s, for a given infinite (resp. finite) specified clock c, (or equivalently, an accepting condition in generalized Buchi mode), there has to be a transition labeled as "c ticks" (resp. "c terminates"), with source and destination states (t, u) in OU . Now these states can also be chosen in U , otherwise the last state s ′ in the path from s would not remain labeled by at least one of the clock during the second step, a contradiction. Now we can reiterate this process from state u and any clock since U ultimately ends in (one or several) SCC, and progressively build an accepting infinite path. This concludes the proof.

Conclusion and future work

In this paper we use a translation from logical time constraint language to automata expansion to study schedulability. This work bears some similarity with previous attempts by Alur and Weiss [START_REF] Alur | Regular specifications of resource requirements for embedded control software[END_REF][START_REF] Alur | Rtcomposer:a framework for real-time components with scheduling interfaces[END_REF], which define schedules as infinite words expressed in regular expressions and then construct corresponding Büchi automatas, and have counterparts in more classical real-time and timed automata in [START_REF] Abdeddaim | Real-time scheduling of energy harvesting embedded systems with timed automata[END_REF][START_REF] Behrmann | Optimal scheduling using priced timed automata[END_REF][START_REF] Fersman | Schedulability analysis of fixedpriority systems using timed automata[END_REF][START_REF] Fersman | Task automata: Schedulability, decidability and undecidability[END_REF]. Comparing to them, we use ccsl as the bridge to connect scheduling and automata theories. This allows to take full of advantages of ccsl in modeling, such as the flexibility provided by logical time, multi-form time bases, unified expressing of both functional requirements and extra performance constraints. And we can construct automata naturally from ccsl specifications, avoiding the high costing of translation from regular expressions to automata.

We propose an algorithm for emptiness-checking and canonical trimming of translated automata, tGBA. As far as we know, it is a new way of handling generalized acceptance conditions, comparing to the existing two classes of emptiness-check algorithms: nested depth-first searches (NDFSs) based [START_REF] Courcoubetis | Memory efficient algorithms for the verification of temporal properties[END_REF][START_REF] Tauriainen | Nested emptiness search for generalized buchi automata[END_REF][START_REF] Tauriainen | On translating linear temporal logic into alternating and nondeterministic automata[END_REF], and strongly connected components (SCCs) computation based [START_REF] Geldenhuys | More efficient on-the-fly ltl verification with tarjan's algorithm[END_REF][START_REF] Moritz | Truly on-the-fly ltl model checking[END_REF][START_REF] Couvreur | On-the-fly emptiness checks for generalized büchi automata[END_REF]. And we believe it has relevance larger than our specific problem, being extended for general model-checking for transition or state GBAs.

In the future we want to inspect how finer relations between logical clocks can be used to detect unschedulability at construction time, in states even prior the full model is built. And we also want to make use of "latency-insensitive" property embedded in the interactions among constraints, which is distinct valid schedules depart from one another by shifting the exact timing placement they each assign to otherwise independent clock ticks.

Inria

Figure 1 :

 1 Figure 1: Execution of specification in Example 1

Inria

  (c) = (c 1 ) = Inf

  (c) = F in, (c 1 ) = Inf

Figure 2 :

 2 Figure 2: c = c 1 filteredBy u(v) ω

  0, and γ = 1 if the next bit of w is 0. If c is specified to be finite ( (c) = F in, (c 1 ) = Inf ), then v must be 0. Corresponding transitions and states are added, c 1 alone ticks infinitely after the termination of c, as Figure 2(b) shows. If c 1 is also specified to be finite ( (c) = (c 1 ) = F in), transitions labeled with T erm(c 1 ) are added. As Figure 2(c) shows, if c 1 terminates, c terminates directly. While c 1 can still tick until the bit of w is 1 after the termination of c.

Figure 3 :

 3 Figure 3: ccAutomata encodings of c = c 1 upto c 2

Constraint c = c 1 union c 2

 2 creates a new clock c which ticks whenever c 1 or c 2 ticks. c is finite if and only if both c 1 and c 2 is finite. If only c 1 (resp. c 2 ) terminates, c can still tick together with c 2 (resp. c 1 ) until the termination of c 2 (resp. c 1 ). On the opposite, c = c 1 intersection c 2 creates c which ticks iff both of c 1 and c 2 tick. If one of c 1 or c 2 terminates, c terminates, and the other clock can tick till its termination; if c terminates, c 1 and c 2 can tick infinitely in the future but not together.• sampledOnThere are two versions of the sampling, either strict (c = c 1 strictSamp c 2 ), or non-strict (c = c 1 nstrictSamp c 2 ). In both versions, clock c 1 is considered as a trigger and c 2 is a time base, they do not affect each other, and clock c ticks in coincidence with the tick of c 2 immediately following a tick of the trigger c 1 . Their difference can be clearly described by the following expressions:

Figure 4 : 2 123425Figure 5 :

 425 Figure 4: ccAutomata encodings of c = c 1 concat c 2

Figure 6 :

 6 Figure 6: ccAutomata encodings of c = c 1 intersect c 2

Figure 8 :

 8 Figure 8: ccAutomata encoding of sampledOn( (c 1 ) = F in, (c 2 ) = Inf )

Figure 4 . 1 .

 41 1 for (c 1 ) = (c 2 ) = Inf case; Figure 4.1.1 for (c 1 ) = F in, (c 2 ) = Inf case; Figure 4.1.1 for (c 1 ) = Inf, (c 2 ) = F in case; and Figure 4.1.1 for (c 1 ) = (c 2 ) = F in case. • inf, sup Constraint c = c 1 sup c 2 defines a new clock c which is slower than both c 1 and c 2 . The k th tick of c is coincident with the later of the i th tick of c 1 and c 2 . The index difference δ between c 1 and c 2 is recorded to determine if c should tick together with which clock, δ = i at s i . When δ = 0, if c 1 (resp. c 2 ) terminates, then c terminates; if c terminates, then either only c 1 ticks or only c 2 ticks in the future. When δ = 0, the slower one c 2 (resp. c 1 ) terminates iff c terminates; if c 1 (resp. c 2 ) terminates, then c can be alive till δ becomes 0, c 2 (resp. c 1 ) can still tick after the terminations of c and c 1 (resp. c 2 ). Constraint c = c 1 inf c 2 is the dual of c = c 1 sup c 2 . The k th tick of c is coincident with the faster of the i th tick of c 1 and c 2 . When δ = 0, if c terminates, then both c 1 and c 2 terminates; if c 1 (resp.c 2 ) terminates, then c 2 (resp. c 1 ) and c ticks together till their terminations. When δ > 0 (resp. δ = 0): if the slower one terminates then c and the faster one tick together in the future up to their terminations; if c terminates, then the faster one terminates immediately, the slower one can tick till δ becomes 0; if the faster one terminates, the slower one ticks alone till δ becomes 0, then ticks together with c before terminate.

Inria

  c = c 1 strictSamp c 2

Figure 9 :

 9 Figure 9: ccAutomata encoding of sampledOn( (c 1 ) = Inf, (c 2 ) = F in)

  c = c 1 strictSamp c 2

Figure 10 :

 10 Figure 10: ccAutomata encoding of sampledOn( (c 1 ) = (c 2 ) = F in)

Figure 14 :

 14 Figure 14: ccAutomata encoding of BoundedDiff

alternate c 1 1↓ 2 7 4 Figure 15 : c 1 alternate c 2 0 ≤ c 1 -c 2 ≤ 1 .

 14152121 Figure 15: c 1 alternate c 2

Figure 16 :

 16 Figure 16: Expand state s 1 in Figure 13(a) into s 11 and s 12

Definition 1 .

 1 A transition-based Generalized Büchi ccAutomata (tGBA) is a tuple (A, , F ) where A is a ccAutomata, is the finiteness specifying function, F is a finite set of elements called acceptance conditions, F = {c|c ∈ Clocks, (c) = Inf } ∪ {T erm(c)|c ∈ Clocks, (c) = F in}. A run ρ is accepted if transitions are labeled by each acceptance condition infinitely often. Let Inf tick(ρ) = {c|s C,T m ----→ s ′ ∈ Inf (ρ), c ∈ C} and Inf term(ρ) = {T erm(c)|s

algorithm Trimming 1 3 if

 13 Input: tgba = ((S, Clocks, T erms, I, T ), , F ) Output: set of states U begin 2 set of states PU,OU,U,Ufront,Uprev; S=∅ then return;

4 PU

 4 :=markState(tgba,S);

5 OUAlgorithm 1 :

 51 tgba = ((S, Clocks, T erms, I, T ), , F ), scope state set SS ⊆ S Output: set of states intersect begin 16 set of states cLabeled,cLabeledfront; 17 intersect := S; 18 for every specified clock c ∈ Clocks do 19 cLabeled :=EnableSS(c); cLabeledfront :=EnableSS(c); 20 while cLabeledfront = ∅ do 21 cLabeledfront:=P reSS(cLabeledfront)\cLabeled; 22 cLabeled:=cLabeled ∪ cLabeledfront; The Trimming algorithm

RR n°8102

Actually this restriction is overstrict, because of the interactions among constraints, the state space of a specification may be finite even if it has unbounded constraints. However, since the semantics is built constructively, it is not easy to check the finiteness from the syntactic constructors directly. So right now, we restrict to bounded ccsl[START_REF] Gascon | Logical time and temporal logics: comparing UML MARTE/CCSL and PSL[END_REF] instead of finiteness checking RR n°8102

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Laurent Albera

Anne Ferréol

Pierre Comon

ICAR, a tool for Blind Source Separation using Fourth Order Statistics only

Keywords: Independent Component Analysis, Blind Source Separation, Overdetermined Mixtures, Fourth Order Statistics

The problem of blind separation of overdetermined mixtures of sources, that is, with fewer sources than (or as many sources as) sensors, is addressed in this paper. A new method, named ICAR (Independent Component Analysis using Redundancies in the quadricovariance), is proposed in order to process complex data. This method, without any whitening operation, only exploits some redundancies of a particular quadricovariance matrix of the data. Computer simulations demonstrate that ICAR offers in general good results and even outperforms classical methods in several situations: ICAR (i) succeeds in separating sources with low signal to noise ratios, (ii) does not require sources with different SO or/and FO spectral densities, (iii) is asymptotically not affected by the presence of a Gaussian noise with unknown spatial correlation, (iv) is not sensitive to an over estimation of the number of sources.

I. INTRODUCTION

I NDEPENDENT Component Analysis (ICA) plays an im- portant role in various application areas, including radiocommunications, radar, sonar, seismology, radio astronomy, data analysis, speech and medical diagnosis [START_REF] Asano | Combined approach of array processing and independent component analysis for blind separation of acoustic signals[END_REF] [START_REF] De Lathauwer | Fetal electrocardiogram extraction by blind source subspace separation[END_REF]. In digital radiocommunications contexts for instance, if some sources are received by an array of sensors, and if the channel delay spread associated with the different sensors is significantly smaller than the symbol durations for each source, a static mixture of complex sources is observed on the sensors. On the other hand in Electrocardiography (ECG), it is possible to record the electrical activity of a fetal heart from ECGrecordings measured on the mother's skin. These cutaneous recordings can also be considered, in a first approximation, as instantaneous linear mixtures of potential signals generated by underlying bioelectric phenomena [START_REF] De Lathauwer | Fetal electrocardiogram extraction by blind source subspace separation[END_REF], hence again the static model considered.

The goal of Blind Source Separation (BSS) is to restore transmitted sources from the sole observation of sensor data. In some applications however, sources are not sought, and it is sufficient to identify the (static) mixture. Direction Of Arrival (DOA) estimation problems belong to this class [START_REF] Viberg | Sensor array processing based on subspace fitting[END_REF], since the column vectors of the mixture contain all the information necessary to determine the location of transmitters. The column vectors of the mixture are the so-called source steering vectors. It is thus legitimate to distinguish between blind identification of source mixtures and blind extraction of sources; we shall go back to this in section II.

Some algorithms utilize Second Order statistics (SO), as classically Principal Component Analysis (PCA) in Factor Analysis. In contrast, ICA attempts to restore the independence of outputs using higher order statistics. The consequence is that the indeterminacy is reduced, so that ICA allows to blindly identify the static mixture, and transmitted sources can eventually be extracted. More precisely, the ICA concept relies on the core assumption that (i) sources should be independent in some way. Additionally, when a contrast functional is sought to be maximized, (ii) the mixture has to be overdetermined, which means that there should be at most as many sources as sensors [START_REF] Zhang | Natural gradient algorithm to blind separation of overdetermined mixture with additive noises[END_REF]. In fact, there must exist a linear source separator [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF] in the latter framework.

On the other hand, the more general case where there may be more sources than sensors is often referred to as Blind Identification of underdetermined mixtures, and is not considered in this paper but is addressed elsewhere; see [START_REF] Cardoso | Super-symetric decomposition of the fourth-order cumulant tensor. Blind identification of more sources than sensors[END_REF] [19] [START_REF] Diamantaras | Blind separation of multiple binary sources using a single linear mixture[END_REF] [START_REF] Taleb | An algorithm for the blind identification of N independent signal with 2 sensors[END_REF] [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF] [3] [START_REF]Blind identification and source separation in 2x3 underdetermined mixtures[END_REF] and references therein.

Since the first paper related to Higher Order (HO) BSS, published in 1985 [START_REF] Herault | Détection de grandeurs primitives dans un message composite par une architecture de calcul neuromimétique en apprentissage non supervisé[END_REF], many concepts and algorithms have come out. For instance, the ICA concept was proposed a few years later, as well as the maximization of a Fourth Order (FO) contrast criterion (subsequently referred to as COM2) [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF]. At the same time, a matrix approach was developed in [START_REF] Cardoso | Blind beamforming for nongaussian signals[END_REF] and gave rise to the joint diagonalization algorithm (JADE). A few years later, Hyvarinen et alterae developed the FastICA method, first for signals with values in the real field [START_REF] Hyvarinen | A fast fixed-point algorithm for independent component analysis[END_REF], and later for complex signals [START_REF] Bingham | A fast fixed-point algorithm for independent component analysis of complex valued signals[END_REF], using the fixed-point algorithm to maximize a FO contrast. This algorithm is of deflation type, as that of Delfosse et alterae [START_REF] Delfosse | Adaptive blind separation of independent sources: a deflation approach[END_REF], and must extract one source at a time, although some versions of FastICA extract all sources simultaneously. In addition, Comon proposed a simple solution [START_REF]From source separation to blind equalization, contrast-based approaches[END_REF], named COM1 in this paper, to the maximization of another FO contrast function previously published in [START_REF] Moreau | New self-adaptive algorithms for source separation based on contrast functions[END_REF] [18] [START_REF] Thirion | New criteria for blind signal separation[END_REF]. Another algorithm of interest is SOBI, based only on SO statistics, developed independently by several authors in the nineties, and addressed in depth later in [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF].

Each of these methods suffers from limitations. To start with, the SOBI algorithm is unable to restore components that have similar spectral densities. Moreover, the JADE method is very sensitive to an over estimation of the number of sources as shown in the simulation section of this paper and in [START_REF] Albera | ICAR: Independent Component Analysis using Redundancies[END_REF]. Note that in electronic war fare contexts, the number of sources needs to be estimated and may be over estimated, especially for low signal to noise ratios [START_REF] Wax | Detection of signals by information theoretic criteria[END_REF] [START_REF] Zhao | Remarks on certain criteria for detection of number of signals[END_REF]. On the other hand, though the previous methods [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF] [START_REF]From source separation to blind equalization, contrast-based approaches[END_REF] [8] [START_REF] Hyvarinen | A fast fixed-point algorithm for independent component analysis[END_REF] [6] perform under some reasonable assumptions, they may be strongly affected by a Gaussian noise with unknown spatial correlation as shown in the simulation section of this paper. Such a noise appears for instance in some radiocommunications applications. It is in particular the case for ionospheric radiocommunications in the HF band where the external noise, composed of multiple emitting sources (industrial noise, atmospheric noise...), is much stronger than the thermal noise generated by the receivers. In order to deal with the correlated noise problem, Ferréol et alterae [START_REF] Ferreol | Higher order blind source separation using cyclostationarity property of the signals[END_REF] and Abed-Meraim et alterae [START_REF] Abed-Meraim | Blind source separation using second-order cyclostationary statistics[END_REF] have proposed a new family of BSS methods respectively exploiting the potential cyclostationarity of the received sources. In fact, the latter family of algorithms uses cyclic statistics of the data. Note that a cyclic covariance matrix associated with a stationary noise is null for non zero cyclic frequencies. Consequently, these cyclic methods allow the optimal separation of independent cyclostationary sources even in the presence of a stationary noise with unknown spatial correlation. However, the use of cyclic methods is more complex because of the estimation of cyclic frequencies and time delays. To overcome this drawback, Ferréol et alterae have recently introduced the FOBIUM algorithm [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF], which, without SO whitening step, performs the blind source separation even in the presence of a Gaussian noise with unknown spatial correlation. Nevertheless, since FOBIUM is an extension of the SOBI method to FO statistics, it requires sources with different FO spectral densities. FOBIUM also allows to address the underdetermined case, but this is out of the scope of the present paper.

In order to overcome the limitations of the previous algorithms, the method named ICAR (Independent Component Analysis using Redundancies in the quadricovariance) shortly presented in [START_REF] Albera | ICAR: Independent Component Analysis using Redundancies[END_REF] is proposed in this paper and addresses the case of complex mixture and sources, in the presence of additive (possibly spatially correlated) Gaussian noise. Only based on FO statistics, ICAR skips the SO whitening step in contrast to classical methods [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF] [15] [START_REF]From source separation to blind equalization, contrast-based approaches[END_REF] [8] [START_REF] Hyvarinen | A fast fixed-point algorithm for independent component analysis[END_REF] [6], and consequently is asymptotically not affected by the presence of a Gaussian noise with unknown spatial correlation. Actually, ICAR exploits redundancies in a particular FO statistical matrix of the data, called quadricovariance. The latter algorithm assumes sources to have non zero FO marginal cumulants with the same sign, assumption that is verified in most radiocommunications contexts. Indeed, the kurtosis of most of radiocommunications signals is negative. For example M -PSK constellations have a kurtosis equal to -2 for M = 2 and to -1 for M ≥ 4. Continuous Phase Modulations (CPM), among which we find the GMSK modulation (GSM standard), are such that their kurtosis is smaller than or equal to -1, due to their constant modulus. Furthermore, the performance of ICAR is also analyzed in this paper, in different practical situations through computer simulations, and compared to those of classical algorithms, namely SOBI, COM1, COM2, JADE, FastICA and FOBIUM. It appears that ICAR exhibits good results in most cases even when classical methods fail.

The paper is organized as follows. Section II introduces the BSS problem, and assumptions needed in ICAR. Section III defines the SO and FO statistics considered in the paper, and section IV describes in detail the ICAR concept. Computer results are reported in section V. Section VI eventually concludes.

II. ASSUMPTIONS AND PROBLEM FORMULATION

A noisy mixture of P statistically independent Narrow-Band (NB) sources s p (k) is assumed to be received by an array of N sensors. In accordance with the usual practice [START_REF] Proakis | Digital Communications[END_REF], only complex envelopes of NB signals are considered. The vector of complex envelopes of the signals at the sensor outputs, x(k), is thus given by

x(k) = A s(k) + ν(k) (1)
where A, s(k), ν(k) are the N×P constant mixing matrix, the P×1 source with components s p (k) and N×1 noise random vectors, respectively. In addition, for any fixed index k, s(k) and ν(k) are statistically independent. We further assume the following hypotheses: A1) Vector s(k) is stationary, ergodic 1 with components a priori in the complex field and mutually uncorrelated at order 4; A2) Noise vector ν(k) is stationary, ergodic and Gaussian with components a priori in the complex field too; A3) FO marginal source cumulants, called kurtosis (if normalized) and defined in section III-B, are not null and have all the same sign; A4) The mixture matrix A does not contain any null entry; A5) A is a full column rank matrix. Note that sources with null kurtosis are tolerated but cannot be seen and processed by method ICAR. Such sources will be considered as noise. Moreover, the second part of (A3) will be discussed in section IV-C.1. Assumption (A4) is not a strong assumption, in particular in digital radiocommunications contexts, since it is more than just reasonable to assume the array of sensors in good repair. On the other hand, if the n-th sensor is defective, the n-th row of A will be null. It is then necessary to erase the contribution of this sensor and to assume that we have N -1 sensor outputs instead of N . As far as the masking phenomenon is concerned, it is more rare and may produce at most one null component in each column of A for arrays with space diversity. Forthcoming works will consist in studying the ICAR robustness with respect to this pathological phenomenon. As far as (A5) is concerned, it implies necessarily P ≤ N . Under the previous assumptions, the problems addressed in the paper are both the blind identification and the blind extraction of the sources using solely the FO statistics of the data. The goal of Blind Mixture Identification (BMI) is to blindly identify the mixing matrix A, to within a trivial matrix T ; recall that a trivial matrix is of the form ΛΠ where Λ is invertible diagonal and Π a permutation. On the other hand, the goal of Blind Source Extraction (BSE) or Separation (BSS) is to blindly find a matrix W , yielding a P×1 output vector y(k) = W H x(k) corresponding to the best estimate, ŝ(k), of the vector s(k), up to a multiplicative trivial matrix. Superscript ( H ) denotes the complex conjugate transpose of a matrix.

III. SO AND FO DATA STATISTICS

A. SO statistics

The SO statistics considered in the paper are given by

C i2 i1, x (k) = Cum{x i1 (k), x i2 (k) * } (2)
Function ( 2) is well-known as the SO cumulant of x(k).

Consequently, the SO marginal cumulant of source s p (k) is defined by

C p p, s (k) = Cum{s p (k), s p (k) * } (3)

B. FO statistics

The FO statistics considered in the paper are given by

C i3, i4 i1, i2, x (k) = Cum{x i1 (k), x i2 (k), x i3 (k) * , x i4 (k) * } (4)
where two terms x i (k) are not conjugate and two terms are conjugate. Function ( 4) is well-known as the FO cumulant of x(k). Consequently, the FO marginal cumulant of source s p (k) is defined by

C p, p p, p, s (k) = Cum{s p (k), s p (k), s p (k) * , s p (k) * } (5)
Likewise, the kurtosis of source s p (k) is given by

κ p,p p,p,s (k) = C p, p p, p, s (k)/(γ p (0)) 2 (6)
where γ p (0) is the variance of source s p (k). Note that in the presence of stationary sources, SO (2) and FO (4) statistics do not depend on time k, so that they can be denoted by C i2 i1, x and C i3, i4 i1, i2, x , respectively.

C. Matrix arrangement

1) SO and FO statistical matrices: SO and FO statistics computed according to (2) and ( 4) may be arranged in two Hermitian statistical matrices, R x = C 2, x and Q x = C 4, x , of size N×N and N 2 ×N 2 , respectively. These matrices are called the covariance and the quadricovariance of x(k), respectively. We limit ourselves to arrangements of SO and FO statistics that give different results, in terms of maximum number of processed sources at the output of the BSS methods. The impact of the chosen way to arrange statistics in a matrix is analyzed in [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF]. It is shown in [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF], through extensions of the Virtual Array concept initially introduced in [START_REF] Dogan | Applications of cumulants to array processing -Part I: Aperture extension and array calibration[END_REF] and [START_REF] Chevalier | On the virtual array concept for the fourth-order direction finding problem[END_REF] for the FO data statistics, that there exists an optimal arrangement of the FO cumulants in a quadricovariance matrix with respect to the maximal number of statistically independent sources to be processed by a method exploiting the algebraic structure of this quadricovariance. As far as SO statistics are concerned, there is a unique non redundant way to store them in a matrix R x under constraints of hermicity. Consider indeed the following arrangement

R x (i 1 , i 2 ) = C i2 i1, x (7)
where

R x (i 1 , i 2 ) is the (i 1 , i 2 )-th entry of matrix R x ; the other possible arrangement R x ′ (i 2 , i 1 ) = C i2 i1, x just leads to R H
x and hence to the same result in terms of maximum number of processed sources. On the other hand, there are two distinct non redundant ways associated with FO statistics under constraints of hermicity, which can be indexed by the integer ℓ (ℓ ∈ {0, 1}). Each way yields a statistical matrix C ℓ 4, x such that its

I ℓ 1 ,I ℓ 2 -th entry (1 ≤ I ℓ 1 ,I ℓ 2 ≤ N 2 ) is given by C ℓ 4, x I ℓ 1 ,I ℓ 2 = C i3, i4 i1, i2, x (8)
where for any ℓ belonging to {0, 1} and for all i 1 , i

2 , i 3 , i 4 (1 ≤ i 1 , i 2 , i 3 , i 4 ≤ N ), I ℓ 1 = i 2 + N (i 1 -1) if ℓ = 0 i 4 + N (i 1 -1) if ℓ = 1 (9)
and

I ℓ 2 = i 4 + N (i 3 -1) if ℓ = 0 i 2 + N (i 3 -1) if ℓ = 1 (10)
Note that the optimal arrangement is shown in [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF] to correspond to ℓ=1 and for this reason, we consider this arrangement in the following sections. So matrices C 1 4, s , C 1 4, x will be denoted by Q s and Q x respectively. Remark 1 Another way, perhaps more intuitive (especially for readers familiar with Matlab), to present the construction of Q x is the following: first, construct a 4-dimensional tensor T , whose elements are given by

T(i 4 , i 1 , i 2 , i 3 ) = C i3, i4 i1, i2, x
The matrix Q x is then given by a simple Matlab reshape operation as follows

Q x = reshape(T , N 2 , N 2 )
2) Multilinearity property: The SO and FO statistical matrices of the data, R x and Q x , have a special structure, due to the multilinearity property under change of coordinate systems, which is enjoyed by all moments and cumulants. Since sources and noise are independent, this property can be expressed, for SO statistical matrices and according to [START_REF] Cardoso | Super-symetric decomposition of the fourth-order cumulant tensor. Blind identification of more sources than sensors[END_REF], by

R x = A R s A H + R ν (11)
Similarly, according to (8), ( 9) and [START_REF] Chevalier | Performances des séparateurs aveugles de sources aux ordres supérieurs[END_REF], and since noise is Gaussian and independent of sources, the FO cumulant matrix can be expressed as follows, using the multilinearity property associated with ℓ = 1:

Q x = [A⊗A * ] Q s [A⊗A * ] H (12)
The P ×P matrix R s and the P 2 ×P 2 matrix Q s are the SO and FO statistical matrices of s(k) respectively. R ν denotes the N ×N SO statistical matrix of ν(k).

D. Statistical estimation

In practical situations, SO and FO statistics have to be estimated from components of x(k). If components are stationary and ergodic, sample statistics may be used to estimate (2) and (4). Nevertheless, if sources are cyclostationary, cycloergodic, potentially non zero-mean, SO and FO continuous time average statistics have to be used instead of ( 2) and ( 4), such as

C i2 i1, x = C i2 i1, x (k) c (13)
and

C i3, i4 i1, i2, x = C i3, i4 i1, i2, x (k) c ( 14
)
where • c is the continuous time average operation defined by

∀ f, f: t -→ f (t), f (t) c = lim T →+∞ 1 T T /2 -T /2 f (t)dt (15)
These continuous-time temporal mean statistics need some knowledge on cyclic frequencies of the received signal and are thus computed using, for instance, the unbiased and consistent estimators described in [START_REF]On the behavior of current second and higher order blind source separation methods for cyclostationary sources[END_REF], [START_REF]Second order blind separation of first and second order cyclostationary sources -Application to AM, FSK, CPFSK and deterministic sources[END_REF] and [START_REF] Ferreol | Higher order blind separation of non zero-mean cyclostationary sources[END_REF]. Moreover, ordering these continuous-time temporal mean statistics in matrices R x and Q x by means of ( 7), ( 8), ( 9) and ( 10) respectively, expressions [START_REF]Optimal separation of independent narrow-band sources : Concept and Performances[END_REF] and ( 12) remain valid.

IV. THE ICAR METHOD We present in this section a new method of BSS named ICAR, which exploits the algebraic structure of an alternative expression of matrix Q x .

A. Matrix notation

Define a columnwise Kronecker product, denoted ⊘ and referred to as the Khatri-Rao product [START_REF] Eriksson | Identifiability, separability and uniquness of linear ICA models[END_REF] [START_REF] Yeredor | Non-orthogonal joint diagonalization in the leastsquares sense with application in blind source separation[END_REF]. For any N×P rectangular matrices G and H, the columns of the N2 ×P matrix G ⊘ H are defined as g j ⊗ h j , where ⊗ denotes the usual Kronecker product, if g j and h j denote the columns of G and H respectively.

B. The core equation

The ICAR method exploits several redundancies present in the quadricovariance matrix of the data, Q x . Although most BSS algorithms, such as JADE, exploit expression [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF], the ICAR method uses an alternative form, described by

Q x = [A⊘A * ] ζ s [A⊘A * ] H (16)
where the P×P diagonal matrix [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF], and where the N 2 ×P matrix A⊘A * is defined by

ζ s = Diag C 1,1 1,1,s , C 2,2 2,2,s , • • •, C P,P P,P,s (i.e. ∀ (p 1 , p 2 ), 1 ≤ p 1 , p 2 ≤ P , ζ s (p 1 , p 2 ) = C p1,p1 p1,p1,s if p 1 =p 2 , 0 otherwise) is full rank, in contrast to Q s in
A⊘A * = [ a 1 ⊗a * 1 a 2 ⊗a * 2 • • • a P ⊗a * P ] (17)
and can be written as

A⊘A * = [ [A * Φ 1 ] T [A * Φ 2 ] T • • • [A * Φ N ] T ] T (18)
with

Φ n = Diag[ A(n, 1), A(n, 2), • • •, A(n, P ) ] (19)
In other words, the non zero elements of the P×P diagonal matrix Φ n are the components of the n-th row of matrix A.

In addition, note that equation ( 16) can be easily derived from equation [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF]. Indeed, the latter equation straight implies, in view of the structure of the diagonal non inversible matrix Q s , that the only column vectors of matrix A⊗A * which generate matrix Q x are a p ⊗a * p (1 ≤ p ≤ P ), hence result [START_REF]From source separation to blind equalization, contrast-based approaches[END_REF].

C. The ICAR concept

The algorithm proposed proceeds in three stages. Firstly, a unitary matrix V is estimated in the Least Square (LS) sense, and allows the estimation of A ⊘ A * from Q x [START_REF]From source separation to blind equalization, contrast-based approaches[END_REF]. In a second stage, several algorithms may be thought of in order to compute an estimate of A from A⊘A * . Finally, an estimation of sources s(k) is computed using the estimate of A.

1) Identification of A⊘A * : Matrix A⊘A * is an unobservable square root of Q x to within a diagonal matrix, as shown by [START_REF]From source separation to blind equalization, contrast-based approaches[END_REF]. In this context, the idea is to built an observable square root,

Q x 1/2 , of Q x , differing from [A⊘A * ] ζ s 1/2
only by a unitary matrix V and then to identify the latter from the exploitation of the algebraic structure of Q x . So consider the following proposition:

Proposition 1 If A is of full column rank (A5) and contains no null entries (A4), then the

N 2 × P matrix A ⊘ A * is full column rank.
The proof is given in appendix I. So proposition 1 and assumption (A3) allow together to prove that matrix Q x , given by ( 16), is of rank P . Moreover, assumption (A3) and equation [START_REF]From source separation to blind equalization, contrast-based approaches[END_REF] imply that Q x is positive if the FO marginal source cumulants are positive, which we assume in the following. Thus, a square root of Q x , denoted Q x 1/2 , and defined such that

Q x = Q x 1/2 [Q x 1/2
] H , may be computed. If the FO marginal source cumulants are negative, matrix -Q x can be considered instead for computing the square root. In the case where there are terms with a different sign, our derivation can be reformulated in terms of an unknown J -unitary matrix 2 V , instead of unitary. Then we deduce from ( 16) that matrix

[A ⊘ A * ] ζ s 1/2 is a natural square root of Q x .
Yet, another possibility is to compute this square root via the Eigen Value Decomposition (EVD) of Q x given by

Q x = E s L s E s H (20)
where L s is the real-valued diagonal matrix of the non zero eigenvalues of Q x . Since matrix Q x is of rank P , L s is of size P×P . Besides, E s is the N 2 ×P matrix of the associated orthonormalized eigenvectors. Consequently, a square root of Q x can be computed as

Q x 1/2 = E s L 1/2 s ( 21
)
where

L 1/2
s denotes a square root of L s . Proposition 2 For a full rank matrix A⊘A * , source kurtoses are not null and have all the same sign (A3) if and only if the diagonal elements of L s are not null and have also the same sign, corresponding to that of the FO marginal source cumulants.

The proof is given in appendix II. In addition, equation [START_REF] Proakis | Digital Communications[END_REF] can be rewritten as

Q x 1/2 = E s L 1/2 s = [A⊘A * ] ζ s 1/2 V H , (22)
showing the link between Q x 1/2 and A⊘A * . Plugging ( 18) into ( 22), matrix Q x 1/2 can be eventually rewritten as

Q x 1/2 = [A * Φ 1 ζ s 1/2 V H ] T ...[A * Φ N ζ s 1/2 V H ] T T =[Γ 1 T ... Γ N T ] T (23)
where the N matrix blocks Γ n of size N×P are given by

∀ n, 1 ≤ n ≤ N, Γ n = A * Φ n ζ s 1/2 V H (24)
Proposition 3 For any n (1 ≤ n ≤ N ), matrix Γ n is full column rank.

The proof is given in appendix III. Using proposition 3, the pseudo-inverse Γ ♯ n of the N×P matrix Γ n is defined by

∀ n, 1 ≤ n ≤ N, Γ ♯ n = (Γ n H Γ n ) -1 Γ n H (25)
Then, consider the N(N -1) matrices Θ n1,n2 below

∀ (n 1 , n 2 ) , 1 ≤ n 1 = n 2 ≤ N, Θ n1,n2 = Γ ♯ n1 Γ n2 (26)
which can be rewritten, from ( 24) and ( 25), as

Θ n1,n2 = V ζ s -1/2 Φ -1 n1 Φ n2 ζ s 1/2 V H (27) = V Φ -1 n1 ζ s -1/2 ζ s 1/2 Φ n2 V H = V Φ -1 n1 Φ n2 V H
where ζ s 1/2 , Φ n1 , Φ n2 and D n1,n2 = Φ -1 n1 Φ n2 , are P×P diagonal full rank matrices (the full rank character of matrices Φ n1 , Φ n2 and D n1,n2 is due to assumption (A4)). It appears from ( 27) that matrix V jointly diagonalizes the N(N -1) matrices Θ n1,n2 .

Proposition 4 If A is of full column rank (A5) and contains no null entries (A4), then, for all pairs

(p 1 , p 2 ), 1 ≤ p 1 = p 2 ≤ P , at least one pair (n 1 , n 2 ), 1 ≤ n 1 = n 2 ≤ N exists such that D n1,n2 (p 1 , p 1 ) = D n1,n2 (p 2 , p 2 ).
The proof is given in appendix IV. Under proposition 4, paper [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF] allows to assert that if V sol jointly diagonalizes matrices Θ n1,n2 , then V sol and V are related through V sol = V T where T is a trivial unitary matrix. So matrix V sol allows, in accordance with [START_REF] Diamantaras | Blind separation of multiple binary sources using a single linear mixture[END_REF], to recover A ⊘ A * up to a multiplicative trivial matrix:

Q x 1/2 V sol = [A⊘A * ] ζ s 1/2 T (28)
2) Identification of mixture A: Three algorithms are proposed in this section to identify A from the estimate, Q x 1/2 V sol , of A⊘A * . These algorithms optimize differently the compromise between performance and complexity.

Note that equation ( 28) can be rewritten from [START_REF] Comon | Improved contrast dedicated to blind separation in communications[END_REF] in the form of N matrix blocks

Σ n = A * Φ n ζ s 1/2 T of size N×P as Q x 1/2 V sol = [Σ 1 T Σ 2 T • • • Σ N T ] T ( 29
)
So a first approach to estimate A up to a trivial matrix, called ICAR1 in the sequel, consists of merely keeping the matrix block Σ 1 * made up of the N first rows of Q x 1/2 V sol such that

Σ 1 = A * Φ 1 ζ s 1/2 T ( 30
)
where Φ 1 and ζ s 1/2 are diagonal matrices, and where T is a unitary trivial matrix.

It is also possible to take into account all the matrix blocks Σ n * and to compute their average. This yields a second algorithm, named ICAR2, of higher complexity.

A third algorithm, called ICAR3, is now described, and yields a more accurate solution to the BMI problem: since matrix A ⊘ A * , given by ( 17), has been identified from the previous section by Q x 1/2 V sol to within a trivial matrix, ICAR3 consists first of mapping each N 2 × 1 column vector b p of Q x 1/2 V sol into a N×N matrix B p (the n-th column of B p is made up from the N consecutive entries of vector b p , between [N(n -1) + 1] and N n), and secondly of diagonalizing each matrix B p * .

Proposition 5 For any matrix B p (1 ≤ p ≤ P ) built from Q x 1/2 V sol , there exists a unique column vector a q (1 ≤ q ≤ P ) of A such that the eigenvector of B p * associated with the largest eigenvalue corresponds, up to a scale factor, to a q .

The proof is given in appendix V. In addition, the indeterminacy of the norms of columns of A is related to matrices Φ p , ζ s 1/2 , a unitary diagonal matrix (whose product by a permutation matrix gives T ) and the way to identify A from matrix Q x 1/2 V sol . As far as the permutation indeterminacy is concerned, it is related to matrix T .

3) Extraction of the P independent components: Finally, to estimate the signal vector s(k) for any value k, it is sufficient, under (A5), to apply a linear filter built from the identified matrix A : such a filter may be the Spatial Matched Filter (SMF) given by [START_REF]Optimal separation of independent narrow-band sources : Concept and Performances[END_REF]

W = R -1
x A, which is optimal in the presence of decorrelated signals. In practical situations, since matrix A is estimated up to a trivial matrix according to section (IV-C.2), neither the order of sources s(k) nor their amplitude can be identified.

D. Implementation of the ICAR methods

The different steps of the ICAR method are summarized hereafter when K samples of the observations, x(k) (1 ≤ k ≤ K), are available.

Step1 Compute an estimate of FO statistics C i3, i4 i1, i2, x from the K samples x(k) and store them, using the (ℓ = 1)arrangement, into matrix Q x , which is an estimate of Q x .

Step2 Compute the EVD of the Hermitian matrix Q x , estimate P , the number of sources, from this EVD. Restrict Q x to the P principal components :

Q x = E s L s E s
H , where L s is the diagonal matrix of the P eigenvalues of largest modulus and E s is the matrix of the associated eigenvectors.

Step3 Estimate the sign, ǫ, of the diagonal elements of L s .

Step4 Compute a square root matrix

[ǫ Q x ] 1/2 of ǫ Q x : [ǫ Q x ] 1/2 = E s | L s | 1/2 , where |•| denotes the absolute value operator. Step5 Compute from [ǫ Q x ] 1/2 the N matrices Γ n , construct matrices Θ n1,n2 = [ Γ ♯ n1 Γ n2 ] for all (n 1 , n 2 ), 1 ≤ n 1 = n 2 ≤
N , and compute V sol , an estimate of V sol , from the joint diagonalization of the N(N-1) matrices Θ n1,n2 ; one possible joint diagonalization algorithm may be found in [START_REF]Jacobi angles for simultaneous diagonalization[END_REF].

Step6 Compute an estimate A of the mixture A from the

N 2 ×P matrix [[ǫ Q x ] 1/2 V sol ]
by either one of the following: 1) (ICAR1) taking the matrix block made up of the N first rows of

[[ǫ Q x ] 1/2 V sol ] * ;
2) (ICAR2) taking the average of the N matrix blocks, of size N × P , made up of the successive rows of

[[ǫ Q x ] 1/2 V sol ] * ; 3) (ICAR3) taking each column vector b p of [[ǫ Q x ] 1/2 V sol ]
remodeling them into N×N matrices B p , and building the matrix whose p-th column vector is the eigenvector of matrix B p * associated with the largest eigenvalue.

Step7 Estimate the signal vector s(k) for any value k by applying to x(k) a linear filter built from A , like for example the SMF defined by

W = R x -1
A .

V. COMPUTER RESULTS

In this section, a comparative performance analysis of seven BSS methods (SOBI, COM1, COM2, JADE, FastICA, FOBIUM and ICAR) in various scenarios is presented. For this purpose, we consider a Uniform Linear Array (ULA) of N = 4 sensors, except for figure 2 where N = 2, equispaced half a wavelength apart [START_REF] Monzingo | Introduction to adaptive arrays[END_REF]. P = 2 QPSK sources are linearly modulated with a pulse shape filter corresponding to a 1/2-Nyquist filters with a roll off equal to 0.3 [START_REF] Proakis | Digital Communications[END_REF]. In addition, the P = 2 sources have the same symbol period T and the same Signal to Noise ratio (SNR) equal to 15 dB, except for figures 4(a) and 4(b). The sources are assumed to be well angularly separated except for figure 6, where the other cases are also considered. The source carrier residuals are such that f c1 T e = 0, f c2 T e = 0.65, except for figure 2 where f c2 T e = 0. The sample period T e corresponds to the symbol period T . As a result, the used SO and FO statistics are time invariant, so that classical sample estimators may be employed. As far as the background noise is concerned, it is temporally and spatially white except for section V-.1. Eventually, the simulation results are averaged over 200 realizations. Note that we resample the sources and the noise between these 200 experiments. On the other hand, the mixing matrix does not change except for figure 6 where its influence on the BSS methods performance is pointed out.

Moreover, the criterion used in this paper, in order to evaluate performance of BSS algorithms, is the well-known SINRM (Signal to Interference plus Noise Ratio Maximum) criterion defined in [11, section 3]. In other words, for each source s p (k) (1 ≤ p ≤ P ), the Signal to Interference plus Noise Ratio for the source p at the output of a spatial filter w i is defined by

SIN R p [w i ] = γ p (0) |w i H a p | 2 w i H R νp w i (31)
where γ p (0) is the variance of the p-th source. Moreover, R νp is the total noise covariance matrix for source p, corresponding to matrix R x in the absence of source p. In these conditions, the restitution quality of source p at the output of separator W , whose columns are the w i , can be evaluated by the maximum value of SINR p [w i ] when i varies from 1 to P , and may be denoted SINRM p .

1) The white noise case: The performance of ICAR at the output of the considered source separator is firstly illustrated in the presence of a Gaussian noise, spatially and temporally white, and compared with some well-known BSS algorithms. Figures 1(a show the good performance of the ICAR algorithm, especially ICAR3 (the third method given in section IV-C.2), facing the well-known SOBI, COM1, COM2, JADE, FastICA and FOBIUM methods. As for the SOBI method, it requires about 450 snapshots to obtain good results, due to a mild difference between the spectral densities of the sources. Note that similar results have been obtained for the other source. In addition, since the best results between the three ICAR methods are obtained in particular for ICAR3, we report in the following sections the comparison results only for this third method. Contrary to the other figures, figure 2 shows performance results when the two QPSK are chosen in baseband, i.e. taking f c1 T e = f c2 T e = 0, which implies that the two source signals have identical trispectra. Consequently, the SOBI and FOBIUM algorithms are unable to separate them correctly. However, we note that the FOBIUM method seems to be more robust than SOBI with respect to a spectrum difference of the sources. Moreover, other simulations have shown that the FOBIUM results are better as quotient N P increases, even if they remain suboptimal.

Figure 3 shows, for a number of 400 samples, the variations of SINRM 2 at the output of the previous methods as a function of the input SNR, identical for the two sources. All the BSS methods have approximately the same behavior. First, when the SNR is very small, they do not succeed perfectly in extracting the third source. On the contrary, for signal to noise ratios beteen -4 and 20 dB, the source separation is optimal. Finally, although the variations of SINRM 2 for signal to noise ratios greater than 20 dB are somewhat surprising, this result has already been observed by Monzingo and Miller in [START_REF] Monzingo | Introduction to adaptive arrays[END_REF] for optimal separators when mixture A is known. Note that similar results have been obtained for the other source.

2) The colored noise case: Then, the ICAR3 method is compared with the other algorithms in the presence of a Gaussian noise with unknown spatial correlation. Figure 4(a) and 4(b) show the variations of SINRM 2 at the output of the previous methods as a function of the noise spatial correlation factor ρ. SNR of the two sources is taken equal to 5 dB (figure 4(a)), next 0 dB (figure 4(b)). In addition, 400 samples are used to extract the two sources. Note that the Gaussian noise model employed in this simulation is the sum of an internal noise ν in (k) and an external noise ν out (k), of covariance matrices R in ν and R out ν respectively such that

R in ν (r, q) def = σ 2 δ(r-q)/2 R out ν (r, q) def = σ 2 ρ |r-q| /2 (32)
where σ 2 , ρ are the total noise variance per sensor and the noise spatial correlation factor respectively. Note that R ν (r, q) def = R in ν (r, q) + R out ν (r, q) is the (r, q)-th component of the total noise covariance matrix. It appears in figure 4(a) that FOBIUM and ICAR3 are insensitive to a Gaussian noise with unknown spatial correlation, whereas ICAR3 seems to be a bit more robust than FOBIUM. On the other hand, the wellknown COM1, COM2, JADE and SOBI methods are strongly affected as soon as the noise spatial correlation increases beyond 0.5. In fact, the classical BSS methods require a prior spatial whitening based on second order moments. This stage theoretically needs the perfect knowledge of the noise covariance. If this is not the case, a whitening of the observed data is performed instead, which is biased. ICAR does not suffer from this drawback, since it uses only FO cumulants, which are (asymptotically) insensitive to Gaussian noise, regardless of its space/time color. Note the poor performance of FastICA due to the presence of weak sources. Besides, similar results have been observed for source 1. As far as figure 4(b) is concerned, it confirms the fact that the performance differences between ICAR3 and the classical BSS methods increases as the source SNR decreases.

3) Over estimation of the number of sources: On the other hand, in operational contexts, the number of sources may be over estimated. It is then interesting to compare the ICAR method with other algorithms in such situations. To this aim, we assume that the estimated number of sources is equal to P = 3. Figure 5 shows the variations of SINRM 2 (source 2 performance) at the output of the previous methods as a function of the number of samples while the input SNR of the two sources is assumed to be equal to 15 dB. Similar results have been observed for source 1. More particularly, it appears that the FastICA and ICAR3 methods are robust with respect to an over estimation of the number of sources whereas, in this simulation configuration, the JADE algorithm looses 15 dB, for less than 1000 samples, with respect to the case where P = 2. As for the other methods, such as the FOBIUM algorithm, they are also affected by this over estimation, but less than the JADE algorithm since they lose on average 3 dB. The explanation of this surprising phenomenon is not easy and is beyond the scope of this paper. However, a similar behavior had been observed in [START_REF] Chevalier | Performance of HO blind source separation methods: Experimental results on ionospheric HF links[END_REF] [10] when comparing JADE and COM algorithms. The lack of robustness of JADE stems from the fact that only a subset of cross-cumulants are minimized, which means that some cross-cumulants are implicitly maximized along with marginal ones.

4) The mixing matrix influence: Finally, the performance of the seven BSS methods (SOBI, COM1, COM2, JADE, FastICA, FOBIUM and ICAR3) are compared for different mixing matrices. Indeed, figure 6 shows the variations of SINRM 2 at the output of the previous methods as a function of the source spatial correlation c 1,2 , which is defined as the normalized modulus of the scalar product between the two

c 1,2 = 1 N |a 1 H a 2 | (33)
The input SNR of the two sources is assumed to be equal to 15 dB. Similar results have been observed for source 1. More particularly, it appears that ICAR3 presents results generally close to the optimum SMF, except for some isolated values. In addition, FastICA seems to be more sensitive to sources that are not enough angularly separated. On the other hand, this simulation section allows to evaluate the robutness of the previous methods with respect to assumption A5, which is a basic but needed assumption in blind source separation, as shown in figure 6.

VI. CONCLUSION

The ICAR algorithm, exploiting the information contained in the data statistics at fourth order only, has been proposed in this paper. This algorithm allows to process overdetermined (including square) mixtures of sources, provided the latter have marginal FO cumulants with the same sign, which is generally the case in radio communications contexts. Three conclusions can be drawn: first, in the presence of a Gaussian noise spatially and temporally white, the proposed method yields satisfactory results. Second, contrary to most BSS algorithms, the ICAR method is not sensitive to a Gaussian colored noise whose spatial coherence is unknown. Last, the ICAR algorithm is robust with respect to an over estimation of the number of sources, which is not the case for some methods such as JADE. Forthcoming works include the search for a contrast criterion associated with ICAR in order to analyse accurately its performance.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Behavior of BSS methods in the presence of a white noise

Fig. 3 .Fig. 4 .

 34 Fig. 3. Behavior of BSS methods for different signal to noise ratios

Fig. 5 .

 5 Fig. 5. Behavior of BSS methods for an over estimated number of sources

Fig. 6 .

 6 Fig. 6. Behavior of BSS methods for different source spatial correlations

the cyclostationary and cycloergodic case is addressed in section III-D

a J-unitary matrix V is such that V J V H = J where J is a sign diagonal matrix

PROOF OF PROPOSITION 1

The proof of proposition 1 follows immediately from equations [START_REF] Comon | Improved contrast dedicated to blind separation in communications[END_REF], [START_REF] De Lathauwer | ICA algorithms for 3 sources and 2 sensors[END_REF] and assumption (A4), i.e. matrix A does not contain any null entry. In fact, suppose that A⊘A * is not full column rank. Then there exists some P ×1 vector β = 0 such that [A⊘A * ] β = 0, which, due to the structure of A⊘A * [START_REF] Comon | Improved contrast dedicated to blind separation in communications[END_REF] implies that for all 1 ≤ n ≤ N , A * Φ n β = 0. So it implies that A cannot be full column rank (since matrices Φ n are P ×P diagonal with nonzero entries, due to equation [START_REF] De Lathauwer | ICA algorithms for 3 sources and 2 sensors[END_REF] and assumption (A4)), which contradicts the fact that A is of full column rank (A5).

APPENDIX II PROOF OF PROPOSITION 2

The proof is straightforward. In fact, two square roots of a matrix (here Q x ) are always equal to within a unitary matrix, which yields

for some P× P unitary matrix V . Equation [START_REF] Proakis | Digital Communications[END_REF] shows that the right-hand side is the SVD of the left-hand side, hence the proposition 2 result, since

APPENDIX III PROOF OF PROPOSITION 3

A * is a full column rank matrix according to (A5). The diagonal matrices ζ s 1/2 and Φ n (note that the diagonal elements of the latter are components of A), are invertible according to (A3) and (A4) respectively, in other words, because source kurtoses are not null and because matrix A does not contain any null entry. As far as the square matrix V is concerned, it is invertible because of its unitary structure. So matrices Γ n are the product of a full column rank matrix, A * , and an invertible matrix

The fact that this particular product is of full column rank remains to be proved. In fact, suppose that A * F n is not of full column rank. Then there exists some P ×1 vector β = 0 such that A * F n β = 0. So it implies that A * cannot be full column rank (since matrices F n are P ×P invertible), which contradicts the first sentence of this section.

APPENDIX IV PROOF OF PROPOSITION 4

If assumptions (A4) and (A5) are equivalent to assume A with no null entries and of full column rank respectivelly then, proposition 4 may be rewritten as

where s.t. means such that. To prove assertion [START_REF] Taleb | An algorithm for the blind identification of N independent signal with 2 sensors[END_REF], assume the contrary:

which is equivalent, according to [START_REF] De Lathauwer | ICA algorithms for 3 sources and 2 sensors[END_REF], to

This means

In other words, assuming (36) under (A4) and (A5) implies that at least two column vectors of A are collinear, which contradicts (A5). Consequently, assertion [START_REF] Taleb | An algorithm for the blind identification of N independent signal with 2 sensors[END_REF] and hence proposition 4 are true.

Note that B p is a rank one matrix. Consequently, a simple diagonalization of each matrix B * p indeed allows to extract, in a unique way up to a scale and permutation factor, each column vector of A.