Mohammed Nabil

E L Korso

Rémy Boyer

Alexandre Renaux

Sylvie Marcos

Statistical Resolution Limit for the Multidimensional Harmonic Retrieval Model: Hypothesis Test and Cramér-Rao Bound Approaches

Keywords: Statistical resolution limit, multidimensional harmonic retrieval, performance analysis, hypothesis test, Cramér-Rao bound, parameter estimation, multidimensional signal processing

The Statistical Resolution Limit (SRL), which is defined as the minimal separation between parameters to allow a correct resolvability, is an important statistical tool to quantify the ultimate performance for parametric estimation problems. In this paper we generalize the concept of the SRL to the Multidimensional SRL (MSRL) applied to the multidimensional harmonic retrieval model. In this paper, we derive the SRL for the so-called multidimensional harmonic retrieval model by using a generalization of the previously introduced SRL concepts that we call Multidimensional SRL (MSRL). We first derive the MSRL using an hypothesis test approach. This statistical test is shown to be asymptotically an uniformly most powerful test which is the strongest optimality statement that one could expect to obtain. Second, we link the proposed asymptotic MSRL based on the hypothesis test approach to a new extension of the SRL based on the Cramér-Rao Bound approach. Thus, a closed-form expression of the asymptotic MSRL is given and analyzed in the framework of the multidimensional harmonic retrieval model. Particularly, it is proved that the optimal MSRL is obtained for equi-powered sources and/or an equi-distributed number of sensors on each multi-way array.

I. INTRODUCTION

The multidimensional harmonic retrieval problem is an important topic which arises in several applications [START_REF] Jiang | Almost-sure identifiability of multidimensional harmonic retrieval[END_REF]. The main reason is that the multidimensional harmonic retrieval model is able to handle a large class of applications. For instance the joint angle and carrier estimation in surveillance radar system [START_REF] Vanpoucke | Algorithms and Architectures for Adaptive Array Signal Processing[END_REF], [START_REF] Haardt | 3-D unitary ESPRIT for joint 2-D angle and carrier estimation[END_REF], the underwater acoustic multisource azimuth and elevation direction finding [START_REF] Wong | Uni-vector-sensor ESPRIT for multisource azimuth, elevation, and polarization estimation[END_REF], the 3-D harmonic retrieval problem for wireless channel sounding [START_REF] Mokios | On 3-D harmonic retrieval for wireless channel sounding[END_REF], [START_REF] Schneider | Performance verification of MIMO concepts using multi-dimensional channel sounding[END_REF] or the detection and localization of multiple targets in a MIMO radar system [START_REF] Nion | A PARAFAC-based technique for detection and localization of multiple targets in a MIMO radar system[END_REF], [START_REF] Nion | Tensor algebra and multi-dimensional harmonic retrieval in signal processing for MIMO radar[END_REF].

One can find many estimation schemes adapted to the multidimensional harmonic retrieval estimation problem, see, e.g., [START_REF] Jiang | Almost-sure identifiability of multidimensional harmonic retrieval[END_REF], [START_REF] Vanpoucke | Algorithms and Architectures for Adaptive Array Signal Processing[END_REF], [START_REF] Wong | Uni-vector-sensor ESPRIT for multisource azimuth, elevation, and polarization estimation[END_REF]- [START_REF] Nion | A PARAFAC-based technique for detection and localization of multiple targets in a MIMO radar system[END_REF], [START_REF] Gershman | Space-time processing for MIMO communications[END_REF], [START_REF] Boyer | Decoupled root-MUSIC algorithm for multidimensional harmonic retrieval[END_REF]. However, to the best of our knowledge, no work has been done on the resolvability of such a multidimensional model.

The resolvability of closely spaced signals, in terms of parameter of interest, for a given scenario (e.g., for a given Signal-to-Noise Ratio (SNR), for a given number of snapshots and/or for a given number of sensors) is a former and challenging problem which was recently updated by Smith [START_REF] Smith | Statistical resolution limits and the complexified Cramér Rao bound[END_REF], Shahram and Milanfar [START_REF] Shahram | On the resolvability of sinusoids with nearby frequencies in the presence of noise[END_REF], Liu and Nehorai [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF], and Amar and Weiss [START_REF] Amar | Fundamental limitations on the resolution of deterministic signals[END_REF]. More precisely, the concept of Statistical Resolution Limit (SRL), i.e., the minimum distance

The authors are with Laboratoire des Signaux et Systèmes (L2S), Université Paris-Sud XI (UPS), CNRS, SUPELEC, 3 rue Joliot Curie, Gif-Sur-Yvette, 91192, France, phone: +331 6985 1763, fax: +331 6985 1765, {elkorso, remy.boyer, alexandre.renaux, marcos}@lss.supelec.fr. This project is funded by region Île de France and Digiteo Research Park. This work has been partially presented in communication [START_REF] El Korso | Statistical resolution limit for multiple signals and parameters of interest[END_REF].

between two closely spaced signals 1 embedded in an additive noise that allows a correct resolvability/parameter estimation, is rising in several applications (especially in problems such as radar, sonar and spectral analysis [START_REF] Vantrees | Detection, Estimation and Modulation Theory[END_REF].)

The concept of the SRL was defined/used in several manners [START_REF] Smith | Statistical resolution limits and the complexified Cramér Rao bound[END_REF]- [START_REF] Amar | Fundamental limitations on the resolution of deterministic signals[END_REF], [START_REF] Cox | Resolving power and sensitivity to mismatch of optimum array processors[END_REF]- [START_REF] Delmas | Statistical resolution limits of DOA for discrete sources[END_REF], which could turn in it to a confusing concept. There exist essentially three approaches to define/obtain the SRL. (i) The first is based on the concept of mean null spectrum: assuming, for example, that two signals are parameterized by the frequencies f1 and f2, the Cox criterion [START_REF] Cox | Resolving power and sensitivity to mismatch of optimum array processors[END_REF] states that these sources are resolved, w.r.t. a given high-resolution estimation algorithm, if the mean null spectrum at each frequency f1 and f2 is lower than the mean of the null spectrum at the midpoint f 1 +f 2

2

. Another commonly used criterion, also based on the concept of the mean null spectrum, is the Sharman and Durrani criterion [START_REF] Sharman | Resolving power of signal subspace methods for fnite data lengths[END_REF], which states that two sources are resolved if the second derivative of the mean of the null spectrum at the midpoint f 1 +f 2 2 is negative. It is clear that the SRL based on the mean null spectrum is relevant to a specific highresolution algorithm (for some applications of these criteria one can see [START_REF] Cox | Resolving power and sensitivity to mismatch of optimum array processors[END_REF]- [START_REF] Abeida | Statistical performance of MUSIC-like algorithms in resolving noncircular sources[END_REF] and references therein.) (ii) The second approach is based on detection theory: the main idea is to use a hypothesis test to decide if one or two closely spaced signals are present in the set of the observations. Then, the challenge herein is to link the minimum separation, between two sources (for example, in terms of frequencies) that is detectable at a given SNR, to the probability of false alarm, P f a and/or to the probability of detection P d . In this spirit, Sharman and Milanfar [START_REF] Shahram | On the resolvability of sinusoids with nearby frequencies in the presence of noise[END_REF] have considered the problem of distinguishing whether the observed signal contains one or two frequencies at a given SNR using the Generalized Likelihood Ratio Test (GLRT). The authors have derived the SRL expressions w.r.t. P f a and P d in the case of real received signals, and unequal and unknown amplitudes and phases. In [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF], Liu and Nehorai have defined a statistical angular resolution limit using the asymptotic equivalence (in terms of number of observations) of the GLRT. The challenge was to determine the minimum angular separation, in the case of complex received signals, which allows to resolve two sources knowing the DOAs (Direction Of Arrival) of one of them for a given P f a and a given P d . Recently, Amar and Weiss [START_REF] Amar | Fundamental limitations on the resolution of deterministic signals[END_REF] have proposed to determine the SRL of complex sinusoids with nearby frequencies using the Bayesian approach for a given correct decision probability. (iii) The third approach is based on a estimation accuracy criteria independent of the estimation algorithm. Since the Cramér-Rao Bound (CRB) expresses a lower bound on the covariance matrix of any unbiased estimator, then it expresses also the ultimate estimation accuracy [START_REF] Liu | Cramér-Rao lower bounds for lowrank decomposition of multidimensional arrays[END_REF], [START_REF] Boyer | Deterministic asymptotic Cramér-Rao bound for the multidimensional harmonic model[END_REF]. Consequently, it could be used to describe/obtain the SRL. In this context, one distinguishes two main criteria for the SRL based on the CRB: 1) the first one was introduced by Lee in [START_REF] Lee | The Cramér-Rao bound on frequency estimates of signals closely spaced in frequency[END_REF] and states that : two signals are said to be resolvable w.r.t. the frequencies if the maximum standard deviation is less than twice the difference between f1 and f2. Assuming that the CRB is a tight bound (under mild/weak conditions), the standard deviation, σ f1 and σ f2 , of an unbiased estimator f = [ f1 f2]

T is given by CRB(f1)

and CRB(f2), respectively. Consequently, the SRL is defined, in the Lee criterion sense, as 2max CRB(f1), CRB(f2) . One can find some results and applications in [START_REF] Lee | The Cramér-Rao bound on frequency estimates of signals closely spaced in frequency[END_REF], [START_REF]The Cramér-Rao bound on frequency estimates of signals closely spaced in frequency (unconditional case)[END_REF] where this criterion is used to derive a matrix-based expression (i.e., without analytic inversion of the Fisher information matrix) of the SRL for the frequency estimates in the case of the conditional and unconditional signal source models. On the other hand, Dilaveroglu [START_REF] Dilaveroglu | Nonmatrix Cramér-Rao bound expressions for highresolution frequency estimators[END_REF] has derived a closed-form expression of the frequency resolution for the real and complex conditional signal source models. However, one can note that the coupling between the parameters, CRB(f1, f2) (i.e., the CRB for the cross parameters f1and f2), is ignored by this latter criterion.

2) To extend this, Smith [START_REF] Smith | Statistical resolution limits and the complexified Cramér Rao bound[END_REF] has proposed the following criterion: two signals are resolvable w.r.t. the frequencies if the difference between the frequencies, δ f , is greater than the standard deviation of the DOA difference estimation. Since, the standard deviation can be approximated by the CRB, then, the SRL, in the Smith criterion sense, is defined as the limit of δ f for which δ f < CRB (δ f ) is achieved. This means that, the SRL is obtained by solving the following implicit equation

δ 2 f = CRB (δ f ) = CRB(f1) + CRB(f2) -2CRB(f1, f2).
In [START_REF] Smith | Statistical resolution limits and the complexified Cramér Rao bound[END_REF] and [START_REF] Smith | Accuracy and resolution bounds for adaptive sensor array processing[END_REF], Smith has derived the SRL for two closely spaced sources in terms of DOA, each one modeled by one complex pole. In [START_REF] Delmas | Statistical resolution limits of DOA for discrete sources[END_REF], Delmas and Abeida have derived the SRL based on the Smith criterion for DOA of discrete sources under QPSK, BPSK and MSK model assumptions. More recently, Kusuma and Goyal [START_REF] Kusuma | On the accuracy and resolution of powersumbased sampling methods[END_REF] have derived the SRL based on the Smith criterion in sampling estimation problems involving a powersum series.

It is important to note that all the criteria listed before take into account only one parameter of interest per signal. Consequently, all the criteria listed before can not be applied to the aforementioned the multidimensional harmonic model. To the best of our knowledge, no results are available on the SRL for multiple parameters of interest per signal. The goal of this paper is to fill this lack by proposing and deriving the so-called Multidimensional SRL (MSRL) for the multidimensional harmonic retrieval model.

More precisely, in this paper, the MSRL for multiple parameters of interest per signal using a hypothesis test is derived. This choice is motivated by the following arguments: i) the hypothesis test approach is not specific to a certain high-resolution algorithm (unlike the mean null spectrum approach), ii) in this paper, we link the asymptotic MSRL based on the hypothesis test approach to a new extension of the MSRL based on the CRB approach. Furthermore, we show that the MSRL based on the CRB approach is equivalent to the MSRL based on the hypothesis test approach for a fixed couple (P f a , P d ), and iii) the hypothesis test is shown to be asymptotically an uniformly most powerful test which is the strongest statement of optimality that one could expect to obtain [START_REF] Scharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF].

The paper is organized as follows. We first begin by introducing the multidimensional harmonic model, in Section II. Then, based on this model, we obtain the MSRL based on the hypothesis test and on the CRB approach. The link between theses two MSRLs is also described in Section III. Section IV is devoted to the derivation of the MSRL closed-form expression, where, as a by product the exact closed-form expressions of the CRB for the multidimensional retrieval model is derived (note that to the best of our knowledge, no exact closed-form expressions of the CRB for such model is available in the literature). Furthermore, theoretical and numerical analysis are given in the same section. Finally, conclusions are given in Section V.

Glossary of notation

The following notations are used through the paper. Column vectors, matrices and multi-way arrays are represented by lowercase bold letters (a, . . . ), upper-case bold letters (A, . . . ) and bold calligraphic letters (A, . . . ), whereas

• R and C denote the body of real and complex values, respectively,

• R D 1 ×D 2 ו••×D I and C D 1 ×D 2 ו••×D I denote
the real and complex multi-way arrays (also called tensors) body of dimension

D1 × D2 × • • • × DI , respectively, • j = the complex number √ -1, • IQ = the identity matrix of dimension Q, • 0Q 1 ×Q 2 = the Q1 × Q2 matrix filled by zeros, • [a] i = the i th element of the vector a, • [A] i 1 ,i 2 = the i th
1 row and the i th 2 column element of the matrix A,

• [A] i 1 ,i 2 ,...,i N = the (i1, i2, . . . , iN ) th entry of the multi-way array A, • [A] i,p:q = the row vector containing the (q -p + 1) elements

[A] i,k , where k = p, . . . , q, • [A] p:q,k = the column vector containing the (q -p+1) elements

[A] i,k , where i = p, . . . , q, • the derivative of vector a w.r.t. to vector b is defined as follows:

∂a ∂b i,j = ∂[a] i ∂[b] j , • A T = the transpose of the matrix A, • A * = the complex conjugate of the matrix A, • A H = (A * ) T ,
• tr {A} = the trace of the matrix A,

• det {A} = the determinant of the matrix A,

• {a} = the real part of the complex number a,

• E {a} = the expectation of the random variable a,

• a 2 = 1 L L t=1 [a] 2
t denotes the normalized norm of the vector a (in which L is the size of a),

• sgn (a) = 1 if a ≥ 0 and -1 otherwise.

• diag(a) is the diagonal operator which forms a diagonal matrix containing the vector a on its diagonal, • vec(.) is the vec-operator stacking the columns of a matrix on top of each other,

• stands for the Hadamard product, • ⊗ stands for the Kronecker product, • • denotes the multi-way array outer-product (recall that for a given multi-way arrays

A ∈ C A 1 ×A 2 ו••×A I and B ∈ C B 1 ×B 2 ו••×B J ,
the result of the outerproduct of A and B denoted by

C A 1 ו••×A I ×B 1 ו••×B J is given by [C] a 1 ,...,a I ,b 1 ,...,b J = [A • B] a 1 ,...,a I ,b 1 ,...,b J = [A] a 1 ,...,a I [B] b 1 ,...,b J ).

II. MODEL SETUP

In this section we introduce the multidimensional harmonic retrieval model in the multi-way array form (also known as tensor form [START_REF] Westin | A Tensor Framework For Multidimensional Signal Processing[END_REF]). Then, we use the PARAFAC (PARallel FACtor) decomposition to obtain a vector form of the observation model. This vector form will be used to derive the closed-form expression of the MSRL.

Let us consider a multidimensional harmonic model consisting of the superposition of two harmonics each one of dimension P contaminated by an additive noise. Thus, the observation model is given as follows [START_REF] Nion | Tensor algebra and multi-dimensional harmonic retrieval in signal processing for MIMO radar[END_REF], [START_REF] Gershman | Space-time processing for MIMO communications[END_REF], [START_REF] Boyer | Deterministic asymptotic Cramér-Rao bound for the multidimensional harmonic model[END_REF], [START_REF] Haardt | Simultaneous Schur decomposition of several nonsymmetric matrices to achieve automatic pairing in multidimensional harmonic retrieval problems[END_REF]- [START_REF] Pesavento | Multidimensional rank reduction estimator for parametric MIMO channel models[END_REF]:

[Y(t)] n 1 ,...,n P = [X (t)] n 1 ,...,n P + [N (t)] n 1 ,...,n P , (1)
t = 1, . . . , L, and np = 0, . . . , Np -1,

where Y(t), X (t) and N (t) denote the noisy observation, the noiseless observation and the noise multi-way array at the t th snapshot. The number of snapshots and the number of sensors on each array are denoted by L and (N1, . . . , NP ), respectively. The noiseless observation multi-way array can be written as follows 2 [27],

[31]- [START_REF] Pesavento | Multidimensional rank reduction estimator for parametric MIMO channel models[END_REF]:

[X (t)] n 1 ,...,n P = 2 m=1 sm(t) P p=1 e jω (p) m np , (2)
where ω

(p)
m and sm(t) denote the m th frequency viewed along the p th dimension or array and the m th complex signal source, respectively. Furthermore, the signal source is given by sm(t) = αm(t)e jφ m (t) where αm(t) and φ m (t) denote the real positive amplitude and the phase for the m th signal source at the t th snapshot, respectively.

Since,

P p=1 e jω (p) m np = a(ω (1)
m ) • a(ω (2)
m ) • • • • • a(ω (P ) m ) n 1 ,n 2 ,...,n P
, where a(.) is a Vandermonde vector defined as

a(ω (p) m ) = 1 e jω (p) m . . . e j(Np-1)ω (p) m T ,
then, the multi-way array X (t) follows a PARAFAC decomposition [START_REF] Nion | A PARAFAC-based technique for detection and localization of multiple targets in a MIMO radar system[END_REF], [START_REF] Harshman | Foundations of the PARAFAC procedure: Models and conditions for an "explanatory" multi-modal factor analysis[END_REF]. Consequently, the noiseless observation multi-way array can be rewritten as follows:

X (t) = 2 m=1 sm(t) a(ω (1) m ) • a(ω (2) m ) • • • • • a(ω (P ) m ) . (3)
First, let us vectorize the noiseless observation as follows:

vec (X (t)) = [[X (t)] 0, 0,...,0 . . . [X (t)] N 1 -1, 0,...,0 (4)
[X (t)] 0, 1,...,0 . . . [X (t)] N 1 -1, N 2 -1,...,N P -1 ] T .
Thus, the full noise-free observation vector is given by

x = vec T (X (1)) vec T (X (2)) . . . vec T (X (L)) T .
Second, and in the same way, we define y, the noisy observation vector, and n, the noise vector, by the concatenation of the proper multi-way array's entries, i.e.,

y = vec T (Y(1)) vec T (Y(2)) . . . vec T (Y(L)) T = x + n.
(5) Consequently, in the following, we will consider the observation model in [START_REF] Wong | Uni-vector-sensor ESPRIT for multisource azimuth, elevation, and polarization estimation[END_REF]. Furthermore, the unknown parameter vector is given by

ξ = ω T ρ T T , (6)
where ω denotes the unknown parameter vector of interest, i.e., containing all the unknown frequencies

ω = ω (1) T . . . ω (P ) T T , in which
ω (p) = ω (p) 1 ω (p) 2 T . (7)
Whereas ρ contains the unknown nuisance/unwanted parameters vector, i.e., characterizing the noise covariance matrix and/or amplitude and phase of each source (e.g., in the case of a covariance noise matrix equal to σ 2 ILN 1 ...N P and unknown deterministic amplitudes and phases, the unknown nuisance/unwanted parameters vector ρ is given by ρ = α1(1) . . . α2(L) φ 1 (1) . . . φ 2 (L) σ 2 T ).

In the following we conduct a hypothesis test formulation on the observation model ( 5) to derive our MSRL expression in the case of two sources.

III. DETERMINATION OF THE MSRL FOR TWO SOURCES

A. Hypothesis test formulation

Resolving two closely spaced sources, with respect to their parameters of interest, can be formulated as a binary hypothesis test [START_REF] Shahram | On the resolvability of sinusoids with nearby frequencies in the presence of noise[END_REF]- [START_REF] Amar | Fundamental limitations on the resolution of deterministic signals[END_REF] (for the special case of P = 1). To determine the MSRL (i.e., P ≥ 1), let us consider the hypothesis H0 which represents the case where the two emitted signal sources are combined into one signal, i.e., the two sources have the same parameters (this hypothesis is described by ∀p ∈ [1 . . . P ], ω

(p) 1 = ω (p)
2 ), whereas the hypothesis H1 embodies the situation where the two signals are resolvable (the latter hypothesis is described by ∃p ∈ [1 . . . P ], such that ω

(p) 1 = ω (p)
2 ). Consequently, one can formulate the hypothesis test, as a simple one-sided binary hypothesis test as follows:

H0 : δ = 0, H1 : δ > 0, (8)
where the parameter δ is the so-called MSRL which indicates us in which hypothesis our observation model belongs. Thus, the question addressed below is how can we define the MSRL δ such that all the P parameters of interest are taken into account? A natural idea is that δ reflects a distance between the P parameters of interest. Let the MSRL denotes the l1 norm3 between two sets containing the parameters of interest of each source (which is the naturally used norm, since in the mono-parameter frequency case that we extend here, the SRL is defined as δ = f1 -f2 [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF], [START_REF] Amar | Fundamental limitations on the resolution of deterministic signals[END_REF], [START_REF] Stoica | Spectral Analysis of Signals[END_REF]). Meaning that, if we denote these sets as C1 and C2 where Cm = ω

(1)

m , ω (2)
m , . . . , ω

(P ) m
, m = 1, 2, thus, δ can be defined as

δ P p=1 ω (p) 2 -ω (p) 1 . (9)
First, note that the proposed MSRL describes well the hypothesis test (8) (i.e., δ = 0 means that the two emitted signal sources are combined into one signal and δ = 0 the two signals are resolvable). Second, since the MSRL δ is unknown, it is impossible to design an optimal detector in the Neyman-Pearson sense. Alternatively, the Generalized Likelihood Ratio Test (GLRT) [START_REF] Scharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF], [START_REF] Kay | Fundamentals of Statistical Signal Processing : Detection Theory[END_REF] is a well known approach appropriate to solve such a problem. To conduct the GLRT on [START_REF] Nion | A PARAFAC-based technique for detection and localization of multiple targets in a MIMO radar system[END_REF], one has to express the probability density function (pdf) of ( 5) w.r.t. δ. Assuming (without loss of generality) that ω

(1) 1 > ω (1)
2 , one can notice that ξ is known if and only if δ and ϑ ω

(1) 2 ω (2) T . . . ω (P ) T T are fixed (i.e., there is a one to one mapping between δ, ϑ and ξ). Consequently, the pdf of (5) can be described as p(y|δ, ϑ). Now, we are ready to conduct the GLRT for this problem:

LG(y) = max δ,ϑ 1 p(y|δ, ϑ1, H1) max ϑ 0 p(y|ϑ0, H0) = p(y| δ, θ1, H1) p(y| θ0, H0)

H 1 ≷ H 0 ς , (10)
where δ, θ1 and θ0 denote the Maximum Likelihood Estimates (MLE) of δ under H1, the MLE of ϑ under H1 and the MLE of ϑ under H0, respectively, and where ς denotes the test threshold. From [START_REF] Gershman | Space-time processing for MIMO communications[END_REF], one obtains

TG(y) = Ln LG(y) H 1 ≷ H 0 ς = Lnς , (11)
in which Ln denotes the natural logarithm.

B. Asymptotic equivalence of the MSRL

Finding the analytical expression of TG(y) in ( 11) is not tractable. This is mainly due to the fact that the derivation of δ is impossible since from (2) one obtains a multimodal likelihood function [START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF]. Consequently, in the following, and as in 4 [14], we consider the asymptotic case (in terms of the number of snapshots). In [38, eq (6C.1)] it has been proven that, for a large number of snapshots, the statistic TG(y) follows a chi-square pdf under H0 and H1 given by

TG(y) ∼ χ 2 1 under H0, χ 2 1 (κ (P f a , P d )) under H1, (12)
where χ 2 1 and χ 2 1 (κ (P f a , P d )) denote the central chi-square and the noncentral chi-square pdf with one degree of freedom, respectively. P f a and P d are, respectively, the probability of false alarm and the probability of detection of the test [START_REF] Nion | A PARAFAC-based technique for detection and localization of multiple targets in a MIMO radar system[END_REF]. In the following, CRB(δ) denotes the CRB for the parameter δ where the unknown vector parameter is given by δ ϑ T T . Consequently, assuming that CRB(δ) exists (under H0 and H1), is well defined (see Section IV for the necessary 5 and sufficient conditions) and is a tight bound (i.e., achievable under quite general/weak conditions [START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF], [START_REF] Renaux | On the high SNR conditional maximum-likelihood estimator full statistical characterization[END_REF]), thus the noncentral parameter κ (P f a , P d ) is given by [38, p.239]

κ (P f a , P d ) = δ 2 (CRB(δ)) -1 . (13)
On the other hand, one can notice that the noncentral parameter κ (P f a , P d ) can be determined numerically by the choice of P f a and P d [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF], [START_REF] Scharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF] as the solution of

Q -1 χ 2 1 (P f a ) = Q -1 χ 2 1 (κ (Pfa,Pd)) (P d ), (14)
in which

Q -1 χ 2 1 ( ) and Q -1 χ 2
1 (κ (P f a ,P d )) ( ) are the inverse of the right tail of the χ 2 1 and χ 2 1 (κ (P f a , P d )) pdf starting at the value . Finally, from ( 13) and ( 14) one obtains 6

δ = κ(P f a , P d ) CRB(δ), ( 15
)
where κ(P f a , P d ) = κ (P f a , P d ) is the so-called translation factor [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF] which is determined for a given probability of false alarm and 4 Note that, due to the specific definition of the SRL in [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF] (i.e., using the same notation as in [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF], δ = cos(u T 1 u 2 )) and the restrictive assumption in [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF] (u 1 and u 2 belong to the same plan), the SRL as defined in [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF] cannot be used in the multidimensional harmonic context. 5 One of the necessary conditions regardless the noise pdf is that ω One can notice that increasing P d or decreasing P f a has the effect to increase the value of the translation factor κ. This is expected since increasing P d or decreasing P f a leads to a more selective decision [START_REF] Scharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF], [START_REF] Kay | Fundamentals of Statistical Signal Processing : Detection Theory[END_REF].

probability of detection (see Fig. 1 for the behavior of the translation factor versus P f a and P d ).

Result 1: The asymptotic MSRL for model [START_REF] Wong | Uni-vector-sensor ESPRIT for multisource azimuth, elevation, and polarization estimation[END_REF] in the case of P parameters of interest per signal (P ≥ 1) is given by δ which is the solution of the following equation:

δ 2 -κ 2 (P f a , P d ) (A direct + Across) = 0, (16)
where A direct denotes the contribution of the parameters of interest belonging to the same dimension as follows

A direct = P p=1 CRB(ω (p) 1 ) + CRB(ω (p) 2 ) -2CRB(ω (p) 1 , ω (p) 2 )
) + CRB(ω

(p) 2 , ω (p ) 2
)

-2CRB(ω (p) 1 , ω (p ) 2
)),

in which gp = sgn ω (p) 1 -ω (p) 2 . Proof: see Appendix A.
Remark 1: It is worth noting that, the hypothesis test ( 8) is a binary one-sided test and that the MLE used is an unconstrained estimator. Thus, one can deduce that the GLRT, used to derive the asymptotic MSRL [START_REF] Liu | Statistical angular resolution limit for point sources[END_REF], [START_REF] Kay | Fundamentals of Statistical Signal Processing : Detection Theory[END_REF]: i) is the asymptotically uniformly most powerful test among all invariant statistical tests, and ii) has an asymptotic Constant False-Alarm Rate (CFAR). Which is, in the asymptotic case, considered as the strongest statement of optimality that one could expect to obtain [START_REF] Scharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF].

Remark 2: • Existence of the MSRL : It is natural to assume that the CRB is a non-increasing (i.e., decreasing or constant) function on R + w.r.t. δ since it is more difficult to estimate two closely-spaced signals than two largely-spaced ones. In the same time the left hand side of ( 15) is a monotonically increasing function w.r.t. δ on R + . Thus for a fixed couple (P f a , P d ), the solution of the implicit equation given by (15) always exists. However, theoretically, there is no assurance that the solution of equation ( 15) is unique. • Note that, in practical situation, the case where CRB(δ) is not a function of δ is important since in this case, CRB(δ) is constant w.r.t. δ and thus the solution of ( 15) exists and is unique (see Section IV).

In the following subsection, we study the explicit effect of this so-called translation factor.

C. The relationship between the MSRL based on the CRB and the Hypothesis Test approaches

In this subsection, we link the asymptotic MSRL (derived using the hypothesis test approach, see Result 1) to a new proposed extension of the SRL based on the Smith criterion [START_REF] Smith | Statistical resolution limits and the complexified Cramér Rao bound[END_REF]. First, we recall that the Smith criterion defines the SRL in the case of P = 1 only. Then, we extend this criterion to P 1 (i.e., the case of the multidimensional harmonic model). Finally, we link the MSRL based on the hypothesis test approach (see Result 1) to the MSRL based on the CRB approach (i.e., the extended SRL based on the Smith Criterion).

1) The Smith criterion: Since the CRB expresses a lower bound on the covariance matrix of any unbiased estimator, then it expresses also the ultimate estimation accuracy. In this context, Smith proposed the following criterion for the case of two source signals parameterized each one by only one frequency [START_REF] Smith | Statistical resolution limits and the complexified Cramér Rao bound[END_REF] : two signals are resolvable if the difference between their frequency, δ ω

(1) = ω (1) 2 -ω (1)
1 , is greater than the standard deviation of the frequency difference estimation. Since, the standard deviation can be approximated by the CRB, then, the SRL, in the Smith criterion sense, is defined as the limit of δ ω [START_REF] El Korso | Statistical resolution limit for multiple signals and parameters of interest[END_REF] for which δ ω (1) < CRB δ ω (1) is achieved. This means that, the SRL is the solution of the following implicit equation

δ 2 ω (1) = CRB δ ω (1) .
2) The extension of the Smith criterion to the case of P 1: Based on the above framework, a straightforward extension of the Smith criterion to the case of P 1 for the multidimensional harmonic model is as follows : two multidimensional harmonic retrieval signals are resolvable if the distance between C1 and C2, is greater than the standard deviation of the δCRB estimation. Consequently, assuming that the CRB exists and is well defined, the MSRL δCRB is given as the solution of the following implicit equation

δ 2 CRB = CRB (δCRB) s.t. δCRB = P p=1 ω (p) 2 -ω (p) 1 . ( 17
)
3) Comparison and link between the MSRL based on the CRB approach and the MSRL based on the hypothesis test approach: The MSRL based on the hypothesis test approach is given as the solution of

δ = κ(P f a , P d ) CRB (δ), s.t. δ = P p=1 ω (p) 2 -ω (p) 1
, whereas the MSRL based on the CRB approach is given as the solution of ( 17). Consequently, one has the following result: Result 2: Upon to a translation factor, the asymptotic MSRL based on the hypothesis test approach (i.e., using the binary one-sided hypothesis test given in [START_REF] Nion | A PARAFAC-based technique for detection and localization of multiple targets in a MIMO radar system[END_REF]) is equivalent to the proposed MSRL based on the CRB approach (i.e., using the extension of the Smith criterion). Consequently, the criterion given in ( 17) is equivalent to an asymptotically uniformly most powerful test among all invariant statistical tests for κ (P f a , P d ) = 1 (see Fig. 2 for the values of (P f a , P d ) such that κ(P f a , P d ) = 1).

The following section is dedicated to the analytical computation of closed-form expression of the MSRL. In Subsection A we introduce the assumptions used to compute the MSRL in the case of a Gaussian random noise and orthogonal waveforms. Then we derive non matrix closed-form expressions of the CRB in Subsection B (note that to the best of our knowledge, no closed-form expressions of the CRB for such model is available in the literature). In Subsection C and thanks to these expressions, the MSRL wil be deduced using [START_REF] Vantrees | Detection, Estimation and Modulation Theory[END_REF]. Finally, the MSRL analysis is given in Subsection D.

IV. MSRL CLOSED FORM EXPRESSION

In Section III we have defined the general model of the multidimensional harmonic model. To derive a closed form expression of the MSRL, we need more assumptions on the covariance noise matrix and/or on the signal sources.

A. Assumptions

• The noise is assumed to be a complex circular white Gaussian random process i.i.d. with zero-mean and unknown variance σ 2 ILN 1 ...N P . • We consider a multidimensional harmonic model due to the superposition of two harmonics each of them of dimension P ≥ 1. Furthermore, for sake of simplicity and clarity, the sources have been assumed known and orthogonal (e.g., [START_REF] Nion | A PARAFAC-based technique for detection and localization of multiple targets in a MIMO radar system[END_REF], [START_REF] Li | Maximum likelihood angle estimation for signals with known waveforms[END_REF]). In this case, the unknown parameter vector is fixed and does not grow with the number of snapshots. Consequently, the CRB is an achievable bound [START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF]. • Each parameter of interest w.r.t. to the first signal, ω (p) 1 p = 1 . . . P , can be as close as possible to the parameter of interest w.r.t. to the second signal ω (p) 2 p = 1 . . . P , but not equal. This is not really a restrictive assumption, since in most applications, having two or more identical parameters of interest is a zero probability event [10, p. 53].

Under these assumptions, the joint probability density function of the noisy observations y for a given unknown deterministic parameter vector ξ is as follows:

p( y| ξ) = L t=1 p( vec (Y(t))| ξ) = 1 (πσ 2 ) LN e -1 σ 2 (y-x) H (y-x) ,
where N = P p=1 Np. The multidimensional harmonic retrieval model with known sources is considered herein, and thus, the parameter vector is given by

ξ = ω T σ 2 T , (18)
where

ω = ω (1) T . . . ω (P ) T T , in which ω (p) = ω (p) 1 ω (p) 2 T . (19)

B. CRB for the multidimensional harmonic model with orthogonal known signal sources

The Fisher Information Matrix (FIM) of the noisy observations y w.r.t. a parameter vector ξ is given by [START_REF] Cramér | Mathematical Methods of Statistics[END_REF]

FIM(ξ) = E ∂lnp( y| ξ) ∂ξ ∂lnp( y| ξ) ∂ξ H .
For a complex circular Gaussian observation model, the i th , k th element of the FIM for the parameter vector ξ is given by [START_REF] Stoica | Spectral Analysis of Signals[END_REF] [

FIM(ξ)] i,k = LN σ 4 ∂σ 2 ∂ [ξ] i ∂σ 2 ∂ [ξ] k + 2 σ 2 ∂x H ∂ [ξ] i ∂x ∂ [ξ] k (i, k) = {1, . . . , 2P + 1} 2 .
(20) Consequently, one can state the following Lemma. Lemma 1: The FIM for the sum of two P -order harmonic models with orthogonal known sources, has a block diagonal structure and is given by

FIM(ξ) = 2 σ 2 Fω 02P ×1 01×2P × , (21)
where, the (2P ) × (2P ) matrix Fω is also a block diagonal matrix given by

Fω = LN (∆ ⊗ G), (22)
in which ∆ = diag{ α1 2 , α2 2 } where αm = αm(1) . . . αm(L)

T for m ∈ {1, 2} , (23)
and

[G] k,l = (2N k -1)(N k -1) 6
for k = l,

(N k -1)(N l -1) 2
for k = l.

Proof: see Appendix B. After some calculation and using Lemma 1, one can state the following result.

Result 3: The closed-form expressions of the CRB for the sum of two P -order harmonic models with orthogonal known signal sources are given by

CRB(ω (p) m ) = 6 LN SNRm Cp, m ∈ {1, 2} , (24)
where SNRm = αm 2 σ 2

denotes the Signal-to-Noise Ratio (SNR) of the m th source and where

Cp = Np(1 -3VP ) + 3VP + 1 (Np + 1)(N 2 p -1)
in which VP = 1

1 + 3 P p=1 Np-1 Np+1
.

Furthermore, the cross-terms are given by

CRB(ω (p) m , ω (p ) m ) = 0 for m = m , - 6
LN SNRm Cp,p for m = m and p = p ,

where Cp,p = 3VP (Np + 1)(N p + 1) .

Proof: see Appendix C.

C. MSRL derivation

Using the previous result, one obtains the unique solution of ( 16), thus, the MSRL for model ( 1) is given by the following result:

Result 4: The MSRL for the sum of P -order harmonic models with orthogonal known signal sources, is given by

δ = 6 LN ESNR     P p=1 Cp - P p, p = 1 p = p gpg p Cp,p     , (26)
where the so-called Extended SNR is given by ESNR

= SNR 1 SNR 2 SNR 1 +SNR 2 . Proof: see Appendix D.

D. Numerical Analysis

Taking advantage of the latter result, one can analyze the MSRL given by (26):

• First, from Fig. 3 note that the numerical solution of the MSRL based on ( 12) is in good agreement with the analytical expression of the MSRL [START_REF] Dilaveroglu | Nonmatrix Cramér-Rao bound expressions for highresolution frequency estimators[END_REF], which validate the closed-form expression given in [START_REF] Dilaveroglu | Nonmatrix Cramér-Rao bound expressions for highresolution frequency estimators[END_REF]. On the other hand, one can notice that, for P d = 0.37 and P f a = 0.1 the MSRL based on the CRB is exactly equal to the MSRL based on hypothesis test approach derived in the asymptotic case. From the case P d = 0.49 and P f a = 0.3 or/and P d = 0.32 and P f a = 0.1, one can notice the influence of the translation factor κ(P f a , P d ) on the MSRL.

• The MSRL 7 is O 1 ESNR which is consistent with some previous results for the case P = 1 (e.g., [START_REF] Shahram | On the resolvability of sinusoids with nearby frequencies in the presence of noise[END_REF], [START_REF] Amar | Fundamental limitations on the resolution of deterministic signals[END_REF], [START_REF] Delmas | Statistical resolution limits of DOA for discrete sources[END_REF]).

• From ( 26) and for a large number of sensors

N1 = N2 = • • • = NP = N 1, one obtains a simple expression δ = 12 LN P +1 ESNR P 1 + 3P , meaning that, the SRL is O 1 N P +1 . • Furthermore, since P ≥ 1, one has (P + 1) (3P + 1) P (3P + 4) < 1,
and consequently, the ratio between the MSRL of a multidimensional harmonic retrieval with P parameters of interest, denoted by δP and the MSRL of a multidimensional harmonic retrieval with P + 1 parameters of interest, denoted by δP +1, is given by δP +1 δP = (P + 1) (3P + 1) meaning that the MSRL for P + 1 parameters of interest is less than the one for P parameters of interest (see Fig. 4). This, can be explained by the estimation additional parameter and also by an increase of the received noisy data thanks to the additional dimension. One should note that this property is proved theoretically thanks to (27) using the assumption of an equal and large number of sensors. However, from Fig. 4 we notice that, in practice, this can be verified even for a small number of sensors (e.g., in Fig. 4 one has 3 ≤ Np ≤ 5 for p = 3, . . . , 6). • Furthermore, since

N P (3P + 4) , (27)
4 LN P +1 ESNR ≤ δP < δP -1 < • • • < δ1
one can note that, the SRL is lower bounded by 4 LN P +1 ESNR . • One can address the problem of finding the optimal distribution of power sources making the SRL the smallest as possible (s.t. the constraint of constant total source power). In this issue, one can state the following corollary: Corollary 1: The optimal power's source distribution that ensures the smallest MSRL is obtained only for the equi-powered sources case. Proof: see Appendix E. This result was observed numerically for P = 1 in [START_REF] Shahram | On the resolvability of sinusoids with nearby frequencies in the presence of noise[END_REF] (see Fig. 5 for the multidimensional harmonic model). Moreover, it has been shown also by simulation for the case P = 1 that the so-called maximum likelihood breakdown (i.e., when the mean square error of the MLE increases rapidly) occurs at higher SNR in the case of different power signal sources than in the case of equi-powered signal sources [START_REF] Abramovich | Statistical nonidentifiability of close emitters: Maximum-likelihood estimation breakdown[END_REF]. The authors explained it by the fact that one source grabs most of the total power, then, this latter will be estimated more accurately, whereas the second one, will take an arbitrary parameter estimation which represents an outlier. • In the same way, let us consider the problem of the optimal placement of the sensors8 N1, . . . , NP , making the minimum MSRL s.t. the constraint that the total number of sensors is constant (i.e., N total = P p=1 Np in which we suppose that N total is a multiple of P ). Corollary 2: If the total number of sensors N total , is a multiple of P , then an optimal placement of the sensors that ensure the lowest MSRL is (see Fig. 6 and Fig. 7)

N1 = • • • = NP = N total P . (28)
Proof: see Appendix F. Remark 3: Note that, in the case where N total is not a multiple of P , one expects that the optimal MSRL is given in the case where the sensors distribution approaches the equi-sensors distribution situation given in corollary 3. Fig. 7 confirms that (in the case of P = 3, N1 = 8 and a total number of sensors N = 22). From Fig. 7, one can notice that the optimal distribution of the number of sensors corresponds to N2 = N3 = 7 and N1 = 8 which is the nearest situation to the equi-sensors distribution.

V. CONCLUSION

In this paper, we have derived the Multidimensional Statistical Resolution Limit (MSRL) for the multidimensional harmonic retrieval model. Toward this end, we have extended the concept of Statistical Resolution Limit (SRL) to multiple parameters of interest per signal. First, we have used a hypothesis test approach. The applied test is shown to be asymptotically an uniformly most powerful test which is the strongest statement of optimality that one could hope to obtain. Second, we have linked the asymptotic MSRL based on the hypothesis test approach to a new extension of the SRL based Fig. 6. The MSRL vs. N 1 and N 2 in the case of P = 3 and a total number of sensors N total = 21. One can notice that the optimal distribution of the number of sensors (which corresponds to the lowest SLR) corresponds to

N 1 = N 2 = N 3 = N total 3
as predicted by [START_REF] Kusuma | On the accuracy and resolution of powersumbased sampling methods[END_REF]. on the Cramér-Rao bound approach. Using the Cramér-Rao bound and a proper change of variable formula, closed-form expression of the MSRL are given. Finally, note that the concept of the MSRL can be used to optimize, for example, the waveform and/or the array geometry for a specific problem.

APPENDIX A:

The proof of Result 1 Appendix A.1 : In this appendix, we derive the MSRL using the l1 norm.

From CRB(ξ)

where ξ = [ω T ρ T ] T in which ω = ω (1) 1 ω (1) 2 ω (2) 1 ω (2) 2
. . . ω Consequently, after some calculus, one obtains CRB(δ)

∆ = CRB ξ 1,1 = h T CRB(ω)h = 2P p=1 2P p =1 [h] p [h] p [CRB(ω)] p,p = P p=1 P p =1 gpg p [CRB(ξ)] 2p,2p + [CRB(ξ)] 2p-1,2p -1 - [CRB(ξ)] 2p,2p -1 -[CRB(ξ)] 2p-1,2p ∆ = A direct + Across, (30)
where A direct = P p=1 CRB(ω ) . Finally, using (30) one obtains [START_REF] Vantrees | Detection, Estimation and Modulation Theory[END_REF].

Appendix A.2 : In this part, we derive the MSRL using the l k norm for a given integer k 1. The aim of this part is to support the footnote number 1, which stays that using the l1 norm computing the MSRL using the l1 norm is for the calculation convenience.

Once again, from CRB(ξ), one can deduce CRB( ξk ) where ξk = g k (ξ) = [δ(k) ϑ T ] T in which the distance between C1 and C2 using the l k norm is given by δ(k)

gp(k) = ∂ P p =1 ω (p ) 1 -ω (p ) 2 2k 1/k ∂ω (p) 1 = 1 k p i=1 ω (i) 1 -ω (i) 2 2k 1 k -1 ∂ ω (i) 1 -ω (i) 2 2k ∂ω (i) 1 = sgn(ω (p) 1 -ω (p) 2 )   P p =1 ω (p ) 1 -ω (p ) 2 2k   1 k -1 ω (p) 1 -ω (p) 2 2(k-1) (31)
= sgn(ω

(p) 1 -ω (p) 2 )δ 1-k δ k-1 p .
Again, by using the change of variable formula [START_REF] Scharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF], one has

CRB( ξk ) = h T k CRB(ω)h k × × I .
On the other hand, one has

a H (ω (p) m )a(ω (p) m ) = Np, (34)
whereas a H (ω (p) m )a (ω (p) m ) = Np (Np -1) 2 (35)
and

a H (ω (p) m )a (ω (p) m ) = Np (2Np -1) (Np -1) 6
Finally, assuming known orthogonal wavefronts [START_REF] Li | Maximum likelihood angle estimation for signals with known waveforms[END_REF] (i.e., s H m sm = 0) and replacing [START_REF] Vantrees | Detection, Estimation and Modulation theory: Optimum Array Processing[END_REF] and ( 34 where αm = αm(1) . . . αm(L) for m ∈ {1, 2} . Consequently, using [START_REF] Golub | Matrix Computations[END_REF], Fω can be expressed as a block diagonal matrix

Fω = J1 0 0 J2 , (37)
where each P × P block Jm is defined by

Jm = L αm 2 N G, (38)
where

G =       (N 1 -1)(2N 1 -1) 6 (N 1 -1)(N 2 -1) 4
. . .

(N 1 -1)(N P -1) 4 (N 2 -1)(N 1 -1) 4 (N 2 -1)(2N 2 -1) 6
. . .

(N 2 -1)(N P -1) 4
. . . . . . . . . . . .

(N P -1)(N 1 -1) 4
(N 2 P -1)(N 2 -1) 4

. . . (N P -1)(2N P -1)

6       .
Consequently, from [START_REF] Stoica | Spectral Analysis of Signals[END_REF] and [START_REF] Kay | Fundamentals of Statistical Signal Processing : Detection Theory[END_REF] one obtains [START_REF]The Cramér-Rao bound on frequency estimates of signals closely spaced in frequency (unconditional case)[END_REF].

APPENDIX C Proof of result 3

Using (22) one obtains

CRB(ω) = σ 2 2 F -1 ω = σ 2 2LN ∆ -1 ⊗ G -1 (39)
where

∆ -1 = diag 1 α 1 2 , 1 α 2 2
. In the following, we give a closed-form expression of G -1 . One can notice that the matrix G has a particular structure such that it can be rewritten as the sum of a diagonal matrix and of a rank-one matrix:

G = Q + γγ T where Q = 1 12 diag{N 2 1 -1, . . . , N 2 P -1} and γ = 1 2 [N1-1, . . . , NP -1]
T . Thanks to this particular structure, an analytical inverse of G can easly be obtained. Indeed, using the matrix inversion Lemma

G -1 = (Q + γγ T ) -1 = Q -1 - Q -1 γγ T Q -1 1 + γ T Q -1 γ . (40)
A straightforward calculus leads to the following results,

Q -1 γγ T Q -1 = (41) 36       1 (N 1 +1) 2 1 (N 1 +1)(N 2 +1) . . . 1 (N 1 +1)(N P +1) 1 (N 2 +1)(N 1 +1) 1 (N 2 +1) 2 . . . 1 (N 2 +1)(N P +1) . . . . . . . . . . . . 1 (N P +1)(N 1 +1) 1 (N P +1)(N 2 +1) . . . 1 (N P +1) 2       , and γ T Q -1 γ = 3 P p=1 Np -1 Np + 1 . (42)
Consequently, replacing ( 41) and ( 42) into [START_REF] Renaux | On the high SNR conditional maximum-likelihood estimator full statistical characterization[END_REF], one obtains

G -1 k,l =    12 Np(1-3V P )+3V P +1 (Np+1)(N 2 p -1)
for k = l,

- 36V P (Np+1)(N p +1) for k = l, (43)
where VP = 1 + 3 P p=1 Np-1 Np+1 -1

. Finally, replacing (43) into (39) one finishes the proof.

) + CRB(ω

Consequently, replacing [START_REF] Abramovich | Statistical nonidentifiability of close emitters: Maximum-likelihood estimation breakdown[END_REF] and ( 45) into [START_REF] Vantrees | Detection, Estimation and Modulation Theory[END_REF], one finishes the proof.

APPENDIX E Proof of corollary 1

In this Appendix we minimize the MSRL under the constraint SNR1 + SNR2 = SNR total (where SNR total is a real fixed value).

Since, the term

2 .Fig. 1 .

 21 Fig.1. The translation factor κ vs. the probability of detection P d and P f a . One can notice that increasing P d or decreasing P f a has the effect to increase the value of the translation factor κ. This is expected since increasing P d or decreasing P f a leads to a more selective decision[START_REF] Scharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF],[START_REF] Kay | Fundamentals of Statistical Signal Processing : Detection Theory[END_REF].

,

  and where Across is the contribution of the cross terms between distinct dimension given by

Fig. 2 .

 2 Fig. 2. All values of (P f a , P d ) such that κ(P f a , P d ) = 1.

Fig. 3 .

 3 Fig. 3. MSRL vs. σ 2 for L = 100.

Fig. 4 .

 4 Fig. 4. The SRL for multidimensional harmonic retrieval with orthogonal known sources for M equally powered sources, where P = 3, 4, 5, 6, L = 100, and the numbers of sensors are given by N 1 = 3, N 2 = 5, N 3 = 4, N 4 = 4, N 5 = 4 and N 6 = 3.

Fig. 5 .

 5 Fig. 5. MSRL vs. SNR 1 , the SNR of the first source, and SNR 2 , the SNR of the second source. One can notice that the optimal distribution of the SNR (which corresponds to the lowest MSLR) corresponds to SNR 1 = SNR 2 = SNR total 2

Fig. 7 .

 7 Fig. 7. The plot of the MSRL vs. N 2 in the case of P = 3, N 1 = 8 and a total number of sensors N = 22.

2 T 1 ) 2 ω 2 )

 2122 , one can deduce CRB( ξ) where ξ = g(ξ) = [δ ϑ T ] T in which ϑ ω ((2) T . . . ω (P ) T T . Thanks to the Jacobian matrix given by ∂g(ξ) ∂ξ = h = g1 g2 . . . gP T ⊗ 1 -1 T , in which gp = and A = 0 I . Using the change of variable formula CRB( ξ) = ∂g(

2 and

 2 h k = 1 -1 T ⊗ g1(k) g2(k) . . . gP (k) T , in which gp(k) = ∂δ(k)∂ω A = 0 I . Since |x| k can be written as √ x 2k . Thus, for x = 0, one has

4 for

 4 m = m and p = p , L αm 2 N (2Np-1)(Np-1) 6 for m = m and p = p ,

  (Np + 1)(N p + 1) .

Pp=1minSNR 1 , 2 .

 12 Cp -P p, p = 1 p = p gpg p Cp,p is independent from SNR1 and SNR2, minimizing δ is equivalent to minimize G (SNR1, SNR2) whereG (SNR1, SNR2) = δ 2 LN6Using the method of Lagrange multipliers, the problem is as follows:SNR 2 G (SNR1, SNR2) s.t. SNR1 + SNR2 = SNR totalThus, the Lagrange function is given byF ( SNR1, SNR2, λ) = G (SNR1, SNR2) + λ (SNR1 + SNR2 -SNR total )where λ denotes the so-called Lagrange multiplier. A simple derivation leads to, ∂λ = SNR1 + SNR2 -SNR total = 0.(48)Consequently, from[START_REF] Petersen | The matrix cookbook[END_REF] and (47), one obtains SNR1 = SNR1. Using (48), one obtains SNR1 = SNR2 = SNR total Using the constraint SNR1 + SNR2 = SNR total one deduces corollary 1. APPENDIX F Proof of corollary 2 Minimizing δ w.r.t. N1, . . . , NP is equivalent to minimizing the function f (N ) = P p=1 Cp -P p, p = 1 p = p gpg p Cp,p , where N = [N1 . . . NP ] T . However, since the numbers of sensors on each array, N1, . . . , NP , are integers, the derivation of f (N ) w.r.t. N is meaningless. Consequently, let us define the function f (.) exactly as f (.) where the set of definition is R P instead of N P . Consequently, f ( N )| N =N = f (N ), where N = [ N1 . . . NP ] T , in which N1, . . . , NP are real (continuous) variables. Using the method of Lagrange multipliers, the problem is as follows: min N f ( N ) P p=1 Np = Ntotal where Ntotal is a real positive constant value. Thus, the Lagrange function is given by Λ N , λ = f ( N ) + λ P p=1 Np -Ntotal where λ denotes the Lagrange multiplier. For a sufficient number of sensors, the Lagrange function can be approximated by Λ N , λ ≈ This system of equations seems hard to solve. However, an obvious solution is given by N1 = • • • = NP = N and λ = 3V +1 N 4 -3 V (P ν-1)+V -1 N 3 in which ν = P p, p = 1 p = p gpg p . Since, P p=1 Np = Ntotal , thus the trivial solution is given by N1 = • • • = NP = Ntotal P . Consequently, if Ntotal is a multiple of P then, the solution of minimizing the function f ( N ) in R P coincides the solution of minimizing the function f (N ) in N P . Thus, the optimal placement minimizing the MSRL is N1 = • • • = NP = Ntotal P . This conclude the proof.

The notion of distance and closely spaced signals used in the following, is w.r.t. to the metric space (d, C), where d : C × C → R in which d and C denote a metric and the set of the parameters of interest, respectively.

See[START_REF] Vanpoucke | Algorithms and Architectures for Adaptive Array Signal Processing[END_REF]-[START_REF] Gershman | Space-time processing for MIMO communications[END_REF] for some pratical examples for the multidimensional harmonic retrieval model.

This study can be straightforwardly extended to other norms. The choice of the l 1 is motivated by its calculation convenience (see the derivation of Result 1 and Appendix A). Furthermore, since the MSRL is considered to be small (this assumption can be argued by the fact that the high resolution algorithms have asymptotically an infinite resolving power[START_REF] Vantrees | Detection, Estimation and Modulation theory: Optimum Array Processing[END_REF]), thus all continuous p-norms are similar to (i.e., looks like) the l 1 norm. More importantly, in a finite dimensional vector space, all continuous p-norms are equivalent [36, p. 53], thus the choice of a specific norm is free.

, but not equal. This is not really a restrictive assumptions, since in most applications, having two or more identical parameters of interest is a zero probability event[10, p. 53].[START_REF] Mokios | On 3-D harmonic retrieval for wireless channel sounding[END_REF] Note that applying[START_REF] Amar | Fundamental limitations on the resolution of deterministic signals[END_REF] for P = 1 and for κ(P f a , P d ) = 1, one obtains the Smith criterion[START_REF] Smith | Statistical resolution limits and the complexified Cramér Rao bound[END_REF].

Where O(.) denotes the Landau notation[START_REF] Cormen | Introduction to algorithms[END_REF].

One should note, that we assumed a uniform linear multi-array, and the problem is to find the optimal distribution of the number of sensors on each array. The more general case, i.e., where the optimization problem considers the non linearity of the multi-way array, is beyond the scope of the problem addressed herein.

Consequently, after some calculus, one obtains CRB(δ(k))

where

and where Across(k) =

2 )sgn(ω

) CRB(ω

)+CRB(ω

)-2CRB(ω

) . Consequently, note that resolving analytically the implicit equation (32) w.r.t. δ(k) is intractable (aside from some special cases). Whereas, resolving analytically the implicit equation ( 30) can be tedious but feasible (see Section IV).

Furthermore, denoting gp(1) = gp, Across(1) Across and A direct (1) A direct and using (32) one obtains [START_REF] Vantrees | Detection, Estimation and Modulation Theory[END_REF].

APPENDIX B Proof of Lemma 1

From (20) one can note the well known property that the model signal parameters are decoupled from the noise variance [START_REF] Kay | Fundamentals of Statistical Signal Processing[END_REF]. Consequently, the block-diagonal structure in ( 21) is self-evident. Now, let us prove [START_REF]The Cramér-Rao bound on frequency estimates of signals closely spaced in frequency (unconditional case)[END_REF]. From (4), one obtains


Matthijs Douze

Hervé Jégou

Cordelia Schmid

Patrick Pérez

Compact video description for copy detection with precise temporal alignment

This paper introduces a very compact yet discriminative video description, which allows example-based search in a large number of frames corresponding to thousands of hours of video. Our description extracts one descriptor per indexed video frame by aggregating a set of local descriptors. These frame descriptors are encoded using a time-aware hierarchical indexing structure. A modified temporal Hough voting scheme is used to rank the retrieved database videos and estimate segments in them that match the query. If we use a dense temporal description of the videos, matched video segments are localized with excellent precision. Experimental results on the Trecvid 2008 copy detection task and a set of 38000 videos from YouTube show that our method offers an excellent trade-off between search accuracy, efficiency and memory usage.

Introduction

We consider the problem of searching a transformed query video, or part of this query, in a large database of videos. This is important, in particular, for detecting video copies that may be illegally delivered on peerto-peer networks and user generated content sites such as YouTube. The most common transformations observed in practice are camcording and re-encoding, though sophisticated video post-processing is also encountered.

In recent evaluations [START_REF] Over | Trecvid 2008-goals, tasks, data, evaluation mechanisms and metrics[END_REF][START_REF] Law-To | Video copy detection: a comparative study[END_REF], the use of local descriptors [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF][START_REF] Mikolajczyk | A performance evaluation of local descriptors[END_REF][START_REF] Joly | New local descriptors based on dissociated dipoles[END_REF][START_REF] Douze | An image-based approach to video copy detection with spatio-temporal post-filtering[END_REF] combined with a frame voting system appeared to be the most successful architecture for video copy detection. These state-of-the-art systems search individually for each local descriptor of the query video in a structure indexing all local descriptors of the video database. The typical memory requirement associated with representing the set of local descriptors of a video frame ranges from 1 to 10 Kbytes. This seriously limits the number of video frames that can be indexed in practice. Therefore, the video frames are subsampled, which reduces the capability to find very Fig. 1. Overview of our video copy detection method. Local descriptors of a video frame are aggregated into a single vector and the dimension of this vector is reduced. The videos to be indexed are encoded using a temporal-aware indexing scheme. No encoding is applied to the query frame descriptors. A weighted temporal Hough transform provides the final ranking of the database videos w.r.t. the query.

short clips and to determine the precise localization of the query in the database videos. Furthermore, even with subsampling, very large video datasets (several thousands hours of videos) cannot be handled efficiently.

The objective of this paper is to address these scalability and localization issues, while maintaining a high recognition rate. Figure 1 gives an overview of our approach for querying and matching video segments. The individual steps are:

1. Local descriptors are extracted from video frames (either query or database video) and subsequently aggregated into a single vector. This aggregation process is similar to recent approaches for descriptor aggregation [START_REF] Perronnin | Fisher kernels on visual vocabularies for image categorization[END_REF][START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF] which outperform the popular bag-of-features (BOF) representation [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] with descriptors of lower dimension. 2. The dimensionality of this frame descriptor is reduced with either a technique based on a multiplication with a sparse matrix or principal component analysis (PCA). 3. On the database side, the reduced descriptors are encoded within an indexing structure that takes into account the temporal regularity. The video is split in segments. A first description is computed for a segment by minimizing a fidelity criterion for frames of this segment.

In the spirit of [START_REF] Jégou | Product quantization for nearest neighbor search[END_REF], this segment descriptor is refined by a code based on a product quantizer. 4. Each frame's approximate description is refined by encoding the difference between the frame descriptor and the vector describing the segment it belongs to. 5. A modified temporal Hough voting scheme is used to fuse the votes obtained at the frame level. Its main difference with the conventional method is that the votes are weighted so that their contribution is penalized if 1) the query frame has received a large amount of votes and 2) the database frame has voted several times.

The paper is organized as follows. The frame description method is introduced in Section 2. Section 3 describes how frame descriptors are indexed and retrieved when a query frame descriptor is submitted. The voting scheme is presented in Section 4. The contributions of the different steps are evaluated in Section 5. Furthermore, we compare to the state of the art on the Trecvid 2008 benchmark, and obtain top results in terms of localization accuracy. The scalability of the approach is demonstrated by experiments on 38000 YouTube videos represented by more than 200 million frames. We show that these videos are indexed in less than 5GB of memory.

Video description

Local description

The first step of our video description extracts a set of local features for each frame. The same approach is used to extract descriptors for the query and database videos. Here, regions are obtained with the scale-invariant Hessian detector [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF] and described by CSLBP [START_REF] Heikkila | Description of interest regions with local binary patterns[END_REF]. Similar to SURF [START_REF] Bay | Surf: Speeded up robust features[END_REF] and DAISY [START_REF] Winder | Picking the best Daisy[END_REF], this descriptor is a variant of the SIFT descriptor which provides comparable results and reduces the computation time significantly: extracting local descriptors from a frame takes about 35 ms on one 2.4GHz processor core. Note that for large databases, the time for feature extraction is not the critical operation at query time, because it only depends on the query length, not the database size.

Local descriptor aggregation: non probabilistic Fisher Kernel

Given a set of local descriptors {x 1 , . . . , x i , . . . } for each video frame, it is impossible to store all of them in memory in the case of large scale video search, even if only a few bytes are required for each of them [START_REF] Chandrasekhar | Chog: Compressed histogram of gradients: A low bit-rate feature descriptor[END_REF][START_REF] Calonder | Compact signatures for high-speed interest point description and matching[END_REF].

In general, several hundreds of local descriptors are extracted for each frame.

We, therefore, aggregate the local descriptors into a single vector. We adopt a variant [START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF] of the Fisher Kernel image representation introduced by Perronnin et al. [START_REF] Perronnin | Fisher kernels on visual vocabularies for image categorization[END_REF]. The resulting vector, called vector of locally aggregated descriptors (VLAD), provides a compact yet effective representation of images. Assuming that a k-means codebook with k centroids {c 1 , . . . , c j , . . . , c k } has been learned, we obtain the VLAD descriptor for a frame, denoted by µ, as follows:

1. As for the bag-of-features representation, each local descriptor x i of the frame is assigned to the closest centroid in the codebook, i.e., to the quantization index NN(x i ) = arg min j ||x i -c j ||. 2. Given the set of descriptors assigned to a centroid c j , the vector µ j is obtained by summing the differences between these descriptors and the centroid:

µ j = i:NN(x i )=j x i -c j . (1)
3. The VLAD descriptor associated with a frame is obtained by concatenating the vectors µ j into a single vector. 4. As proposed in [START_REF] Perronnin | Large-scale image retrieval with compressed Fisher vectors[END_REF] for the Fisher Kernel image representation, we apply a power-law normalization to the components to reduce the contribution of the most energetic ones. Here, we use the signed square root of each component. The vector is subsequently L2-normalized and is denoted µ in the following.

The resulting vector is of dimension k times the dimensionality of the input vector, e.g., k × 128 for the CSLBP descriptor. For the same codebook size, the dimensionality of the descriptor is significantly larger than for the bag-of-features representation [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF]. However, the VLAD description is already discriminant for low values of k in contrast to BOF, which requires very large codebooks (up to 1 million) to provide the best results [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF][START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF]. Therefore, the dimensionality of the VLAD vector is typically lower than for BOF. It is worth noting that this representation can be seen as a non-probabilistic version of the Fisher kernel. In the latter, a Gaussian mixture model and soft assignment are used instead of k-means centroids, and additional information (variance and count) are used to obtain a richer (but longer) representation.

Dimensionality reduction of frame descriptors

Local descriptor aggregation results in one VLAD descriptor per video frame. This descriptor is highly dimensional: for a typical value of k = 64, the vector µ has D = 128×k = 8192 components. Such a vector is difficult to index due to its dimensionality [START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF]. We, therefore, use and compare two different methods to reduce the dimensionality:

1. Principal component analysis (PCA) allows to reduce the dimension D of the VLAD descriptor to a smaller dimension d [START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF]. The vector µ is multiplied with a projection matrix M formed by the first principal eigenvectors of an empirical covariance matrix. The PCA matrix is pre-multiplied with a random rotation to "whiten" the output; 2. Alternatively, we define M as a d × D sparse matrix obtained as M = P σ, where σ is a D × D random permutation matrix and P is a d×D aggregation matrix that sums several consecutive components. For example with D = 6 and d = 2, a possible matrix M is:

0 1 1 0 1 0 1 0 0 1 0 1 M = 1 1 1 0 0 0 0 0 0 1 1 1 P ×       0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1       σ . (2)
The two techniques are evaluated in Section 5. The advantage of using a structured matrix is that there is no need for a training stage. The dimensionality reduction is also cheaper to compute, because the multiplication is more efficient with the sparse matrix M than with the full matrix obtained by PCA. However, during the search, the dimensionality reduction typically has a low computing cost compared to the other steps.

The descriptor f ∈ d of reduced dimensionality is obtained by multiplying the matrix M with the VLAD descriptor µ and by L2-normalizing it. The resulting vector f is used in the following as the frame descriptor. The dot product is used as a similarity measure.

Indexing frame descriptors with temporal integration

The objective of this section is to provide a compact representation and an efficient search mechanism for the frames of a database video. Let us consider a video to be indexed 4 , for which we have extracted a set f 1 , . . . , f t , . . . , f T of d-dimensional descriptors using the method introduced in Section 2. The individual descriptors will be approximated, indexed and searched using three successive refinement levels:

1. joint description of a group of contiguous frames: all the frames associated with the same time segment have the same coarse approximation; 2. refinement of this video segment descriptor; 3. refinement of the individual frame descriptors.

Fig
Each of these levels provides an improved approximation of an indexed frame descriptor, as illustrated by Figure 2. This procedure is not used for the frames of the query, i.e., their VLAD descriptors are not approximated. Our approach is, to some extent, similar in spirit to the method proposed in [START_REF] Jégou | Product quantization for nearest neighbor search[END_REF]. However, a major difference is the integration of the temporal aspect into the indexing scheme.

Level 1: Coarse segment description

Due to the temporal dependency between frames, contiguous frame descriptors of a video shot are similar. We exploit this property by providing a shared description for consecutive frames. Let us define a segment {t b , . . . , t e } as an interval of consecutive frames, for which the same level 1+2 approximation of the descriptor is used.

If segments are of fixed size, we have t e = t b + 1/r -1, where r is the ratio between the number of segments and the total number of frames in the video. The first approximation is given by a coarse vector quantizer q c (.), for which the codebook C c = {c 1 , . . . , c L } is learned using a spherical k-means algorithm. The coarse quantization index i c (t b : t e ) associated with the segment {t b , . . . , t e } aims at best representing the set of frame descriptors {f t b , . . . , f te } with respect to the total square error, i.e.,

i c (t b : t e ) = arg min i∈Cc t=t b :te ||f t -c i || 2 , (3)
which is equivalent to

c ic(t b :te) = q c r t=t b :te f t (4)
The vector c ic(t b :te) is the level-1 approximation of the frame descriptors in the segment {t b , . . . , t e }. When searching the nearest neighbors of a given query descriptor, only the database descriptors associated with the k c closest elements in C c are explored.

We have tested both a fixed and adaptive number of frames per segment. Several variants for selecting keyframes have been tested in the adaptive case. Best results were obtained when constructing the segments based on the consistency of the k-nearest neighbors in C c for the frame descriptors. However, experimental results showed that no variant is better than uniform subsampling. We, therefore, only use segments of fixed size in the following.

Level 2: segment descriptor refinement

The index associated with a given video segment is not precise, as an approximation with a centroid in C c introduces a large quantization error. Therefore, similar to [START_REF] Jégou | Product quantization for nearest neighbor search[END_REF], we refine this first approximation by using a product quantizer q f , whose codebook5 is denoted by C f . The total number of centroids implicitly defined by a product quantizer composed of m f subquantizers having L f centroids each is equal to (L f ) m f . This quantizer aims at reducing, over the set of frames associated with a given segment, the average energy of the error vector f t -c ic(t b :te) made by the first approximation q c (f t ) = c ic(t b :te) . The new approximation of a frame descriptor f t associated with the segment {t b , . . . , t e } is, therefore, of the form

f t ≈ c ic(t b :te) + c i f (t b :te) , (5)
where the centroid c i f (t b :te) ∈ C f is obtained by the minimization

c i f (t b :te) = arg min c i ∈C f t=t b :te ||f t -c ic(t b :te) -c i || 2 . ( 6
)
The minimization is efficiently done using the decomposition associated with the product quantizer. Note that this quantizer q f is more precise than the coarse quantizer used in the first stage, because the set of centroids C f that is implicitly defined by the product quantizer is large: it is 2 8×64 for the typical 64-byte codes we use (m f = 64, L f = 256). The product quantizer decomposition is used to obtain a complexity comparable to that of a quantizer comprising L f elements.

Level 3: refinement of individual frame descriptors

So far, the frames of a segment are described by the same approximation.

We now refine the description of each individual frame f t by using another refinement product quantizer q r induced by m r subquantizers with L r centroids each. This quantizer encodes the error resulting from the two previous approximations by minimizing the quantization error of f t . For a time instant t such that t ∈ {t b , . . . , t e }, this is done by quantizing the residual error vector f t -c ic(t b :te) -c i f (t b :te) . The frame descriptor f t is therefore approximated by

ft = c ic(t b :te) + c i f (t b :te) + q r (f t -c ic(t b :te) -c i f (t b :te) ). (7)

Search procedure

Searching a query frame vector y in a database of frame descriptors B = {f 1 , . . . , f T } proceeds in a hierarchical manner.

1. The k c nearest neighbors of y in C c identify the segments to be considered: only those associated with one of the selected k c indexes are explored. 2. For each vector f i in the set of selected lists, the distance approximation

l 2 (f i , y) = l 2 (f i -q c (f i ), y -q c (f i )) ≈ l 2 (q f (f i -q c (f i )) , y -q c (f i )) (8)
is efficiently obtained from the quantizer indexes q f (f i -q c (f i )) by exploiting ADC method of [START_REF] Jégou | Product quantization for nearest neighbor search[END_REF]. The best segments corresponding to a query vector are found based on the approximation of the square distance of Equation 8. This step returns the set of the k f nearest segment descriptors. 3. The query frame descriptor y is now compared to all the approximated ft frame descriptors associated with the k f segments found in the previous stage. This step returns a set of k r nearest frame descriptors.

Complexity

The cost of searching a frame descriptor in a database containing T frames is expressed in terms of the number C dist of regular d-dimensional vector comparisons and the amount C mem of memory scanned in the indexing structure. These are given by

C dist = L + k c L f + k f r (9)
and

C mem = α k c L r T m f log 2 L f + k f r m r log 2 L r , (10)
where log 2 L f = log 2 L r = 8 bits = 1 byte in our case. The factor α ≥ 1 accounts for the fact that the probability to assign a frame descriptor to an index is not uniform over C c . As observed in [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF] in the context of the BOF representation, this increases the expectation of the number of elements that are processed. We measured that α ≈ 8.4 with our parameters. Note that our calculation of C dist assumes that the level-2 search is optimized by using look-up tables computed on-the-fly.

Temporal alignment: improved Hough transform

Once each query frame has been matched to k r putative frames of the database, the video search matches a sequence of query descriptors to sequence(s) from the database. This sequence matching can be cast in terms of temporal sequence alignment and addressed by dynamic programming techniques, such as dynamic time warping. However, this type of approaches requires a complete frame-to-frame matching cost matrix, which is not feasible at this stage of the detection system. Furthermore, they require a good initialization of the starting and endpoint of the matching sequences. Simplified approaches can be used instead, e.g., partial alignment [START_REF] Yeh | Video copy detection by fast sequence matching[END_REF] or classic voting techniques, such as Hough transform or RANSAC. We adopt a temporal Hough transform [START_REF] Douze | An image-based approach to video copy detection with spatio-temporal post-filtering[END_REF], as it allows the efficient estimation of 1D temporal shifts. Individual frame votes are re-weighted according to suitable normalizations. In particular, re-weighting is used to reduce the harmful influence of temporal "burstiness", i.e., the presence of query frames that match strongly with multiple, possibly unrelated frames in the database. This is similar to the burstiness of local descriptors in images [START_REF] Jégou | On the burstiness of visual elements[END_REF].

Hough transform

As output of the indexing structure, we have a set of matches, each of them represented by a score s(τ, t) > 0 between the query frame timestamp τ and the database timestamp t. This score is given by the inner product between the query frame descriptor and the approximation of the database descriptor in Equation 7. We set to 0 the scores of frames that are not matched:

s(τ, t) = y τ , ft if t is retrieved by the index 0 otherwise. ( 11
)
The temporal Hough transform consists in computing a histogram h(δ) to accumulate these scores for all δ = t -τ hypotheses. Denoting by Y and B = {1, . . . , T } the sets of time instants on the query side and the database side, respectively, the score is formally obtained as

h(δ) = τ ∈Y s(τ, τ + δ), (12)
where s(τ, t) = 0 if t / ∈ B. Peaks (maximum values) are then searched in the histogram h. We typically select 100 peaks and then apply non maximum suppression, i.e., peaks that are closer than 1 minute to a stronger one are discarded. For each peak identified by δ, the boundaries of the matching segments are then identified by collecting all the matches (τ, t) associated with a hypothesis δ such that |τ -t -δ| < 10. The final score is the sum of scores of these matches.

Handling the temporal burstiness

As mentioned in Section 3, consecutive video frames are often very similar, and so are their descriptors. This temporal dependency biases the scores returned by the Hough histogram, as bursts of votes occur for some frames, both on the query and database. This emphasizes them, i.e., they gather an abnormally large amount of scores.

We address this problem by modifying the scoring strategy in a way that mitigates this effect, in the spirit of the re-weighting scheme proposed in [START_REF] Douze | An image-based approach to video copy detection with spatio-temporal post-filtering[END_REF]. This is done by updating the score, prior to the Hough histogram computation in two steps:

s 1 (τ, t) = s(τ, t)/ τ ∈Y s(τ, t) and s 2 (τ, t) = s 1 (τ, t)/ t∈B s 1 (τ, t), (13)
where the computation is done efficiently by considering only the non-zero score values in the summations. The updated score s 2 is used instead of the original scores in Equation 12. We will show in Section 5 that this procedure significantly improves the quality of the Hough estimation.

Experiments

Datasets and evaluation protocol

Trecvid'08. This dataset contains 200 hours of Dutch television broadcasts. It was used for the copy detection pilot task in the Trecvid'08 evaluation campaign. A set of 134 query clips was extracted from the dataset and 67 clips from other videos were added as negatives, i.e., with no corresponding videos in the database. Some clips were embedded in a distractor video and all were then transformed with 10 different transformations, see Table 2. As a result, 2010 queries with varying degrees of difficulty are used to evaluate a system. The performance measure used to evaluate the search quality in the Trecvid competition is the Normalized Detection Cost Ratio (NDCR), which integrates the cost of missing a true positive with the cost of retrieving an incorrect video. It is equal to 0 if all the true positives are returned before the false positives (no matter how many there are) and lower values of the NDCR correspond to better results. A result video segment is considered as a true positive if it overlaps with the ground-truth. We have used this measure6 to compare our results with those obtained by the participants of the Trecvid'08 evaluation, see Subsection 5.3.

STV (Small Trecvid). In order to evaluate the parameters of our approach, we created a reduced version of the Trecvid dataset, referred to in the following as STV. It uses a subset of 21 h of the videos. From these videos we extracted a set of 100 clips not used by in the Trecvid'08 queries, embedded them in an independent distractor set and transformed them with some of the most challenging transformations.

This dataset is smaller than the Trecvid'08 dataset. However, the transformations are more challenging on average. We, therefore, obtain comparable conclusions with reduced runtime. Furthermore, using this dataset for parameter evaluation avoids optimizing parameters on the Trecvid query set, and provides a fair comparison with the state of the art on this dataset.

YouTube. In order to evaluate the scalability of our system, we have collected a dataset of YouTube videos. We downloaded a total of 38,000 videos from YouTube, corresponding to 189 million frames (or 2100 h). Most of the videos have a similar resolution as the Trecvid ones (about 352*288) and the number of interest points extracted per frame is similar (about 300 per frame on average). These videos are used as distractors in our large scale experiment in Section 5.4.

Evaluation measures. In addition to the NDCR measure used for Trecvid, we have used two additional measures to evaluate performance, localization accuracy and average precision. The localization accuracy for a result segment is measured as the overlap between the ground-truth and the returned result: Ω = |T gt ∩ T found |/|T gt ∪ T found |. If the match is incorrect, Ω = 0, and if the localization of a match is perfect, Ω = 1. Better matches have a higher overlap.

We have used the overall Average Precision (AP) as a quality measure. The results returned for all queries are evaluated together ranked by their scores. A match is considered a true positive if the overlap measure is above 0.5. A precision-recall curve evaluates the results as a function of the score. The area under this curve is the AP measure. In the following, we measure the impact of these parameters. The performance is reported for the STV dataset.

The dimensionality reduction is evaluated on the two first levels of the method, i.e., without frame grouping (r = 1) nor refinement (Levels 1+2 only). Our aggregation method is compared with PCA to reduce the D = 8192 dimensions of the frame descriptor to d = 2048. For this operating point, dimensionality reduction with an aggregator (see Section 2.3) gives AP=87.7, and the PCA-based dimensionality reduction achieves AP=90.1. In the following, we use the PCA-based dimensionality reduction. Indexing: impact of the quantization and of the refinement step.

Table 1 shows the influence of the descriptor quantization and frame grouping on the search accuracy, memory usage and search time. Bruteforce search gives an upper bound on the performance that can be achieved by using our frame descriptor, but is unreasonably expensive. We denote by Levels 1+2 the methods that do not refine the segment level representation on the frame level. If r = 1, then frames are directly considered as segments, while Level 1+2, no refinement has the same effect as a subsampling, except that the average frame descriptor over a segment is used instead of a particular frame of this segment. This variant provides lower search quality, but is interesting to index very large datasets. One can observe that subsampling the video and indexing subsampled frames strongly degrades the performance, by 16 points of AP. The refinement improves the results. Short codes (m r = 16 bytes) are sufficient to capture most of the possible improvement. Note that, for large databases, this last refinement stage (Level 3) is the limiting factor in terms of memory usage, even with the setting m r = 16.

Burstiness handling. We have evaluated the impact of our vote regularization procedure which addresses the problem of burstiness (cf. section 4.2). This method significantly improves the results, as shown below:

bursts regularization AP none 83.8 database-side (using s 1 in Equation13) 84.3 full regularization (s 2 ) 91.2 We also provide the rank associated with our score. One can see that our system is very competitive, in particular for the most interesting transformations encountered on a peer-to-peer network: camcording and re-encoding. Note that the best score for these transformations is obtained with the approach of [START_REF] Douze | An image-based approach to video copy detection with spatio-temporal post-filtering[END_REF]. This approach requires 10 GB of RAM against 300 MB used here, i.e., it is difficult to scale to very large video sets. Furthermore, it is more than 5 times slower than the approach presented in this paper (13× "real time" against 2.47× here on a single computing core). We also compare our localization accuracy with the four best runs (from different participants) of the Trecvid'08 competition. We measure the accuracy on the 195 videos that were correctly retrieved by all 4 runs as well as by our method (i.e. the resulting segment has an overlap with the the ground-truth above 0.5). The measure used is average overlap. The results are: rank of the participant 1 2 3 4 ours mean overlap 0.952 0.858 0.846 0.884 0.973

We can observe that our approach localizes the segments very precisely, i.e., with a better precision than the competing approaches.

Large scale experiment

The large scale experiments are performed on the YouTube dataset merged with our STV dataset. Figure 3 shows the AP obtained as the function of the growing database size (up to 2316 hours). When performing the experiments on the whole set, the index requires 4.6 GB of RAM to index 208 million frames. The search is slower on that scale: 23.5 real time for a single processor core. One can observe that the AP measure decreases as to be expected, but that results are still good, i.e., we obtain an AP=0.53 on the entire set. Typical retrieval results of our system are shown in Figure 4.

5. 2

 2 Impact of parameters Unless stated otherwise, for the parameters introduced in Sections 2 and 3, we have used the following values: Descriptor D = 128 × k = 8192 d = 2048 r = 1/10 Coarse quantizer L = 2048 k c = 64 Fine quantizer L f = 256 m f = 64 k f = 128 Refinement L r = 256 m r = 16 k r = 32

Fig. 3 .

 3 Fig. 3. Retrieval performance on the STV dataset combined with a varying number of videos from the remaining Trecvid videos and YouTube.

Fig. 4 .

 4 Fig. 4. Example results for video retrieval in our large scale dataset, (left) query and (right) best retrieved video. Lefts pairs: correct retrieval results. Right pairs: incorrect retrieval. Note the visual similarity between the queries and the retrieved videos.

  . 2. Hierarchical representation of the frame descriptors. At levels 1 and 2, segments are represented by a frame descriptor. At level 1, the descriptor is quantized to a single index, that is refined at Level 2 by a product quantizer code. At level 3, the individual frames of the segment are represented as a refinement of their segment descriptor. Typically, a segment descriptor is encoded in 4+64 bytes, while each frame is refined by a 16-byte code.

Table 1 .

 1 Evaluation of memory usage, search accuracy and timings on the STV dataset. The method "Levels 1+2, r = 1" indicates that frames are directly considered as segments in our method without any further refinement. The timings are given as a slowdown factor w.r.t. the "real" video time for one single 2.4GHz processor core (not including the frame description time).

	structure/algorithm	total mem	Cmem	C dist	AP	time
	none/brute-force search	62 GB	62 GB	19 M 96.7	89 ×
	Levels 1+2, r = 1	122 MB	25 MB	18432 90.1 2.52 ×
	Levels 1+2, no refinement	12 MB	2.71 MB 18432 74.1 1.10 ×
	Levels 1+2+3 (mr = 16)	43 MB	2.73 MB 19710 91.2 1.33 ×
	Levels 1+2+3 (mr = 32)	73 MB	2.75 MB 19710 90.5 1.42 ×
	Levels 1+2+3 (mr = 64)	134 MB	2.79 MB 19710 91.7 1.43 ×

Table 2 .

 2 Evaluation of our method relative to other Trecvid'08 participants. The score is NDCR (the lower the better). The rank is obtained by taking the best run from each of the 22 participants. Comparison with state of the art Table2compares the NDCR scores of our system with the best and second best run (from different participants) of the Trecvid'08 competition 7 .

	no transformation	best	second ours	rank (/23)
	1 camcording	0.079	0.363	0.224	2
	2 picture in picture	0.015	0.148	0.321	4
	3 insertion of patterns	0.015	0.076	0.079	3
	4 strong re-encoding	0.023	0.095	0.064	2
	5 change of gamma	0.000	0.000	0.023	3
	6 photometric attacks	0.038	0.192	0.064	2
	7 geometric attacks	0.065	0.436	0.140	2
	8 3 random transformations from 6/7 0.045	0.076	0.437	5
	9 5 random transformations from 6/7 0.038	0.173	0.693	5
	10 5 random transformations	0.201	0.558	0.537	2
	5.3

We treat several videos as a single long video, except that we constrain a segment of frames not to cross video boundaries.

A product quantizer decomposes the space into a Cartesian product of low dimensional subspaces and quantizes each space separately. As a result, learning codebooks and searching the quantization index have a low complexity even for very large codebooks. The codebook C f has not to be stored explicitly.

NIST (the institute organizing Trecvid) provided the software to compute this measure as well as the results of the other participants.

The competitors can not be identified by name due to the non disclosure agreement of Trecvid.

Acknowledgements

This work was partly realized as part of the Quaero Programme, funded by OSEO, French State agency for innovation.