Bernard Mourrain
email: bernard.mourrai@inria.fr

Nelly Villamizar
email: nellyyv@cma.uio.no

Galaad

HOMOLOGICAL TECHNIQUES FOR THE ANALYSIS OF THE DIMENSION OF TRIANGULAR SPLINE SPACES

The spline space C r k (∆) attached to a subdivided domain ∆ of R d is the vector space of functions of class C r which are polynomials of degree ≤ k on each piece of this subdivision. Classical splines on planar rectangular grids play an important role in Computer Aided Geometric Design, and spline spaces over arbitrary subdivisions of planar domains are now considered for isogeometric analysis applications. We address the problem of determining the dimension of the space of bivariate splines C r k (∆) for a triangulated region ∆ in the plane. Using the homological introduced by Billera (1988), we number the vertices and establish a formula for an upper bound on the dimension. There is no restriction on the ordering and we obtain more accurate approximations to the dimension than previous methods and furthermore, in certain cases even an exact value can be found. The construction makes also possible to get a short proof for the dimension formula when k ≥ 4r + 1, and the same method we use in this proof yields the dimension straightaway for many other cases.

Introduction

Let ∆ be a connected, finite two dimensional simplicial complex, supported on |∆| ⊂ R 2 , with |∆| homotopy equivalent to a disk. We denote by C r k (∆) the vector space of all C r functions on ∆ that, restricted to any simplex in ∆, are given by polynomials of degree less or equal than k. These functions are called splines, and they have many practical applications, including the finite element method for solving differential equations [START_REF] Fix | An analysis of the finite element method[END_REF]). Recently they have also been considered for isogeometric analysis applications [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF]), in computer-aided design for modeling surfaces of arbitrary topology [START_REF] Farin | and surfaces for computer aided geometric design, 3rd Edition[END_REF]). In these application areas, spline functions of low degree are of particular interest, and the degree of smoothness attainable is an important design consideration.

Date: 10 October 2012.

A fundamental problem is to determine the dimension of this vector space as a function of known information about the subdivision ∆. A serious difficulty for solving this problem is that the dimension of the space C r k (∆) can depend not only on the combinatorics of the subdivision, but also on the geometry of the triangulation i.e., how |∆| is embedded in R 2 . In [START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF], the author presented a lower and an upper bound on the dimension of spline spaces of arbitrary degree and smoothness for general triangulations; using Bernstein-Bézier methods, a formula for the dimension for k ≥ 4r + 1 was obtained by [START_REF] Alfeld | The dimension of bivariate spline spaces of smoothness r for degree d ≥ 4r + 1[END_REF]. The result was extended to k ≥ 3r + 2 in [START_REF] Hong | Spaces of bivariate spline functions over triangulation[END_REF], and in [START_REF] Alfeld | On the dimension of bivariate spline spaces of smoothness r and degree d = 3r+1[END_REF] a dimension formula is proved for almost all triangulations for k = 3r + 1. However, there are no explicit formulas for the dimension of the spline spaces C r k (∆) for degree k < 3r + 2 for general triangulations. The use of homological algebra in spline theory was introduced by [START_REF] Billera | Homology of smooth splines: generic triangulations and a conjecture of Strang[END_REF]. He obtained the dimension of C 1 k (∆) for all k for generic triangulations. Schenck and Stillman (1997a), introduced a chain complex different from that used by Billera; this complex was studied by Schenck and Stillman (1997a), Schenck and Stillman (1997b), [START_REF] Schenck | A spectral sequence for splines[END_REF], [START_REF] Geramita | Fat points, inverse systems, and piecewise polynomial functions[END_REF]. The lower homology modules of the chain complex in this construction, differ from the one introduced by Billera, and they have nicer properties. The connection between fat points and the spline space defined on ∆, allows to give a complete characterization of the free resolutions of the ideals generated by power of linear forms appearing in the chain complex, and so to prove the dimension formula for C r k (∆) for sufficiently high degree (Schenck and Stillman (1997b)). The main contribution of the paper is a new formula for an upper bound on the dimension of the spline space. The formula applies to any ordering established on the interior vertices of the partition, contrarily to the upper bound formulas in [START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF], [START_REF] Lai | Encyclopedia of Mathematics and its Applications[END_REF]. Having no restriction on the ordering makes it possible to obtain accurate approximation to the dimension and even exact value in many cases. As a consequence, we give a simple proof for dimension formula when k ≥ 4r + 1.

The paper is structured as follows. In Section 2 we recall the construction of this chain complex and some of the properties of the homology modules. We describe the dimension as a sum of a combinatorial part and the dimension of an homology module. This latter term which happens to be zero in high degree is always non-negative, so the formula gives a lower bound for the dimension of the spline space for any degree k. We present the main result in Section 3, which is an upper bound on the dimension of the spline space. In Section 4 we compare the formulas of the lower and the upper bound with those appearing in [START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF]. The result about the exact dimension for k ≥ 4r + 1 is proved in Section 5. This latter result and some other examples that we present in Section 6, illustrate the interest of the homology construction for proving exact dimension formulas.

Construction of the chain complex

We reproduce some notations and definitions presented in Schenck and Stillman (1997a), restricting them to the case where ∆ is a planar simplicial complex supported on a disk.

Denote by ∆ 0 the set of interior faces of ∆, and by ∆ 0 i (i = 0, 1, 2) the set of i-dimensional interior faces whose support is not contained in the boundary of |∆|. We denote by f 0 i the cardinality of these sets, and by ∂∆ the complex consisting of all 1-faces lying on just one 2face as well as all subsets of them. As it will be convenient to study the dimension of the vector space C r k (∆), we embed ∆ in the plane {z = 1} ⊆ R 3 and form the cone ∆ over ∆ with vertex at the origin. Denote by C r k ( ∆) the set of splines on ∆ of C r -smoothness and degree exactly k. Then C r ( ∆) := ⊕ ≥0 C r k ( ∆) is a graded R-algebra and there is an isomorphism of R-vector spaces between C r k (∆) and the elements in C r ( ∆) of degree exactly k [START_REF] Billera | A dimension series for multivariate splines[END_REF]), in particular , y, z]. For an edge τ ∈ ∆ 0 1 , let ℓ τ denote a non-zero homogeneous linear form vanishing on τ , and define the ideal J (β) of R for each simplex β ∈ ∆ 0 as follows:

dim C r k (∆) = dim C r k ( ∆). Define R := R[x
J (σ) = 0 for each σ ∈ ∆ 0 2 J (τ ) = ℓ r+1 τ for each τ ∈ ∆ 0 1 J (γ) = ℓ r+1 τ τ ∋γ
for each γ ∈ ∆ 0 0 , τ ∈ ∆ 0 1 . Consider the chain complex R defined on ∆ 0 as R i = R f 0 i , with the usual simplicial boundary maps ∂ i used to compute the relative (modulo ∂∆) homology with coefficients in R as in [START_REF] Billera | Homology of smooth splines: generic triangulations and a conjecture of Strang[END_REF]. Let R/J be the chain complex obtained as the quotient of R by J ,

0 -→ σ∈∆ 2 R ∂ 2 -→ τ ∈∆ 0 1 R/J (τ ) ∂ 1 -→ γ∈∆ 0 0 R/J (γ) -→ 0
where the maps ∂ i are the induced by the simplicial boundary maps ∂ i . The complex R/J was introduced by Schenck and Stillman (1997a), and agrees with the complex studied by Billera except at the vertices.

It was shown in [START_REF] Billera | Homology of smooth splines: generic triangulations and a conjecture of Strang[END_REF] that C r ( ∆) is isomorphic to the top homology module of R/J H 2 (R/J ) := ker(∂ 2 ).

Let us consider the short exact sequence of complexes 0 -→ J -→ R -→ R/J -→ 0 that gives rise to the long exact sequence of homology modules

0 → H 2 (R) →H 2 (R/J ) → H 1 (J ) → H 1 (R) → H 1 (R/J ) → H 0 (J ) → H 0 (R) → H 0 (R/J ) → 0 (1)
Since ∆ is supported on a disk, both H 0 (R) and H 1 (R) are zero. Hence, by the long exact sequence, it implies that H 0 (R/J ) is also zero and H 1 (R/J ) is isomorphic to H 0 (J ). Applying the Euler characteristic equation to the chain complex R/J χ(H(R/J )) = χ(R/J ), and considering the modules in degree exactly k, leads to the formula

(2) dim C r k (∆) = 2 i=0 (-1) i β∈∆ 0 2-i dim R/J (β) k + dim H 0 (J ) k .
The aim is to determine the modules in the previous formula as functions of known information about the subdivision ∆. We know that

σ∈∆ 0 2 dim R k = f 0 2 k + 2 2 (3) τ ∈∆ 0 1 dim R/J (τ ) k = f 0 1 k + 2 2 - k + 2 -(r + 1) 2 . (4)
For the computation of dim R/J (γ) k , Schenck and Stillman (1997b) proposed the following resolution for R/J (γ) k . Without loss of generality we translate γ to the origin and assume that the linear forms in J (γ) involve only variables x, y. Letting ℓ r+1 1 , . . . , ℓ r+1 t be a minimal generating set for J (γ), then a free resolution for R/J (γ) k is given by:

(5) 0 → R(-Ω -1) a ⊕ R(-Ω) b → ⊕ t j=1 R(-r -1)-→R → R/J (γ) → 0
where Ω -1 is the socle degree of R/J (γ); Ω and the multiplicities a and b are given by

(6) Ω = t r t -1 + 1, a = t (r + 1) + (1 -t) Ω, b = t -1 -a.
In the case t = 1, we take A = b = Ω = 0 so that R(-Ω -1) a ⊕ R(-Ω) b = 0. Applying this to each vertex γ i ∈ ∆ 0 0 , and letting t i , Ω i , a i and b i denote the values for t, Ω, a and b at γ i respectively, leads to dim

γ i ∈∆ 0 0 R/J (γ i ) k = f 0 0 i=1 k + 2 2 -t i k + 2 -(r + 1) 2 + b i k + 2 -Ω i 2 + a i k + 2 -(Ω i + 1) 2 . ( 7
)
Let us notice that the previous formula does not change when we take t i as the number of different slopes of the edges containing the vertex γ i . Then (3), (4), and ( 7), and the fact that dim H 1 (R/J ) k ≥ 0, yields the following theorem:

Theorem 2.1. The dimension of C r k (∆) is bounded below by dim C r k (∆) ≥ k + 2 2 + f 0 1 k + 2 -(r + 1) 2 - f 0 0 i=1 t i k + 2 -(r + 1) 2 -b i k + 2 -Ω i 2 -a i k + 2 -(Ω i + 1) 2 ,
where t i is the number of different slopes of the edges containing the vertex γ i , and

Ω i = t i r t i -1 + 1, a i = t i (r + 1) + (1 -t i ) Ω i , b i = t i -1 -a i .
In [START_REF] Schenck | A spectral sequence for splines[END_REF], it was proved that the homology module H 1 (R/J ) k vanishes for sufficiently high degree, thus the lower bound in the latter theorem is actually the exact dimension formula for C r k (∆) when k ≫ 0 (Geramita and Schenck, 1998, Theorem 4.2).

3. An upper bound on the dimension of C r k (∆) Let us fix an ordering γ 1 , . . . , γ f 0 0 for the vertices in ∆ 0 0 . For each vertex γ i , denote by N(γ i ) the set edges that contain this vertex, and define ti as the number of different slopes of the edges connecting γ i to one of the first i -1 vertices in the list or to a vertex on the boundary.

Theorem 3.1. The dimension of C r k (∆) is bounded by dim C r k (∆) ≤ k + 2 2 + f 0 1 k + 2 -(r + 1) 2 - i ti k + 2 -(r + 1) 2 -bi k + 2 -Ωi 2 -ãi k + 2 -( Ωi + 1) 2
with ti as we have defined above and

Ωi = ti r ti -1 + 1, ãi = ti (r + 1) + (1 -ti ) Ωi , bi = ti -1 -ãi . if ti > 1 and ãi = bi = Ωi = 0 if ti = 1 .
Proof. By the long exact sequence (1), and the fact that H 1 (R) = 0, we have the short exact sequence

0 -→ H 2 (R) -→ H 2 (R/J ) -→ H 1 (J ) -→ 0
The Euler characteristic equation applied to the complex R implies that

H 2 (R) k = R k .
Hence the isomorphism between C r ( ∆) and the homology module

H 2 (R/J ) implies that dim C r k (∆) = dim R k + dim H 1 (J ) k
where the complex of ideals J (as defined above) is given by

0 -→ τ ∈∆ 0 1 J (τ ) ∂ 1 -→ γ∈∆ 0 0 J (γ) -→ 0 and where H 1 (J ) = ker ∂ 1 = K 1 . Define W 1 := Im (∂ 1 ). By the exact sequence 0 -→ K 1 -→ τ ∈∆ 0 1 J (τ ) -→ W 1 -→ 0 we get dim C r k (∆) = dim R k + τ ∈∆ 0 1 dim J (τ ) k -dim(W 1 ) k
Therefore, to find an upper bound on dim C r k (∆) it is enough to find a lower bound on the dimension of W 1 = Im ∂ 1 in degree k. Let us consider the map

δ 1 : τ =(γ,γ ′ )∈∆ 0 1 J (τ ) [τ ] → γ i ∈∆ 0 0 τ ∈N (γ) R [τ |γ] such that δ 1 ([τ ]) = [τ |γ] -[τ |γ ′ ] for τ = (γ, γ ′ ) ∈ ∆ 0 1
, and the map ϕ 1 :

γ i ∈∆ 0 0 τ ∈N (γ) R [τ |γ] → γ∈∆ 0 0 R [γ] with ϕ 1 ([τ |γ]) = [γ] if γ ∈ ∆ 0 0 , = 0 if γ ∈ ∆ 0 0 .
Then, we have ∂ 1 = ϕ 1 • δ 1 . We consider now the map π 1 :

γ∈∆ 0 0 τ ∈N (γ) R [τ |γ] → γ∈∆ 0 0 τ ∈N (γ) R [τ |γ] with π 1 ([τ |γ]) = 0 if γ is the end point of biggest index of τ and π 1 ([τ |γ]) = [τ |γ] otherwise. We denote by ∂1 = ϕ 1 • π 1 • δ 1 .
For γ ∈ ∆ 0 0 , we denote by Ñ (γ) the set of interior edges τ connecting γ to another vertex which is not of bigger index. Let J (γ) =

τ ∈ Ñ (γ) R ℓ r+1 τ ⊂ J (γ). By construction, we have Im ∂1 = γ∈∆ 0 0 J (γ) [γ] and dim(W 1 ) k = dim(Im ∂ 1 ) k ≥ dim(Im ∂1 ) k .
Using formulas (3), ( 4) and the resolution (5) applied to J (γ), we obtain the upper bound

dim C r k (∆) ≤ k + 2 2 + f 0 1 k + 2 -(r + 1) 2 - f 0 0 i=1 ti k + 2 -(r + 1) 2 -bi k + 2 -Ωi 2 -ãi k + 2 -( Ωi + 1) 2 ,
with ti = | Ñ (γ i )| and Ωi , ãi , bi defined as in (6).

As an immediate consequence of this theorem we mention the following result.

Corollary 3.2. If for a numbering on the vertices in ∆ 0 0 , either t i = ti or ti ≥ r + 1 for every vertex γ i , then the upper bound we get, equals the lower bound, obtaining so the exact dimension formula for the spline space.

Proof. We compare the terms corresponding to each interior vertex in both formulas. If t i = ti then they are trivially the same. By definition (6), for t ≥ r + 1, Ω and a are both constant equal to r + 1. Hence if r + 1 ≤ ti (≤ t i ), the terms in the binomials are the same when we compare them in the formulas in Theorem 2.1 and 3.1. Since the respective terms in t i and ti cancel out then we get the equality of the bounds.

4. The bounds on dim C r k (∆) given by [START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF] In this section we compare the bounds on dim C r k (∆) found by [START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF], with the lower and upper bounds given in the previous two sections. With the notation as before, the upper bound presented by Schumaker can be stated as follows.

Theorem 4.1. [START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF], Theorem 2.1) Suppose that the vertices γ i ∈ ∆ 0 0 of the partition are numbered in such a way that each pair of consecutive vertices in the list are corners of a common triangle in ∆. For each γ i define ti as the number of edges with different slopes joining the vertex γ i to a vertex in the boundary of ∆ or to one of the first i -1 vertices. Then

dim C r k (∆) ≤ k + 2 2 + f 0 1 k -r + 1 2 -f 0 0 k + 2 2 - r + 2 2 + f 0 0 i=1 k-r j=1 (r + j + 1 -j • ti ) + . (8)
In the same article, the author also presented a lower bound on dim C r k (∆). The formula can be obtained by replacing ti by t i in (8). We shall prove that that formula for the lower bound is the same as the formula we presented in Theorem 2.1. Essentially the same proof shows that the upper bound formulas, the one in the previous theorem and the one we presented in Theorem (3.1), coincide. With the exception that the upper bound presented by Schumaker can only be applied for certain numberings on the vertices. That restriction makes sometimes not possible to find an upper bound. In Section 6 we will include some examples of those situations, and also some cases when not having that restriction leads to find the exact dimension of the space.

Remark. The above cited article defined ti as the number of slopes of the edges containing γ i but not containing any of the i -1 vertices in the list. This coincides with the definition of ti we consider here, except that the reverse counting on the vertices in ∆ 0 0 is used. Lemma 4.2. Let t ≥ 2 and Ω defined as in (6) above. The notations in the upper bound formulas (and the lower bounds, respectively) are linked in the following way: if Ω = r + ℓ then r + j + 1 -j • t > 0 only when j ≤ ℓ -1.

Proof. For any value of t ≥ 2 we have

1 < t t -1 ≤ 2 .
Thus, r + 1 ≤ Ω ≤ 2r + 1, and we can write Ω = r + ℓ for some integer ℓ between 1 and r + 1. The interval of t where the value of Ω is r + ℓ is given as follows:

(9)

     r+ℓ ℓ < t ≤ r+(ℓ-1) ℓ-1 when ℓ ≥ 2 t > r + 1 when ℓ = 1.
On the other hand, for fixed r and t, the condition

(10) r + j + 1 -j • t ≥ 1
is satisfied if and only if t ≤ (r + j)/j. Thus, the biggest number j subject to condition (10) must satisfy

r + (j + 1) j + 1 < t ≤ r + j j .
From ( 9), the previous relation holds if and only if Ω = r +(j +1). Let t ≥ 2 and Ω, a and b be defined as in ( 6). We consider the following two formulas:

(11) t k + 2 -(r + 1) 2 -b k + 2 -Ω i 2 -a k + 2 -(Ω + 1) 2 (12) k + 2 2 - r + 2 2 - k-r j=1 (r + j + 1 -j • t) + .
Lemma 4.3. If t > r + 1 and k ≥ r then the formulas (11) and ( 12) are equal.

Proof. By (9) we know that if t > r + 1 one has Ω = r + 1. From the definition in (6), we have a = r + 1 and b = t -(r + 2). Thus, (11) reduces to

       0 if k + 1 -r < 2 (k = r) r + 2 if k -r = 1 k+2 2 -r+2 2 otherwise.
This can be simplified to k+2 2 -r+2 2 for any value of k and r. By Lemma 4.2, the relation (10) in this case is not satisfied by any value of j, then also (12) simplifies to that expression. Lemma 4.4. Let t be an integer ≥ 2. If t ≤ r + 1 and k ≥ r then the formulas (11) and ( 12) are equal.

Proof. By Lemma 4.2 we know that Ω = r + ℓ for some integer ℓ ≥ 2. Depending on the value of ℓ the binomials in (11) could be automatically zero. We need to consider three possible situations k -r ≥ ℓ + 1, k -r = 1 and k -r < ℓ. The formula reduces to a different expression in each case, but it is not difficult to check that they are respectively equivalent to the expressions we get from (12). Proof. Since any interior vertex of ∆ is contained in at least two edges with different slopes, then t i ≥ 2 for every i = 1, 2, . . . f 0 0 . Collecting the vertices γ i with the same value t i , the statement follows directly applying Lemmas 4.3 and 4.4 for those values of t i .

Proposition 4.6. If the vertices in ∆ 0 are numbered as in Theorem (4.1), then the upper bound on dim C r k (∆) given in Theorem 3.1 coincides with the bound (8) in [START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF][Theorem 2.1 p. 252].

Proof. This follows similarly as in the latter proposition. By collecting the vertices in ∆ 0 with the same value ti , we apply Lemmas 4.3 and 4.4 with the values ti . The only case that remains to be considered is ti = 1. It corresponds to the third term in the formula in Theorem (3.1). Since

k + 2 -(r + 1) 2 = k + 2 2 - r + 2 2 - k-r j=1 (r + 1)
the statement follows.

We summarize the results in the following theorem. Let us denote by lbh, ubh and lbs, ubs the respective lower and upper bounds obtained by Theorems 2.1 and 3.1 and the ones obtained in [START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF].

Theorem 4.7. Let ∆ be a connected, finite two dimensional simplicial complex, supported on a disk. Then,

lbh(dim C r k (∆)) = lbs(dim C r k (∆)) ubh(dim C r k (∆)) ≤ ubs(dim C r k (∆))
. Proof. This follows from Propositions 4.5 and 4.6, and the fact that the formula in Theorem 3.1 can be applied to any numbering on the interior vertices of ∆.

Dimension formula for degree k ≥ 4r + 1

In this section we present an alternative proof for the dimension formula of the spline space C r k (∆) for k ≥ 4r + 1. The proof is substantially shorter than the one presented in [START_REF] Alfeld | The dimension of bivariate spline spaces of smoothness r for degree d ≥ 4r + 1[END_REF].

We include the following notation: let ∆ i and ∆ ∂ i be (respectively) the sets of i-dimensional faces and i-dimensional boundary faces. Denote by f i (∆) and f ∂ i (∆) the cardinality of the preceding sets, respectively. We begin by stating the following lemma which we will need later on.

Lemma 5.1. (Schenck and Stillman, 1997a, Lemma 3.3) If ∆ is a triangulated region in R 2 , then there exists a total order on ∆ 0 such that for every γ in ∆ 0 0 , there exist vertices γ ′ , γ ′′ adjacent to γ, with γ ≻ γ ′ , γ ′′ , and such that γγ ′ , γγ ′′ have distinct slopes.

For an ordering on ∆ 0 as in the previous lemma, we assign indices to the vertices in such a way that γ i ≻ γ j when i > j. The first indices 1, 2, . . . , f ∂ 0 are assigned to the vertices lying on the boundary. To those interior vertices which are joined to the boundary by two or more edges of distinct slope we assign the indices f ∂ 0 + 1, . . . , n. Let us recall the following notation and remarks from Schenck and Stillman (1997a). For each interior vertex γ, and each f ∈ J (γ), let f [γ] denote the corresponding element in H 0 (J ). Then H 0 (J ) is generated by

{f [γ] | f ∈ J (γ)}. By definition of J (γ), we know that (13) f [γ] = τ ∋γ ℓ r+1 τ f τ [γ]
for some polynomials f τ ∈ R. Notice that if τ is an edge whose vertices are γ and γ ′ then the ( 14)

ℓ r+1 τ f [γ] = ℓ r+1 τ [γ ′
] in H 0 (J ); in particular ℓ r+1 τ f [γ] = 0 when τ is an edge connecting γ to the boundary.

Here is another lemma that we need to prove that H 0 (J k ) = 0 for degree k ≥ 4r + 1: Lemma 5.2. Let ℓ 1 , ℓ 2 and ℓ 3 be three equations of distinct lines through a point p and L the equation of a line that does not contain the point p. Then for any polynomial g of degree d ≥ r + 1, there exist u, v ∈ R of degree d and w ∈ R of degree r -1 such that

ℓ r+1 3 g = ℓ r+1 1 u + ℓ r+1 2 v + L d-r+1 ℓ r+1 3 w,
with w = 0 in the case r = 0.

Proof. The case r = 0 is direct from the linear dependency of ℓ 1 , ℓ 2 and ℓ 3 . Suppose that r > 0. Since ℓ 1 , ℓ 2 and L are linearly independent, we can make the following change of coordinates:

ℓ 1 = x ; ℓ 3 = x + ay for some a = 0 ; ℓ 2 = y ; L = z .
Let g = x i y j z k be a monomial such that i

+ j + k = d. If k ≤ d -r then i + j ≥ r and the polynomial ℓ r+1 2 g is in the ideal ∈ x r+1 , y r+1 . If k ≥ d -r + 1 then g is a multiple of z d-r+1 = L d-r+1 . Thus we can write ℓ r+1 3 g = ℓ r+1 1 u ′ + ℓ r+1 2 v ′ + L d-r+1 w ′ , for some polynomials u ′ , v ′ ∈ R of degree d and w ′ ∈ R of degree 2 r. As L d-r+1 w ′ is in the ideal ℓ r+1 1 , ℓ r+1 2 , ℓ r+1
3 and as L is a non-zero divisor modulo this ideal, we deduce that w ′ = ℓ r+1

1 u ′′ + ℓ r+1 2 v ′′ + ℓ r+1 3 w ′′ , with u ′′ , v ′′ , w ′′ ∈ R of degree r -1. Collecting the coefficients of ℓ r+1
1 , ℓ r+1 2 , we obtain the desired decomposition:

ℓ r+1 3 g = ℓ r+1 1 (u ′ + L d-r+1 u ′′ ) + ℓ r+1 2 (v ′ + L d-r+1 v ′′ ) + L d-r+1 ℓ r+1 3 w ′′ .
We use this lemma to prove the following result:

Theorem 5.3. The dimension of C r k (∆) when k ≥ 4r + 1 is given by the lower bound formula of Theorem 2.1.

Proof. From (2), and the formulas for the dimension of the modules (3), ( 4) and ( 7), it suffices to show that H 0 (J ) k = 0 for k ≥ 4r + 1. Equivalently, we need to show that f [γ] = 0 in H 0 (J ) k for all γ and f ∈ J (γ) k . Ordering the vertices as in Lemma 5.1, we consider the first interior vertex in the ordering, and denote it by γ. Let τ 1 , τ 2 , . . . , τ t be the edges (not necessarily with different slopes) containing γ and ω 1 , ω 2 , ..., ω t be respectively the end points of τ i , which are distinct from γ. We number first the edges τ 1 and τ 2 that connect γ to the boundary, and the remaining edges counterclockwise, starting from τ 2 . Let us denote by ℓ 1 , ℓ 2 , . . . , ℓ t (respectively) equations of these edges and by L i equations of the edges connecting ω i to ω i-1 for i = 2, . . . t.

ℓ 1 ℓ 2 ℓ 3 ℓ i-1 ℓ i L i-1 L i L 3 γ ω 3 ω i-1 ω i ω 1 ω 2 Figure 1. Notation in Theorem 5.3.
We are going to prove by induction on j > 2 that for any homogeneous polynomial f of degree ≥ 3r, we have

• ℓ r+1 j f [ω j ] = ℓ r+1 j f [γ] = 0, • L r+1 j f [ω j ] = 0 in H 0 (J ).
Let us prove it first for j = 3. Let f ∈ R k with k ≥ 3r. By construction ω 2 is on the boundary, thus we have

L r+1 3 f [ω 3 ] = L r+1 3 f [ω 2 ] = 0. By Lemma 5.2, ℓ r+1 3 f = ℓ r+1 1 u + ℓ r+1 2 v + L 2r+1 3 ℓ r+1
3 w for some polynomials u, v, w ∈ R. Then we have

ℓ r+1 3 f [γ] = ℓ r+1 1 u [γ] + ℓ r+1 2 v [γ] + L 2 r+1 3 ℓ r+1 3 w [γ] = L 2 r+1 3 ℓ r+1 3 w [ω 3 ] = 0.
Let us take now i > 2 and assume that the induction hypothesis is true for 3 ≤ j < i. We consider first ℓ r+1 i f [ω i ] with f homogeneous of degree ≥ 3r. By Lemma 5.2 applied to ℓ 1 , ℓ 2 , ℓ 3 and L i , we have

ℓ r+1 i f = ℓ r+1 1 u + ℓ r+1 2 v + L 2r+1 i ℓ r+1 i w,
for some polynomials u, v, w ∈ R. Then we have

ℓ r+1 i f [γ] = ℓ r+1 1 u [γ] + ℓ r+1 2 v [γ] + L 2r+1 i ℓ r+1 i w [γ] = L 2r+1 i ℓ r+1 i w [ω i ] = L 2r+1 i ℓ r+1 i w [ω i-1 ].
As L 2r+1 i is in the ideal generated by ℓ r+1 i-1 , L r+1 i-1 , by the induction hypothesis, we have

L 2r+1 i ℓ r+1 i w [ω i-1 ] = 0. This proves our claim for ℓ r+1 i f [ω i ] = 0. Let us consider now L r+1 i f [ω i ] = L r+1 i f [ω i-1 ] with f homogeneous of degree ≥ 3r. By Lemma 5.2, L r+1 i f = ℓ r+1 i-1 u + L r+1 i-1 v + ℓ r+1 i L r+1 i w,
for some polynomials u, v, w ∈ R. We deduce that

L r+1 i f [ω i-1 ] = ℓ r+1 i-1 u [ω i-1 ] + L r+1 i-1 v [ω i-1 ] + ℓ r+1 i L r+1 i w [ω i-1 ] = ℓ r+1 i L r+1 i w [ω i ] = ℓ r+1 i L r+1 i w [γ] = 0
by the induction hypothesis and the previous computation.

This concludes the induction proof and shows that for any f ∈ J (γ) 4r+1 , we have

f [γ] = t i=1 ℓ r+1 i f i [γ] = t i=3 ℓ r+1 i f i [γ] = 0.
Therefore H 0 (J ) = 0 and the dimension of C r k (∆) is given by the lower bound when k ≤ 4r + 1.

Examples and remarks

In this section we compute the dimension formula for some complexes. We prove that the dimension formula is in fact given by the corresponding lower bound given in Theorem 2.1. In the first example we show directly that H 0 (J ) k = 0. In the second example we number the vertices in such a way that the upper bound we obtain agrees with the lower one. The subdivisions we consider are two of the so-called Powell-Sabin subdivisions [START_REF] Powell | Piecewise quadratic approximations on triangles[END_REF]).

Example 6.1. Powell-Sabin 12-split.

Let ∆ be the simplicial complex supported on a triangle |∆| subdivided into twelve smaller triangles, as in Figure 2. Let us note that by using any numbering on ∆ 0 0 , the upper bound we get from Theorem 3.1 equals the lower bound, leading easily to the dimension formula. On the other hand, it is not possible to find a numbering on the interior vertices of ∆ with the condition required to apply Theorem 4.1. Hence it is not possible to find an upper bound for this spline space by using the formula of Schumaker.

By the remarks in Section 5, it is easy to check that ℓ r+1 j f [γ i ] = 0 for every γ i and ℓ r+1 j f ∈ J (γ i ). Thus, for any k and any r ≤ k, the homology module H 0 (J ) k is zero and the dimension of the spline space C r k (∆) is given by the formula in Theorem 2.1. This example is an instance of a triangulation where all edges are pseudoboundaries (in the terminology of Schenck and Stillman (1997b)), and the computation also follows by Lemma 2.5 of Schenck and Stillman (1997b).

Example 6.2. Powell-Sabin 6-split.

Let ∆ be the simplicial complex supported on a simply connected triangulated region |∆| in R 2 . Assume the triangulation consists of f 0 vertices (interior and exterior) and f 0 i i-dimensional interior faces (i = 0, 1, 2), where f 0 2 is taken as the number of all the triangles σ i in the subdivision. We refine this triangulation by subdividing each triangle into 6 triangles to get ∆, as follows.

(1) In each triangle σ i ∈ ∆ 0 2 choose an interior point ν i , in such a way that if two triangles σ i , σ j have an edge τ in common, then the line joining ν i and ν j intersects τ at an interior point µ ij .

(2) Join each new point ν i to the vertices of the triangle σ i , and to the points µ ij on the edges of σ i . (3) For the triangles having an edge (or more) on the boundary, choose a point in the interior of such edge and join it with ν i .

∆Figure 3. Powell-Sabin 6-split of ∆.

We want to apply Corollary 3.2 and give a formula for dim C 1 k ( ∆). Let us consider the following numbering on the vertices in ∆0 0 . We take first a triangle σ i ∈ ∆ with (at least) one of its edges on the boundary, denote by γ i the vertex in σ i that is in the interior of ∆. With the notation as above, we assign the index 1 to the vertex ν i , the indices 2 and 3 to the vertices µ ij on the edges of σ i , and finally the index 4 to the vertex γ i . After this, we consider the edges of σ i as if they were all part of the boundary and iterate the previous process. We choose another triangle in ∆ having (at least) one edge on the "new" boundary. We index its vertices following the same order as we did for the ones in σ i . We continue this numbering until we have considered all the triangles in ∆ and hence have index all the vertices in ∆0 0 . Figure 4 shows the triangles that could be considered first, and the different colors refer to the values of ti for the interior vertices in each case. Let us notice that for this numbering, the values for ti corresponding to these vertices are either ≥ 3, as is the case of the vertices ν i and γ i , or ti = t i = 2, as it happens for the vertices µ ij . Thus, we arrive to the following result. Proposition 6.3. For any k ≥ 3:

dim C 1 k ( ∆) = k + 2 2 + f 0 1 k -2 2 + 2f 0 0 k -1 2 + 2k -1 2 (f 0 -2).
For k = 2: dim C 1 2 ( ∆) = 3f 0 . Proof. It follows directly from the remarks we made above about ti for the counting, and Corollary 3.2.

For the numbering we defined on ∆0 0 , it is not possible to apply Theorem 4.1 to find an exact value for the dimension. In fact, when the initial subdivision consists of more than one triangle, by using the numberings allowed in that theorem, we get upper bounds which are strictly bigger that the actual dimension of the space.

Remark 1. In this paper we have confined our attention to triangular partitions but it is easy to check that our proofs of the lower and the upper bounds can be extended to rectilinear partitions. The construction also allows to consider regions subdivided by curved boundaries, but then the ideals must be considered individually. The boundary of |∆| itself can, of course, be curved.

Remark 2. The resolution for ideals generated by power of linear forms presented in the paper by [START_REF] Geramita | Fat points, inverse systems, and piecewise polynomial functions[END_REF], applied also when the power of the linear forms is different. It makes possible to apply the ideas we present in this paper to mixed splines, i.e., splines where the order of smoothness may differ on the various edges of the subdivision.

Remark 3. The proof of the exact dimension formula when the degree of the spline space k is ≥ 4r + 1 cannot be directly extended to spaces of polynomials of degree k ≥ 3r + 1. With some restrictions associated to the number of different slopes of the edges containing each vertex it is possible to prove the result for this degree.

Remark 4. We hope that the two methods we used in the examples to find the exact dimension, namely, by showing that the lower and the upper bounds coincide for certain numbering on the vertices, or that H 0 (J ) = 0 directly by considering the equations of the edges, illustrate the way to easily find the dimension of the spline space for many particular triangulations.

Remark 5. There is a version of the upper bound given by [START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF] in [START_REF] Lai | Encyclopedia of Mathematics and its Applications[END_REF], where the constraint on the numbering of the vertices is relaxed. However, for this new version of the theorem, there are still examples for which the upper bound formula does not give the correct result, see for instance the example p. 243 of [START_REF] Lai | Encyclopedia of Mathematics and its Applications[END_REF]. Remark 6. For the three dimensional version of the problem, it is necessary to analyse ideals generated by powers of linear forms in three variables. The formula of the dimension, analogous to formula (2) in Section 2, have two homology modules, and in order to approximate the dimensions the dimension of these two modules must be bound. Some issues and positive results about this problem will be in an upcoming paper.

Proposition 4. 5 .

 5 The lower bound on dim C r k (∆) given in Theorem 2.1 coincides with the lower bound in[START_REF] Schumaker | Bounds on the dimension of spaces of multivariate piecewise polynomials[END_REF][Theorem 3.1 p. 256].

Figure 2 .

 2 Figure 2. Powell-Sabin 12-split.

Figure 4 .

 4 Figure 4. Numbering on the vertices in ∆0 0 .


Jonathan Le Roux
email: jonathan.le-roux@normalesup.org

E Vincent

come

Consistent Wiener filtering for audio source separation

Consistent Wiener Filtering for Audio Source Separation

Jonathan Le Roux, Member, IEEE, and Emmanuel Vincent, Senior Member, IEEE Abstract-Wiener filtering is one of the most ubiquitous tools in signal processing, in particular for signal denoising and source separation. In the context of audio, it is typically applied in the time-frequency domain by means of the short-time Fourier transform (STFT). Such processing does generally not take into account the relationship between STFT coefficients in different time-frequency bins due to the redundancy of the STFT, which we refer to as consistency. We propose to enforce this relationship in the design of the Wiener filter, either as a hard constraint or as a soft penalty. We derive two conjugate gradient algorithms for the computation of the filter coefficients and show improved audio source separation performance compared to the classical Wiener filter both in oracle and in blind conditions. Index Terms-Wiener filtering, Short-time Fourier transform, Spectrogram consistency, Source separation, Conjugate gradient EDICS Category: SAS-ICAB

I. INTRODUCTION

W IENER filtering is one of the most widely used tools in signal processing, in particular for signal denoising and source separation. In the context of audio, where signals are not stationary but short-term stationary, it is typically applied in the timefrequency domain via the short-time Fourier transform (STFT). The filter coefficients, which are nonnegative, are independently computed in each time-frequency bin from the variances of the signals to be separated without accounting for the consistency of the obtained set of STFT coefficients, i.e., the fact that they actually correspond to the STFT of a time-domain signal [START_REF] Le Roux | Explicit consistency constraints for STFT spectrograms and their application to phase reconstruction[END_REF]. Together with inaccurate estimation of the signal variances, this independent processing is one of the reasons behind the presence of residual interference or musical noise in the separated signals [START_REF] Vincent | The signal separation evaluation campaign (2007-2010): Achievements and remaining challenges[END_REF].

One way to promote consistency of the separated STFT coefficients is to employ phase reconstruction methods, which estimate consistent phases given magnitude spectra [START_REF] Le Roux | Explicit consistency constraints for STFT spectrograms and their application to phase reconstruction[END_REF], [START_REF] Griffin | Signal estimation from modified short-time Fourier transform[END_REF]. Multiple input spectrogram inversion (MISI) [START_REF] Gunawan | Iterative phase estimation for the synthesis of separated sources from single-channel mixtures[END_REF] and partitioned phase retrieval (PPR) [START_REF] Sturmel | Iterative phase reconstruction of Wiener filtered signals[END_REF] methods were specifically designed to handle multiple sources. Because they modify phase only, these methods implicitly assume that the input magnitude spectra are close to the true source spectra. This assumption is valid in the context of informed source separation [START_REF]Informed source separation using iterative reconstruction[END_REF] but not in the usual context of blind or semi-blind source separation of interest here, where both the estimated source variances and the magnitude spectra resulting from Wiener filtering may differ from the true source spectra. As a result, these methods provide limited improvement as we shall see in the following.

In earlier work [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF], we proposed a complex-valued variant of the Wiener filter promoting consistency of the separated STFT coefficients via a soft penalty term. In contrast to the above methods, the magnitude and phase of the coefficients are jointly optimized and the resulting filter retains all properties of the classical Wiener filter, namely minimum mean square error (MMSE) optimality under Gaussian assumptions and exact reconstruction of the mixture from the separated sources. Although the problem can be formulated as the minimization of a quadratic function, exact inversion of the associated linear system is typically intractable. In [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF], we proposed an iterative auxiliary function-based algorithm for mixtures of two sources which converged slowly for large values of the penalty weight parameter, so that a dynamical update scheme for this parameter had to be introduced. Another algorithm for the joint modification of magnitude and phase was recently proposed in [START_REF]Informed source separation using iterative reconstruction[END_REF] under the name of informed source separation using iterative reconstruction (ISSIR).

We present here two new optimization algorithms based on the conjugate gradient method [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF] with well-chosen preconditioners, which enforce consistency either as a hard constraint or as a soft penalty for any number of sources. We show that, for a suitable choice of the penalty weight, the latter provides better separation performance than the former. It also outperforms the algorithms in [START_REF] Gunawan | Iterative phase estimation for the synthesis of separated sources from single-channel mixtures[END_REF]- [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF] both in oracle conditions, when the variances of the sources are known, and in blind conditions, when they are estimated from the noisy input.

After presenting a general maximum-likelihood formulation of the problem in Section II, we describe the proposed algorithms in Sections III and IV, and evaluate their performance through audio source separation experiments in Section V.

II. WIENER FILTERING AND CONSISTENCY

Let us consider an observed single-channel audio signal x(t) that is a mixture of J sound sources sj(t),

x(t) = J j=1 sj(t), (1)
where J is known and t denotes discrete time. This is equivalently written in the time-frequency domain as

X nf = j S jnf . (2)
Assuming that the STFT coefficients S jnf of the sources are independent zero-mean Gaussian random variables with positive variance v jnf in each time frame n and frequency bin f , the STFT coefficients X nf of the mixture also follow independent zero-mean Gaussian distributions with variance

v x nf = J j=1 v jnf . (3)
We are interested in the conditional distribution of the sources S jnf given the mixture X nf in each time-frequency bin (n, f ). Due to the constraint (2), this conditional distribution lies on a subspace of dimension J -1 so we focus on a subset of free variables. Without loss of generality, we consider the first J = J -1 sources as free variables given the mixture and denote them as S nf = [S 1nf , . . . , S J nf ] T . Since S nf and X nf are jointly Gaussian with zero mean, their conditional distribution p(S nf |X nf ) is also Gaussian and its mean μnf and covariance Λ -1 nf are given by (see for example [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF])

μnf = Σ SX Σ -1 XX X nf , (4)
Λ -1 nf = Σ SS -Σ SX Σ -1 XX Σ XS . (5)
Here,

Σ XS = Σ T SX = [v 1nf , . . . , v J nf ], ΣXX = v x nf
, and Σ SS is the diagonal matrix with entries v 1nf , . . . , v J nf . The conditional mean simplifies in particular to

μnf = v 1nf v x nf , . . . , v J nf v x nf T X nf , (6)
which is none other than the classical Wiener filter output, i.e., the MMSE estimate of the sources. Using Woodbury's identity on Eq. ( 5), we obtain a simple expression for the precision Λ nf :

Λ nf =     1 v 1nf 0 . . . 0 1 v J nf     + 1 v Jnf    1 • • • 1 . . . . . . . . . 1 • • • 1    . ( 7
)
We shall consider in the following the negative log-likelihood of the conditional distribution -log P (S|X), which is equal up to a constant to the quadratic loss function

ψ(S) = n,f (S nf -μnf ) H Λ nf (S nf -μnf ), (8)
where S and X respectively denote the sets of all variables S nf and X nf . Without further constraint on S, the maximum-likelihood solution is obviously the classical Wiener filter output S = μ. However, we need to keep in mind that the STFT, when classically computed using overlapping windows, is a redundant representation which results in a certain relationship between the coefficients in different time-frequency bins. Denoting by N the number of time frames and F the number of frequency bins, let us consider the set S of arrays S of C J ×N ×F such that S jn,F -f = Sjnf for all j, n and f . The STFT of J real-valued time-domain signals is an element of S, which we call consistent, but not all elements of S can be obtained as such [START_REF] Le Roux | Explicit consistency constraints for STFT spectrograms and their application to phase reconstruction[END_REF], [START_REF] Griffin | Signal estimation from modified short-time Fourier transform[END_REF]. If the inverse STFT (iSTFT) is defined in such a way that time-domain signals can be exactly reconstructed from their STFT through iSTFT without additional time-varying normalization, a necessary and sufficient condition for an element S of S to be consistent is that it is equal to the STFT of its iSTFT [START_REF] Le Roux | Explicit consistency constraints for STFT spectrograms and their application to phase reconstruction[END_REF] or, in other words, that it belongs to the null space Ker F of the R-linear operator F from S to itself defined by

F = Id -STFT • iSTFT (9)
where Id denotes the identity operator and the STFT and iSTFT operators are separately applied to each signal. The classical Wiener filter output S = μ does generally not satisfy this constraint, so that the STFT of the resynthesized signals iSTFT( μ) differs from μ and no longer maximizes the conditional log-likelihood.

III. CONSISTENCY AS A HARD CONSTRAINT

A first natural approach to enforcing consistency in the Wiener filter is to minimize ψ(S) under the constraint that F(S) = 0. This is equivalent to finding time-domain source signals s(t) = [s1(t), . . . , s J (t)] T minimizing ψ(STFT(s)).

Setting the gradient of ψ(STFT(s)) with respect to s to 0, we see that its global minimum verifies

STFT * • Λ • STFT(s) = STFT * • Λ( μ) (10)
where STFT * is the adjoint of the STFT and the operator Λ denotes multiplication by the matrix Λ nf in each time-frequency bin (n, f ), i.e., Λ(S) nf = Λ nf S nf . Under the common assumption that the synthesis window of the iSTFT is equal to the analysis window of the STFT up to a multiplicative constant, it can be shown that STFT * = M iSTFT where M is the frame length [START_REF] Yang | A study of inverse short-time Fourier transform[END_REF]. The global minimum is therefore equivalently given by iSTFT

• Λ • STFT(s) = iSTFT • Λ( μ). ( 11
)
Algorithm 1 Conjugate gradient with hard constraint Inputs:

Λ, μ, > 0 s0 ← iSTFT( μ) r0 ← iSTFT • Λ( μ) -iSTFT • Λ • STFT(s0) p 0 ← iSTFT • Λ -1 • STFT(r0) δnew ← r0, p 0 k ← 0 repeat q k ← iSTFT • Λ • STFT(p k ) α k ← δnew p k ,q k s k+1 ← s k + α k p k r k+1 ← r k -α k q k z k+1 ← iSTFT • Λ -1 • STFT(r k+1 ) δold ← δnew δnew ← r k+1 , z k+1 β k ← δnew δ old p k+1 ← z k+1 + β k p k k ← k + 1 until α 2 k-1 p k-1 2 < s k 2 return s k
Although the invertibility of the Hermitian operator iSTFT • Λ • STFT seems difficult to prove, we believe that it always holds in practice. For it to be singular, there would need to exist signals whose STFTs, once multiplied in each bin by the real-valued matrix Λ nf , would become perfectly inconsistent despite the lack of phase modification. We believe this situation to be unlikely to exist, and indeed always observed finite condition numbers in our experiments (see Fig. 1).

A solution to the linear system (11) can be iteratively estimated through the preconditioned conjugate gradient method [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]. Because the operator Λ is typically badly conditioned, iSTFT • Λ • STFT is itself badly conditioned and the choice of a suitable preconditioner M is crucial for fast convergence. We propose to choose the preconditioner as

M -1 = iSTFT • Λ -1 • STFT ( 12
)
where Λ -1 denotes multiplication by the matrix Λ -1 nf in each timefrequency bin (n, f ), such that M -1 approximates the inverse of the operator to be inverted.

The resulting algorithm is described in Algorithm 1, with ., . denoting the dot product. The computational complexity of each iteration is dominated by four STFT or iSTFT operations.

IV. CONSISTENCY AS A SOFT PENALTY

As an alternative approach, we consider relaxing the hard consistency constraint F(S) = 0, which may be inadequate when the estimated source variances are unreliable, by introducing the L 2 norm of F(S) as a soft penalty term in [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF] with some weight γ. We thus obtain a new objective function:

ψγ(S) = ψ(S) + γ n,f F(S) nf 2 . ( 13
)
If the weight of the penalty is chosen sufficiently large, the estimated spectrograms should finally both be consistent and minimize ψ among all consistent spectrograms. Setting the gradient of ψγ(S) with respect to S to 0, we see that its global minimum verifies

(Λ + γF * • F )(S) = Λ( μ). (14)
It can easily be shown that F is a projector, i.e., F • F = F, and that, under the same assumption on synthesis and analysis windows Algorithm 2 Conjugate gradient with soft penalty Inputs:

Λ, μ, γ > 0, > 0 S0 ← μ R0 ← -γF(S0) P 0 ← (Λ + γ F N -T F N Id) -1 (R0) δnew ← R0, P 0 k ← 0 repeat Q k ← Λ(P k ) + γF(P k ) α k ← δnew P k ,Q k S k+1 ← S k + α k P k R k+1 ← R k -α k Q k Z k+1 ← (Λ + γ F N -T F N Id) -1 (r k+1 ) δold ← δnew δnew ← R k+1 , Z k+1 β k ← δnew δ old P k+1 ← Z k+1 + β k P k k ← k + 1 until α 2 k-1 P k-1 2 < S k 2 return S k
as in Section III, it is Hermitian, i.e., F * = F. As a consequence, F * • F = F, and ( 14) becomes

(Λ + γF)(S) = Λ( μ). ( 15
)
We use the preconditioned conjugate gradient method to invert (15). As the operator Λ + γF is ill-conditioned, the choice of a good preconditioner is crucial again. Denoting by T the signal length, the projector F can be diagonalized on two eigenspaces: its null space Ker F for the eigenvalue 0 with multiplicity J T and its image space Im F for the eigenvalue 1 with multiplicity J (F N -T ). On average, the eigenvalues of F are thus equal to F N -T F N . A suitable choice for the preconditioner M is then

M -1 = Λ + γ F N -T F N Id -1 , ( 16
)
where M -1 amounts to binwise matrix multiplication. The resulting algorithm is described in Algorithm 2. The computational complexity of each iteration is dominated by two STFT or iSTFT operations, which is comparable to the complexity of the algorithm in [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF] or half that of Algorithm 1.

V. EXPERIMENTAL EVALUATION

A. Setup

We now evaluate the separation performance of the proposed algorithms on mixtures of single-channel speech and real-world nonstationary noise signals taken from the development set of the 2010 Signal Separation Evaluation Campaign (SiSEC) for the task "Source separation in the presence of real-world background noise" [START_REF] Vincent | The signal separation evaluation campaign (2007-2010): Achievements and remaining challenges[END_REF]. The noise signals were recorded in a subway, on a square and in a cafeteria. For each of nine pairs of speech and noise signals in the SiSEC data, we created mixtures at three input Signal to Noise Ratios (SNR): -10, 0 and +10 decibels (dB). In order to avoid dependency of the algorithm behavior on the scale of Λ, each mixture signal was normalized so that its root mean square amplitude was equal to 0.063. All signals were sampled at 16 kHz and had a duration of 10 s. The STFT and iSTFT were implemented as in [START_REF] Vincent | BSS Oracle Toolbox Version 2.1[END_REF] and computed using half-overlapping sine windows of length M = 1024.

The soft penalty-based algorithm was implemented with a fixed weight γ, for γ ∈ {10 0 , 10 1 , 10 2 , 10 3 , 10 4 , 10 5 , 10 6 }. The algorithm using a hard constraint will be formally denoted by γ = +∞. The stopping criterion was chosen as = 10 -6 , so that the separation performance is similar to that obtained by exact inversion of ( 11) and ( 15), as verified by preliminary experiments with much smaller . Note that is the sole parameter of both algorithms, while the algorithm in [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF] required additional parameters to be set for the dynamic update scheme of the penalty weight. For comparison, we also evaluated the classical Wiener filter output μ, which corresponds to γ = 0, our previous method [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF], PPR [START_REF] Sturmel | Iterative phase reconstruction of Wiener filtered signals[END_REF], and the result of applying MISI [START_REF] Gunawan | Iterative phase estimation for the synthesis of separated sources from single-channel mixtures[END_REF] and single-resolution ISSIR [START_REF]Informed source separation using iterative reconstruction[END_REF] to the magnitude of μ, using the same stopping criterion as above and default values for the other parameters as given in [START_REF] Gunawan | Iterative phase estimation for the synthesis of separated sources from single-channel mixtures[END_REF]- [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF]. Note that stopping some of the methods prior to convergence may lead to a shorter computation time with similar performance [START_REF] Sturmel | Iterative phase reconstruction of Wiener filtered signals[END_REF], but we preferred to use the same convergence criterion for all methods to ease comparison. We assessed all algorithms in two different conditions: an oracle condition yielding an upper bound on performance, in which the source variances are set to the true power STFTs of the speech and noise signals, and a blind condition, in which only the stationary noise spectrum is assumed to be known and fixed to the time average of the power STFT of the noise signal and the speech variance is estimated from the mixture by means of spectral subtraction [START_REF] Boll | Suppression of acoustic noise in speech using spectral subtraction[END_REF]. For each algorithm, we evaluated the separation performance via the overall Signal-to-Distortion ratio (SDR) and further analyzed it by means of the Signal-to-Interference Ratio (SIR) and Signal-to-Artifacts Ratio (SAR) metrics [START_REF] Vincent | The signal separation evaluation campaign (2007-2010): Achievements and remaining challenges[END_REF]. We also report the computation time using Matlab on a 2.93 GHz quad core CPU.

B. Results

In Fig. 1, we show the average condition numbers of the linear systems [START_REF] Vincent | BSS Oracle Toolbox Version 2.1[END_REF] and (15) with and without preconditioning in the oracle case (the blind case leads to similar numbers). This quantity is equal to the ratio of the largest and smallest eigenvalues of the operators to be inverted, which can be computed via the well-known power iteration algorithm. For all values of γ, the proposed preconditioners reduce the condition numbers by several orders of magnitude and hence lead to a considerable convergence speedup [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF].

In Fig. 2, we show the SDR achieved by the proposed constraintbased and penalty-based algorithms as a function of γ. We see that

TABLE I

Comparison of the SDR, SIR, SAR (dB) and computation time (s), averaged over all mixtures for each input SNR, for the classical Wiener filter, MISI [START_REF] Gunawan | Iterative phase estimation for the synthesis of separated sources from single-channel mixtures[END_REF], ISSIR [START_REF]Informed source separation using iterative reconstruction[END_REF], PPR [START_REF] Sturmel | Iterative phase reconstruction of Wiener filtered signals[END_REF], our previous method [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF], and the proposed penalty-based algorithm with γ chosen by cross-validation for each input SNR, in oracle and blind conditions.

Input the soft penalty-based approach with γ = 10 5 is generally the best overall, especially in blind conditions when the noise and speech variances are inaccurately estimated. Numerical results for the classical Wiener filter, MISI, ISSIR, PPR, our previous method and the proposed penalty-based algorithm are given in Table I. The weight γ in the proposed algorithm is set for each file and each input SNR by cross-validation, using the weight that led to highest average SDR on all other files (note that similar results are obtained with γ fixed to 10 5 ). These confirm that the proposed method leads to a significant improvement over the classical Wiener filter and all previous algorithms in terms of SDR for both oracle and blind conditions, and to similar or better performance compared to our previous algorithm [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF]. In the oracle case, the SIR and SAR also improve, while in the blind case the increase of the SDR can be further analyzed as a strong improvement of the SIR, offset by a small deterioration of the SAR. Informal listening tests support this conclusion in terms of the perceptual salience of interference and artifacts. The perceptual quality improvement is particularly noticeable for the SNR of -10 dB in the oracle case and for the SNR of 0 dB in the blind case. MISI and PPR did not lead to improvement in blind conditions, while ISSIR led to a smaller improvement than the proposed method. The computation time is in almost all cases shorter than the signal duration. The proposed method is faster than [START_REF] Le Roux | Consistent Wiener filtering: Generalized time-frequency masking respecting spectrogram consistency[END_REF] for better or similar performance, and runs comparably to other methods but with better performance.

We also tested MISI and ISSIR on the magnitude spectrograms obtained by spectral subtraction in blind conditions, leading to worse performance by MISI and slightly better for ISSIR while still lower than our method (results not shown in the table).

Fig. 1 .

 1 Fig. 1. Condition numbers of the linear systems (11) and (15) in oracle conditions with (+P ) and without (-P ) preconditioning, for the classical Wiener filter (γ = 0) and for the proposed constraint-based (γ = +∞) and penalty-based algorithms, averaged over all mixtures for each input SNR.

Fig. 2 .

 2 Fig. 2. Output SDR of the classical Wiener filter (γ = 0) and the proposed constraint-based (γ = +∞) and penalty-based algorithms in oracle and blind conditions, averaged over all mixtures for each input SNR.

VI. CONCLUSION

We showed how to account for the redundancy of the STFT in Wiener filtering and presented efficient optimization algorithms based on the conjugate gradient method with suitable preconditioning. The proposed algorithms lead to significant improvements of source separation performance under both oracle and blind conditions. Future work will deal with the extension of these algorithms to the multichannel case.


Jonathan Le Roux

come

Consistent Wiener filtering for audio source separation


Iacopo ; Galli
email: iacopo.galli@ino.it

Pastor, Gianfranco Pablo Cancio

Di Lonardo

Luciano Fusina

Giovanni Giusfredi

Davide Mazzotti

Filippo Tamassia

Paolo De Natale

Pablo Cancio Pastor

Gianfranco Di Lonardo

The ν3 band of 14C16O2 molecule measured

Keywords: radiocarbon dioxide, infrared spectroscopy, v3 fundamental band, optical frequency comb, cavity ring-down 14 C 16 O 2, infrared spectroscopy, ν 3 fundamental band, optical frequency comb, cavity ring-down

Introduction

Infrared spectroscopy of radiocarbon dioxide ( 14 C 16 O 2 ) has been the subject of very few papers [1][2][3][4][5][6]. The most recent systematic works on the ν 3 fundamental band were performed by Sams and De Voe in 1988 [5], and by Dobos et al. [6] in 1989. In Ref. [5], they measured twelve P and R-branch transitions by diode laser spectroscopy. The lines were assigned an uncertainty ranging between 0.001 and 0.005 cm -1 . In Ref. [6] a medium resolution (0.18 cm -1 ) FTIR spectrum and a high-resolution diode laser spectrum were recorded. In the latter case, the precision estimated by the authors was better than 0.001 cm -1 . Some of the papers published on this subject dealt with CO 2 lasers [3,4], because the use of 14 C 16 O 2 can broaden the working range of spectrometers using CO 2 lasers as source of radiation. However, this species is of particular interest for several other reasons. A spectroscopic study can provide physical information on a still poorly investigated isotopologue of CO 2 . This molecule is produced in nuclear power reactors and it is detected in order to monitor its environmental concentration. This is usually done through accelerated mass spectrometry [7]. Labrie and Reid [8]showed that radiocarbon dating can be performed by infrared spectroscopy. Also, 14 C measurements have been carried out on planetary atmospheres [9]. In this paper we report on the high-resolution investigation of the fundamental ν 3 band of 14 C 16 O 2 by comb-referenced cavity ring-down spectroscopy. The main advantages of this technique are: high sensitivity (CO 2 absorptions with line-strength as low as 5 × 10 -32 cm can be detected), high resolution (1 kHz) and absolute frequency measurements, which in turn implies that no calibration of the spectra is needed.

Experimental details

Cavity ring-down (CRD) spectroscopy at 4.5 µm was realized to detect thirtythree ro-vibrational transitions belonging to the ν 3 band of 14 C 16 O 2 molecule. The experimental set-up, described elsewhere [10][11][START_REF] Galli | Molecular gas sensing down to parts per quadrillion: Radiocarbon-dioxide optical detection[END_REF], uses a cw IR coherent radiation at 4.5 µm generated by non-linear frequency difference (DFG) on a periodically-poled MgO:LiNbO 3 crystal, placed inside a cw ring Ti:sapphire laser cavity [13]. Such a laser, single-mode operated by optical injection of an external cavity diode laser (ECDL), is the DFG "pump" source, whereas the DFG "signal" light is provided by an external Nd:YAG laser at 1064 nm. A continuous IR tuning between 3850 ÷ 4540 nm (2597 ÷ 2203 cm -1 ) is allowed thanks to the tuning range of the ECDL-Ti:sapphire laser between 834÷862 nm, with a maximum IR power of 25 mW. The IR source is frequency stabilized against an octave-spanning (500 ÷ 1100 nm) fs Ti:sapphire Optical Frequency Comb Synthesizer (OFCS), as "pump" and "signal" frequencies fall within the OFCS spectral operation range [14]. As a consequence, the absolute frequency ν i of the generated "idler" radiation is given by the following equation:

ν i = ν p -ν s = N p N s -1 ν s + ν µ (1)
where N p and N s are the integer orders of the closest comb modes to the "pump" and "signal" absolute frequencies, ν p and ν s respectively. As the Nd:YAG laser frequency drift is controlled against the nearest OFCS mode, the absolute frequency of the IR radiation is stable to about 6 × 10 -13 at 1 s and accurate to 10 -12 , both limited by the performance of the Rb/GPS-disciplined 10-MHz quartz reference oscillator of our OFCS. In addition, synthesized IR frequency scans with a Cs-frequency standard traceability can be performed by changing the microwave synthesizer frequency ν µ in Eq. 1. With this procedure a maximum scan of 450 MHz (a bit less than half the repetition rate of the OFCS) can be performed, wide enough to detect the whole absorption line profiles of the observed transitions for the molecular gas conditions of the present measurements. The generated "idler" IR radiation, sent out of the Ti:sapphire cavity thanks to dispersion at the λ p -Brewster-cut output facet of the non-linear crystal, is efficiently coupled into a 1-m-long high-finesse (F ≈ 11600) Fabry-Perot cavity (FPC) for CRD spectroscopy. When the FPC is loaded up to a threshold level, the IR light is quickly switched off-resonance by using an AOM placed on the Nd:YAG laser, and transmitted photons leaking out of the cavity are detected by a InSb photodiode. The CRD decay signal is acquired by a 18-bit digitizing oscilloscope with a sampling rate of 10 MS/s. The measured decay time, τ c = 1/γ c for empty-cavity conditions is about 12 µs (γ c ≈ 86 ms -1 ). Such a set-up is suitable for linear-absorption CRD spectroscopy [15] and saturated-absorption cavity ring-down (SCAR) spectroscopy [10], with sensitivity figures, expressed in terms of the minimum detectable absorption coefficient α min , of 10 -9 and 10 -10 cm -1 Hz -1/2 , respectively. Although line-strengths of the ro-vibrational transitions of of magnitude as other CO 2 isotopologues (about 10 -18 cm for the most intense), the very low natural abundance of such species (≈ 10 -12 ) needs high-sensitivity spectroscopy techniques, as CRD, to be detected. For the measurements presented here a pure 14 C 16 O 2 sample would be preferred, to get the highest possible signal without absorption interferences from other CO 2 lines. Unfortunately, it is very difficult to handle such a sample, taking into account also radioactivity safety reasons.

Only spectroscopy of a few orders of magnitude 14 C-enriched samples (with respect to natural abundance) is feasible, which has partly prevented the spectroscopic measurements of such molecule. Even in these enriched samples, ultra-sensitive techniques such as CRD or SCAR are required. As a price to pay, other CO 2 absorptions lying in the same spectral region are also detected, which are added as an interference background to be accurately evaluated (see fitting procedure in the "results and discussion" section for other details). For the measurements presented here, we used a 2200 times 14 C-enriched CO 2 gas. Fig. 1 shows the IR spectra of most part of the ro-vibrational transitions between the P(40) and R(30) lines of 14 C 16 O 2 ν 3 fundamental band. Lines from the P(30) to R(30) were recorded by linear-absorption CRD spectroscopy with a total CO 2 pressure of 12 mbar. The FPC was dry-ice cooled in this case to get a gas temperature of 195 K, thus minimizing the interference effects from hot-bands of other CO 2 isotopologues. Nevertheless, strong interferences are detected, as it can be seen in the R-side of Fig. 1.

Lines from the P(40) to P(32) were recorded by SCAR spectroscopy with a total CO 2 pressure of 18 mbar1 . For these transitions, room temperature (296 K) was used to increase the ground state population, and hence the absorption signal. At the same time, the interfering background is increased, but it is still not larger than the 14 C 16 O 2 transition to be measured. In addition, the SCAR spectroscopic technique independently measures, in the same decay event, both γ g and γ c . The former decay rate γ g = cα g contains only the saturable part of the gas absorption. The latter includes all other contributions: cavity mirror losses, technical fringes, absorption from non-saturable transitions interfering with the target line.

We plot in Fig. 1 α g values (left vertical scale) for these SCAR lines, all having a negligible background except P(40), at which frequency the achievable saturation parameter was too low. On the contrary, for the CRD spectra, α = (γ g + γ c )/c was measured (right vertical scale in Fig. 1, where a constant α c = 2.87 × 10 -6 cm -1 due to the empty FPC losses was removed), and residual interfering background is present. In fact, the P(10), R(24) and R(28) transitions were not recorded because of the absorbing background.

For all the lines shown in Fig. 1, the acquisition procedure was the same independently of the CRD or SCAR technique used. For each IR absolute frequency, 128 consecutive decay events are averaged, and fitted to the expected CRD or SCAR decay function to obtain as a result α or α g , respectively. Ten consecutive α (or α g ) values are measured before a ±10-MHz IR frequency step is applied and the FPC length adjusted to retrieve the resonance condition. Each full 400 MHz scan is performed by stepping the frequency once forward (+ sign) and once backward (-sign). After each scan, the α (or α g ) plotted value at each scan frequency is the weighted mean of the 20 measured values (ten + and ten -). This procedure takes about 5 min for each complete scan. In the present gas conditions, linewidths are limited by Gaussian thermal broadening (FWHM 98 MHz and 120 MHz for 195 K and 296 K respectively) and Lorentzian pressure broadening [16]. Hence, a total 400 MHz scan is adequate to record a whole line profile. We note finally that in Fig. 1 all spectra were arranged together, one close to the other in a relative frequency scale, for the sake of clarity. For the analyzed spectra, the frequency scale is absolute and reproducible against the OFCS, as described above.

Results and Discussion

Only thirty-three transitions, plotted in Fig. 1, were detected in our experiments. The main limitations to the detection of transitions with higher values of J arose from the operating frequency range of the spectrometer, the weakness of the high-J transitions, even at room temperature, as well as the overlap with lines of more abundant isotopologues. The central frequency of each transition was measured by fitting the line profile to the expected Voigt profile, plus absorbing background. The fit residuals, both for CRD and SCAR measurements, are flat within noise, with no evidence of collisional narrowing [17]. Initial parameters for frequencies and intensities of the interfering transitions were obtained from the HITRAN database [16]. Unfortunately, the uncertainty of these data, if compared with the frequency precision and sensitivity of our measurements, and the presence of CO 2 transitions not listed by the database made difficult, and sometimes impossible, to evaluate the absorbing background. This effect is more important for CRD spectra, and the P(22), P( 16), R(0) and R(10) transitions are an example of that; in addition, background due to technical fringes was also considered. For the Voigt profile of the 14 C 16 O 2 transitions, Gaussian linewidth was always fixed to the required value by each temperature condition. The other line parameters, area, Lorentzian linewidth, and center frequency were left free. In Fig. 1, the agreement between experimental data and fit function (solid line trace) is evident in both CRD and SCAR experimental methodologies.

For each transition, several fit attempts were done with different absorbing background parameters, either fixing or not the Lorentzian linewidth, in order to get the influence of such change in the final line center uncertainty. We have also checked for possible systematic errors due to use of two different recording methods (CRD and SCAR spectroscopy) by using the P(26) transition as a trial test. The good agreement between different recordings is evident in Fig. 1, where this line is plotted in the two cases. In particular, the P(26) measured by CRD is affected only by technical fringes and background offset (no interfering lines), while the SCAR measurement exhibits no background at all. It is worth noting that this is one of the main advantages intrinsic to the SCAR technique. The different spectral areas are simply due to different temperature conditions. The results, including final standard errors are listed in Table 1. The line center uncertainty from the fit is several times lower (one order of magnitude for some transitions) than the final error reported in Table 1. However, we believe that the background parameters' uncertainties derived from the line profile fits were underestimated. We thought that a more realistic approach was to correct the fits uncertainties by individual factors, based on the goodness of the profile and the influence of contributions such as overlaps with other transitions and interference fringes.

Systematic corrections and uncertainties were considered in the transition frequencies reported in Table 1. Concerning the pressure shift, there are no literature data for transitions of this molecule, and we considered the contribution to be equal to that of the transitions of the ν 3 band of 13 C 16 O 2 [16]. In addition, we did not consider the J dependence of such shift, but we preferred to apply a more conservative procedure. A mean pressure shift correction value was applied for all measured transitions, as a result of the average of the tabulated values for the correspond- ing transitions of 13 C 16 O 2 . The uncertainties of these corrections are given by the standard deviation of such averages and it was added to the estimated line center fit error. As a result, the CRD measurements were blue-shifted by 0.9 MHz and their uncertainty increased by 0.5 MHz, while the SCAR measurements were blueshifted by 1.3 MHz and their uncertainty increased by 0.7 MHz. Finally, frequency calibration was not required, thanks to the primary Cs-standard traceability of our measurements. Also, the uncertainty due to the OFCS accuracy was not considered because, even in the worst case (67 Hz for R(30) transition), it is negligible with respect to other contributions. This ν 3 fundamental infrared band of CO 2 is a Σ -Σ parallel band with the usual P-and R-branch structure and lacking a Q-branch. The data were analysed using the basic Hamiltonian of a linear molecule, with centrifugal distortion corrections up to the quartic term for the ground state and to the octic term for the upper vibrational state. The ro-vibrational energy term values for non-degenerate vibrational modes of a linear molecule are given by

T (v, J) = G v + B v J(J + 1) -D v [J(J + 1)] 2 + H v [J(J + 1)] 3 + L v [J(J + 1)] 4 (2)
The band centre is defined as

ν 0 = G v -G v (3)
where G v and G v represent the purely vibrational energy of the upper and lower state, respectively. The parameters for the upper state v 3 = 1, B 1 , D 1 , H 1 and L 1 , were held fixed at the very accurate values of Ref. [4] while the ground state constants B 0 and D 0 were fitted, together with the band centre ν 0 . The results of the fit are shown in Table 2. In the fitting procedure the lines were given a weight according to the uncertainties reported in Table 1. The tolerance for rejection was chosen line by line and fixed at twice the standard error associated to the lines. Only four lines were rejected from the fit. The parameters obtained from the fit are well determined and are considerably improved in terms of both accuracy and precision with respect to the best constants available in literature [6]. The great precision obtained in this work for the band centre, together with the very precise constants of Ref. [4] allowed the determination of the vibrational energy of the levels of the Fermi dyad involved in the CO 2 laser excitation mechanism. In Ref. [4], they reported ν(001 -I) = 2.5965917618(13) × 10 7 MHz and ν(001 -II) = 2.94600023911(49) × 10 7 MHz. We eventually derived ν(I) = 4.076192962(51) × 10 7 MHz and ν(I) = 3.726784485(51) × 10 7 MHz.

Conclusions

The detection of the ν 3 fundamental band of 14 C 16 O 2 has been performed with a OFCS-referenced cavity ring-down spectrometer. Several ro-vibrational transitions ranging from P(40) to R(30) have been detected and analyzed. The band origin ν 0 and the ground state parameters B 0 and D 0 have been determined with accuracy improved by about one order of magnitude. This experiment shows the potential of this technique in terms of accuracy and sensitivity and provides important information on a rare isotopologue of CO

P 3 0 R 3 0 R 1 4 R 1 2 R 1 0 R 8 R 4 R 2 R 2 6 R 2 2 R 1 8 R 2 0 R 1 6 R 0 α -2 8 . 7 ( x 1 0 -7 c m -1 ) R 6 0 1 2 3 4 1 0 0 M H z / d i v P 2 6
Figure 1. Overview of the infrared spectrum of 14 C 16 O 2 molecule in the ν 3 band region. P(40)-P(32) lines were recorded by SCAR spectroscopy, while P(30)-R(30) lines were recorded by CRD spectroscopy. The fit of all measured transitions to a Voigt profile, that takes into account also the background from other CO 2 lines, is also shown. The spectra are plotted in a relative horizontal frequency scale, and frequency gaps between transitions are not shown for the sake of clarity. The right side vertical scale applies to CRD spectra and the left side vertical scale to SCAR spectra. An empty cavity loss rate of αc = 2.87×10 -6 cm -1 was removed from the measured rates in the case of CRD spectra to get a comparable scale with the SCAR ones. The CO 2 pressure and temperature are 18 mbar and 296 K for the SCAR spectra, 12 mbar and 195 K for the CRD ones. The acquisition conditions and fit procedure are described in the text. ) P

P 3 0 R 3 0 R 1 4 R 1 2 R 1 0 R 8 R 4 R 2 R 2 6 R 2 2

The ν 3

 3 band of 14 C 16 O 2 molecule

  14 C 16 O 2 ν 3 band are expected to be of the same order The ν 3 band of 14 C 16 O 2 molecule

	July 28, 2011 15:33 Molecular Physics 14CO2pap˙rev
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 1 .

 1 2 .The ν 3 band of 14 C 16 O 2 molecule 7 Transitions frequencies, uncertainties and residuals from the fit of the lines measured for the ν3 fundamental band of 14 C 16 O2. All numbers are in MHz.

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

a See text for a discussion about the uncertainty of the lines.

b Lines marked with an asterisk exceeded the tolerance for rejection and were excluded from the fit.

Table 2 .

 2 Spectroscopic parameters (in cm -1 ) obtained for the ν3 fundamental band of 14 C 16 O2. The quantity in parentheses is one standard deviation in units of the last quoted digit.The ν 3 band of 14 C 16 O 2 molecule 9

	Parameter	This work a	Ref. [6]
	ν 0	2225.801399(17)	2225.80239(16)
	B 0	0.390253082(58)	0.39025488(18)
	D 0 × 10 -7	1.33144(37)	1.3372(20)
	B 1	0.38739025628 b
	D 1 × 10 -7	1.327523 b
	H 1 × 10 -14	-3.08 b
	L 1 × 10 -17	1.244 b
	σ F IT × 10 -5	3.1
	b Fixed at the values of Ref. [4]

a

With respect to Ref.[10], valid in the inhomogeneous broadening regime, the theoretical model to fit the SCAR signals was adjusted to work at these pressure conditions[START_REF] Galli | Molecular gas sensing down to parts per quadrillion: Radiocarbon-dioxide optical detection[END_REF].

Acknowledgements

This work was partially supported by Ente Cassa di Risparmio di Firenze and by Regione Toscana through the projects CTOTUS and SIMPAS, in the framework of POR-CReO FESR 2007-2013.