Pierre-Vincent Koseleff
email: pierre-vincent.koseleff@imj-prg.fr

Fabrice Rouillier
email: fabrice.rouillier@inria.fr

Cuong Tran
email: cuong.tran@imj-prg.fr

On the sign of a trigonometric expression

Keywords: Chebyshev polynomial, minimal polynomial, cyclotomic polynomial, root isolation

We propose a set of simple and fast algorithms for evaluating and using trigonometric expressions in the form

for a fixed n ∈ Z>0: computing the sign of such an expression, evaluating it numerically and computing its minimal polynomial in Q[x]. As critical byproducts, we propose simple and efficient algorithms for performing arithmetic operations (multiplication, division, gcd) on polynomials expressed in a Chebyshev basis (with the same bit-complexity than in the monomial basis) and for computing the minimal polynomial of 2 cos π n in O(n 2 0 ) bit operations with n0 < n is the odd squarefree part of n. Within such a framework, we can decide if F = 0 in O(d(τ +d)) bit operations , compute the sign of F in O(d 2 τ ) bit operations and compute the minimal polynomial of F in O(n 3 τ ) bit operations, where τ denotes the maximum bitsize of the f k 's.

INTRODUCTION

The present contribution was initially motivated by the computation of Chebyshev knots with the objective of decreasing the complexity bound of the algorithm proposed in [START_REF] Koseleff | The first rational Chebyshev knots[END_REF]. The bottleneck for this algorithm was the computation of the sign of a multivariate polynomial at the zeroes of a zero-dimensional system. It has been shown later ( [START_REF] Koseleff | Computing chebyshev knot diagrams[END_REF]) that this problem finally reduces to the computation of the sign of a polynomial with rational coefficients at some particular real points in the form cos π n , n ∈ Z>0, which can easily be turn into the computation of the sign of an expression in the form F = d k=0 f k cos k π n , f k ∈ Z for a fixed n in Z>0. Evaluating numerically such an expression can be done straightforwardly, but getting its sign (in particular deciding if it is null or not) or computing its minimal polynomial to follow using it for further exact computations need some efforts to be done efficiently in practice.

As

F = f (2 cos π n ) with f (x) = f0 + d k=1 f k T k (x)
, where T k (x) ∈ Z[x] is the k-th monic Chebyshev polynomial, evaluating F can be seen as evaluating f ∈ Q[x] at the real value 2 cos π n , testing if F = 0 can resume to testing if the minimal polynomial Mn ∈ Q[x] of 2 cos π n divides f and, finally, computing the minimal polynomial MF ∈ Q[x] of F . remains to computing the minimal polynomial M f of f in Q[x]/(Mn).

Our primary objective was thus to find a way for computing efficiently Mn ∈ Q[x], the minimal polynomial of 2 cos π n in order, for example, to use state of the art algorithms for testing if F = 0 (which then resumes to a single polynomial division) or computing MF (for example using [START_REF] Shoup | Efficient computation of minimal polynomials in algebraic extensions of finite fields[END_REF]).

The algorithm we propose in section 4 for computing Mn is an adaptation of the algorithm from [START_REF] Arnold | A high-performance algorithm for calculating cyclotomic polynomials[END_REF] that basically computes Φn, the minimal polynomial of the n-th primitive root of unity. The relation between Φn and Mn is simple: Φ2n = D(Mn), where D is the transformation Q[x] defined by D(P ) = x deg P P (x + 1

x ). This simple transformation is also used in Section 2 to set simple and fast algorithms for multiplying, dividing, computing greatest common divisors for polynomials expressed in the Chebyshev basis. In particular, we provide, for the multiplication, a very simple alternative to [START_REF] Giorgi | On polynomial multiplication in chebyshev basis[END_REF] with an equivalent bit complexity. More generally, the complexity of these operations is the same as for polynomials in the monomial basis (see [START_REF] Von Jur Gathen | Modern Computer Algebra[END_REF]).

It turns out that expressing the polynomials in the Chebyshev basis instead of in the usual monomial basis simplifies some parts of the algorithms and make easier the computation of complexity bounds. For example, we show that the computation of Mn can be done in O(n 2 0 ) bit operations where n0 < n is the odd squarefree part of n.

The efficiency of this algorithm depends also on the way the transformations from the Chebyshev basis to the monomial basis and vice-versa are computed. In order not to reinvent the wheel, we have considered these transformations as a particular case of composition method developed in [START_REF] Hart | Practical divide-and-conquer algorithms for polynomial arithmetic[END_REF] with additional refinements linked to the properties of Chebyshev polynomials.

In section 3, we show that passing from the monomial basis to the Chebyshev basis consists essentially to inject in the algorithm from [START_REF] Hart | Practical divide-and-conquer algorithms for polynomial arithmetic[END_REF] the fast operations we have introduced in Section 2. For the reverse operation, we make use of a the linear recurrence introduced in [START_REF] Bostan | Fast conversion algorithms for orthogonal polynomials[END_REF] and essentially keep the structure of the composition algorithm from [START_REF] Hart | Practical divide-and-conquer algorithms for polynomial arithmetic[END_REF], thanks to some remarkable properties of Chebyshev polynomials. We finally show that the two transformations we propose perform O(d) arithmetic operations in Z and O(d(τ + d)) bit operations if τ is the maximal bitsize of the f k 's.

In section 5, we first come back to our initial problem and propose some algorithms together with their complexity analysis. We show that we can test whether F = 0 in O(n 2 + nτ ) binary operations and can decide the sign of F in O(n 2 τ ) binary operations, using an algorithm of [START_REF] Brent | Fast multiple precision evaluation of elementary functions[END_REF] for a fast evaluation of cos x.

As F could also be used inside other algebraic expressions, we propose in section 6 an algorithm for computing efficiently the minimal polynomial of F in O(n 3 τ ) binary operations.

CHEBYSHEV POLYNOMIALS

As seen in the introduction, Chebyshev polynomials and their algebraic properties play a central role in the present contribution. We consider the Chebyshev monic polynomials defined by Tn(x) = 2 cos nt, Un(x) = sin nt sin t where x = 2 cos t. Both of them satisfy the linear recurrence relation:

Pn+1(x) = xPn(x) -Pn-1(x), (1)
but their first terms differ:

{T0(x) = 2, T1(x) = x}, {U0(x) = 0, U1(x) = 1}.
The roots of Tn and Un are real roots and we have:

Tn(x) = n-1 k=0 (x -2 cos (2k+1)π 2n
), (2a)

Un(x) = n-1 k=1 (x -2 cos kπ n ). ( 2b
)
By derivating the relation Tn(

2 cos t) = 2 cos nt, we obtain sin t T ′ n (2 cos t) = n sin nt that is T ′ n = nUn. (3)
By derivating once again we get

4(sin t) 2 T ′′ n (2 cos t) + 2 cos t T ′ n (2 cos t) = 2n 2 cos nt, or (4 -x 2 )T ′′ n (x) -xT ′ n (x) + n 2 Tn(x) = 0. ( 4
)
It means that Tn is an hypergeometric function. By identifying the coefficients, we get

Tn(x) = ⌊ n 2 ⌋ k=0 (-1) k n n-k n-k k x n-2k . (5)
On the other hand, from that x = 2 cos t we get x n = exp it + exp -it n whose expansion gives

x n = ⌊ n 2 ⌋-1 k=0 n k T n-2k (x) + 1+(-1) n 2 , (6)
From Equation (3) and Formula (5), we also obtain

Un+1(x) = ⌊ n 2 ⌋ k=0 (-1) k n-k k x n-2k . (7)
The family T = (1, Tj, j ≥ 1) is therefore a basis of Z[x].

We will also denote by X the monomial basis (x n , n ≥ 0). Given a basis B = (Bi) i≥0 of Q[x] and a polynomial P of degree d, we denote by τB(P ) the maximum bitsize of its coefficients [B]iP in the basis B.

Proposition 1. We have τX (Tn) ≤ n, τX (Un) ≤ n -1, τT (x n ) ≤ n. Proof. From n-k k ≤ n n-k n-k k
, we obtain first that τX (Un+1) ≤ τX (Tn), using Formulas (5) and [START_REF] Giorgi | On polynomial multiplication in chebyshev basis[END_REF]. The sum of the absolute values of the coefficients of Tn is k

n n-k n-k k = φ n + (-1/φ) n ,
where φ is the golden ratio. We thus deduce that [X ] k Tn ≤ ϕ n + 1 ≤ 2 n and τX (Tn) ≤ n. We get

τ (x n ) ≤ n because that k n k = 2 n .
We thus deduce Corrolary 2. Let P be a polynomial of degree d. Then we have τX (P ) ≤ τT (P ) + d + log 2 d, and τT (P ) ≤ τX (P ) + d + log 2 d.

Chebyshev forms

Let f ∈ Q[x]
. We will say that f is in Chebyshev form when it is expressed in the Chebyshev basis: [T ]f = f0 + d j=1 fjTj. The set of Chebyshev forms with rational coefficients is an Euclidean ring and we focus on computing efficiently in this ring, starting with elementary operations.

Let us introduce a simple shortcut that will straightforwardly offer a simple link between basic operations in the Chebyshev basis and in the usual monomial basis:

Definition 3. D is the transformation from Q[x]
to the space of even degree reversible polynomials given by

D(P ) = x deg P P (x + 1 x ).
Note that D satisfies the following relation:

D(f • g) = D(f ) • D(g). (8)
so that, as Tn(x

+ 1 x ) = x n + 1 x n , then D(a0 + n j=1 ajTj) = a0x n + n j=1 aj(x n-j + x n+j ). (9)
Note that Equation [START_REF] Hart | Practical divide-and-conquer algorithms for polynomial arithmetic[END_REF] shows that passing from [X ]D(P ) to [T ]P and vice versa is just a way of changing the numbering of the coefficients, having in mind that [X ]D(P ) is a reversible polynomial. In particular, τX (D(P )) = τT (P ) and deg D(P ) = 2 deg P . Also, passing from [X ]D(P ) to [T ]P and vice-versa can be done in O(d) bit operations.

From this remark, it is then sufficient to apply the classical fast multiplication for polynomials in Q[x] after applying D on the operands. We will denote by M (δ, τ ) = O(δτ ) the number of bit operations for multiplying two polynomials of Z[x] with degrees at most δ and coeffitients of bitsize at most t with respect to the usual monomial basis X and by M (τ ) = O(τ ) the cost of the multiplication of two integers of bitsize dominated by τ . We obtain

Proposition 4. Let f, g ∈ Z[x] given by their Chebyshev forms [T ]f and [T ]g. Then we have τT (f • g) ≤ τT (f ) + τT (g) + log 2 min(deg f, deg g) + 1.
Suppose that τT (f ), τT (g) ≤ τ and deg f,

deg g ≤ d then computing [T ](f • g) claims M (2d, τ ) + O(d) bit operations.
Note that this very simple algorithm might be replaced by the one from Giorgi ( [START_REF] Giorgi | On polynomial multiplication in chebyshev basis[END_REF]) which claims the equivalent of two multiplications between polynomials of degree d with coefficients of bitsize in O(τ ).

Let us now extend these results to the division. Applying the transformation D on both side of the relation f = g•q+r, we then get:

D(f ) = D(g)D(q) + x deg f -deg r D(r). (10)
In particular, Formulas ( 8) and [START_REF] Hsiao | On factorization of chebyshev's polynomial of the first kind[END_REF] imply that g | f iff D(g) | D(f ). Also, applying the same strategy as for the product of Chebyshev forms, we immediately get:

Corrolary 5. Let f, g ∈ Z[x] given by their Chebyshev forms [T ]f and [T ]g. Suppose that deg f, deg g ≤ d and τT (f ), τT (g) ≤ τ . Then, 1. (a) deciding if g divides f knowing [T ]f and [T ]g can be done in O(d 2 + dτ ) bit operations; (b) if g divides f then computing [T ]Quo(f, g) from [T ]f and [T ]g can be done in O(d 2 + dτ ) bit op- erations; (c) computing [T ]gcd(f, g) from [T ]f and [T ]g can be done in O(d 2 τ ) bit operations. 2. (a) if g|f then τT (g) = O(τT (f ) + d). (b) τT (gcd(f, g)) ≤ O(τ + d).
Proof. The first items are direct consequences of [16, exercice 9.14] while the third item is a direct consequence of [START_REF] Von Jur Gathen | Modern Computer Algebra[END_REF]Corollary 11.20]. The two last items come from the Mignotte bound ([16, p. 166]).

BASES CHANGE

It is a classical result that T is an orthogonal basis. Bostan et al. (see [START_REF] Bostan | Fast conversion algorithms for orthogonal polynomials[END_REF]) give a fast general conversion method from the monomial basis to orthogonal bases. The alternatives we propose are based on the general method of Hart & Novocin for polynomial composition (see [START_REF] Hart | Practical divide-and-conquer algorithms for polynomial arithmetic[END_REF]) which lead to simplier descriptions and complexity analysis. Given f = d i=0 aix i in the Chebyshev basis and f = f0 + d i=1 fiTi in the monomial basis, in both cases, the method is based on the same divide and conquer principle as in [START_REF] Hart | Practical divide-and-conquer algorithms for polynomial arithmetic[END_REF]: splitting f into ⌈ d 2 ⌉ sub-polynomials of length 2.

Suppose that we have to evaluate the sum S = d j=0 ujv j . We first define:

d0 = d, v0 = v, S0,j = uj, j = 0, . . . , d.
Let c = ⌈log 2 d⌉ -1. Then we define for k = 0, . . . , c:

d k+1 = ⌈ d k 2 ⌉, v k+1 = v 2 k , S k+1,j = S k,2j + S k,2j+1 • v k , j = 0, . . . , d k+1 . Then, for k = 0, . . . , c, S = d k j=0 S k,j • v j k .
At the end, we get S = Sc,0.

1: From X to T Input: f = n i=0 aix i Output: [T ]f = f0 + n i=1 fiTi 1 begin 2 d ← ⌊ n 2 ⌋; c ← 0; v ← T1 ; 3 for j ← 0 to d do 4 Sc,j ← a2j + a2j+1 • T1 5 v ← ([T ]v) 2 ; 6 while d > 0 do 7 d ← ⌊ d 2 ⌋; 8 for j ← 0 to d do 9 Sc+1,j ← [T ]Sc,2j + [T ]Sc,2j+1 • [T ]v 10 c ← c + 1; v ← ([T ]v) 2 ; 11 return Sc,0 Proof. We use Algorithm 1. Let c = ⌈log 2 d⌉. We start with v = [T ]x = T1, and compute successively v = x 2 k , k = 0, . . . , c. We know that x 2 k has size O(2 k ) in the Chebyshev basis. From Proposition 4, one can compute [T ]x 2 k+1 from [T ]x 2 k in O(2 2k ) bit operations.
At each step, the bitsize τ k+1 of the Chebyshev form

S k+1,j = [T ]S k,2j + [T ]S k,2j+1 • [T ]x 2 k satisfies τ k+1 ≤ τ k + 2 k + k + 1. We therefore obtain τ k = τ + O(2 k ). One can compute the Chebyshev form S k+1,j = S k,2j + S k,2j+1 • x 2 k in O(2 k τ k ) binary operations from Proposition 4.
The total computing time is bounded by

c k=1 τ k 2 k = O(d 2 + dτ ).
Conversely, given a Chebyshev form f = f0 + d j=1 fjTj, we have to find the expression f = d j=0 aix i . We use a key idea from the general method of [START_REF] Bostan | Fast conversion algorithms for orthogonal polynomials[END_REF], considering the characteristic matrix of the Chebyshev polynomials:

X = 0 1 -1 x .
Using the recurrence (1), we have:

Tj Tj+1 = X j T0 T1 , j ≥ 0.
Algorithm 2 computes the matrix polynomial

M (f ) = ⌊ d 2 ⌋ i=0 f2i f2i+1 X 2i
using the composition method in [START_REF] Hart | Practical divide-and-conquer algorithms for polynomial arithmetic[END_REF]. Then, following [START_REF] Bostan | Fast conversion algorithms for orthogonal polynomials[END_REF], we get (because T0 = 2) Proof. An induction shows that X n = -Un-1 Un -Un

f = -f0 + M (f ) • T0 T1 . Proposition 7. Let f = f0 + d i=1 fiTi
Un+1 .

Thus the bitsize of X n is O(n), using Proposition 1.

Algorithm 2: From T to X Input: f = n i=0 fiTi Output: [X ]f = n i=0 aix i 1 begin 2 d ← ⌊ n 2 ⌋; c ← 0; V ← X 2 ; 3 for j ← 0 to d do 4 Sc,j ← [f2j f2j+1] 5 while d > 0 do 6 d ← ⌊ d 2 ⌋; 7 for j ← 0 to d do 8 Sc+1,j ← Sc,2j + Sc,2j+1 • V 9 c ← c + 1; V ← [X ]V 2 ; 10 return Sc,0 • 2 x Let c = ⌈log 2 d⌉ and 0 < k ≤ c. We compute X 2 k from X 2 k-1 in O(2 2k
) binary operations. We thus compute X 2 k for k = 0, . . . , c in O(d 2 ) binary operations.

The bitsize τ k+1,j of the 1

× 2 matrix S k+1,j = S k,2j + S k,2j+1 • X 2 k is bounded by τ k,j + O(2 k ). It follows that τ k,j ≤ τ k ≤ τ + O(2 k ). At the end the bitsize of [X ]f is bounded by τ + O(d).
S k+1 is deduced from S k by evaluating the S k+1,j for j = 0, . . . ,

⌊ d 2 k ⌋. This claims d/2 k × O(τ k 2 k ) = O(d(τ + 2 k )) binary operations.
At the end we compute Sc,0 in O(d 2 + dτ ) binary operations.

THE MINIMAL POLYNOMIAL

The degree of the minimal polynomial of cos 2π n is 1 2 ϕ(n) (see Rivlin, [14, Chapter 5] or Watkin & Zeitlin , [START_REF] Watkins | The minimal polynomial of cos 2π n[END_REF]). This minimal polynomial appears as a factor of the Chebyshev monic polynomials in [START_REF] Rivlin | Chebyshev Methods in Numerical Approximation[END_REF]. It is a factor of

T ⌊ n 2 ⌋+1 -T ⌊ n 2 ⌋
in [START_REF] Watkins | The minimal polynomial of cos 2π n[END_REF]. Bayard & Cangul ( [START_REF] Bayad | The minimal polynomial of 2 cos π q and dickson polynomials[END_REF]) use this formula to get an induction formula for the minimal polynomial of 2 cos π n with the help of the Möbius inversion formula.

Here, we will show that the minimal polynomial Mn of 2 cos π n satisfies Φ2n = D(Mn). We will adapt the algorithm given in [START_REF] Arnold | A high-performance algorithm for calculating cyclotomic polynomials[END_REF] for the computation of Φn to the computation of Mn in the Chebyshev basis.

Definition 8. Let n > 2, Mn denotes the minimal poly- nomial of 2 cos π n in Q[x]. We set M1 = 1.
Let ζn be a primitive n-th root of the unity, its minimal polynomial is the cyclotomic polynomial Φn. It is monic and of degree ϕ(n) where ϕ is the Euler function. The roots of Φn are all the primitive n-th roots and we deduce that Φn(1/ζn) = 0 and thus Φn is reversible. We have Q(cos 2π 2n ) = Q(ζn) ∩ R and the minimal polynomial over Q of cos 2π 2n has degree 1 2 ϕ(2n), when n > 2. Its roots are the cos kπ n where k and 2n are coprime. Note that the minimal polynomial of cos 2 π n is ±Mn(-t).

lemma 9. If m ≥ 3 is odd, then M 2 k m = Mm(T 2 k ). Proof. As Mm • T 2 k (cos π 2 k m ) = 0 and (2 k , m) = 1, then M 2 k m |Mm(T 2 k ). We conclude since M 2 k m and Mm(T 2 k ) have same leading term.
The number of factors of Tn is given in [START_REF] Hsiao | On factorization of chebyshev's polynomial of the first kind[END_REF]. We give here the complete factorization of the Chebyshev polynomials Tn and Un as products of polynomials Mn, see also [START_REF] Koseleff | Computing chebyshev knot diagrams[END_REF]:

Proposition 10. We have

U 2 k (2m+1) = (-1) m d|2m+1 M d (-t) k i=0 M d (T 2 i ) T 2 k (2m+1) = d|2m+1 M d (T 2 k+1 )
where Mn is the minimal polynomial of cos π n . Proof. The factorization of Un is obtained by comparing its roots with those of M d (±t), when d|2m + 1 and using the Euler identity n = d|n ϕ(d).

Let d be an odd divisor of n. We write n = 2 k •d1 •d, where d1 is odd. cos d 1 π 2n = cos π 2 k+1 d is a root of Tn so M 2 k+1 d |Tn. We conclude by comparing the leading terms.

From Mn(ζ2n +1/ζ2n) = 0 when n > 2, we deduce a relation between the Chebyshev form of Mn and Φ2n, see also [14, Chapter 5]: Proposition 11. We have Φ2n = D(Mn), for n > 2.

Let x = 2 cos t and Vn(x

) = sin(n + 1 2 )t cos t 2 then Vn = d|2n+1 M d .
This is an analogous formula to X n -1 = d|n Φ d , see [START_REF] Koseleff | The first rational Chebyshev knots[END_REF].

Computation of Mn

The relation Tn •Tm = Tn •Tm = Tmn makes the Chebyshev basis very convenient for the computation of the minimal polynomial Mn. Lemma 12 is a generalization of Lemma 9:

lemma 12. Let n = 2 α p α 1 1 . . . p α k k
where the pi are distinct odd primes. Let n0 = p1 • • • p k . Then we have

Mn = Mn 0 (T n/n 0 ), if k ≥ 1 M2α = T 2 α-1 .
Proof. If n = 2 α > 1 then T n/2 (2 cos π/n) = 2 cos π/2 = 0. On the other hand T n/2 and Mn are monic with the same degree. Thus they are equal.

If n0 > 1 then Mn 0 (T n/n 0 (2 cos π/n)) = Mn 0 (2 cos π/n0) = 0. We thus have Mn|Mn 0 (T n/n 0 ). On the other hand, Mn and Mn 0 (T n/n 0 ) are monic with the same degree. Thus they are equal.

Lemma 12 shows that the Chebyshev basis is particularly adapted for the computation of Mn.

lemma 13. n = 2 α p α 1 1 . . . p α k k
where the pi are distinct odd primes. Let n0 = p1 • • • p k . Then the coefficients of Mn (in the monomial basis, as well as in the Chebyshev basis) have bitsize O(n0).

Proof. As Mn 0 |Un 0 the conclusion is a consequence of Proposition 5, because we deduce from

Un = ⌊ n 2 ⌋ k=1 T n-2k-1 + 1 2 (1 -(-1) n ),
that the Chebyshev form Un has maximum bitsize 1.

When n is an odd prime, we get : lemma 14. Let p a prime odd number, then

Mp = T k -T k-1 + • • • + (-1) k-1 T1 + (-1) k ,
where k = 1/2(p -1).

Proof. Let p = 2k + 1 and h

= T k -T k-1 + • • • + (-1) k-1 T1 + (-1) k . We have sin π 2k+1 • h(2 cos π 2k+1 ) = 0. But deg h = k = 1
2 ϕ(p) and thus h = Mp. lemma 15. Let p be an odd prime and n such that p ∤ n, then Mnp = Mn(Tp) Mn .

Proof. We just need to show that both Mnp and Mn divide Mn(Tp) because these polynomials are monic and their degrees satisfy 2(deg Mn

+ deg Mnp) = ϕ(2n) + ϕ(2pn) = pϕ(n) = 2 deg Mn(Tp).
If t = 2 cos π np then Tp(t) = 2 cos π n and Mn(Tp)(t) = 0 so that Mnp|Mn(Tp). If t = 2 cos π n then Tp(t) = 2 cos pπ n is a root of Mn.

Following Arnold and Monagan in [START_REF] Arnold | A high-performance algorithm for calculating cyclotomic polynomials[END_REF] for the computation of Φn, we deduce Algorithm 3 to compute the Chebyshev form of Mn: We have to compute successively

Algorithm 3: Compute [T ]Mn Input: n Output: [T ]Mn 1 begin 2 Factorize n = 2 α 0 p α 1 1 • • • p α k k ; n0 ← p1 • • • p k ; 3 if k = 0 then 4 if α0 = 0 then 5 return 1 6 else 7 return T 2 α 0 -1 8 i ← 1 2 (p1 -1); 9 m ← Ti -Ti-1 + • • • + (-1) i-1 T1 + (-1) i ;
Mq i+1 = (Mq i •Tp i+1 )/Mq i , where qi = p1 • • • pi.
We know that τ

(Mq i ) = O(qi), τ (Mq i (Tp i+1 )) = τ (Mq i ) and deg(Mq i (Tp i+1 )) = ϕ(2qi)pi+1/2 ≤ qi+1. By Proposi- tion 4, we compute Mq i+1 from Mq i in O(qi+1) arithmetic operations or O(q 2
j+1 ) binary operations. Finally we compute Mn 0 in O(q2

+ • • • + q k ) = O(n0) arithmetic operations and O(q 2 1 + • • • + q 2 k ) = O(n 2 0 ) binary operations.
The analogy between Mn and the cyclotomic polynomial Φn allows us to deduce the following complexity bound for the computation of Φn from the algorithm proposed by Arnold and Monagan in [START_REF] Arnold | A high-performance algorithm for calculating cyclotomic polynomials[END_REF] :

Corrolary 17. One can compute Φn in the complexity O(n0) and in the running time O(n 2 0 ), where n0 is the odd squarefree part of n.

Example

The minimal polynomial M1260 is

M105(T12) = T288 -T276 + T264 + T228 -T216 + 2 T204 -T192 + T180 + T144 -T132 + T120 -T108 +T96 -T84 -T48 -T24 -1.
The maximum bitsize of [X ]M1260 is τ ([X]128M1260) = 197 while the total bitsize of its coefficients is 20329.

The minimal polynomial M936 is

M39(T24) = T288 + T264 -T216 -T192 +T144 + T120 -T72 -T48 + 1.
Its maximum bitsize is also 197 while the total bitsize of its coefficients is 20338.

SUM OF COSINES

The aim of this section is to give methods for determining the sign of F = f (γ) where γ is one root of Mn, for evaluating F and for getting the minimal polynomial of F .

We consider here the expression

F = f0 + 2 d j=1 fj cos jπ n , (11)
where d ≤ n and fi ∈ Z. It is convenient to write F as f (γ) where γ is the algebraic number 2 cos π n and f = f0

+ d k=1 f k T k ∈ Z[x].

Evaluating the sum

The first step is the approximation of γ = 2 cos k π n , using Brent's methods ( [START_REF] Brent | Multipleprecision zero-finding methods and the complexity of elementary function evaluation[END_REF][START_REF] Brent | Fast multiple precision evaluation of elementary functions[END_REF])

lemma 18. Let 0 ≤ k ≤ n and γ = 2 cos k π n . Let ℓ ∈ Z>0. One can compute c ∈ Q, of bitsize τ (c) ≤ ℓ such that |c -γ| ≤ 2 -ℓ in O(ℓ + log n) bit operations. Proof. We first compute r ∈ Q of bitsize 2(ℓ + 2 log n) such that r -k π n < 2 -ℓ-1 .
It claims O(l + log n) bit operations, using [START_REF] Brent | Multipleprecision zero-finding methods and the complexity of elementary function evaluation[END_REF].

Then we compute c such that |c -cos r| < 2 -ℓ-1 , using [START_REF] Brent | Fast multiple precision evaluation of elementary functions[END_REF], in O(ℓ + log n) binary operations.

We thus deduce

Corrolary 19. Let 0 ≤ k ≤ n and γ = 2 cos k π n . Let ℓ ∈ Z>0. Let f be the Chebyshev form f0 + n-1 i=1 fiTi. One computes F ∈ Q of bitsize O(nτ (f ) + ℓ) such that F -f (γ) ≤ 2 -ℓ in O(nℓ + nτ ) bit operations. One computes F -and F + of bitsize O(nτ (f ) + ℓ) such that F -≤ F ≤ F + and F + -F -≤ 2 -ℓ in O(nℓ + nτ ) bit operations.
Proof. Using Lemma 18, one can compute the rational numbers ci, i = 1, . . . , n -1 such that ci -cos ik π n ≤ 2 -ℓ-τ /n. They are of bitsize O(τ + ℓ + log n) and we can compute them in O(nτ + nℓ) bit operations.

If F = f0 + 2 n-1 i=1 fici, then F -F ≤ 2 n-1 i=1 |fi| ci -cos ik π n ≤ 2 -ℓ .
Here, the fici have bitsize τ + ℓ + log n and each can be computed in O(τ + ℓ + log n) bit operations, using the fast multiplication described in [START_REF] Von Jur Gathen | Modern Computer Algebra[END_REF]Chap. 8]. F has bitsize O(τ + ℓ + log n) and can be computed in O(nτ + nℓ) bit operations.

The proof of the second part of the Lemma is similar.

Evaluating the sign of the sum

Here we want to evaluate the sign of f (γ) where γ = 2 cos k π n . and f is the Chebyshev form f0 + n-1 i=1 fiTi. Without loss of generality, we can suppose that (k, 2n) = 1, that is to say that γ is a root of Mn. If (k, n) = d > 1, then we can change n by n/d. If k is even, then γ will be a root of Mn(-t).

As Mn is the minimal polynomial of γ then f (γ) = 0 iff Mn | f . If Mn does not divide f then we have |Res (Mn, f )| ≥ 1. As Mn is monic, we have also

Res (Mn, f ) = Mn(γ)=0 f (γ).
Let γ be a root of Mn, then we have

|f (γ)| ≤ ||f || 1 , where ||f || 1 = |f0| + 2 d i=1 |fi|.
We thus obtain for every root γ of Mn:

|f (γ)| ≥ (||f || 1 ) 1-deg Mn = δ. (12)
Then, it is then sufficient to compute an approximation F of f (γ) such that F -f (γ) < δ 2 in order to ensure that this approximation has the same sign than f (γ). We thus deduce Proposition 20. Let f of degree d < n with τT (f ) ≤ τ . Let γ = 2 cos k π n where (k, 2n) = 1.

1. we can decide whether f (γ) = 0 in O(n 2 + nτ ) bit operations, 2. we can compute sign f (γ) in O(n 2 τ ) bit operations, 3. we can evaluate f (γ) with precision 2 -ℓ in the complexity O(nℓ + nτ ).

Proof. We first decide if f (γ) = 0 by testing if Mn|f . It claims O(n 2 + nτ ) bit operations by Propositions 16 and 5.

If f (γ) = 0, then we compute an approximation F of f (γ) with accuracy 1 2 δ where δ = (||f || 1 ) 1-deg Mn . But -log 2 δ ≤ (deg Mn -1)(log 2 (2n + 1) + τ ) = O(nτ ). We deduce, using Corollary 19, that F can be computed in O(n 2 τ ) bit operations.

The last point of the Proposition is a consequence of Corollary 19.

The estimate of δ in Equation ( 12) may be too sharp. The alternative we propose is to explicitly compute [T ]Mn, then test if [T ]Mn divides F using Corollary 5, and then compute F using interval arithmetic, starting in low precision, and doubling the precision if the resulting interval contains 0. According to the above result, Algorithm 4 will end, performing O(n 2 τ ) bit operations in the worst case, but will use much less precision in generic situations.

Algorithm 4: Determine the sign of f (γ) Indeed, min Mn(x)=0 |f (x) -3| ≥ 2 -29 and Algorithm 4 provides the value and the sign of f (γ) very quickly.

Input: f = f0 + n i fiTi, k such that (k, 2n) = 1 Output: sign f (γ) where γ = 2 cos k π n ) 1 begin 2 Compute [T ]Mn ; 3 if Mn | f then 4 return f (2 cos k π n ) = 0 5 ℓ ← 1, F -← -1, F + ← 1; 6 while F + • F -< 0 do 7 Compute [F -, F + ] of length 2 -ℓ , containing f (γ);

THE MINIMAL POLYNOMIAL OF A SUM

We give an algorithm for the minimal polynomial of F = f (γ). It is classical to consider the monic polynomial (because Mn is monic)

P f (z) = Res x(z -f (x), Mn(x)) = Mn(γ)=0
(z -f (γ)). [START_REF] Myerson | Rational products of cosine of rational angles[END_REF] The roots of PF are

z k = f (2 cos k π n ) = f (T k (2 cos π n ))
where (k, 2n) = 1. We thus deduce that Proposition 21. There exists ν ∈ N * such that P f = M ν f , where M f is the minimal polynomial of f (2 cos π n ). Proof. Let γ k = 2 cos k π n , be a root of Mn. The minimal polynomial µ k of f (γ k ) satisfies µ k • f (γ k ) = 0. Thus Mn divides µ k • f and µ k (f (γj)) = 0 for every root γj of Mn. We thus deduce that µj divides µ k and that µ k = µj = M f . P f factorizes into P f = M ν f .

We will compute M f by computing first P f and then M f = P f / gcd(P f , P ′ f ). On might directly compute P f using a general algorithm, following [START_REF] Von Jur Gathen | Modern Computer Algebra[END_REF]Section 11.2]. It would then require O(δτs) bit operations where δ = deg(f ) + deg(Mn) and τs is a bound on the total bitsize required to store any principal subresultant coefficient which are all polynomials of degree in O(δ) in z with coefficients of bitsizes in O(δτ ) (see Proposition 8.50 in [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]). Such an algorithm would then perform, in our case, O(n 3 τ ) bit operations, and we do not pretend to give a better theoretical upperbound. However, in pratice, this classical strategy is based on the Half-Gcd algorithm and makes use if fast operations on univariate polynomials, both starting to be being efficient for rather high degrees (several hundred). The algorithm we propose in the sequel runs efficiently even for small degrees.

It is natural to consider P f as a polynomial whose coefficients are in Q[x]/(Mn).

As Tn+i(γ) = 2 cos(n + i)k π n = (-1) k Ti(γ) we have Tn+i ≡ Tn-i ≡ -Ti (mod Mn).

The Chebyshev basis (1, T1, . . . T ⌊ n-1 2 ⌋ ) is then particularly adapted as generating family of Q[x]/(Mn):

lemma 22. f = f0+ ⌊ n-1 2 ⌋ i=1 fiTi and g = g0+ ⌊ n-1 2 ⌋ i=1 giTi be Chebyshev forms. Then one can compute h = h0 + ⌊ n-1 2 ⌋ i=1 hiTi where h ≡ f • g (mod Mn) in O(nτ ) bit op- erations and τ (h) ≤ τ (f ) + τ (g) + log 2 n + 1. Proof. We compute h = [T ](f • g) in O(nτ ) bit opera- tions. We obtain h = h0 + n-1 i=1 hiTi = h0 + ⌊ n-1 2 ⌋ i=1 (hi - hn-i)Ti.
Let N = 1 2 ϕ(2n) be the degree of P f . The N roots of PF are the f (T k (x)) where (k, 2n) = 1.

We shall consider the i-th Newton sum:

Si(P f ) = 1≤k≤n (k,2n)=1 f (T k (x)) i . Consider first f i = f (x) i = fi,0+ ⌊ n-1 2 ⌋ j=1 fi,jTj in Q[x]/(Mn), then f (T k (x)) i ≡ fi,0 + ⌊ n-1 2 ⌋ j=1
fi,jT jk (mod Mn). We thus deduce that

Si(P f ) ≡ N fi,0 + ⌊ n-1 2 ⌋ j=1 fi,j 1≤k≤n (k,2n)=1
T jk (mod Mn) lemma 23. Let j and n be nonnegative integers, then 1≤k≤n (k,2n)=1

T jk ≡ 2Sj(Φ2n) (mod Mn)

Proof. Let γ = 2 cos π n then 1≤k≤n (k,2n)=1 T jk (γ) = 2 Re 1≤k≤n (k,2n)=1 e jki π n = 2Sj(Φ2n)
because Sj(Φ2n) is real.

Note that the sums Sj(Φn) are known as Ramanujan sums, see Remark 25. We thus deduce that one can obtain Sj(P f ) from coefficients fi,j of

f j (x) in Q[x]/(Mn) with Si(P f ) = fi,0N + 2 ⌊ n-1 2 ⌋ j=1 fi,jSj(Φ2n) (14)
Equation ( 14) gives a linear relation between the Newton sums of Φ2n and the Newton sums of P f . The associated matrix is (fi,j) = [I]j(f i ). Note that (fi,j) are not unique but the Formula ( 14) does not depend on them. We thus deduce Algorithm 5.

Algorithm 5:

Compute P f Input: [T ]f = f0 + n-1 j=1 fiTi Output: P f = Res x(z -f (x), Mn(x)) 1 begin 2 Compute [T ]Mn = σ0 + N i=1 σiTi; 3 for i ← 2 to N do 4 Compute the Newton sums Si = iσi -i-1 j=1 σjSi-j 5 for i ← 2 to N do 6 Compute f i = fi,0 + ⌊ n-1 2 ⌋ j=1 fi,jTj; 7 Compute the Newton sums Si(P f ) = fi,0 • N + 2 ⌊ n-1 2 ⌋ j=1 fi,jSj; 8 Compute σi(P f ) = - 1 i Si(P f ) + i-1 j=1 σi-j(P f )Sj(P f ) 9 return z N + N i=1 σi(P f )z N -i
Proposition 24. We can compute P f in O(n 3 τ ) bit operations.

Proof. We first compute [T ]Mn in O(n 2 ) bit operations. We obtain [T ]Mn = σ0 + N i=1 σiTi. As Φ2n = D(Mn) we deduce that Φ2n = 2N i=0 bix i where bi = b2N-i = σi. We then obtain the Si(Φ2n) using the Newton relations:

S k + σ1S k-1 + • • • + σ k-1 S1 + kσ k = 0, 1 ≤ k ≤ N. (15)
Here the σi's have bitsize O(n). We thus deduce by induction that all the S k have bitsize O(n). S k is computed from S1, . . . , S k-1 in n O(n) bit operations. In conclusion S1, . . . , SN may be computed in O(n 3 ) bit operations.

By induction f i (mod Mn) has bitsize O(iτ ) and can be computed from f i-1 (mod Mn) in O(inτ ) bit operations, using Lemma 22. We can compute all fi,j, 1 ≤ i ≤ N , 1 ≤ j ≤ ⌊ n-1 2 ⌋ in O(n 3 τ ) bit operations. We then compute each Sj(P f ) in O(n 2 τ ) bit operations, using Formula [START_REF] Rivlin | Chebyshev Methods in Numerical Approximation[END_REF]. Il claims O(n 3 τ ) bit operations.

We then use the Newton formula [START_REF] Shoup | Efficient computation of minimal polynomials in algebraic extensions of finite fields[END_REF] to get the coefficients σ1(P f ), . . . , σN (P f ) of P f :

kσ k (P f ) = - k-1 i=1
Si(P f )σ k-i (P f ). Also, all kσ k have bitsize O(nτ ) and may be computed in O(n 3 τ ) bit operations. As σ k ∈ Z they can all be obtained in O(n 3 τ ) bit operations Remark 25. The Ramanujan sums are known (see [START_REF] Hardy | An Introduction to the Theory of numbers[END_REF]Th. 272]) to satisfy Sm(Φn) = µ(n/ gcd(n, m)) ϕ(n) ϕ(n/ gcd(n, m)) ,

where µ is the Möbius function and ϕ is the Euler function.

It means that Sm(Φn) may be computed in O(n) binary operations (see [START_REF] Von Jur Gathen | Modern Computer Algebra[END_REF]Chap. 19]). But this does not allow to compute Φn with the Newton formulas, which would require O(n 3 ) binary operations.

Once P f is computed, we can compute M f as P f / gcd(P f , P ′ f ).

Here P f has degree N and bitsize O(nτ ). Then h = gcd(P f , P ′ f ) is calculated in time O(n 3 τ ), using Proposition 5.

Using Proposition 5, we know that τT (h) = O(nτ ) so that the last operation M f = P f /h claims O(n 2 τ ) bit operations.

Finally we deduce Proposition 26. We can compute the minimal polynomial M f of f (2 cos π n ) in O(n 3 τ ) binary operations.

Example

Consider F = 2a cos π 7 + 2b cos π 5 . We write F = f (2 cos π 35 ) where f = aT5 + bT7.

The minimal polynomial of 2 cos π 35 is M35 = T12 + T11 -T7 -T6 -T5 -T4 + T2 + T1 + 1 = t 12

Algorithm 1 Proposition 6 .

 16 is a direct application of this strategy for the evaluation of f = d j=1 ajx j in the Chebyshev basis T . Let f = d i=0 aix i be a polynomial of bitsize τ . One can compute the Chebyshev form of f in O(d) arithmetic operations, or in binary complexity O(d 2 + dτ ). The Chebyshev form [T ]f has bitsize τ + O(d).

Algorithm

  be a Chebyshev form of bitsize τ . One can compute f in the monomial basis in O(d) arithmetic operations, or a binary complexity O(d 2 + dτ ). [X ]f has bitsize τ + O(d).

10 for j ← 2 to k do 11 m

 1011 ← [T ] m(Tp j ) m 12 return m(T n/n 0 ) Proposition 16. One computes Mn in the Chebyshev basis in complexity O(n0) or O(n 2 0 ) bit operations. Proof. We will use Algorithm 3. The factorization of n requires O(log n) binary operations, see [16, Chap. 19].

1 =

 1 to 3. We write a as f (γ) where f = -T771 -T459 -T451 -T277 + T251 + T243 + T69 -T43, and γ = 2 cos π 936 is the biggest root of M936. Equation (12) asserts that|f (γ) -3| ≥ ||f -3|| -288 19 -288 ≃ 2 -1224 .

Computing 4 + 5 2 - 3

 4523 P f = Res (z -aT5 -bT7, M35) we obtain P f = M 2f whereM f = z 6 -(3 b + 2 c) z 5 + c (-3 c + 5 b) z b 3 + 6 bc 2 + 6 c 3 z 3 -c 5 b 3 + 7 b 2 c + 9 bc 2 -2 c 3 z b 5 -4 b 3 c 2 -7 b 2 c 3 + 2 bc 4 + 4 c 5 z -b 6 + b 5 c + 7 c 2 b 4 -2 b 3 c 3 -7 b 2 c 4 + 2 bc 5 + c 6 .

  + t 11 -12 t 10 -11 t 9 + 54 t 8 + 43 t 7 -113 t 6 -71 t 5 + 110 t 4 + 46 t 3 -40 t 2 -8 t + 1.


3D Head-Mounted Antenna Array Architecture Optimization based on the Fisher Information Matrix

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Abstract-This paper introduces a method to optimize the configuration of a 3D helmet-mounted antenna array carried by emergency rescuers expected to locate natural disaster survivors. Configuration optimization is based on metrics extracted from the 3D Fisher Information Matrix (FIM) relative to the considered array architecture for single source Direction Of Arrival (DOA) estimation. Comprehensive simulations illustrate the optimization metrics behaviour with regard to survivor's emitter 3D DOA and array configuration. Design constraints related to helmetmounted antenna requirements reduce the usually exhausting FIM based design optimization to a 2 degrees of freedom exploration of the array DOA estimation performances.

Index Terms-Direction Of Arrival Estimation, Fisher Information Matrix, Helmet-mounted Antenna Array

I. INTRODUCTION

To efficiently rescue survivors in the aftermath of earthquakes or tsunamis, it is of paramount importance to know where to aim the research once the emergency team is deployed on the field. In these situations, one may reasonably assume that most survivors carry a mobile phone, even one without GPS capability in the worst-case scenario. In such conditions, the rescue protocol requires a local activation of these mobile phones to allow Direction Of Arrival (DOA) estimation process. These estimations are based on the electromagnetic measurements collected through the antenna array carried by the rescuers. To individually reach and assist as many survivors as possible, rescuers would need to be equipped with handfree antenna arrays that could estimate their DOAs on the go [START_REF] Fischer | Location and navigation support for emergency responders: A survey[END_REF].

While very well known high resolution DOA estimation methods (MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], ESPRIT [START_REF] Roy | Esprit estimation of signal parameters via rotational invariance techniques[END_REF], WSF [START_REF] Viberg | Subspace Fitting Concepts in Sensor Array Processing[END_REF], etc.) exist in the literature along with extensive studies concerning 2D and 3D antenna array optimization ( [START_REF] Manikas | Comparison of the ultimate direction-finding capabilities of a number of planar array geometries[END_REF], [START_REF] Verhaevert | Uniform spherical distributions for adaptive array applications[END_REF]), design optimization of helmet-mounted devices remains a challenge. The involved dimensions and geometry constraints obviously raise several issues compared to usually considered arrays. This paper is part of a broad project that assumes a network operator substitution on the local wireless network, which ensures the ability to activate one handset at a time therefore reducing the considered scenario to a single source 3D DOA estimation in a noised environment. A few papers tackle the 3D antenna array optimization problem based on the Fisher Information Matrix (FIM) in the single source case ( [START_REF] Ferréol | High resolution direction finding: From performance toward antenna array optimization-the mono-source case[END_REF], [START_REF] Gazzah | Cramer-rao bounds for antenna array design[END_REF]), but for general array configurations not adapted to hereafter considered architecture. This paper thus introduces a dedicated optimization method for 3D helmet-mounted antenna arrays based on metrics derived from the 3D FIM.

The remainder of the paper is organized as follows. Section 2 presents the geometrical configuration of the helmetmounted antenna array and introduces the considered signal model. Section 3 details the geometrically oriented form of the 3D FIM for the considered model while its derivation is given in Appendix A. Section 4 suggests antenna array architecture optimization metrics extracted from the 3D FIM. Finally, Section 5 presents simulation results and analysis on suggested metrics with relation to DOA angle and array configuration while Section 6 holds the conclusion of this paper. We consider the 3 elements spherical array depicted in Figure 1. Each supposedly isotropic and identical sensor M i for i P 1, 3 evolves on a sphere of radius r centered at O. OM 1 and OM 3 (resp. OM 2 ) belongs to the cartesian plane Oxy (resp. Ozx), their angular positions being defined by angles ˘β (resp. η). In keeping with Figure 1 and according to the activation scenario, the single source S propagates in free space at radio wavelength λ from angle coordinates pθ, ϕq while the two angular degrees of freedom of the helmetmounted antenna array are pβ, ηq. For notational convenience, pθ, ϕq and pβ, ηq will also be denoted as Θ and Ω, or even omitted in thereafter derivations whenever there is no possible confusions. Consequently, the vector measurement used to perform the DOA estimation is yrts " e Ω pΘqsrts `nrts

Notations

where srts P C is one of the learning sequences of the considered wireless protocol, at this stage assumed to be a stochastic zero mean white Gaussian noise, while nrts P C M ˆ1 is supposed to be a complex, circular, temporally and spatially white Gaussian noise. Respective signal and noise variances are denoted σ 2 s and σ 2 n . The next section is focused on the associated Fisher Information Matrix, subsequently used in the array configuration optimization. A complete study on the 3D FIM derivation can be found in [START_REF] Stoica | The stochastic crb for array processing: a textbook derivation[END_REF] for general array architectures, not necessarily formulated to apprehend the simple optimization method adopted in this paper.

III. FISHER INFORMATION MATRIX AND ASSOCIATED STOCHASTIC CRAMER-RAO BOUND

As a reminder, the Cramer-Rao Bound (CRB) is a lower bound of the covariance matrix of any unbiased estimator, first introduced by R.A. Fisher in [START_REF] Fisher | The logic of inductive inference[END_REF], defined as the inverse of the Fisher Information Matrix. Various expressions and applications of the CRB may be found in the literature, e.g. for MUSIC (MUltiple SIgnal Classification) performances analysis in [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] or [START_REF] Porat | Analysis of the asymptotic relative efficiency of the music algorithm[END_REF].

Defining the Fisher Information Matrix as F pΘq " ´Er B 2 ΛpΘq BΘ 2 s, the Cramer-Rao theorem states that Γp Θq ´F pΘq ´1 is a positive definite matrix, where Γp Θq is the covariance matrix of Θ, unbiased estimator of Θ. Due to the exponential nature of the steering vector elements in e Ω pΘq " re j 2πr λ unpΩ,Θq s, the aforementioned derivative leads to define the steering derivative vector v ξ " δ ξ pUqe Ω pΘq for ξ P tθ, ϕu, where U = diag(u n pΩ, Θq). After calculations detailed in Appendix A, under the assumptions made in Section II and using the simplified notation e for e Ω pΘq, the FIM relevant to our problem takes the following form

F pΘq " γ N ˆxv θ , v θ y Ke xv θ , v ϕ y Ke xv ϕ , v θ y Ke xv ϕ , v ϕ y Ke ˙(3) with γ N " 2αT 1 p1 `1 αN q ˆ2πr λ ˙2 ( 4
)
where T is the number of measured samples, α " σ 2 s {σ 2 n the linear SNR and N the number of array elements. x., .y Ke denotes the scalar product weighted by the orthonormal projector on the MUSIC defined noise subspace Π Ke " I N ´ẽẽ : , i.e.

xv θ , v ϕ y Ke " v : θ Π Ke v ϕ with ẽ " e ?
e : e Note that except the explicit form of γ N in (4), equation ( 3) is valid for any spherical array configuration in the single source DOA estimation scenario.

IV. ARRAY OPTIMISATION METRICS

The inversion of (3) being straightforward, let us consider only the upper left and lower right terms of F pΘq ´1, lower bounds for the variance of respectively any θ and ϕ estimation, designated thereafter by CRB ξ,ξ for ξ P tθ, ϕu. Defining the hereafter two correlation coefficients

ρ 2 " |xv ϕ , v θ y Ke | 2 xv θ , v θ y Ke xv ϕ , v ϕ y Ke and r 2 ξ " |ẽ : v ξ | 2 ||v ξ || 2 (5)
the Cramer-Rao Bound for parameter ξ adopts the simplified form

CRB ξ,ξ " 1 γ N 1 p1 ´ρ2 q 1 ||v ξ || 2 p1 ´r2 ξ q (6)
As shown in Appendix A, one can see that the coupling correlation coefficient ρ 2 is common to both CRB, while the expression of r 2 ξ depends upon the estimated angle. The interest of this formulation is twofold. First hand, it splits each scalar CRB into three distinct angle related components, easing further analysis of the array DOA estimation performances regarding its configuration. Secondly, these two metrics explicitly expose the numerical stability of both CRB ξ,ξ (always positive according to the Schwarz inequality).

The optimization criterion presented in next section consists in minimizing ρ 2 and r 2 ξ in [START_REF] Verhaevert | Uniform spherical distributions for adaptive array applications[END_REF]. It is closer to an iterative adaptation of the sensors positions for considered DOA areas to extract inherent geometric principles than a strictly speaking mathematical optimization. The aim is to highlight the role of ρ 2 and r 2 ξ played in the CRB variations with respect to array configurations rather than to deal with the analytic formulations of these metrics.

V. ARRAY CONFIGURATION ANALYSIS BASED ON CRB METRICS

In this section we first confirm both CRB θ,θ and CRB ϕ,ϕ expressions through MUSIC Monte-Carlo simulations, since this algorithm has been proved to be asymptotically efficient in the single source case [START_REF] Porat | Analysis of the asymptotic relative efficiency of the music algorithm[END_REF]. We next address the problem of array configuration optimization regarding source DOA based on metrics defined in Section IV.

As seen in ( 6), scalar CRB both depend on SNR (through N ), array configuration parameters pΩ, rq and DOA vector Θ. Hereafter simulations thus concern SNR and Ω variations for nearly random values of Θ. We deliberately indeed avoided DOA values that would raise trivial ambiguities due to array symmetry [START_REF] Proukakis | Study of ambiguities of linear arrays[END_REF], not considered in this paper.

In what follows the sphere radius is assumed fixed at a quarter wavelength to ensure spatial sampling conditions for non-exotic values of Ω. Each Mean Square Error (MSE) value is obtained through a 400 point Monte-Carlo simulation of T " 100 samples. Figure 2 shows the evolution of both MUSIC estimation of θ and ϕ with respect to SNR compared to CRB ξ,ξ for DOA Θ " p0, 60q and array configuration Ω " p80, 10q. Accordingly we will assume for next simulations that the MUSIC estimator here adopts its asymptotic behaviour above a threshold SNR of approximately 10dB.

Achieving CRB validation, Figure 3 displays the MUSIC estimation MSE of DOA Θ " p0, 60q compared to CRB ξ,ξ with respect to β and η for a 15dB SNR.

Trivial behaviours already stand out, as the dependency of an azimuthal estimation on the Oxy array aperture (related to β), or the drop of performances in any estimation of ϕ while M 2 comes closer to the horizontal plane (i.e. for large values of η). Next part extends these observations through 3D plots and analysis of ρ 2 and r 2 ξ variations with respect to sensors locations.

An attentive examination of the variations of ρ 2 with respect to array configurations Ω shows that this metric is related to the orthogonality of the considered source DOA regarding the array plane. As seen in Figure 4, the nulls of p1 ´ρ2 q are located on DOAs coplanar to M 1 M 2 M 3 , while it approaches 1 for orthogonal values of Θ.

On the other hand, coefficients r 2 θ and r 2 ϕ are to be considered as cost metrics, binding the considered array configuration to its would have been configuration for a 2D DOA estimation. To be more specific, it is the whole term ||v ξ || 2 p1 ´r2 ξ q that has to be considered to pursue the analogy. These coefficients (for ξ P tθ, ϕu) concern the DOA estimation efficiency once projected on the azimuthal or elevation planes. The analysis is simpler for ξ " θ, since a projection on the azimuthal plane is very much alike one on the Oxy plane (while the elevational plane depends on θ). Thus Figure 5 represents two spherical plots of ||v θ || 2 p1´r 2 θ q for Ω " p90, 20q and p10, 20q, showing that as β drops the main lobe progressively shifts form θ " 0 to π 2 . The shift occurs when the projection of tM 1 , M 2 , M 3 u in the azimuthal plane is an equilateral triangle, i.e. when M 1 M 2 ¨Ox " ||M 1 M 2 ||cosp π 6 q, which would correspond to the isotropic condition for a 2D estimation on θ. Observing each colour bar, one may also observe the awaited drop of ||v θ || 2 p1 ´r2 θ q maxima while the array characteristic length goes from half to a quarter wavelength.

The same analysis eventually holds for ξ " ϕ. Since the projection of tM 1 , M 2 , M 3 u on the elevational plane is less intuitive, Figure 6 mimics Figure 5 for the term ||v ϕ || 2 p1 ŕ2 ϕ q. As for θ, the general orientation of the antenna array toward Oyz or Ozx defines the main lobe position, partially conditioning DOA estimation acuteness regarding ϕ.

These results are in accordance with those of Gazzah [START_REF] Gazzah | Cramer-rao bounds for antenna array design[END_REF] and Dogandzic [START_REF] Dogandzic | Cramer-rao bounds for estimating range, velocity, and direction with an active array[END_REF] , who observed that least and most desirable DOAs are perpendicular to each other and correspond to the antenna array "main" axes of symmetry (termed principal axes of inertia in [START_REF] Dogandzic | Cramer-rao bounds for estimating range, velocity, and direction with an active array[END_REF]).

Finally, Figures 7 and8 present both products ||v ξ || 2 p1 ŕ2 ξ qp1 ´ρ2 q, i.e. the terms of interest for each considered CRB ξ,ξ in ( 6) regarding angular dependency of the DOA estimation acuteness relative to sensors positions. Reflecting the CRBs behaviours for various array configurations, these DOA is fixed at Θ " p0, 60q for a 15dB SNR.

representations are to be afterwards exploited in the real time design optimization process.

Based on these results, one may iteratively adapt the antenna array geometry to various research scenarios, assuming rough DOA and distance measurements feedback. From afar, the victim's elevation with respect to the man held antenna array approaches zero, while it may obviously be positioned anywhere in the horizontal plane. Preferred configurations thus would be those granting a low and constant CRB θ,θ for ϕ " 0 and θ P r0, 2πs. While the array geometry constraints and its reduced number of elements prevent isotropic estimations on θ, the above-mentioned equilateral condition would here represent the most favourable configuration, giving a pattern close to the one depicted in Figure 7

VI. CONCLUSION

In this paper, we adapted the stochastic Cramer-Rao Bound analytic expression to a helmet-mounted antenna array configuration for the single source DOA estimation case. We decomposed the obtained CRBs into metrics to give physical insight on the DOA estimation acuteness with respect to array configurations. Array elements positioning being reduced to a 2 Degree-of-Freedom problem, a few spherical plots along with proposed metrics alleviate the FIM based design optimization. Observing simple geometrical rules relative to these metrics, one may control each CRB ξ,ξ main lobe position while avoiding critical precision scenarios for a given set of DOAs. Further analysis based on 2D parametric representations of ρ 2 against r 2 ξ with respect to Θ for a given Ω might facilitate zones of interests (i.e. areas where both ρ 2 and r 2 ξ approach zero) characterization. Along with network logistic and human-machine interface considerations not discussed in this paper, the presented helmet-mounted antenna array and its optimization tool are supposed to enhance rescuers efficiency to reach natural disaster survivors. General expressions herein proposed hold for any spherical array configuration in the single source case, and would gain in being extended to any array geometry to pursue metrics behaviour analysis.

VII. ACKNOWLEDGMENTS

The work described in this paper was carried out as part of an ANRT-CIFRE industrial thesis funded by EADS Astrium Services.

APPENDIX

This appendix holds an exhaustive demonstration of the FIM formulation of equation (3) used to establish CRB ξ,ξ scalar expressions in [START_REF] Verhaevert | Uniform spherical distributions for adaptive array applications[END_REF].

According to hypothesis stated in Section II, the probability density functions (pdf) of received snapshots yrts along the T measurements are independents. Their joint pdf thus takes the form ϕ q on a sphere spanning pθ, ϕq for (a) Ω " p66, 12q and (b) p9, 12q. The main lobe inversion occurs when the elevational projections of tM 1 , M 2 , M 3 u form an equilateral triangle.

T ź t"1 ppyrts|Θ, σ 2 s , σ 2 n q " 1 π N T 1 |Γ y | T expp´T ÿ t"1 y : rtsΓ ´1 y yrtsq (7)
Where the assumptions made on srts and nrts imply (using the steering vector simplified notation e " e Ω pΘq)

Γ y " σ 2 s ee : `σ2 n I N (8)
Γ y being the spatial covariance matrix of vector y. The exponential term in [START_REF] Ferréol | High resolution direction finding: From performance toward antenna array optimization-the mono-source case[END_REF] being a scalar, we may use the circular The basic formulation of the log-likelihood function is then obtained by replacing ( 9) in [START_REF] Ferréol | High resolution direction finding: From performance toward antenna array optimization-the mono-source case[END_REF] and taking its natural logarithm

Λ " ´N T logpπq ´T log |Γ y | ´T trpΓ ´1 y Γy q (10)
Λ thus needs to be shaped for derivation regarding Θ. Using [START_REF] Gazzah | Cramer-rao bounds for antenna array design[END_REF]

Since the considered derivation concerns Θ dependent elements, only the last term of ( 11) is of interest. For the sake of clarity, derivatives like δ ξ peq or δ 2 ξ,ζ pe : q will thereafter be denoted as δ ξ e and δ 2 ξ,ζ e : , the derivative affecting its adjacent term only. The first derivative of [START_REF] Porat | Analysis of the asymptotic relative efficiency of the music algorithm[END_REF] with respect to ξ P tθ, ϕu is then ξ pI N ´ẽẽ : qv ξ " ||v ξ ||p1 ´r2 ξ q for ξ P tθ, ϕu.

δ ξ Λ " T σ 2

3D

  Head-Mounted Antenna Array Architecture Optimization based on the Fisher Information Matrix D. Poirier-Quinot :πε , P. Duvaut ε , L. Girardeau : , B.F.G. Katz π : EADS Astrium Services, 92150 Suresnes, France ε ETIS / ENSEA -Université Cergy-Pontoise -CNRS, 95000 France π LIMSI-CNRS, Université Paris Sud, 91403 Orsay, France

Fig. 1 .

 1 Fig. 1. Helmet-mounted antenna for single source S location

Fig. 2 .

 2 Fig. 2.Variations of CRB θ,θ , CRBϕ,ϕ and MSE of MUSIC DOA estimation of θ and ϕ with respect to SNR for fixed array configuration Ω " p80, 10q and DOA Θ " p0, 60q.

Fig. 3 .

 3 Fig. 3. Variations of CRB θ,θ , CRBϕ,ϕ and MSE of MUSIC estimation of θ and ϕ regarding (a) β and (b) η for fixed η " 0 ˝and β " 60 ˝respectively.DOA is fixed at Θ " p0, 60q for a 15dB SNR.

  Fig. 4. Variations of ρ 2 along θ and ϕ for Ω " p80, 45q. The least desirable DOAs are for ρ 2 =1, i.e for a source in the antenna array plane M 1 M 2 M 3 .

Fig. 5 .

 5 Fig. 5. Variations of ||v θ || 2 p1 ´r2θ q on a sphere spanning pθ, ϕq for (a) Ω " p66, 12q and (b) p9, 12q. The main lobe inversion occurs when the azimuthal projections of tM 1 , M 2 , M 3 u (i.e. on Oxy) form an equilateral triangle.

Fig. 6 .

 6 Fig. 6. Variations of ||vϕ|| 2 p1 ´r2ϕ q on a sphere spanning pθ, ϕq for (a) Ω " p66, 12q and (b) p9, 12q. The main lobe inversion occurs when the elevational projections of tM 1 , M 2 , M 3 u form an equilateral triangle.

Fig. 7 ." 1 y 1 yrtsyrtsFig. 8 .

 7118 Fig. 7. Variations of p1 ´ρ2 q||v θ || 2 p1 ´r2 θ q for (a) Ω " p87, 51q and (b) p45, 0q. The general array orientation favors θ estimations on larger areas for DOAs orthogonal to the array plane.

  of (3) reveals the determinant p1 ´ρ2 q " xv θ , v θ y Ke .xv ϕ , v ϕ y Ke ´xv θ , v ϕ y Ke .xv ϕ , v θ y Ke and leads to the CRB matrix, which holds both CRB ξ,ξ diagonal terms introduced in (6)F pΘq ´1 " 1 γ N 1 p1 ´ρ2 q ˆxv θ , v θ y ´1 Ke ´ρ2 xv θ , v ϕ y ´1 Ke ´ρ2 xv ϕ , v θ y ´1 Ke xv ϕ , v ϕ y ´1 Ke ẇhere xv ξ , v ξ y Ke " v :

  to replace |Γ y | |Γ

	trpΓ ´1 y Γy q "	1 σ 2 n	trprαee : `IN s ´1 Γy q
	"	1 σ 2 n	trprI N	´αee : 1 `αN	s Γy q
	in (10), the log-likelihood becomes
	Λ " ´N T plogpπq `logpσ 2 n qq ´T σ 2 n	trp Γy q
	´T logp1 `αN q `T σ 2 n	α p1 `αN q		e :	Γy e
						n " σ 2N n pαN `1q pσ 2 s ee : `σ2 n q and Γ ´1 y (according to Woodbury inversion formula)

y | " |σ 2 s ee : `σ2 n I N | " σ 2pN ´1q

  ℜtδ ξ e : δ ζ eu `αℜtδ ξ e : ee : δ ζ eus we get the first and second order derivatives δ ξ e : " ´j 2πr λ e : δ ξ U δ 2 ξ,ζ e : " ´j 2πr λ e : δ 2 ξ,ζ U ´p 2πr λ q 2 e : δ ξ Uδ ζ U used to develop (13) into ´Erδ 2 ξ,ζ Λs " γ N N pe : ee : δ ζ Uδ ξ Ue ´e: δ ζ Uee : δ ξ Ueq (14) Remark here that Erδ 2 ξ,ζ Λs " Erδ 2 ζ,ξ Λs, since U is diagonal with real coefficients and ee : is a hermitian matrix. Noticing the weighted scalar product in (14) pe : ee : δ ζ Uδ ξ Ue ´e: δ ζ Uee : δ ξ Ueq " e : epv : ζ Π Ke v ξ q " N xv ξ , v ζ y Ke where Π Ke is the orthonormal projector on the noise space defined in the MUSIC subspace decomposition Π Ke " I N ´ẽẽ : with ẽ " Λs " γ N xv ξ , v ζ y Ke (15) which proves the proposed FIM in (3): F pΘq " γ N ˆxv θ , v θ y Ke xv θ , v ϕ y Ke xv ϕ , v θ y Ke xv ϕ , v ϕ y Ke

										with
										γ N " 2αT	1 p1 `1 αN q	p	2πr λ	q 2
				n	α p1 `αN q	pδ ξ e :	Γy e `e:	Γy δ ξ eq
	the second derivative related to ζ P tθ, ϕu naturally reveals cross conjugated terms, hereafter gathered in ℜt¨u
	δ 2 ξ,ζ Λ " 2	T σ 2 n	α p1 `αN q	pℜtδ 2 ξ,ζ e :	Γy eu `ℜtδ ξ e :	Γy δ ζ euq (12)	equation (14) becomes	?	e e : e
	According to the FIM definition stated in Section III, we
	need to evaluate the negative expectation of the second order derivative log-likelihood. Replacing the estimate Γy by its true value Γ y in (12) gives Er´δ 2 ξ,ζ Λs " ´2 T σ 2 n α p1 `αN q pℜtδ 2 ξ,ζ e Er´δ 2 ξ,ζ Λs " ´2αT p1 `αN q rp1 `αN qℜtδ 2 `(13) ξ,ζ e : eu	´Erδ 2 ξ,ζ ȯr its equivalent vectorial form that highlights the MUSIC noise space projection F pΘq " γ N ˆv:
	We then develop the differential terms δ ξ e : and δ 2 ξ,ζ e : ,
	introducing the notations
	e "	¨ej 2πr λ u1 e j 2πr λ u2 . . .	‹ ‹ ‹ ‹ '	and U "	¨u1 0 . . . 0 0 u 2 . . . 0 . . . . . . . . . . . .	‹ ‹ ‹ '
		e j 2πr λ u N					0	0 . . . u N

:

Γ y eu`ℜtδ ξ e : Γ y δ ζ euq which once developed according to

[START_REF] Gazzah | Cramer-rao bounds for antenna array design[END_REF]

gives


Nathan Souviraà-Labastie

Emmanuel Vincent

Frédéric Bimbot

MUSIC SEPARATION GUIDED BY COVER TRACKS: DESIGNING THE JOINT NMF MODEL

Keywords: Music separation, Joint-NMF, Cover song

In audio source separation, reference guided approaches are a class of methods that use reference signals to guide the separation. In prior work, we proposed a general framework to model the deformation between the sources and the references. In this paper, we investigate a specific scenario within this framework: music separation guided by the multitrack recording of a cover interpretation of the song to be processed. We report a series of experiments highlighting the relevance of joint Non-negative Matrix Factorization (NMF), dictionary transformation, and specific transformation models for different types of sources. A signal-to-distortion ratio improvement (SDRI) of almost 11 decibels (dB) is achieved, improving by 2 dB compared to previous study on the same data set. These observations contribute to validate the relevance of the theoretical general framework and can be useful in practice for designing models for other reference guided source separation problems.

INTRODUCTION

In signal processing, audio source separation is the task of recovering the different audio sources that compose an observed audio mixture. In the case of music, this task aims to provide signals for each instrument or voice. As original songs are rarely released in multitrack formats, this step is compulsory to open new possibilities in music post-production, e.g., respatialization, upmixing and more widely in audio edition.

The complexity of the mixing process (not necessarily linear) as well as the fact that there are more sources than input channels make the demixing of professionally produced music difficult in the blind case, i.e., without any prior information. Thus, blind separation shows certain limitations for professional music applications that require high audio quality [START_REF] Vincent | From blind to guided audio source separation: How models and side information can improve the separation of sound[END_REF]. Many approaches have taken several kinds of additional information into account with the objective of overcoming these limitations [START_REF] Liutkus | An overview of informed audio source sep-aration[END_REF]. For instance, spatial and spectral information about the sources [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF], information about the recording/mixing conditions [START_REF] Ngoc | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], musical scores [START_REF] Umut Simsekli | Score guided audio restoration via generalised coupled tensor factorisation[END_REF][START_REF] Ewert | Score-informed source separation for musical audio recordings: An overview[END_REF], or even selection of spectrogram areas [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for user-guided audio source separation[END_REF][START_REF] Durrieu | Musical audio source separation based on user-selected F0 track[END_REF][START_REF] Fitzgerald | User assisted separation using tensor factorisations[END_REF], potentially in an interactive way [START_REF] Fuentes | Blind harmonic adaptive decomposition applied to supervised source separation[END_REF][START_REF] Bryan | Interactive refinement of supervised and semi-supervised sound source separation estimates[END_REF][START_REF] Ngoc | An interactive audio source separation framework based on non-negative matrix factorization[END_REF] have been proposed in the literature. It is also possible to consider reference signals [START_REF] Souviraà-Labastie | Multi-channel audio source separation using multiple deformed references[END_REF] that are similar to the sources to be separated, for instance uttered by a user [START_REF] Smaragdis | Separation by humming : sound extraction from monophonic mixtures[END_REF][START_REF] Le Magoarou | Text-informed audio source separation using nonnegative matrix partial co-factorization[END_REF][START_REF] Hennequin | Speech-guided source separation using a pitch-adaptive guide signal model[END_REF], or synthesized from symbolic information [START_REF] Le Magoarou | Text-informed audio source separation using nonnegative matrix partial co-factorization[END_REF][START_REF] Fritsch | Score informed audio source separation using constrained nonnegative matrix factorization and score synthesis[END_REF] or even retrieved from a large set [START_REF] Souviraà-Labastie | Audio source separation using multiple deformed references[END_REF][START_REF] Badawy | Onthe-fly audio source separation[END_REF].

In this paper, we focus on source separation guided by reference signals, and more precisely music separation guided by cover multi- * Work supported by Maia Studio and Bretagne Region scholarship track songs [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF]. A cover song is another performance of an original song. It can differ from the original song by its musical interpretation and it is potentially performed by a different singer and with different instruments. Multitrack recordings of such covers are more likely to be found on the market than the original multitrack recording and contrary to expectations, they are (for commercial reasons) musically faithful to the original [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF]. Furthermore, most separation algorithms are sensitive to initialization and using cover multitrack recordings for initialization is an efficient way to sidestep this problem [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF]. All these reasons make cover guided music separation a very promising approach for high quality music separation.

In the following, rather than using the cover multitrack recordings for initialization only, we also use them to constrain the power spectrum of each source. In addition, although the considered covers are musically faithful, deformations exist between the original sources and the covers at the signal level. These deformations are significant enough not to be ignored. Here, different configurations of deformations as formalized in [START_REF] Souviraà-Labastie | Multi-channel audio source separation using multiple deformed references[END_REF] are tested. Finally, the optimal deformation model is selected for each type of source (bass, drums, guitar, vocal ...).

The paper is organized as follows. Section 2 recalls the general model of reference guided source separation proposed in [START_REF] Souviraà-Labastie | Multi-channel audio source separation using multiple deformed references[END_REF]. Section 3 presents a specific model adapted to the task, and its estimation procedure. Section 4 describes the data and the settings of the experiments. Section 5 reports on results illustrating the proposed contributions and provides to the reader useful advice on how to design deformation models for reference guided source separation. Section 6 draws the conclusions and gives some perspectives.

GENERAL FRAMEWORK

In this section, we recall the general framework for multi-channel source separation guided by multiple deformed references [START_REF] Souviraà-Labastie | Multi-channel audio source separation using multiple deformed references[END_REF].

The observations are M multi-channel audio mixtures x m (t) indexed by m and containing I m channels. For instance, one mixture is to be separated, and the other mixtures contain the reference signals used to guide the separation process. Each mixture x m (t) is assumed to be a sum of source spatial images yj(t) indexed by j ∈ Jm. In the Short-Time Fourier Transform (STFT) domain, this can be written as

x m f n = j∈Jm y j,f n with x m f n , y j,f n ∈ C I m , (1)
where f = 1, ..., F and n = 1, ..., N are respectively the frequency and the time indexes of the STFT. We assume that the timefrequency coefficients of the source spatial images y j,f n have a zero-mean Gaussian distribution [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF] :

y j,f n ∼ N C (0, v j,f n R j,f n ) (2)
whose covariance factors into a scalar power spectrum v j,f n ∈ R+ and a spatial covariance matrix R j,f n ∈ C I m ×I m . The spatial covariance matrices model the spatial characteristics of the sources, like phase or intensity difference between channels. The power spectrogram of each source j is denoted as

Vj = [v j,f n ] f n ∈ R F ×N + .
Each Vj is split into the product of an excitation spectrogram V e j and a filter spectrogram V φ j . The excitation spectrogram (resp. the filter spectrogram) is decomposed by a Non-negative Matrix Factorization (NMF) into a matrix of spectral patterns W e j ∈ R F ×D e + (resp.

W φ j ∈ R F ×D φ +
) and a matrix of temporal activations

H e j ∈ R D e ×N + (resp. H φ j ∈ R D φ ×N +
). D e (resp. D φ ) denotes the number of spectral patterns used in the NMF decomposition of the excitation (resp. filter) part. This results in the following decomposition :

Vj = V e j V φ j = W e j H e j W φ j H φ j (3)
where denotes point wise multiplication.

As the different audio mixtures are composed of similar sources, the matrices W and H can be shared (i.e., jointly estimated) between a given source j ∈ J m and one or more related sources j ∈ J m with m = m. A reference signal is by nature deformed : if there were no deformation, the reference would be equal to the true source. These differences can be taken into account by partly sharing the model parameters (e.g., sharing only H e ) and/or by using transformation matrices T jj (between sources j and j ). For instance, the sharing of excitation properties between two sources j and j is modeled by one of the three following configurations, depending on the sharing of W e and H e :

V e j = T f e jj W e j H e j (4)
V e j = W e j H e j T te jj (5)

V e j = T f e jj W e j T de jj H e j T te jj [START_REF] Ewert | Score-informed source separation for musical audio recordings: An overview[END_REF] As it would be redundant, T de jj only appears if the corresponding W e j , H e j are shared.

MODEL AND ALGORITHM

In this section, we introduce the specific framework used to model the cover guided music separation problem in reference to the general framework recalled in Section 2 from which we removed the excitation-filter decomposition and the multi-channel formulation.

Here, the reference mixtures are the different cover tracks and contain a single source. This configuration with isolated references leads to a very efficient initialization (detailed in Section 3.2)

Proposed model

Here, we consider that x 1 (t) is the song to be separated, and x m (t) for m > 1 are the different tracks of the cover version used to guide the separation process. All mixtures are single channel, and the mixture x m (t) for m > 1 are assumed to contain only one source. Each Vj is decomposed by a NMF into a matrix of spectral patterns Wj ∈ R F ×D + and a matrix of temporal activations Hj ∈

R D×N + . Vj = WjHj (7)
D denotes the number of spectral patterns used in the NMF decomposition. Hereafter, we only consider frequency and dictionary transformation matrices that are now denoted T f jj ∈ R F ×F + , and

T d jj ∈ R D×D + .
As each track of the two versions are sufficiently aligned in time, we do not consider T t matrices that induce unwanted smoothing effects. Thus, the related sources are modeled using equation ( 6):

V j = T f jj WjT d jj Hj (8)
It can be noticed that this formulation leaves the possibility to set these transformation matrices either in the reference model (j ∈ J 1 and j ∈ J m , m = 1) or in the source model (j ∈ J 1 and j ∈ J m , m = 1). See Tables 3 and4 in Section 5 for concrete cases. For both T f and T d matrices, we consider two possible initializations:

• Diag : Identity matrix

• Full : The sum of an identity matrix and a random matrix drawn from a rectified Gaussian distribution

Parameter estimation

Here, we present a method for parameter estimation in the maximum likelihood (ML) sense. In the single-channel case, maximizing the log-likelihood is equivalent to minimizing the Itakura-Saito divergence [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF]:

θ = argmin θ M m=1 F,N f,n=1 dIS(X m f n |V m f n ) (9)
where θ is the set of parameters to be estimated, i.e., the matrices W , H and T .

X m = [|x m f n | 2 ]
f n and V m = j∈J m Vj are respectively the observed and estimated power spectrograms, and dIS(a|b) = a/b -log(a/b) -1 is the Itakura-Saito divergence. A common way to estimate the parameters is the use of a multiplicative gradient descent approach [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF] in which each parameter is updated at each iteration without increasing criterion (9) [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF]. The update of one parameter consists in multiplying it by the ratio of the negative and positive parts of the derivative of the criterion with respect to this parameter. Different multiplicative updates (MU) are derived for each parameter. Examples of such updates are given in [START_REF] Souviraà-Labastie | Audio source separation using multiple deformed references[END_REF].

The results of MU depend on initialization. With respect to blind source separation, reference guided separation provides better initial parameter values taking advantage of the provided references. For instance, in our case, we can use MU to minimize the following criterion :

θref = argmin θ ref M m=2 F,N f,n=1 dIS(X m f n | Ṽ m f n ) (10)
where Ṽ m f n = WjHj with j ∈ J m (m > 1), and θref is the set of corresponding Wj and Hj parameters. This is especially efficient when there is a single source in the reference signal, as here. In the experiments, we will distinguish the following successive initialization and algorithmic stages :

• Pre-NMF : MU that try to minimize the criterion (10),

• Joint-NMF : MU that try to minimize the criterion [START_REF] Fitzgerald | User assisted separation using tensor factorisations[END_REF].

In all experiments, 50 iterations of Pre-NMF are beforehand performed and the final source estimates are obtained using an adaptive Wiener filter.

EXPERIMENTS

Data

In order to compare our results with previous study, we use the same data set and equivalent settings as [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF]. Both original and cover multitracks are available in order to evaluate the separation. The number of cover tracks is the same as the number of sources to be separated, as each track is used as reference for the related source. They are also used in the mirror configuration, i.e., the cover is to be separated and the original multitracks are used as references. Experiments are conducted on 30 second examples typically consisting of half of a verse and half of a chorus. The considered tracks of four songs are listed in Table 1. Some examples are available online 1 .

Settings

We here make an exhaustive list of settings that vary from [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF] and refer the reader to [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF] for common details. To simplify the analysis of the results, the two channels are summed into a single one. We used 50 iterations for Pre-NMF and Joint-NMF steps instead of 500. The number of components D is kept to 50.

RESULTS

The quality of the estimated sources is evaluated in terms of signalto-distortion ratio improvement (SDRI) that is the difference between the output SDR [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF] and the input SDR. The input SDR is defined as the power ratio between the source to be estimated and the mixture to be separated. It is given for each type of source in Table 4. The samples that we selected lead to an input SDR mean of -8.44 dB instead of -7.60 dB in [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF].

As we work with MU, we underline that zeros in the parameters remain unchanged over the iterations. Hence, a matrix initialized at identity will stay diagonal over the iterations. Moreover, T matrices are not present during the Pre-NMF.

Comparison with a previous study

In [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF], the cover multitrack signals are only used to initialize the source parameters W and H. Here, we decided to share these pa- rameters between the source and the reference hence the reference signals are used during the global estimation stage too. The results are shown in Table 2. First, we reproduce the experiments in [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF] with the differences of settings presented in Section 4.2. An equivalent SDRI mean is obtained compared to [START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF] for the case where the parameters are not shared (8.74 dB instead of 8.98 dB). The previous configuration leads in fact to an important decrease of the SDRI mean, compared to what is obtained if the sources are reconstructed directly after the Pre-NMF (10.06 dB). This can be explained by the great level of similarity between the covers and the original tracks. Conversely, sharing the NMF parameters during the final estimation (Joint-NMF) guarantees not to move away too much from this relevant starting point while getting closer to a solution that fits better the original tracks. In our case, a marginal improvement is observed (10.27 dB).

These first results show the strong similarity between each original track and its related cover track. In that case, sharing W and H during the Joint-NMF estimation is the most relevant method even without considering any deformations.

Designing the deformation model

Here, we analyze whether the transformation matrices are useful in the reference or the source model. The comparison of different initializations of frequency and dictionary deformations matrices is done as well. Results of exhaustive experiments are displayed in Table 3.

Bold values indicate improvements compared to a fully shared Joint-NMF (10.27 dB). It can be noticed that in those two cases the number of non zero coefficients Z in T is of the same order of magnitude (Z = D 2 = 2500 for a Full T d and Z = F = 1025 for a Diag T f ), while Z vary from D = 50 to F 2 ≈ 10 6 .

We also observe that for almost all cases, SDRIs are always higher when T matrices are positioned in the reference model. This can be explained by the fact that the final signals are reconstructed based on the source model, and T matrices can induce abrupt changes. Conversely, inserting a Full T matrix in the reference model instead of in the source model would distort more the output of the Pre-NMF, as the product of Wj and Hj estimated during the Pre-NMF step try to fit the reference signal. So, it is difficult to distinguish which effect is predominant, especially since the number of non-zero coefficients has also an impact.

One can remember that when inserting transformation matrices it is important not to distort the output of Pre-NMF. For instance in

Specific source model

In this last set of experiments, SDRI for each type of source are given in Table 4 for different configurations. The combination of T d and T f gives interesting results, especially for the drums and the bass. Moreover, we observe for each type of source clear differences between models while the SDRI means are similar. In a final experiment, referred to as "Best", the most proper configurations are chosen for each type of source (indicated in bold in Table 4). Note that values in italic and bold were chosen based on the experiment with T f (Diag) in the source model that shows promising results for these two sources. The result (10.92 dB) shows the relevance of defining specific models depending on the source type. It should also be noted that using transformation matrices for only one source of a song leads to a huge decrease of separation quality for that source. This is certainly due to the estimation algorithm that we used. So, defining specific source models is indeed interesting but the number of parameters should be balanced between the different models.

CONCLUSION

In this paper, we have addressed the problem of cover tracks guided music separation with a previously defined general framework for audio source separation using multiple deformed references. This study leads us to consider specific aspects of the general framework, in particular dictionary deformations. Several deformation models were tested for different types of source (Bass, Drums, Guitar, Vocal ...) which lead to 2 dB improvement compared to a previous study on the same data set. The results obtained in this article may be useful in other scenarios of reference guided audio source separation.

As this paper entirely focused on spectral modeling, considering spatial information would be an asset, for instance using [START_REF] Souviraà-Labastie | Multi-channel audio source separation using multiple deformed references[END_REF]. However, it will be needed to compensate the potential differences of spatialization between each source and its related cover track.

Table 1 .

 1 Cover multitrack data set

	Title		Track names
	I Will Survive	Bass, Brass, Drums, Guitar, Strings, Vocal.
	Pride and Joy		Bass, Drums, Guitar, Vocal.
	Rocket Man	Bass, Choirs, Drums, Others, Piano, Vocal.
	Walk this Way		Bass, Drums, Guitar, Vocal.
			Joint-NMF	SDRI
			Source Reference mean
	Result in [20]		W H		8.98
	Reproduction of [20]	W H		8.74
	Only Pre-NMF			10.06
	Joint-NMF		W H	W H	10.27

Table 2 .

 2 SDRI (dB) compared with a previous study[START_REF] Gerber | Professionally-produced music separation guided by covers[END_REF].

Table 3 .

 3 1 http://speech-demos.gforge.inria.fr/source separation/icassp2015/ SDRI (dB) for different configurations.

	Init		Joint-NMF	SDRI
	T f	T d	Source	Reference mean
			W H	W H	10.27
	Full		W H	T f W H	10.08
	Full		T f W H	W H	9.23
	Diag		T f W H	W H	10.09
	Diag		W H	T f W H	10.35
	Diag W T d H	W H	9.25
	Diag	W H	W T d H	9.88
	Full	W H	W T d H	9.79
	Full W T d H	W H	10.64

Table 3 ,

 3 it is the case for values in bold.

		Number of tracks	4	4	3	4	1	1	1	1	1
		Input SDR (dB)	-7.42	-7.17	-9.98	-4.18	-12.34	-9.75	-12.48 -18.64	-10.55
	Joint-NMF Source Reference mean SDRI Bass Drums Guitar Vocal Choirs Others Piano Brass Strings
			10.06	9.33	9.02	9.71	9.60	13.70	9.79	10.80	16.25	9.80
	W H	W H	10.27	9.26	9.28	9.82	10.24	13.23	10.27	11.11	15.67	10.62
	W H	T f W H	10.35	9.39	9.29	9.69	10.22	13.85	9.98	12.08	15.68	10.71
	W T d H	W H	10.64	9.23	9.94	10.80	10.73	13.01	10.07	11.36	15.80	10.61
	W T d H	T f W H	10.66	9.88	10.48	9.41	10.44	12.74	10.01	12.24	16.66	10.11
		Best	10.92	9.66	10.58	10.35	10.85	13.73	10.74	11.91	14.91	11.71

Table 4 .

 4 SDRI (dB) for each type of source.


J Ér Émie Garcia

Jean Bresson
email: bresson@ircam.fr

Thibaut Carpentier
email: thibaut.carpentier@ircam.fr

Towards Interactive Authoring Tools for Composing Spatialization

Keywords: H, 5, 5 [Information Systems]: Information Interf, and Presentation-Sound and Music Computing;

We present interactive tools designed to help music composers controlling sound spatialization processes in a computer-aided composition environment. We conducted interviews with composers to understand their needs and inform the design of new compositional interfaces. These interfaces support quick input, visualization and edition of three-dimensional trajectories, as well as the control of the temporal dimension in spatial scene descriptions.

INTRODUCTION

Contemporary music composers extensively use sound spatialization techniques as a means to increase the formal or acoustic dimensionality of music [START_REF] Baalman | Spatial composition techniques and sound spatialisation technologies[END_REF]. These techniques create the illusion that sound sources come from various directions in space and organize them in three-dimensional scenes.

Despite the development of advanced spatial audio renderers [START_REF] Jot | Real-time Spatial Processing of Sounds for Music, Multimedia and Interactive Human-Computer Interfaces[END_REF], the control of sound spatialization remains a challenging task for composers [START_REF] Peters | Current technologies and compositional practices for spatialization: A qualitative and quantitative analysis[END_REF]. Spatialization processes require users to define the location of sound sources, as well as a great number of other spatial features such as their orientation, the amount of doppler effect or the reverberating characteristics of the virtual space. As these parameters often vary to produce sound motions or other dynamic effects, composers also need to specify their temporal evolutions. In addition, the compositional focus implies this time determination to be carefully related with the score parts.

Many composers use visual programming environments such as OpenMusic [START_REF] Assayag | Computer Assisted Composition at IRCAM: From PatchWork to OpenMusic[END_REF] or Max to create and run programs that combine algorithmic processing with musical material or notation. Previous works explored tools and formats to generate, represent and manipulate spatial parameters within this kind of environment [START_REF] Bresson | Spatial Structures Programming for Music[END_REF]. However they mostly focused on the integration of complex multidimensional data in compositional processes rather than providing interactive tools to support end user interactions.

Spatialization control requires interactions in a threedimensional environment that represents the acoustic world. It presents interesting challenges for 3D user interfaces because it also involves a temporal dimension, sound perception factors, and the need to provide live aural feedback [START_REF] Thibaut | A design space for three-dimensional curve edition[END_REF].

Several researchers already explored interaction techniques to select, move and interact with objects in space [START_REF] Jankowski | A Survey of Interaction Techniques for Interactive 3D Environments[END_REF]. We focus on a different task: the musical composition of temporal spatial scenes. Our goal is to design interactive systems to support composers' work with spatialization. We aim at providing tools for defining sound sources trajectories in potentially complex scenes that can be connected to temporal structures and to other aspects of compositional processes.

OBSERVATIONS

We conducted interviews with four professional composers experienced in sound spatialization in order to better understand their needs and reflect on new interactive tools to facilitate the monitoring and writing of spatialization processes. We studied the representations and the software they use during their composition process and we focused on the problems they encounter.

We observed that composers primarily work with visual representations of the sound space, for sketches as well as for the final scores. Each composer follows a unique process, combining personal graphical representations with various musical concepts. They explained that they need such representations to explore ideas and articulate them with other musical elements usually described in the score.

All four composers mentioned the temporal specification of spatial processes as a major problem. They also emphasized their need for a global vision of the piece and a more detailed one, generally multidimensional (many sources and many parameters), in a unified context. They were all interested in potential interactive ways to capture, visualize and render trajectories as well as other spatial and temporal parameters. Two composers discussed the need to quickly input and hear trajectories to explore their ideas before developing more complex processes.

PROTOTYPES

Based on these observations, we decided to improve the tools available in the OpenMusic computer-aided composition environment for specifying 3D-trajectories [START_REF] Bresson | Spatial Structures Programming for Music[END_REF] and describing temporal spatial scenes with multiple sound sources.

Figure 1 shows an OpenMusic visual program containing our prototype objects 1 . OpenMusic "objects" are data structures built and used in visual programs. Composers can create them programmatically from lists of parameters or using graphical editors. Objects and functions connected together in visual programs constitute a data-flow graph corresponding to a compositional process.

Editing / Visualizing Trajectories

The 3DC object allows to view a 3D curve and edit it with two synchronized 2D editors (top/front view). Figure 2 shows the graphical user interface of a 3DC object with the 3D and 2D views, as well as the parameters panel providing a number of visualization options. For instance, the user can map a color scale to the length or speed of the trajectory. These features are meant to help composers to easily understand the motion of a trajectory and to quickly recognize fast trajectories from slow ones.

Gestural 3D Input

To support gestural input from external controllers, the 3DC interface embeds an OSC2 server receiving position messages from various input devices. The protocol consists in two main messages: the move message (move a sound source in the 3D space) and the draw message (append points to the trajectory). A spatial rendering engine provides real-time audio feedback and the graphical user interface displays a green or a red sphere in the 3D view when the user respectively moves or "draws" in the trajectory (Figure 3.a).

We experimented this system with a GameTrak controller. This controller provides two 3D positions localized from the length and angle of two wires articulated on a base. The wires are tied to two fingers of one hand. Figures 3.b and 3.c present the two hand postures used to switch between move and draw modes. The posture detection is made using a distance threshold between the points. This simple interaction technique allows the user to explore sound source movements with audio-visual feedback and record the gesture by joining fingers and performing mid-air hand movements.

Temporal Spatial Scenes

SPAT-SCENE is an other OpenMusic object used to create, view and edit spatial sound scenes and their evolution in time. This object is built from a set of sound sources and 3DC objects. Figure 1 contains a SPAT-SCENE object at the bottom-left. Figure 4 shows its graphical editor. The SPAT-SCENE editor is divided in two main areas: the sound scene (top) and the timeline (bottom). The sound scene view controls the sound sources and the speakers positions at a given moment. The timeline contains, for each source, the times of all source position changes. These time positions are equivalent to keyframes in video editors. Users can determine the rendering interpolation rules for the source positions between successive keyframes.

CONCLUSION AND PERSPECTIVES

We presented interactive tools for authoring sound spatialization processes in the OpenMusic computer-aided composition environment. Interviews with four professional composers helped us understanding their needs for input and visualization of complex temporal scenes with different levels of details. Our prototypes address these needs for specific scenarios which can be generalized and adapted to other devices or environments.

The tools can be integrated into composers' workflow to control the production of spatialized audio while exploring alternatives with real-time audio-visual feedback. We believe this approach will enhance composers' creative processes by fostering the expression and evaluation of musical ideas. Our future plan is to collaborate with composers to refine and assess these tools in real music production contexts.

Figure 1 :

 1 Figure 1: Construction of a spatial scene object in an OpenMusic visual program. Each 3DC object contains a trajectory. The SPAT-SCENE displays the initial positions of the three trajectories.

Figure 2 :

 2 Figure 2: Graphical editor for a 3D trajectory object (3DC).

Figure 3 :

 3 Figure 3: Gestural input of a new trajectory using the GameTrak controller: a) 3D view with visual feedback; b) free move hand posture; c) draw mode hand posture.

Figure 4 :

 4 Figure 4: Graphical editor for the SPAT-SCENE object.

The OpenMusic objects are implemented in Common Lisp. We used OpenGL for the 3D rendering.

OpenSoundControl: http://opensoundcontrol.org

ACKNOWLEDGEMENTS

We thank the composers for their time and valuable insights. This work is part of the EFFICACe project (ANR-13-JS02-0004-01).