Nikolaos Stefanakis

Jacques Marchal

Valentin Emiya

Nancy Bertin

Rémi Gribonval

Pierre Cervenka

SPARSE UNDERWATER ACOUSTIC IMAGING: A CASE STUDY

Keywords: Sonar imaging, underwater acoustics, beamforming, orthogonal matching pursuit, sparse reconstruction, compressed sensing

Underwater acoustic imaging is traditionally performed with beamforming: beams are formed at emission to insonify limited angular regions; beams are (synthetically) formed at reception to form the image. We propose to exploit a natural sparsity prior to perform 3D underwater imaging using a newly built flexible-configuration sonar device. The computational challenges raised by the highdimensionality of the problem are highlighted, and we describe a strategy to overcome them. As a proof of concept, the proposed approach is used on real data acquired with the new sonar to obtain an image of an underwater target. We discuss the merits of the obtained image in comparison with standard beamforming, as well as the main challenges lying ahead, and the bottlenecks that will need to be solved before sparse methods can be fully exploited in the context of underwater compressed 3D sonar imaging.

INTRODUCTION

The use of sparsity prior have proven its potential and efficiency, as well as its versatility in a wide range of signal processing contexts. However, scaling up from synthetic experiments to real world applications remains a challenge in many cases. Representing a signal in a redundant dictionary implies handling huge volumes of data, and a discrete model in which the considered analog data fits. This paper addresses this difficulty in the case of underwater acoustic imaging.

In underwater acoustic scenarii, seafloor imaging and bathymetry are usually performed by sidescan sonars or multibeam echosounders [START_REF] Lurton | An introduction to underwater acoustics : principles and applications[END_REF]. These classical sonar systems insonify an across-track area at each emission burst: 3D images are built by juxtaposing successive strips, obtained one after the other while the ship is moving ahead. Real 3D-imaging systems are not available, especially due to the large amount of data to acquire and to process.

In this context, the use of a sparse modelling of the problem and of compressed sensing techniques [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF][START_REF] David | Compressed sensing[END_REF][START_REF] Baraniuk | Compressive sensing[END_REF] opens interesting doors to new devices that could image directly the 3D space. Such sparse models have been proposed for the related radar imaging problem [START_REF] Herman | High-resolution radar via compressed sensing[END_REF], and for over-the-air sonar imaging [START_REF] Boufounos | Compressive sensing for over-the-air ultrasound[END_REF]. However, existing work essentially tested ideas on simple synthetic data. Dealing with real data implies challenges such as: handling a 3D grid with a higher number of points; detecting targets that are not located on the grid points; detecting complex-shape objects rather than a simple This work was supported by Agence Nationale de la Recherche (ANR), project ECHANGE (ANR-08-EMER-006) and project ASAP (ANR-09-EMER-001). pattern like a square; using non-ideal transducers for emission and reception, with directivity patterns and calibration issues; and handling phase issues due to the propagation and to the modulation by a carrier frequency.

In this paper, we perform underwater acoustics imaging with a new flexible-configuration sonar device that we built and used to record real data. It leads to a problem statement, described in Section 2, that includes the formulation of the theoretical estimation problem as well as the issues raised by the actual orders of magnitudes of the considered dimensions. In Section 3, we describe the standard beamforming approach as a linear estimator and introduce an alternate linear estimator where emission beams are replaced by elementary emission patterns. In Section 4, we then propose a nonlinear estimator based on a sparse model and on a modified version of the Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] algorithm, in which algorithmic enhancements makes the estimation tractable. As a proof of concept, the proposed approach provides an image of the target comparable to the basic beamforming technology, as shown in Section 5. Though not yet competing with refined beamformingbased multibeam echosounders, the results show the potential of our approach. We discuss this first applicative and algorithmic milestone towards new compressive-sampling-inspired 3D-imaging sonar systems and draw perspectives in Section 6.

PROBLEM STATEMENT

Specifications of the new sonar prototype

The proposed sonar prototype is made of two 26 cm long linear arrays of 64 transducers each. The antennas work at a carrier frequency Fe = 480 kHz and with a bandwidth FB = 160 kHz. Each elementary transducer can operate either as a receiver or as an emitter; nevertheless, it is not possible to interchange the transducers utility during the imaging process. A computer drives both the emission of NE = 64 amplified channels, and the reception of NR = 64 others channels. The system can emit several (NP ) successive transmission sequences called pings. At each ping, the NE transmitted signals can be different from each other in order to form an emission beam or in order to achieve a complex radiation pattern.

These specifications allow to carry out experiments in a tank or in the sea with a range up to a few tens of meters. Due to its flexibility, sidescan sonar, multibeam systems or original configurations (such as acoustical camera) can be simulated.

Description of the experiment

The underwater imaging experiments take place in a tank of dimensions 2.9 m × 0.9 m × 0.5 m filled with fresh water. The sound speed c is constant and equal to 1484 m/s; thus, the wavelength for Fig. 1. Picture of the target (left) and its relative position to the coordinate system and to the antenna (right). Fe = 480 kHz is λ ≈ 3 mm. The target to be imaged is a round wheel with outer diameter of 52 cm made of plywood. The wheel is partially covered with sand and is lying on the surface of the water in the tank. For the experiment performed in this paper, the transmitting array is placed along the x-axis and the receiving array along the y-axis (Fig 1). The elementary signal that drives the emitters is a frequency-modulated window of short duration T = 20 µs with 10 periods of the carrier frequency Fe:

e(t) = e (j2πFet) e (-π 2 F 2 B t 2 /2) rect(t/T ). (1)

Direct model formulation

The 3D acoustic scene to be sensed is discretized into M small volumes called voxels with centers at locations p k , k = 1, ..., M . The NE emitters are indexed by i and the NR receivers are indexed by j.

Positions of the emitters (resp. receivers) are denoted by vectors pi (resp. pj). The signal emitted by transducer i at ping p is denoted ei,p(t). The signal recorded at the receiver j at ping p is given by

mj,p (t) = N k k=1 b k N E i=1 eip (t -τ ik -τ kj ) (2)
where τ ik

p k -p i 2 c
and τ kj

p k -p j 2 c
are the forward and backward propagation delays from the elements of the two antennas to the k-th voxel and b k is the backscattering coefficient of the voxel. Note that only the times of flight are taken into account into the model, and not the spreading attenuation term. The signal mj,p(t) recorded at the j th receiver at ping p is eventually sampled as mj,p[n] = mj,p(nTs), with Fs = 1/Ts = 2 M Hz the sampling frequency.

Inverse problem formulation

The 3D acoustic imaging problem is the estimation of the relative energy backscattered from each voxel k, given the NP × NR recorded signals: in other words, given the discrete time-series mj,p[n] measured at all receivers for all pings, we wish to estimate the backscattering coefficients b k . This is an inverse problem that can be expressed in more compact form:

m = Φb, (3)
using the following conventions. Discrete time column vectors of length NT associated to continuous time signals x(t) are represented by small bold letters x. For a given delay τ we denote x τ [n]

x (nTs -τ ). The dot product is denoted by x, y x H .y with (•) H the Hermitian adjoint. The following compact notations are used, with (•) T the transpose:

m T j m T j,1 , . . . , m T j,N P , m T m T 1 , . . . , m T N R , (4) b T b T 1 , . . . , b T M .
(5)

a j,p,k N E i=1 e τ ikj i,p , τ ikj τ ik + τ kj (6) a T j,k a T j,1,k , . . . , a T j,N P ,k , (7)
a T k a T 1,k , . . . , a T N R ,k , (8)
Φ [a1, . . . , aM ] .

In other words, m ∈ C N P N R N T is the full measurement, b ∈ C M is the vector with the unknown backscattering coefficients and Φ ∈ C N P N R N T ×M is made with delayed versions of the emitted signals.

STANDARD BEAMFORMING APPROACH

Beamforming (BF) is a standard imaging approach to estimate the backscattering coefficients b as a linear combination of the measurements, for specific configurations of the emission [START_REF] Lurton | An introduction to underwater acoustics : principles and applications[END_REF].

Emission

In the classical approach, at each ping p, all emitters are simultaneously active. Introducing small delays between them allow to concentrate the emitted energy in a specific 2D-slice by forming a beam.

Illuminating only a region of interest shall improve the signal-tonoise ratio (SNR). Moreover, it reduces the computational cost, by dividing the 3D imaging problem into several independent, smaller 2D subproblems. However, by nature, this strategy does not take into account across-slice structure.

We target a more flexible system, able to deal with 3D and to adopt various configurations. At each ping, we choose to have only one active emittor (mono-emission), which insonifies then the whole half-space below the ship. Thanks to linearity, such a configuration can be used to simulate the previous one, but also a wider variety of schemes. As the insonification is global, information about all the scene to image is carried at each ping, opening to a better scalability and versatility in the number of pings and acquisition time, even though at the price of a possible lower SNR.

Imaging

In the classical beamforming approach, images are formed by scanning the space. Reception beams, perpendicular to the emission beams, are synthetically formed for each ping. The imaging process basically amounts to a linear estimate, similar to matched filtering b Φ H m. One of the main drawbacks of this method is the appearance, in the reconstructed image, of energetic artifacts around bright points (e.g. points of the wheel where a specular reflexion occur). These ghost reflectors are due to the sidelobes in the antenna directivity, or, in our formalism, as due to the non-zero off-diagonal terms in the matrix Φ H Φ. Any linear estimator of b is likely to suffer from these sidelobes.

Classical beamforming neglects an important feature of the scene to image: its sparsity in the 3D space. We will show in section 4 how this sparsity can be taken into account, and used to build an image from the mono-emission sonar measurements. Being able to deal with 3D-insonification and a sparse model will then open many doors to new acquisition strategies.

OMP-LIKE RECONSTRUCTION ALGORITHM

We introduce our approach by noticing that the coefficients b k are expected to be zero for most of the voxels (voxel containing only water). Exploiting this sparse nature, we can formulate a non-linear estimator from known greedy algorithms that have been developed for generic sparse recovery problems, such as Orthogonal Matching Pursuit [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] (OMP, see Algorithm 1). We propose algorithmic enhancements for the configuration where only one emittor is active at a time (eip δipe).

Discretization and dimensionality issues

Discretizing the full tank with a step of λ would results in a M = 48.10 6 -point grid, and, in time, the measurement length is about NP NRNT = 13.10 6 samples. This makes the problem obviously intractable. Before detailing an algorithmic proposal, we alleviate this dimensionality burden in two ways.

We focus on the region of interest defined by 0.6 m ≤ r ≤ 1.2 m, -1.4 • ≤ θ ≤ 6.4 • , -4.7 ≤ φ ≤ 34.4 • . This area is depicted in the right part of Fig. 1 using dashed lines. A tradeoff between resolution and complexity must be found. We have chosen a discretization using 61 uniform radial steps (∆r ≈ 9.9 mm), 16 uniform angular steps along θ (∆θ ≈ 0.52 • ), and 72 uniform angular steps along φ (∆φ ≈ 0.54 • ). In the region of interest, each voxel is approximately a cube with edge length of 1 cm ≈ 3λ. As shown below, scanning this region remains a large-size problem.

Signals can be downsampled at F s = Fs/5 to reduce their size according to their bandwidth. Then, they are cropped according to the radial limitation of the region of interest. As a result, the downsampled cropped measurements are NT = 324-sample series for each of the NR = 64 receivers and of the NP = 64 pings.

The problem formulation ( 3) is now realistic, with NP NRNT = 1, 327, 104 measures, and M = 61 × 16 × 72 = 70, 272 voxels.

Reconstruction algorithm

Assuming that b is sparse, it can be estimated by many sparsitybased methods. Here, we propose to use an approach based on the OMP algorithm, which has a relatively affordable computational cost to face a high-size problem. Algorithm 1 details a naive implementation of OMP. Several days of computations and a hundred gigabytes of memory would be necessary to run such a program. The memory issue comes from the size of the dictionary. As a result, each atom must be built when needed and then forgotten. At each iteration, two stages are computationally demanding: the atom selection and the sparse representation update. The complexity of the latter is cubic and can be reduced -as proposed in [START_REF] Blumensath | Gradient pursuits[END_REF] -to a quadratic complexity by using a block matrix inversion of Φ * Ω ΦΩ, where the results from the previous iteration is used. The complexity is approximately reduced by a factor Kmax (about one thousand in our experiments). Reducing the complexity of the atom selection stage remains the bottleneck addressed in the next section.

Fast atom selection

The cost for reconstructing all atoms and the inner products with the residual at each iteration is about O (M NT NRNP ), where M is the total number of atoms. We propose two strategies to reduce this cost. We will denote by τikj [τ ikj /F s ] F s the approximation of τikj using the sampling frequency F s . We denote the part of r related to ping p and receiver j by rp,j.

Algorithm 1 OMP Initialisation: r ← m; Ω ← ∅. for K = 1 to Kmax do Atom selection: k ← arg max k | a k , r | Ω ← Ω ∪ k Sparse representation update: bΩ ← Φ + Ω m Residue update: r ← m -ΦΩ bΩ end for Output: b OMP ← bΩ.
Atom selection based on FFT-filtering: the first strategy relies on the structure of the atoms, which are only formed by delayed versions of the elementary emission pattern e (see eq. ( 6)-( 8)). The strategy is close to [START_REF] Krstulovic | MPTK: Matching pursuit made tractable[END_REF]:

1. first, compute the auxiliary result e * r using an FFT filtering technique, the resulting complexity being O (NT log Ne), where Ne NT is the length of e when sampled at the same sampling frequency as r;

2. then, approximate τ ikj by τikj and according to eq. ( 6), use

a k , r ≈ N R j=1 N E
i=1 [e * rp=i,j](τ ikj )e 2jπFe(τ ikj -τ ikj ) , with a complexity of O (M NRNE).

Approximate atom selection: in the second strategy, we approximate the contribution of rp,j to the inner product by the single, most significant value rp=i,j(τ ikj ) so that a k , r is approximated by

N R j=1 N E i=1 rp=i,j(τ ikj )e 2iπFe(τ ikj -τ ikj ) , (10)
the complexity being O (M NENR).

IMPLEMENTATION AND RESULTS

A significant variation in the level of the recorded signals was observed, which revealed a calibration problem at the arrays. In particular, the acquisition gains of some transducers were below -10 dB compared to an estimation of the average acquisition gain. In order to circumvent this issue, 12 receivers were completely discarded when applying the sparse reconstruction approach. Nevertheless, these transducers were not removed during the classical imaging approach; amplitude mismatch doesn't adversely affect the classical approach, as long as there is no phase mismatch. Both OMP algorithms were performed for Kmax = 1200 iterations. The computation time was about 8 hours for both selection tactics. Each algorithm provided a residual r and a sparse set of coefficients b k , k ∈ Ω. An image was built in a hybrid way representing both the sparse coefficient values and the magnitude of the residual projected on the non-selected atoms:

I k |b| k , k ∈ Ω | a k , r | , k ∈ Ω c . ( 11
)
The resulting figures are shown in grey scale in Fig. 2. Focusing on column (a) of Fig. 2, it can be observed that there are regions with very large amplitudes at x = 0.6 and x = 1.1 m. This is due to specular echos, which result to the appearance of side lobes, i.e. spatial spreading along φ and θ. This can lead to misintepretation of the object shape. While the main shape of the scatter object is more or less apparent with all approaches, it can be seen that the To be noticed is the fact that after 1200 iterations, the ratio between the initial energy of the measurement to the energy of the residual, calculated as 20 log 10 m r 2 , was only +0.6 dB for OMP algorithms. This rather small convergence is believed to be a consequence of the calibration issues mentioned earlier and possibly, of modeling errors concerning the relevant position of the two antennas or the calculation of the sound speed. Given a better modeling or estimation of these parameters, it is believed that the sparse formulation will further enhance the imaging quality compared to the standard approach.

CONCLUSION AND PERSPECTIVES

In this paper, we have applied sparse reconstruction to an underwater acoustic imaging problem. We highlighted the very large amount of data it implies to deal with, the issues it raises, and showed that this burden can be overcome in our case by applying an optimized version of OMP algorithm. Although under conditions of poor signalto-noise ratio and with serious calibration issues concerning the receiving antenna, the sparse reconstruction approach proved to be a viable option to produce images of the target, and even to bring some improvement in certain aspects of the image quality, compared to standard imaging. Though perfectible, the sparse approach shows here its potential on real world data.

Besides possible refinements of the sparse model itself and improvement of the computation speed, some other deeper changes must be considered. One of the major interest of modelling the imaging problem as a sparse inverse problem, is the possibility to exploit the compressive sampling framework. Some major changes in the acquisition schemes can be performed: for instance, introduction of randomness in the emitted sequences, and different distributions of the transducers utilities (for instance, placing receivers and emitters on both antennas, with irregular distribution, which is allowed by our flexible sonar). This framework lets foresee dramatic improvements in terms of acquisition time and data flow, which is a major issue in underwater acoustics.

Fig. 2 .

 2 Fig. 2. Images obtained by classical imaging (a), by OMP with approximate selection (b), by OMP with FFT based selection step (c). The selected voxels, as obtained by the last method, are shown with respective amplitudes (d). Upper row corresponds to the slice θ = -1.4 o and lower row to slice at θ = 3 o . All figures are displayed at the same dynamic range and the colour scale is normalized with respect to the maximum voxel amplitude. The x axis corresponds to range and the y axis to the y coordinate


Jean B Lasserre
email: lasserre@laas.fr

NEW APPROXIMATIONS FOR THE CONE OF COPOSITIVE MATRICES AND ITS DUAL

Keywords: 1991 Mathematics Subject Classification. 15B48 90C22 copositive matrices, completely positive matrices, semidefinite relaxations

We provide convergent hierarchies for the convex cone C of copositive matrices and its dual C * , the cone of completely positive matrices. In both cases the corresponding hierarchy consists of nested spectrahedra and provide outer (resp. inner) approximations for C (resp. for its dual C * ), thus complementing previous inner (resp. outer) approximations for C (for C * ). In particular, both inner and outer approximations have a very simple interpretation. Finally, extension to K-copositivity and K-complete positivity for a closed convex cone K, is straightforward.

Introduction

In recent years the convex cone C of copositive matrices and its dual cone C * of completely positive matrices have attracted a lot of attention, in part because several interesting NP-hard problems can be modelled as convex conic optimization problems over those cones. For a survey of results and discussion on C and its dual, the interested reader is referred to e.g. Anstreicher and Burer [START_REF] Anstreicher | Computable representations for convex hulls of low-dimensional quadratic forms[END_REF], Bomze [START_REF] Bomze | Copositive optimization -recent developments and applications[END_REF], Bomze et al. [START_REF] Bomze | Think co(mpletely) positive! -matrix properties, examples and a clustered bibliography on copositive optimization[END_REF], Burer [START_REF] Burer | Copositive programming, in Handbook of Semidefinite, Conic and Polynomial Optimization[END_REF], Dür [START_REF] Dür | Copositive Programming -a survey[END_REF] and Hiriart-Urruty and Seeger [START_REF] Hiriart-Urruty | A variational approach to copositive matrices[END_REF].. As optimizing over C (or its dual) is in general difficult, a typical approach is to optimize over simpler and more tractable cones. In particular, nested hierarchies of tractable convex cones C k , k ∈ N, that provide inner approximations of C have been proposed, notably by Parrilo [START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF], deKlerk and Pasechnik [START_REF] De Klerk | Approximation of the stability number of a graph via copositive programming[END_REF], Bomze and deKlerk [START_REF] Bomze | Solving standard quadratic optimization problems via linear, semidenite and copositive programming[END_REF], as well as Peña et al. [START_REF] Peña | Computing the stability number of a graph via linear and semidenite programming[END_REF]. For example, denoting by N (resp. S + ) the convex cone of nonnegative (resp. positive semidefinite) matrices, the first cone in the hierarchy of [START_REF] De Klerk | Approximation of the stability number of a graph via copositive programming[END_REF] is N , and N + S + in that of [START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF], whereas the hierarchy of [START_REF] Peña | Computing the stability number of a graph via linear and semidenite programming[END_REF] is in sandwich between that of [START_REF] De Klerk | Approximation of the stability number of a graph via copositive programming[END_REF] and [START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF]. Of course, to each such hierarchy of inner approximations (C k ), k ∈ N, of C, one may associate the hierarchy (C * k ), k ∈ N, of dual cones which provides outer approximations of C * . However, quoting Dür in [START_REF] Dür | Copositive Programming -a survey[END_REF]: "We are not aware of comparable approximation schemes that approximate the completely positive cone (i.e. C * ) from the interior."

Let us also mention Gaddum's characterization of copositive matrices described in Gaddum [START_REF] Gaddum | Linear inequalities and quadratic forms[END_REF] for which the set membership problem "A ∈ C" (with A known) reduces to solving a linear programming (LP) problem (whose size is not polynomially bounded in the input size of the matrix A). However, Gaddum's characterization is not convex in the coefficients of the matrix A, and so cannot be used to provide a hierarchy of outer approximations of C. On the other hand, de Klerk and Pasechnik [START_REF] De Klerk | A linear programming formulation of the standard quadratic optimization problem[END_REF] have used Gaddum's characterization to help solve quadratic optimization problems on the simplex via solving a hierarchy of LP-relaxations of increasing size and with finite convergence of the process.

Contribution. The contribution of this note is precisely to describe an explicit hierarchy of tractable convex cones that provide outer approximations of C which converge monotonically and asymptotically to C. And so, by duality, the corresponding hierarchy of dual cones provides inner approximations of C * converging to C * , answering Dür's question and also showing that C and C * can be sandwiched between two converging hierarchies of tractable convex cones. The present outer approximations are not polyhedral and hence are different from the outer polyhedral approximations O n r of C defined in Yildirim [START_REF] Yildirim | On the accuracy of uniform polyhedral approximations of the copositive cone[END_REF]. O n r are based on a certain discretization ∆(n, r) of the simplex ∆ n := {x ∈ R n + : e T x = 1} whose size is parametrized by r. All matrices associated with quadratic forms nonnegative on ∆(n, r) belong to O r n for all r, and C = ∩ ∞ r=1 O r n ; Combining with the inner approximations of de Klerk and Pasechnik [START_REF] De Klerk | Approximation of the stability number of a graph via copositive programming[END_REF], the author provides tight bounds on the gap between upper and lower bounds for quadratic optimization problems; for more details the interested reader is referred to Yildirim [START_REF] Yildirim | On the accuracy of uniform polyhedral approximations of the copositive cone[END_REF].

In fact, our result is a specialization of a more general result of [START_REF] Lasserre | A new look at nonnegativity and polynomial optimization[END_REF] about nonnegativity of polynomials on a closed set K ⊂ R n , in the specific context of quadratic forms and K = R n + . In such a context, and identifying copositive matrices with quadratic forms nonnegative on K = R n + , this specialization yields inner approximations for C * with a very simple interpretation directly related to the definition of C * , which might be of interest for the community interested in C and C * but might be "lost" in the general case treated in [START_REF] Lasserre | A new look at nonnegativity and polynomial optimization[END_REF], whence the present note. Finally, following Burer [START_REF] Burer | Copositive programming, in Handbook of Semidefinite, Conic and Polynomial Optimization[END_REF], and K ⊂ R n being a closed convex cone, one may also consider the convex cone C K of K-copositive matrices, i.e., real symmetric matrices A such that x T Ax ≥ 0 on K, and its dual cone C * K of K-completely positive matrices. Then the outer approximations previously defined have an immediate and straightforward analogue (as well as the inner approximations of the dual).

Main result

2.1. Notation and definitions. Let R[x] be the ring of polynomials in the variables x = (x 1 , . . . , x n ). Denote by R[x] d ⊂ R[x] the vector space of polynomials of degree at most d, which forms a vector space of dimension s(d) = n+d d , with e.g., the usual canonical basis (x α ) of monomials. Also, let

N n d := {α ∈ N n : i α i ≤ d} and denote by Σ[x] ⊂ R[x] (resp. Σ[x] d ⊂ R[x] 2d
) the space of sums of squares (s.o.s.) polynomials (resp. s.o.s. polynomials of degree at most 2d).

If f ∈ R[x] d , write f (x) = α∈N n d f α x α in
the canonical basis and denote by f = (f α ) ∈ R s(d) its vector of coefficients. Finally, let S n denote the space of n × n real symmetric matrices, with inner product A, B = trace AB, and where the notation A 0 (resp. A ≻ 0) stands for A is positive semidefinite.

Given Moment matrix. With a sequence y = (y α ), α ∈ N n , let L y : R[x] → R be the linear functional

K ⊆ R n , denote by cl K (resp. conv K) the closure (resp. the convex hull) of K. Recall that given a convex cone K ⊆ R n , the convex cone K * = {y ∈ R n : y, x ≥ 0 ∀x ∈ K} is called the dual cone of K, and satisfies (K * ) * = cl K . Moreover, given two convex cones K 1 , K 2 ⊆ R n , K * 1 ∩ K * 2 = (K 1 + K 2 ) * = (K 1 ∪ K 2 ) * (K 1 ∩ K 2 ) * = cl (K * 1 + K * 2 ) = cl (conv (K * 1 ∪ K * 2
h (= α h α x α ) → L y (h) = α h α y α , h ∈ R[x].
With d ∈ N, let M d (y) be the symmetric matrix with rows and columns indexed in N n d , and defined by: (2.1)

M d (y)(α, β) := L y (x α+β ) = y α+β , (α, β) ∈ N n d × N n d .
The matrix M d (y) is called the moment matrix associated with y, and it is straightforward to check that

L y (g 2 ) ≥ 0 ∀g ∈ R[x] ⇔ M d (y) 0, d = 0, 1, . . . . Localizing matrix. Similarly, with y = (y α ), α ∈ N n , and f ∈ R[x] written x → f (x) = γ∈N n f γ x γ ,
let M d (f y) be the symmetric matrix with rows and columns indexed in N n d , and defined by:

(2.2) M d (f y)(α, β) := L y f (x) x α+β = γ f γ y α+β+γ , (α, β) ∈ N n d × N n d .
The matrix M d (f y) is called the localizing matrix associated with y and

f ∈ R[x].
Observe that

(2.3) g, M d (f y)g = L y (g 2 f ), ∀g ∈ R[x] d ,
and so if y has a representing finite Borel measure µ, i.e., if

y α = R n x α dµ, ∀ α ∈ N n , then (2.3) reads (2.4) g, M d (f y)g = L y (g 2 f ) = R n g(x) 2 f (x) dµ(x), ∀g ∈ R[x] d .
Actually, the localizing matrix M d (f y) is nothing less than the moment matrix associated with the sequence z = f y = (z α ), α ∈ N n , with z α = γ f γ y α+γ . In particular, if f is nonnegative on the support supp µ of µ then the localizing matrix M d (f y) is just the moment matrix associated with the finite Borel measure dµ f = f dµ, absolutely continuous with respect to µ (denoted µ f ≪ µ), and with density f . For instance, with d = 1 one may write

M 1 (f y) = R n   1 | x T - - x | xx T   f (x) dµ(x) = R n   1 | x T - - x | xx T   dµ f (x), or, equivalently, (2.5) M 1 (f y) = mass(µ f ) ×   1 | E μf (x) T - - E μf (x) | E μf (xx T )   ,
where E μf (•) denotes the expectation operator associated with the normalization μf of µ f (as a probability measure), and E μf (xx T ) denotes the matrix of noncentral second-order moments of μf (and the covariance matrix of μf if E μf (x) = 0).

Main result. In [12, Theorem 3.2] the author has shown in a general context that a polynomial

f ∈ R[x] is nonnegative on a closed set K ⊆ R n if and only if (2.6) K g 2 f dµ ≥ 0, ∀g ∈ R[x],
where µ is a given finite Borel measure with support supp µ being exactly K; if K is compact then µ is arbitrary whereas if K is not compact then µ has to satisfy a certain growth condition on its moments. If y = (y α ), α ∈ N n , is the sequence of moments of µ then (2.6) is in turn equivalent to

M d (f y) 0, ∀d = 0, 1, . . . ,
where M d (f y) is the localizing matrix associated with f and y, defined in (2.2).

In this section we particularize this result to the case of copositive matrices viewed as homogeneous forms of degree 2, nonnegative on the closed set K = R n + .

So with A = (a ij ) ∈ S n , let denote by f A ∈ R[x] 2 the quadratic form x → x T Ax, and let µ be the joint probability measure associated with n i.i.d. exponential variates (with mean 1), with support supp µ = R n + , and with moments y = (y α ), α ∈ N n , given by:

(2.7) y α = R n + x α dµ(x) = R n + x α exp(- n i=1 x i ) dx = n i=1 α i !, ∀ α ∈ N n .
Recall that a matrix A ∈ S n is copositive if f A (x) ≥ 0 for all x ∈ R n + , and denote by C ⊂ S n the cone of copositive matrices, i.e., (2.8)

C := { A ∈ S n : f A (x) ≥ 0 ∀ x ∈ R n + }.
Its dual cone is the closed convex cone of completely positive, i.e., matrices of S n that can be written as the sum of finitely many rank-one matrices xx T , with x ∈ R n + , i.e., (2.9)

C * = conv { xx T : x ∈ R n + }.
Next, introduce the following sets C d ⊂ S n , d = 0, 1, . . ., defined by: (2.10)

C d := { A ∈ S n : M d (f A y) 0 }, d = 0, 1, . . .
where M d (f A y) is the localizing matrix defined in (2.2), associated with the quadratic form f A and the sequence y in (2.7).

Observe that in view of the definition (2.2) of the localizing matrix, the entries of the matrix M d (f A y) are homogeneous and linear in A. Therefore, the condition

M d (f A y)
0 is a homogeneous Linear Matrix Inequality (LMI) and defines a closed convex cone (in fact a spectrahedron) of S n (or equivalently, of R n(n+1)/2 ). Each C d ⊂ S n is a convex cone defined solely in terms of the entries (a ij ) of A ∈ S n , and the hierarchy of spectrahedra (C d ), d ∈ N, provides a nested sequence of outer approximations of C. Theorem 2.1. Let y be as in (2.7) and let C d ⊂ S n , d = 0, 1, . . ., be the hierarchy of convex cones defined in (2.10)

. Then C 0 ⊃ C 1 • • • ⊃ C d • • • ⊃ C and C = ∞ d=0 C d .
The proof is a direct consequence of [START_REF] Lasserre | A new look at nonnegativity and polynomial optimization[END_REF]Theorem 3.3] with K = R n + and f = f A . Since f A is homogeneous, alternatively one may use the probability measure ν uniformly supported on the n-dimensional simplex ∆ = {x ∈ R n + : i x i ≤ 1} and invoke [START_REF] Lasserre | A new look at nonnegativity and polynomial optimization[END_REF]Theorem 3.2]. The moments ỹ = (ỹ α ) of ν are also quite simple to obtain and read:

(2.11) ỹα = ∆ x α dx = α 1 ! • • • α n ! (n + i α i )! , ∀α ∈ N n .
(See e.g. Grundmann [START_REF] Grundmann | Invariant integration formulas for the n-simplex by combinatorial methods[END_REF].) Observe that the set membership problem "A ∈ C d ", i.e., testing whether a given matrix A ∈ S n belongs to C d , is an eigenvalue problem as one has to check whether the smallest eigenvalue of M d (f A y) is nonnegative. Therefore, instead of using standard packages for Linear Matrix Inequalities, one may use powerful specialized softwares for computing eigenvalues of real symmetric matrices.

We next describe an inner approximation of the convex cone C * via the hierarchy of convex cones (C * d ), d ∈ N, where each C * d is the dual cone of C d in Theorem 2.1. Recall that Σ[x] d is the space of polynomials that are sums of squares of polynomials of degree at most d. A matrix A ∈ S n is also identified with a vector a ∈ R n(n+1)/2 , and conversely, with any vector a ∈ R n(n+1)/2 is associated a matrix A ∈ S n . For instance, with n = 2, (2.12)

A = a b b c ↔ a =   a 2b c   .
So we will not distinguish between a convex cone in R n(n+1)/2 and the corresponding cone in S n .

Theorem 2.2. Let C d ⊂ S n be the convex cone defined in (2.10). Then

(2.13) C * d = cl ( X, M d (x i x j y) ) 1≤i≤j≤n : X ∈ S s(d) +
.

Equivalently:

C * d = cl      R n + xx T σ(x) dµ(x) dµσ(x) : σ ∈ Σ[x] d      = cl mass(µ σ )E μσ (xx T ) : σ ∈ Σ[x] d , (2.14)
with μσ and E μσ (xx T ) as in (2.5).

Proof. Let ∆ d := ( X, M d (x i x j y) ) 1≤i≤j≤n : X ∈ S s(d) + , so that ∆ * d =    a ∈ R n(n+1)/2 : 1≤i≤j≤n a ij X, M d (x i x j y ≥ 0 ∀ X ∈ S s(d) +    , =    a ∈ R n(n+1)/2 : X, M d   ( 1≤i≤j≤n a ij x i x j ) y   ≥ 0 ∀ X ∈ S s(d) +    , = {A ∈ S n : X, M d (f A y) ≥ 0 ∀ X ∈ S s(d) + } [with A, a as in (2.12)] = {A ∈ S n : M d (f A y) 0 } = C d .
And so we obtain the desired result

C * d = (∆ * d ) * = cl (∆ d ).
Next, writing the singular decomposition of X as s k=0 q k q T k for some s ∈ N and some vectors (q k ) ⊂ R s(d) , one obtains that for every 1

≤ i ≤ j ≤ n, X, M d (x i x j y) = s k=0 q k q T k , M d (x i x j y) = s k=0 q k , M d (x i x j y)q k = s k=0 R n + x i x j q k (x) 2 dµ(x) [by (2.4)] = R n + x i x j σ(x) dµ(x) dµσ(x)
,

where σ(x) = s k=0 q k (x) 2 ∈ Σ[x] d ,
C * = cl ∞ d=0 C * d .
Proof. As C * d ⊂ C * d+1 for all d ∈ N, the result follows from

C * = ∞ d=0 C d * = cl conv ∞ d=0 C * d = cl ∞ d=0 C * d .
In other words, C * d approximates C * = conv {xx T : x ∈ R n + } (i.e., the convex hull of second-order moments of Dirac measures with support in R n + ) from inside by second-order moments of measures µ σ ≪ µ whose density is a s.o.s. polynomial σ of degree at most 2d, and better and better approximations are obtained by letting d increase. (2.16)

3a 3 + 15a 2 b + 29a 2 c + 16ab 2 + 50abc + 29ac 2 + 4b 3 + 16b 2 c + 15bc 2 + 3c 3 ≥ 0.
Alternatively, by homogeneity, instead of µ we may take the probability measure ν uniformly supported on the simplex ∆ and with moments ỹ = (ỹ α ) given in (2.11), in which case the corresponding matrix M 1 (f A ỹ) now reads:

(2.17 It turns out that up to a multiplicative factor both matrices (2.15) and (2.17) have same determinant and so define the same convex cone C 1 (the same connected component of (2.16)). Figure 1 below displays the projection on the (a, b)-plane of the sets C 1 and C intersected with the unit ball.

2.

3. An alternative representation of the cone C * d . From its definition (2.13) in Theorem 2.2, the cone C * d ⊂ S n is defined through the matrix variable X ∈ S s(d) which lives in a (lifted) space of dimension s(d)(s(d) + 1)/2 and with the linear matrix inequality (LMI) constraint X 0 of size s(d). In contrast, the convex cone C d ⊂ S n defined in (2.10) is defined solely in terms of the entries of the matrix A ∈ S n , that is, no projection from a higher dimensional space (or no lifting) is needed.

We next provide another explicit description on how C * d can be generated with no LMI constraint and with only s(d) variables, but of course this characterization is not well suited for optimization purposes.

Since every g ∈ Σ[x] d can be written ℓ g 2 ℓ for finitely many polynomials (g ℓ ) ⊂ R[x] d , the convex cone Σ[x] d of s.o.s. polynomials can be written

Σ[x] d = conv {g 2 : g ∈ R[x] d }.
Next, for g ∈ R[x] d , with vector of coefficients g = (g α ) ∈ R s(d) , let g (2) = (g

(2)
α ) ∈ R s(2d) be the vector of coefficients of g 2 , that is,

g(x) = α∈N n d g α x α → g(x) 2 = α∈N n 2d g (2) α x α . Notice that for each α ∈ N n 2d , g (2)
α is quadratic in g. For instance, with n = 2 and d = 1, g(x) = g 00 + g 10 x 1 + g 01 x 2 with g = (g 00 , g 10 , g 01 ) T ∈ R s (1) , and so g (2) = (g 2 00 , 2g 00 g 10 , 2g 00 g 01 , g 2 10 , 2g 10 g 01 , g 2 02 ) T ∈ R s (2) .

Next, for g

∈ R[x] d and every 1 ≤ i ≤ j ≤ n, let G d = (G ij ) ∈ S n be defined by: (2.18) G ij := R n g(x) 2 x i x j dµ(x) = α∈N n 2d g (2) α (α i + 1)!(α j + 1)! k =i,j α k !,
for all 1 ≤ i ≤ j ≤ n. We can now describe C * d .

Corollary 2.4. Let G d ∈ S n be as in (2.18). Then:

(2.19) C * d = cl conv {G d : g ∈ R s(d) } .
Proof. From Theorem 2.2

C * d = cl ( X, M d (x i x j y) ) 1≤i≤j≤n : X ∈ S s(d) + .
As in the proof of Theorem 2.2, writing X 0 as

s k=1 g k g T k for some vectors (g k ) ⊂ R s(d) (and associated polynomials (g k ) ⊂ R[x] d ), X, M d (x i x j y) = s k=1 R n + x i x j g k (x) 2 dµ, = s k=1 G k ij , for all 1 ≤ i ≤ j ≤ n, where G k ij is as in (2.18) (but now associated with g k instead of g). Hence ( X, M d (x i x j y) ) 1≤i≤j≤n = s k=1 G k d ∈ conv {G d : g ∈ R s(d) }. Hence C * ⊆ cl conv {G d : g ∈ R s(d) } .
Conversely, let h ∈ conv {G d : g ∈ R s(d) }, i.e., for some integer s and polynomials (g k ) ⊂ R[x] d , and for all 1 ≤ i ≤ j ≤ n,

h ij = s k=1 λ k >0 G k ij = R n + s k=1 λ k g 2 k (x) x i x j dµ(x) = X, M d (x i x j y) where X = s k=1 λ k g k g T k ∈ S s(d) + . Hence cl conv {G d : g ∈ R s(d) } ⊆ C * .
The 2.4. K-copositive and K-completely positive matrices. Let K ⊂ R n be a closed convex cone and let C K be the convex cone of K-copositive matrices, i.e., matrices A ∈ S n such that x T Ax ≥ 0 on K. Its dual cone C * K ⊂ S n is the cone of K-completely positive matrices.

Then one may define a hierarchy of convex cones (C K ) d and (C * K ) d , d ∈ N, formally exactly as in Theorem 2.1 and Theorem 2.2, but now y is the moment sequence of a finite Borel measure µ with supp µ = K (instead of supp µ = R n + in (2.7)), i.e., y = (y α ), α ∈ N n , with: But of course, for practical implementation, one need to know the sequence y = (y α ), α ∈ N n , which was easy when K = R n + . For instance, if K is a polyhedral cone, by homogeneity of f A , one may equivalently consider a compact base of K (which is a polytope K ′ ), and take for µ the Lebesgue measure on K ′ . Then all moments of µ can be calculated exactly. The same argument works for every convex cone K for which one may compute all moments of a finite Borel measure µ whose support is a compact base of K.

  ) ). See for instance Rockafellar [15, Theorem 3.8; Corollary 16.4.2].

Example 1 .Figure 1 . n = 2 :a

 112 Figure 1. n = 2: Projection on the (a, b)-plane of C versus C 1 , both intersected with the unit ball the exponential measure on R 2 + with moments defined in (2.7). With f A (x) := ax 2 1 + 2bx 1 x 2 + cx 2 2 and d = 1, the condition M d (f A y) 0 which reads

  + b + c) 6(3a + 2b + c) 6(a + 2b + 3c) 6(3a + 2b + c) 12a + 6b + 2c 3a + 4b + 3c 6(a + 2b + 3c) 3a + 4b + 3c 2a

(2. 20 )

 20 y α := K x α dµ, ∀α ∈ N n (and where supp µ = K).And so with y as in (2.20) Theorem 2.1 and 2.2, as well as Corollary 2.3, are still valid. (In (2.14 replace R n + with K ). Therefore,C K = ∞ d=0 (C K ) d and C * K = cl ∞ d=0 (C K ) * d .

  and µ σ (B) = B σ dµ for all Borel sets B. So Theorem 2.2 states that C * d is the closure of the convex cone generated by second-order moments of measures dµ σ = σdµ, absolutely continuous with respect to µ (hence with support on R n + ) and with density being a s.o.s. polynomial σ of degree at most 2d. Of course we immediately have:

	Corollary 2.3. Let C *

d , d ∈ N, be as in

(2.14)

. Then C * d ⊂ C * d+1 for all d ∈ N, and

  characterization (2.19) of C * d should be compared with the characterization conv {xx T : x ∈ R n + } of C * . Example 2. With n = 2 and d = 1, G 1 reads:

	2g 2 00 + 12g10(g00 + g01) + 4g01(g00 + g01) + 24g 2 10 g 2 00 + 4g00(g10 + g01) + 6(g 2 10 + g 2 01 ) + 8g10g01	g 2 00 + 4g00(g10 + g01) + 6(g 2 10 + g 2 01 ) + 8g10g01 2g 2 00 + 4g10(g00 + g10) + 12g01(g00 + g10) + 24g 2 01

Acknowledgement

The author wishes to thank two anonymous referees for their very helpful remarks and suggestions to improve the initial version of this note.


Ulf Bauchinger

Scott R Mcwilliams

Harald Kolb

Vladimir M Popenko

Edwin R Price

Herbert Biebach

Flight muscle shape reliably

Keywords: Calibration, validation, pectoral muscle, body condition, passerine

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

and validation of a "muscle-meter" device designed to accurately measure the shape of the flight muscle. For two species of different size, the European starling (Sturnus vulgaris) and the garden warbler (Sylvia borin), we compared the accuracy and precision of different linear regression models for predicting flight muscle mass. The multi factorial linear regression model with the most support for both species included "muscle-meter score" (mm score ), tarsus length and body mass (m b ), although a simpler model with mm score and m b had as much support for predicting flight muscle mass of European starlings. A validation exercise revealed that flight muscle mass of these two species could be estimated with a relative error of about 3%. The "muscle-meter" is a simple device, easy and quick to handle, that can reliably and non-invasively estimate flight muscle mass of captive and wild birds when used in conjunction with standard measurements of tarsus length and m b . The pectoral muscle in birds amounts to 10-25% of body mass for most actively flying species with some extreme cases of even more than 30% [START_REF] Hartman | Locomotor mechanisms of birds[END_REF]. It is thus the biggest single organ within a passerines" body. Besides its main function for flight locomotion it serves as a protein source during times with high protein demands or low protein intake [START_REF] Jones | Level of reserve protein as proximate factor controlling timing of breeding and clutch-size in red-billed quelea Quelea quelea[END_REF]Piersma and Jukema 1990;[START_REF] Houston | Use of labeled methionine to investigate the contribution of muscle proteins to egg-production in zebra finches[END_REF][START_REF] Cottam | The use of muscle protein for egg production in the Zebra Finch Taeniopygia guttata[END_REF]; for review see: Piersma and Lindstöm 1997;[START_REF] Bauchinger | The role of protein during migration in passerine birds[END_REF][START_REF] Jenni | Fuel supply and metabolic constraints in migrating birds[END_REF]. The pectoral muscle generates the force for the down stroke of the wing, while the supracoracoideus muscle lifts the wing. Both muscles together are often referred to as flight muscle. Because the pectoral muscle roughly contributes to about 90% of the flight muscle mass, with the remainder being associated to the supracoracoideus muscle [START_REF] Hartman | Locomotor mechanisms of birds[END_REF], the term flight and pectoral muscle are often used interchangeably.

Zusammenfassung

Mass of the flight muscle in birds changes during the annual cycle and especially during growth, migration, egg laying, molt and seasonal acclimatization (for review see: [START_REF] Piersma | Rapid reversible changes in organ size as a component of adaptive behaviour[END_REF][START_REF] Bauchinger | The role of protein during migration in passerine birds[END_REF][START_REF] Cottam | The use of muscle protein for egg production in the Zebra Finch Taeniopygia guttata[END_REF][START_REF] Piersma | Phenotypic flexibility and the evolution of organismal design[END_REF][START_REF] Dietz | Intraspecific variation in avian pectoral muscle mass: constraints on maintaining manoeuvrability with increasing body mass[END_REF], Swanson 2010). In general, flight muscle mass is positively related to intraspecific variation in body size [START_REF] Calder | Size, Function, and Life History[END_REF]. Flight muscle mass often tracks body mass (m b ) during various conditions like flight, starvation and re-feeding (Lindström et al. 2000), during the annual cycle [START_REF] Dietz | Intraspecific variation in avian pectoral muscle mass: constraints on maintaining manoeuvrability with increasing body mass[END_REF]) and during migration (Fry et al.1972;[START_REF] Marsh | Adaptations of the gray catbird Dumetella carolinensis to long-distance migration: flight muscle hypertrophy associated with elevated body mass[END_REF][START_REF] Bauchinger | Differential catabolism of muscle protein in garden warblers (Sylvia borin): flight and leg muscle act as a protein source during longdistance migration[END_REF][START_REF] Dietz | Intraspecific variation in avian pectoral muscle mass: constraints on maintaining manoeuvrability with increasing body mass[END_REF]). However, changes in flight muscle mass are not always associated with changes in body mass. Eared grebes (Podiceps nigricollis) increased m b while flight muscle mass decreased when they were flightless at an autumn staging ground, and then m b decreased while flight muscle mass increased just prior to the grebes leaving the staging area [START_REF] Gaunt | Rapid atrophy and hypertrophy of an avian flight muscle[END_REF], Jehl 1997). Migratory garden warblers (Sylvia borin) during migration showed a general positive relationship between flight muscle mass and m b when sampled immediately before and after long flight periods [START_REF] Bauchinger | Differential catabolism of muscle protein in garden warblers (Sylvia borin): flight and leg muscle act as a protein source during longdistance migration[END_REF]. However, garden warblers arriving at stopover sites at the edge of an ecological barrier and prior to the flight across this barrier had relatively heavy flight muscle mass even though m b was relatively low (Bauchinger and Biebach 2005).

Several methods are available for non-invasively estimating flight muscle mass of larger birds, although those used for songbirds provide only qualitative indices of flight muscle mass that are subjectively assessed. Measuring flight muscle thickness by use of ultrasound has proven a useful method to estimate flight muscle mass in different wader species (Lindsström et al. 2000, Dietz et al. 1999;Dietz et al. 1999b). Dietz et al. (1999) report that flight muscle mass of two wader species can be estimated with an individual error of 20 to 25 %. Contouring the shape of the muscle with solder wire over a defined point at the Carina sterni provided a reliable estimate of flight muscle mass in gulls [START_REF] Bolton | An improved technique for estimating pectoral muscle protein condition from body measurements of live gulls[END_REF] and geese [START_REF] Nyeland | Field methods to assess pectoral muscle mass in moulting geese[END_REF]. [START_REF] Selman | A technique for measuring lean pectoral muscle mass in live small birds[END_REF] molded the shape of the flight muscle of passerine birds into a gel used for dental imprints and quantified the volume of the flight muscle from the negative imprint. An index of flight muscle size of geese and songbirds was obtained by scoring the shape of the flight muscle (e.g., concave, straight, convex, bulging) on both sides of the Carina sterni [START_REF] Bairlein | Manual of field methods ESF[END_REF][START_REF] Nyeland | Field methods to assess pectoral muscle mass in moulting geese[END_REF]. The advantage of such an index is that the muscle score is simple, repeatable and fast, and can thus be performed on hundreds of individuals at ringing stations. However, such indices provide an estimate of only relative flight muscle shape and not flight muscle mass (but see for geese [START_REF] Nyeland | Field methods to assess pectoral muscle mass in moulting geese[END_REF]. Other techniques measure flight muscle thickness or shape, but require time and handling effort, more than one person, or costly equipment which limits their usefulness especially for studies of freeliving birds (see [START_REF] Bairlein | Manual of field methods ESF[END_REF].

Herein, we describe a new technique to estimate flight muscle shape that, when combined with other size measures, it can quickly and reliably estimate flight muscle mass of free-living wild-caught birds. The "muscle-meter" device is simple, easy to handle by one person, and provides fast yet accurate measurements of the flight muscle shape. For two species of different size, the European Starling (Sturnus vulgaris) and the garden warbler (Sylvia borin), we compared the accuracy and precision of different linear regression models for predicting flight muscle mass. Our data sets were as heterogeneous as possible, providing maximum variability in terms of structural size, m b and flight muscle mass in order to demonstrate the predictability of flight muscle mass by application of "muscle-meter" measurements.

Materials

Accurate use of any non-invasive technique for estimating body composition of free-living birds requires performance of a validation study [START_REF] Scott | The use of total body electrical conductivity (TOBEC) to determine body composition in vertebrates[END_REF]. One validation approach involves building predictive models for estimating flight muscle mass using a subset of birds (calibration) and then using the models to predict flight muscle mass in another subset of birds that were not used to develop the predictive models (validation). We predicted flight muscle mass (m fm ) of European starlings (Sturnus vulgaris) and garden warblers (Sylvia borin) given their body mass, tarsus length, and "muscle-meter" measures of their flight muscle shape. Body mass of European Starlings ranges between 60 and 90 g (Snow et al. 1998), whereas body mass of garden warblers ranges between 15 and 25 g [START_REF] Bairlein | Nutritional requirements for maintenance of body weight and fat deposition in the long-distance migratory garden warbler, Sylvia borin (Boddaert)[END_REF][START_REF] Bairlein | Body mass of Garden Warblers (Sylvia borin) on migration: a review of field data[END_REF], with extreme values beyond those ranges for both species. All birds used in this contribution were part of investigations on migration and flight physiology, and so spanned the range of documented body mass (and flight muscle mass). We captured 116 starlings during August in Ontario (43°10´N; 81°19´W), Canada, and maintained them for up to four months in captive facilities at the AFAR (Advanced Facility for Avian Research) facility at University of Western Ontario.

A total of 45 garden warblers were used in this study with 14 birds captured in Turkey (36°40´N; 33°05´E) during autumn migration [START_REF] Bauchinger | The role of protein during migration in passerine birds[END_REF][START_REF] Bauchinger | Differential catabolism of muscle protein in garden warblers (Sylvia borin): flight and leg muscle act as a protein source during longdistance migration[END_REF]) and 31 birds captured on Crimea peninsula (45°25´N; 32°32´E) in Ukraine during spring migration. Data sets for both species were decidedly heterogeneous with respect to sex, body size and body mass.

For all birds we measured tarsus length, (Svenson 1992), m b , and the shape of the flight muscle at a defined location (subsequently, mm score ). Figure 1 shows the "musclemeter" device and describes the measurement method in detail. We took three consecutive measures of mm score for each starling and then used the mean of these three measures for subsequent modeling, whereas we took only one measure of mm score for each garden warbler. Subsequently, birds were anesthetized under Isoflurane and killed by cervical dislocation and the left pectoral and supracoracoideus muscles were dissected. Fresh mass of the muscles was determined immediately for starlings and within weeks upon storage in liquid nitrogen for warblers (for methods see [START_REF] Bauchinger | Differential catabolism of muscle protein in garden warblers (Sylvia borin): flight and leg muscle act as a protein source during longdistance migration[END_REF]. For calculation of flight muscle mass (m fm ) of each bird we combined the samples of the pectoral and the supracoracoideus muscles, multiplied those by two, and hereafter refer to this as m fm . We calculated repeatability for the measurements of mm score by use of the "muscle-meter" according to method described by [START_REF] Lessels | Unrepeatable Repeatabilities: A common Mistake[END_REF].

Procedure for developing predictive models

The complete data set for each species was used to build regression models with flight muscle mass as the dependent variable and all possible combinations of the following independent variables: tarsus length, mm score , or m b . We used the r 2 -value and the Akaike Information Criterion (AICc) for low sample size to evaluate which model(s) explained the most variance and had the highest level of support [START_REF] Burnham | Model selection and multi-model inference: a practical information theoretic approach[END_REF]. All data were tested for normality (K-S test) prior to regression analysis. Linear regression models were used to describe the relationship between the independent variables to m fm , because cubic and quadratic models did not explain more of the variance and were more parsimonious.

We estimated the accuracy and precision of the predictive models using a validation approach. We randomly selected 88 % of the birds of each species ( 105European starlings and 39 garden warbler) and built regression models with flight muscle mass as the dependent variable and all possible combinations of the following independent variables: tarsus length mm score , or m b . We then estimated the m fm of the remaining birds (11 European starlings and 6 garden warblers) using these regression models built with the 88 % subset of birds. To ensure that the birds used to build each model spanned the full range of body mass for that species, we assigned individuals with the two most extreme values for tarsus length, mm score and m b of each of the two species to the calibration group used to build the model, and randomly selected the remaining individuals for the calibration and the validation group. We estimate the absolute and relative error terms for the predictive models for the comparison between the estimated m fm and the m fm value determined after dissection by weighing the tissue, because coefficients of determination are not always indicative of error when evaluating the usefulness of regression equations (e.g., Skagen et al. 1993). Absolute error was calculated as predicted-actual , and relative error as [100 ( predictedactual ) / actual]. We used PASW Statistics (version 18) for the analysis.

Results

Distribution and range of measurements

We used untransformed tarsus length, mm score and m b to predict m fm because these data were normally distributed (K-S test p > 0.40 for all parameters and both species). For starlings (n =115), the range for tarsus length was 5.1 mm (min. to max. 26.7 to 31.8 mm), the range for mm score was 3.46 mm (-6.27 to -2.79 mm), the range for m b was 26.2 g (55.3 to 81.5 g), and the range for flight muscle mass was 7.44 g (10.86 to 18.30 g). For garden warblers (n = 45), the range for tarsus length was 2.7 mm (min. to max., 18.9 to 21.6 mm), the range for mm score was 1.75 mm (-3.74 to -1.99 mm), the range for m b was 12.8 g (14.7 to 27.5 g) and range for flight muscle mass was 1.26 g (2.23 to 3.49 g).

In European starlings measurements of mm score showed a repeatability (r) of 0.79 (r = s 2 A / (s 2 + s 2 A ); where s 2 A = 0.39 is the among group variance component and s 2 = 0.10 the within-group variance component).

Calibration and Akaike information criterion

We compared r 2 -values and Akaike information criteria for the seven linear regression models for each species with flight muscle mass as the dependent variable and all possible combinations of the following independent variables: tarsus length, mm score , and m b (Table 1). All starlings (n = 116) and garden warblers (n = 45) were used to build these regression models. For garden warblers, the three-factor model (eq 2) was the only model that was substantially supported, with ∆ AICc values of nine and higher for the remaining models (Table 1). For European starlings, both the three-factor model (eq 1) as well as the twofactor model (eq 3) were substantially supported (Table 1), with ∆ AICc values of zero for the two factorial linear regression model (eq. 3), and a ∆ AICc value of 2.9 for the three factorial model (eq. 2). All other models for the estimation of m fm for European starlings failed to substantially explain variation in these data (Table 1). This most parsimonious two-factor linear regression model for predicting flight muscle mass (m fm ) was as follows:

European Starling: m fm = 5.003+ (0.860 mm score ) + (0.195 m b ) (eq. 3).

Validation

Consistent with the above results from the calibration, when we compared predicted and actual m fm using the validation subset of both bird species, the absolute and relative errors for estimated m fm were lowest for the three-factor linear regression model (Table 2). Error terms for both species increased with increasing ∆ AICc values. For the European starling, mean absolute error increased from 0.39 g for the three-factor model to 1.16 g for the single-factor model that included only tarsus length. For the garden warbler, mean absolute error increased from 0.12 g for the three-factor model to 0.27 g for the single-factor model that included mm score . Mean relative error for the three-factor model was 2.6 % for the European starling and 3.1 for the garden warbler and it increased to 7.8 % and 9.1 %, respectively, for the single-factor models.

Discussion

Our calibration and validation experiments for two passerines show that using the newly developed "muscle-meter" in conjunction with measurements for whole animal mass and structural size results in reliable estimates of m fm in live birds. In order to obtain estimates for m fm with a relative error margin of about 3%, two additional measurements were included in the linear regression model, namely the tarsus length as a measure of structural size and m b . Both these additional measures are routinely collected at ringing stations. For European starling use of the mm score in conjunction with m b already provides reliable estimates of m fm , even without additional measurements of the tarsus length.

Whether or not tarsus length is taken into account for estimating m fm does not alter the relative error margins (Table 2). Besides providing a reliable estimate of flight muscle mass of an individual captured once (as commonly occurs at a ringing station), this technique allows repeated measurements of the same individual either in captive studies or in mark re-capture studies in the field. We used the "muscle-meter" with two different mouth gaps, 10 mm for bigger birds like the European starlings and 6 mm for smaller birds like the garden warbler, and this provided accurate estimates of the m fm with similar relative error (3 and 4 %, respectively) when used in conjunction with measurements for tarsus length and m b . We estimated the m fm of European starlings with 0.39 ± 0.33 g absolute error which was relatively small given mean flight muscle of starlings was 14.97 ± SD 5.16 g. Likewise, we estimated the m fm of garden warblers with 0.12 ± 0.09 g absolute error which was relatively small given mean flight muscle of these warblers was 2.89 g ± SD 0.29.

Although we measured mm score only once for each garden warbler and repeatability for the measurements of mm score is high (0.79 for the starlings), we recommend three repeated measures of the mm score , as done for each European Starling, in order to ensure the accuracy and precision of these estimates.

The accuracy and precision of the muscle meter is sufficient for studying phenotypic flexibility of flight muscle in wild birds. Flight muscle mass of garden warblers decreased on average by 0.51 g during their spring flight across the Sahara desert (Bauchinger et al. 2005) which is more than four times the error that we report for the estimates of m fm from the "muscle-meter" measurements in conjunction with tarsus length and m b . Similarly, m fm of garden warblers increased on average by 0.27 g during molt while in the African wintering area [START_REF] Bauchinger | Transition between moult and migration in a longdistance migratory passerine: organ flexibility in the African wintering area[END_REF], which is more than double the error that we report for the estimates of m fm from the "muscle-meter" measurements. Other studies on phenotypic mass changes in garden warblers indicate similar mass changes that could be assessed using the muscle meter [START_REF] Hume | Digestive tract function in the long-distance migratory garden warbler, Sylvia borin[END_REF]Biebach 1996, Schwilch et al. 2002). In European starlings the maximum difference in fresh m fm between exercised and unexercised birds amounted to 1.8 g (based on mass for pectoral muscle, Swaddle and Biewener 2000) and this is more than four times higher than the error that we report for the estimates of m fm from the muscle meter measurements. These documented differences for m fm in wild birds or in captive birds subjected to ecologically relevant treatment indicate that phenotypic changes in the flight muscle mass frequently occurs in response to exercise, seasonal acclimatization, migration and/or moult and that those phenotypic mass changes could be detected by estimation of m fm using the "musclemeter".

The "muscle-meter" device has been used to measure the change in flight muscle shape in several studies of birds, although none of these studies predicted m fm from measured mm score because no validation study like that presented here had been completed.

For example, [START_REF] Schmidt-Wellenburg | Energy expenditure and wing beat frequency in relation to body mass in free flying Barn Swallows (Hirundo rustica)[END_REF][START_REF] Schmidt-Wellenburg | Energy expenditure during flight in relation to body mass: effects of natural increases in mass and artificial load in Rose Coloured Starlings[END_REF] reported changes in mm score of up to 0.55 mm for rose-colored starlings (Sturnus roseus) before and after a 6-hour flight in a windtunnel [START_REF] Schmidt-Wellenburg | Energy expenditure during flight in relation to body mass: effects of natural increases in mass and artificial load in Rose Coloured Starlings[END_REF]. [START_REF] Bize | Intra-individual varation in body temperature and pectoral muscle size in nestling Alpine swifts Apus melba in response to an episode of inclement weather[END_REF] used the "muscle-meter" (mouth gap of 10 mm) to determine changes in flight muscle shape in alpine swifts (Apus melba) in response to a prolonged period of inclement weather. They documented a rapid decrease in mm score of about 1.4 mm during inclement weather. European starlings that experienced rich, intermediate or poor foraging environments showed marked differences in mm score of as much as 1.08 mm (Wiersma et al. 2005). Unfortunately, Wiersma et al.

(2005) used a mouth gap of 6 mm for their muscle meter whereas we used a wider mouth gap of 10 mm. Therefore, we are not able to the actual change in muscle mass of the starlings from the data provided by Wiersma et al. (2005). We recommend that standard mouth gaps be used for each species so that results from different studies are most comparable. Clearly, the muscle meter detects significant changes in flight muscle shape that, given validation studies such as ours, can be used to predict changes in flight muscle mass. In concert with a calibration like the ones presented here for European "Muscle-meter" and mounted caliper used to measure the shape of the flight muscle. Left top image shows three "muscle-meter" devices of which two have a free gap of 6mm (A1) and one has a free gap of 10mm (A2). The larger gap is used for bigger passerines birds like starlings, while the 6mm gap is used for measurements on smaller sized birds, like the garden warbler. Picture on the top right demonstrates handling of the "muscle-meter", with the front panel sliding freely above the back panel. The front panel is then placed perpendicular on the flight muscle of a bird, which is held horizontally in the alternate hand (see bottom left picture). The location where the "Muscle-meter" is placed on the flight muscle is defined as the middle of the Carina sterni between its posterior and anterior end points. The "muscle-meter" is lowered perpendicularly onto the flight muscle until the back panel rests on the Carina sterni. This position of the muscle meter is now secured by pressing the free thumb of the hand that is holding the "musclemeter" against the front panel thereby locking the position of the "muscle-meter" (see bottom middle picture), and subsequently by closing the screw with a finger from the hand that holds the bird (see bottom right picture). The distance measured between the back panel (= level of the Carina sterni) and the front panel that measures the shape of the muscle (mm score in the muscle sketch) is now fixed at the lower part of the "muscle-meter". This distance is similar to the distance between the two panels at the upper end of the "muscle-meter", because both panels have equal length. The upper part of the "musclemeter" is now pushed gently into the slide at the measuring device B that holds a mounted caliper and the final mm score is read from the caliper display. The caliper is mounted onto a metal block to provide stability and allows the measurement to be performed with one hand only, with the other hand still holding the bird. Before a set of measurements is recorded the calibration tool C is pushed into the slide and the value is set to zero to calibrate the caliper and its position in the device. The caliper is mounted in the measuring device B in a way that the zero point already represents an input tension, i.e. without the calibration tool C the caliper reads a value > zero. We used an input tension of 0.14mm; actual values for the mm score are read as negative values so that increasing values signify increasing muscle bulging. Using three repeated mm score measures to create a mean value reduces the likelihood of odd values to find their way into the data set.

  (Musculus pectoralis) ist das größte Organ im Singvogelkörper. Er ermöglicht das Fliegen, stellt aber auch eine körpereigene Proteinquelle dar, die bei erhöhtem Proteinbedarf oder reduzierter Proteinverfügbarkeit genutzt wird. Die Masse des großen Flugmuskels zeigt ein dynamisches Verhalten und ändert sich beispielsweise während des jugendlichen Wachstums, der Reproduktion, der jahreszeitlichen Anpassung, des Fastens und des Zuges. Ein Hilfsmittel, zu einer exakten und nicht-invasiven Erfassung dieser phänotypischen Flexibilität des Flugmuskels ist deshalb von großem Interesse für Ornithologen. Hierzu stellen wir hier das "Muskel-Meter" sowie seine Kalibrierung und Validierung, zur exakten Messung der Form des Flugmuskels vor. Bei zwei unterschiedlich großen Vogelarten, dem Europäischen Star (Sturnus vulgaris) und der Gartengrasmücke (Sylvia borin), haben wir multifaktorielle lineare Regressionsmodelle auf ihre Genauigkeit hinsichtlich einer präzisen Voraussage der Flugmuskelmasse überprüft. Das Regressionsmodell mit der besten Anpassung an die empirischen Daten beider Vogelarten beinhaltet den "Muskel-Meter Score" (mm score ), die Tarsuslänge und die Körpermasse (m b ). Das einfachere Modell mit mm score und m b erzielt eine ähnliche Übereinstimmung, allerdings nur für den Europäischen Star. Eine Validitätsprüfung zeigt, dass die Flugmuskelmasse der beiden Arten mit einem relativen Fehler von 3% abgeschätzt werden kann. In Kombination mit Standardmessungen, wie der Tarsuslänge und der Körpermasse, ist somit das Muskel-Meter bestens für eine rasche und genaue, nicht-invasive Abschätzung der Flugmuskelmasse von Wildvögeln und Vögeln in Tierhaltung geeignet.

  r 2 -values increased and Akaike information criteria decreased from single factor to multiple factor regression analysis and in both species the linear regression model combining the three factors of tarsus length, m b and mm score had the highest r 2 -values and substantially more support based on Akaike information criteria. For European starling, both, a two-factor linear regression model and the three-factor linear regression model provide good support based on Akaike information criteria (please see below). The three-factor linear regression model for predicting flight muscle mass (m fm ) was as follows for each species: European starling: m fm = 3.797 + (0.886 mm score ) + (0.189 m b ) + (0.059 tarsus length) (eq. 1) Garden warbler: m fm = -1.212 + (0.293 mm score ) + (0.045 m b ) + (0.199 tarsus length) (eq. 2)

  Figure1

Acknowledgements

The "muscle-meter" development started in 1996 when H.B, H.K and U.B developed the first prototype. We thank Franz Antoni and Klaus Pichler for their excellent support constructing the "Muscle-meter" and their help in finding solutions to simplify the measurement envisioning an easy to handle device. Josef Chernichko and other members of the Azov-Black Sea Ornithological Station, Ukraine made the collection of garden warblers possible. We thank Chris Guglielmo, and other "AFAR" members, for their generous support of the European starling project, and Lillie Langlois, Michelle Boyles and Alexander Gerson for their contributions in running the experiment. David Swanson and an anonymous reviewer provided constructive comments that helped to improve the manuscript. The research project on starlings performed in Canada was approved by the University of Western Ontario Animal Use Sub-Committee (Starlings, protocol # 2006-011-04). The field trip to Ukraine was funded by an ESF-travel grant to U.B., the research on European starlings was funded by a NSF-research grant to S.McW.


Hervé Jégou

Jonathan Delhumeau

Jiangbo Yuan

Guillaume Gravier

Patrick Gros

Babaz: a large scale audio search system for video copy detection

Keywords: audio fingerprinting, audio search, copy detection, reciprocal neighbors, TRECVID

INTRODUCTION

Due to the emergence of networks and personal capture devices, large volumes of image and video are now available in digital format, which has enabled new services and easier access to large amount of content. On the downside, the marginal cost of copying a video is extremely low. In conjunction with the generalization of high-speed Internet connection, a pirate can easily and massively diffuse illegal content on Internet, in particular through peer-to-peer networks. This is a threat which has a critical impact on copyright holders and artists, which have less revenues because the content is used without being paid for. Moreover, the videos may be uploaded on websites hosting user-generated content (UGC), such as Youtube or Dailymotion. These companies are likely to host illegal content without knowing it, given that the mount of videos uploaded every day forbids a manual verification.

For these reasons, automatic detection of pirated videos on peerto-peer networks and UGC websites has received a significant attention from the multimedia and signal processing communities in the last few years, in particular using the image channel [START_REF] Law-To | Video copy detection: a comparative study[END_REF]. Research on that topic is supported by the National Institute of Standard and Technology (NIST), which since 2008 organizes a specific task in TRECVID [START_REF] Smeaton | Evaluation campaigns and trecvid[END_REF] to evaluate the different techniques. However, as observed in TRECVID's 2010 evaluation campaign, most of the systems only consider the visual channel, or use off-the-shelf solutions adapted from the audio community [START_REF] Cano | A Review of Audio Fingerprinting[END_REF], such as the one proposed by Wang [START_REF] Wang | An industrial-strength audio search algorithm[END_REF] in the original Shazaam audio fingerprinting system or the Echoprint project 1 . Although the current Shazaam implementation probably departs from that of the original paper, a Matlab implementation by Dan Ellis of this method is available online [START_REF] Ellis | Robust Landmark-Based Audio Fingerprinting[END_REF] and will be used as the reference method for this reason. In this paper, we propose an audio fingerprint system, called BABAZ, that is intended for video copy detection based on audio, unlike most concurrent schemes which were mainly addressing the problem of query-by-example on mobile devices [START_REF] Ellis | Robust Landmark-Based Audio Fingerprinting[END_REF][START_REF] Chandrasekhar | Survey and evaluation of audio fingerprinting schemes for mobile query-by-example applications[END_REF]. It has been successfully used for the copy detection task of the TRECVID'2011 campaign, see [START_REF] Ayari | Inria@trecvid'2011: Copy detection & multimedia event detection[END_REF]. It is available online2 under open-source license.

BABAZ is inspired by state-of-the-art image search systems. The first step extracts descriptors from the waveform, which are treated as sequences of Euclidean descriptors. This departs from [START_REF] Wang | An industrial-strength audio search algorithm[END_REF], based on landmark hashes. The nearest neighbors of the query descriptors are retrieved using an indexing structured shared by all audio tracks, and compute a similarity measure that takes into account the approximate reciprocal neighborhood of the selected nearest neighbors. The temporal consistency is exploited through a temporal Hough transform, similar to [START_REF] Wang | An industrial-strength audio search algorithm[END_REF] and visual video search systems [START_REF] Law-To | Video copy detection: a comparative study[END_REF][START_REF] Douze | An image-based approach to video copy detection with spatio-temporal post-filtering[END_REF]. Finally, we introduce the idea of time-shifting on query side. It consists in submitting to the system several versions of the query, shifted by a few milliseconds. This approach, which approximately multiplies the query time by the number of versions considered, increases the recall by avoiding the cases where the query and database videos are sampled in phase opposition.

This paper is organized as follows. Section 2 describes the main components of BABAZ and Section 3 reports our experiments in two application scenarios: the dataset used in TRECVID'2010 copy detection task and a music landmark detection setup as proposed in [START_REF] Wang | An industrial-strength audio search algorithm[END_REF]. Our experiments demonstrate the better performance of our system for video copy detection.

BABAZ AUDIO SEARCH SYSTEM

This section describes the main components of BABAZ. We will consider two application cases, leading to two variants of our system, namely BABAZv and BABAZm:

• BABAZv considers the application case for which our system was originally designed, i.e., a copy detection setup such as the one considered in TRECVID [START_REF] Smeaton | Evaluation campaigns and trecvid[END_REF], where the signal is the audio track of a video. It typically includes voices, silences and occasionally music. The audio may have been transformed by different kinds of transformations, such as strong pass-band filter, compression, mixing, single-or multi-band companding, etc.

• BABAZm considers an application scenario comparable to that of [START_REF] Wang | An industrial-strength audio search algorithm[END_REF], i.e., music retrieval in a noisy environment. This setup is considered, in particular, for the sake of comparison with [START_REF] Wang | An industrial-strength audio search algorithm[END_REF] in their own application scenario.

BABAZv and BABAZm only differ in the choice of the filter banks, as detailed in the feature extraction step described below.

Pre-processing

The audio tracks extracted from an audio corpus are not necessarily homogeneous. Sample rates as well as encoding quality vary significantly from one track to another. It is in particular the case in the Internet Archive dataset used in the TRECVID's copy detection task, where the videos are mainly amateur videos captured and encoded by different devices and audio codecs. In order to deal with this variability in a consistent way, all the tracks are resampled to 32,000 Hz. We use the right stereo channel only when stereo is available.

Feature extraction: filter banks

Hereafter, we detail how we extract descriptors from the audio waveform. The signal is ran through a pre-emphasis filter to compensate for the spectral slope and divided into overlapping short-term windows of 25 ms taken every 10 ms. In each window, the short-term spectrum is represented by log-energies at the output of overlapping band-pass filters. As mentioned above, two variants are considered:

• BABAZv: For video copy detection, the signal is severely attacked, and may have undergone a strong band-pass filter. In this case, we use 40 filters spread along the [500 Hz,3000 Hz] frequency range on a Mel scale.

• BABAZm: In music retrieval, the audio tracks have a relatively good frequency localization, which is exploited in [START_REF] Wang | An industrial-strength audio search algorithm[END_REF]. In this application scenario, we consider 64 filters spread along the [200 Hz,4000 Hz] frequency range (Mel scale).

As a result, the dimensionality d of the descriptors is either d = 40 or d = 64. The representation based on these filters gives a rough approximation of the signal's spectral shape in the frequency range considered while smoothing out the harmonic structure, if any, and is therefore robust to many spectral distortions. We have used the freely available spro software 3 for the generation of filter banks. This software also includes an efficient implementation of the widely used MFCC descriptors. However, in our experiments, these descriptors are significantly outperformed by the filter banks.

Compound descriptors and energy invariance

The temporal consistency provided by a single filter bank is limited, as their temporal span is limited and only frequencies are considered. This is problematic since the database is large: the filter banks themselves might not be discriminative enough to identify a matching hypothesis with sufficient reliability. In order to increase the discriminative power of the descriptor, the temporal aspect is emphasized by concatenating several filter banks, as done in Serra's thesis [START_REF] Serrà | Identification of versions of the same musical composition by processing audio descriptions[END_REF] in a context of cover detection.

For a given timestamp t, 3 successive filter banks are extracted at timestamps tδ, t and t + δ, producing a compound descriptor of dimensionality 3d (i.e., 120 or 192, depending on the application 3 http://gforge.inria.fr/projects/spro scenario). We set δ = 30 ms in order to avoid overlapping. We have performed a few experiments on a validation dataset to decide on how to take into account this dynamic aspect, e.g., using derivatives of the filter bank with respect to time. Compounding the descriptors appeared a reasonable choice. As illustrated in Figure 1, the resulting span of this descriptors is 85 ms. This approach favors the temporal aspect by taking into account the dynamic behavior of the frequency energies, at the cost of an increased descriptor dimensionality.

Descriptors are compared with the Euclidean distance. For large vector databases it allows for efficient indexing algorithms, as the one described in the next subsection. Note however that we have not performed a thorough evaluation of other possible distances/divergences because the search algorithms are generally less efficient for sophisticated metrics. In particular, we have not considered the metrics specifically developed to compare filter banks, such as the Itakura-Saito divergence or the log-spectral distance,

In order to take into account attacks on the volume (signal energy), the descriptor is finally made invariant by subtracting its mean.

Approximate nearest neighbor search

As the exact search is not efficient enough, BABAZ uses an approximate nearest neighbor search technique. Many methods exist for this task, such as locality sensitive hashing or the FLANN package [START_REF] Muja | Fast approximate nearest neighbors with automatic algorithm configuration[END_REF]. However, this step has a major impact on both efficiency and search quality, and only a few methods are able to search in hundreds of millions of descriptors with reasonable quality, as required by our method to index thousands of hours of audio.

BABAZ uses the IVFADC indexing method of [START_REF] Jégou | Product quantization for nearest neighbor search[END_REF], which is able to index billions of descriptors on a commodity server. It finds the approximate k nearest neighbors using a compression-based approach, and relies on an inverted structure to avoid exhaustive search. This approximate nearest neighbor method implicitly sees multi-dimensional indexing as a vector approximation problem. It is proved [START_REF] Jégou | Product quantization for nearest neighbor search[END_REF] that the square error between the distance and its estimation is bounded, on average, by the quantization error. This ensures, asymptotically, near perfect search results when increasing the number of bits allocated for the quantization indexes:

The main parameters of this method are the number of bytes b used per database vector and the number c of inverted lists associated with the partitioning of the feature space (learned by k-means). In our case, we set b = 24 bytes and use multiple assignment [START_REF] Jégou | Product quantization for nearest neighbor search[END_REF] on query side, leading to visit 16 inverted lists out of c = 16, 384.

Scoring vote and reciprocal nearest neighbors

The search technique returns a list of k (approximate) nearest neighbors. A conventional method to exploit them consists in assigning a vote of 1 to all the corresponding audio tracks, or alternatively, a function of the rank or of the distance. Based on a recent stateof-the-art work [START_REF] Jégou | Exploiting descriptor distances for precise image search[END_REF] in image search, we adopt a different strategy, which is illustrated in Figure 2.

Denoting by d k (q) the distance between the query descriptor and its k-th nearest neighbor, the quantity d k (q) -d(q, i) is shown, based on a mutual information criterion [START_REF] Jégou | Exploiting descriptor distances for precise image search[END_REF] measured on image descriptors, to better reflect the quality of the match. This is also the case for our audio descriptors, so we adopt this weighting scheme.

The distance d k (q) is relative to the query. In order to symmetrize the relationship between the query and database descriptors, it is worth considering the reciprocal nearest neighbors of the database vector, or more specifically the typical distance between the database vector and its own k-nearest neighbors. In practice, computing the reciprocal nearest neighbors is impractical: the audio descriptor database may contains up to billions of vectors. If exact nearest neighbor search is used, then it turns out that each database vector has to be submitted to the system. Although some approximate strategies [START_REF] Dong | Efficient k-nearest neighbor graph construction for generic similarity measures[END_REF] were proposed to compute the nearest neighbor graph, these approaches were only tested on up to 1 million vectors. However, we are not interested in the neighbors themselves, but in the typical distance of a database vector to its neighborhood. This reciprocal typical distance is estimated on a limited subset of 1 million vectors. In this case, the parameter k associated with the database vectors has to be adjusted to account for the smaller size of this subset.

Re-ranking: Finally, in the spirit of [START_REF] Jégou | Searching in one billion vectors: re-rank with source coding[END_REF], the hypotheses are reranked based on exact descriptors to obtain the exact distances, in order to increase the precision of the proposed similarity. The difference with [START_REF] Jégou | Searching in one billion vectors: re-rank with source coding[END_REF] is that we use the original descriptors and not only a compressed version of these.

Energy weighting: Video tracks contain many silences. Those are filtered when the signal and consequently the descriptor is zero. However, there are also many descriptors extracted on audio frames containing almost no energy, but which are not pure silence. Filtering audio segments with low energy may lead to loose some precious information, and reduce the accuracy of the localization. For this reason, we adopt a smoother strategy and multiply the score associated with the match with the energy of the query descriptor.

Hough matching

BABAZ assumes that the transformations do not include any acceleration. Given a query, for each of its audio descriptors we first search for the k approximate nearest neighbors and compute their weighting score based on the strategy exposed above. We then vote for several alignment hypotheses (a b , ∆T ) using the scoring method introduced above. Compared with uniform voting, this brings a slight improvement at almost no cost in efficiency. The video hypotheses with low scores are filtered. On output, this Hough matching system returns a maximum of 40 hypotheses per query. Each database track is associated with a maximum of 3 ∆T hypotheses.

Detection of boundaries

At this point, for each query we may have several alignment hypotheses (id, ∆T ), where id is the database track identifier and ∆T is the difference between the query and the database timestamps. We use the whole set of descriptors and weight their distance to produce a score per time instant. This score is filtered in time and used to detect the boundaries defining a matching segment. Each segment receives a score computed from the individual scores for each time instant.

Shifted query

The audio descriptors are extracted every 10 ms, which leads to reduce the quality of the comparison if the sampling of the database track occurs in phase opposition, i.e., with a shift of 5 ms relative to the query track. To address this problem, we submit several shifted version of the query to the system. For instance, we create shifted versions of the query with shifts of 2, 4, 6 and 8 ms. This, obviously, significantly impacts the efficiency of the search by a significant factor, and should be used when high precision is required only.

EXPERIMENTS

We have conducted some experiments in two applications cases: audio landmark detection, and copy detection. Our system is compared with the approach of [START_REF] Wang | An industrial-strength audio search algorithm[END_REF], using the Matlab implementation by Dan Ellis [START_REF] Ellis | Robust Landmark-Based Audio Fingerprinting[END_REF] available online. In order for this approach to scale to the larger datasets we consider here, we have adapted the size of the hash table to avoid collusions.

Video copy detection based on audio

For audio-based video copy detection, we have considered the dataset provided for the TRECVID'2010 copy detection task, as TRECVID can now be considered as the standard evaluation setup for this application scenario. It comprises more than 11,000 video clips corresponding to about 400 hours of video. The queries are transformed excerpts (3 to 30 seconds) of some database video, and are inserted in a longer distractor track. About 1/3 of the queries have no corresponding video in the database. In this case the system should return no result. Our optimization has been performed on an independent validation dataset.

The comparison of BABAZv with [START_REF] Wang | An industrial-strength audio search algorithm[END_REF] is shown in Figure 4 using a precision-recall curve. Our system is significantly better in this application scenario. Our system is 13 times faster than real time on a 24-core machines without the shift extension. This is slower than the Wang's system. The shift provides a relatively low improvement given the significantly increased complexity. It is however worth using it in an evaluation like TRECVID, where the recognition performance is the most important measure.

Music detection

In this evaluation, we have considered an audio dataset comprising 1,000 tracks (approximately 55 hours) of mixed sound files (music, voice and other) extracted from TRECVID'dataset. We have considered an evaluation protocol similar to that of [START_REF] Wang | An industrial-strength audio search algorithm[END_REF]: 10 musical queries of various genres cropped to individual length of 5, 10 or 15 seconds. The energy of the signal is normalized before being mixed with a distractor consisting of crowd noise from a pub (also normalized). The respective signal to noise ratios are between -15dB an 15dB for the set of 330 queries. The performance of our system BABAZv, and its variant BABAZm dedicated to music retrieval, are first shown in Figure 5 for queries of 10s. Note that our result for the Wang's system [START_REF] Wang | An industrial-strength audio search algorithm[END_REF] are consistent with those reported in his paper. Our system reports comparable performance.

Another set of 330 queries were extracted by applying GSM 6.10 compression and 8kHz subsampling to the first query set. The results are presented in Figure 6 with respect to precision-recall, our system is much better. The discrepancy with Figure 5 is explained by the fact that the scores, in our system, are better normalized across queries, which is taken into account by precision-recall (which interleaves the hypotheses of all queries based on the recognition score) but not by the recognition rate.

CONCLUSION

We have presented the BABAZ audio search system, an open system dedicated to video copy detection based on audio track. It also exhibits competitive search quality for landmark retrieval, although in this context it is slower than the very efficient approaches based on spectral peaks.

  1 http://echoprint.me

Fig. 1 .

 1 Fig. 1. Compound descriptor.

Fig. 2 .

 2 Fig. 2. Reciprocal nearest neighbors and of our voting strategy.

Fig. 3 .

 3 Fig. 3. Illustration of the temporal Hough transform: the audio matches output by the approximate search engine are collected and summed up for each hypothesis (id, ∆T ). This dilutes the scores of false positives over time shift hypotheses.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Precision recall on Trecvid'10.

Fig. 6 .

 6 Fig.6. Precision recall on the music retrieval dataset.

http://babaz.gforge.inria.fr