Souleymen Sahnoun
email: souleymen.sahnoun@gipsa-lab.grenoble-inp.fr

El-Hadi Djermoune
email: el-hadi.djermoune@univ-lorraine.fr

David Brie
email: david.brie@univ-lorraine.fr

Pierre Comon
email: pierre.comon@gipsa-lab.grenoble-inp.fr

A Simultaneous Sparse Approximation Method for Multidimensional Harmonic Retrieval

Keywords: Multidimensional harmonic retrieval, frequency estimation, simultaneous sparse approximation, multigrid dictionary refinement, Cramér-Rao lower bound

In this paper, a new method for the estimation of the parameters of multidimensional (R-D) harmonic and damped complex signals in noise is presented. The problem is formulated as R simultaneous sparse approximations of multiple 1-D signals. To get a method able to handle large size signals while maintaining a sufficient resolution, a multigrid dictionary refinement technique is associated to the simultaneous sparse approximation. The refinement procedure is proved to converge in the single R-D mode case. Then, for the general multiple modes case, the signal tensor model is decomposed in order to handle each mode separately in an iterative scheme. The proposed method does not require an association step since the estimated modes are automatically "paired". We also derive the Cramér-Rao lower bounds of the parameters of modal R-D signals. The expressions are given in compact form in the single tone case. Finally, numerical simulations are conducted to demonstrate the effectiveness of the proposed method.

Introduction

The problem of estimating the parameters of sinusoidal signals from noisy measurements is an important topic in signal processing and several parametric and nonparameteric approaches have been developed for one-dimensional (1-D) signals [START_REF] Stoica | Introduction to Spectral Analysis[END_REF]. Recently, this problem has received a renewed interest thanks to the emergence of multidimensional (R-D) applications. Indeed, parameter estimation from R-D signals is required in numerous applications in signal processing and communications such as nuclear magnetic resonance (NMR) spectroscopy, wireless communication channel estimation [START_REF] Gershman | Space-time processing for MIMO communications[END_REF] and MIMO radar imaging [START_REF] Nion | Tensor algebra and multidimensional retrieval in signal processing for MIMO radar[END_REF]. In all these applications, signals are assumed to be a superposition of R-D sinusoids or, more generally, of exponentially decaying R-D complex exponentials (modal signals). As for the 1-D case, the crucial step is the estimation of the R-D modes (including frequencies and damping factors) because they are nonlinear functions of the data. In this paper, we consider the single snapshot R-D signal model described in [START_REF] Li | A high-resolution technique fro multidimensional NMR spectroscopy[END_REF].

In order to achieve high resolution estimates, parametric approaches are often preferred to nonparamatric ones. Several parametric R-D methods (R ≥ 2) have been proposed.

They include linear prediction-based methods such as 2-D TLS-Prony [START_REF] Sacchini | Two-dimensional Prony modeling and parameter estimation[END_REF], and subspace approaches such as matrix enhancement and matrix pencil (MEMP) [START_REF] Hua | Estimating two-dimensional frequencies by matrix enhancement and matrix pencil[END_REF], 2-D ESPRIT [START_REF] Rouquette | Estimation of frequencies and damping factors by two-dimensional ESPRIT type methods[END_REF], multidimensional folding (MDF) [START_REF] Mokios | On 3-D harmonic retrieval for wireless channel sounding[END_REF], improved multidimensional folding (IMDF) [START_REF] Liu | An eigenvector-based approach for multidimensional frequency estimation with improved identifiability[END_REF][START_REF] Liu | Multidimensional frequency estimation with finite snapshots in the presence of identical frequencies[END_REF], Tensor-ESPRIT [START_REF] Haardt | Higher-order SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems[END_REF], principal-singularvector utilization for modal analysis (PUMA) [START_REF] Sun | Accurate and computationally efficient tensor-based subspace approach for multi-dimensional harmonic retrieval[END_REF][START_REF] So | An efficient approach for twodimensional parameter estimation of a single-tone[END_REF] and the methods proposed in [START_REF] Huang | Multidimensional sinusoidal frequency estimation using subspace and projection separation approaches[END_REF][START_REF] Lin | Efficient multidimensional harmonic retrieval: A hierarchical signal separation framework[END_REF]. All these methods perform at various degrees but it is generally admitted that they yield accurate estimates at high SNR scenarios and/or when the frequencies are well separated. This is obtained at the expense of computational effort. In [START_REF] Sun | Accurate and computationally efficient tensor-based subspace approach for multi-dimensional harmonic retrieval[END_REF], tensor PUMA was proposed as an accurate and computationally efficient multidimensional harmonic retrieval method, which attains the Cramér-Rao lower bound (CRLB) and does not require to build large size matrix or tensor. However its performance degrades rapidly with the increase of the number of components in the R-D signal.

Recently, methods based on sparse approximations have been proposed to address the harmonic or modal retrieval problem [START_REF] Sahnoun | Sparse multigrid modal estimation: Initial grid selection[END_REF][START_REF] Goodwin | Matching pursuit and atomic signal models based on recursive filter banks[END_REF][START_REF] Malioutov | A sparse signal recontruction perspective for source localization with sensor arrays[END_REF][START_REF] Stoica | SPICE: A sparse covariance-based estimation method for array processing[END_REF][START_REF] Sahnoun | Sparse multidimensional modal analysis using a multigrid dictionary refinement[END_REF][START_REF] Sahnoun | Sparse modal estimation of 2-D NMR signals[END_REF][START_REF] Sward | High resolution sparse estimation of exponentially decaying signals[END_REF][START_REF] Adalbjörnsson | High resolution sparse estimation of exponentially decaying two-dimensional signals[END_REF]. For time-data spectral estimation, the dictionary is formed from a set of (normalized) complex exponentials potentially embedded in the data, which allows one to easily include some prior knowledge about the position of some known modes. More generally, the usual choice is a uniform spectral grid obtained by sampling the frequency and damping factor lines. Clearly, a fine grid is required to get a good resolution but, on the other hand, it will result in a huge dictionary [START_REF] Sahnoun | Sparse multigrid modal estimation: Initial grid selection[END_REF]. This complexity is further increased in the case of R-D signals in which we are confronted with 2R-D grids. In order to reduce the computational burden, a multigrid scheme for sparse approximation was proposed in [START_REF] Sahnoun | Sparse multidimensional modal analysis using a multigrid dictionary refinement[END_REF] to iteratively refine the dictionary starting from a coarse one. At each iteration, a sparse approximation is performed and then new grid points (called "atoms") are inserted in the vicinity of active ones leading to a multiresolution-like scheme. This algorithm, which refines jointly R 2-D grids, is efficient but has mainly two drawbacks: 1) it does not have convergence guarantees, 2) the dictionary becomes intractable for large signals when R ≥ 2. Recently, several studies have also focused on gridless sparse recovery methods based on continuous dictionaries [START_REF] Tang | Compressed sensing off the grid, Information Theory[END_REF][START_REF] Tan | Discrete and continuous sparse recovery methods and their applications[END_REF]. However, the proposed algorithms demand a large computational burden even for 1-D signals.

The goal of the present paper is to propose a fast multidimensional modal estimation technique able to handle large signals and yielding a good estimation accuracy.

1. First, the proposed approach, as for some parametric methods for modal retrieval, is based on the idea of estimating the parameters independently along each dimension r = 1, . . . , R. It will be shown that the simultaneous sparse approximation concept [START_REF] Tropp | Algorithms for simultaneous sparse approximation[END_REF][START_REF] Sahnoun | Sparse modal estimation of 2-D NMR signals[END_REF] is well-suited for R-D modal retrieval (R ≥ 2).

2. The second contribution consists in the proposition of a new multigrid scheme which amounts to consider a two-step refinement of 1-D grids, the first step for frequencies and the second one for damping factors. One advantage of this procedure is that it reduces the computational time. The convergence of the proposed multigrid strategy is analyzed in the single tone case (F = 1), and convergence conditions are expressed in terms of atom positions in the initial dictionaries.

3. The extension of this result to the multiple tones case (F > 1) is not trivial because, not only it depends on the selected sparse approximation algorithm, but also on the coherence of the dictionary [START_REF] Tropp | Algorithms for simultaneous sparse approximation[END_REF]. Indeed, due to the multigrid strategy, the columns of the refined dictionary are increasingly correlated, which may prevent convergence even in the noiseless case. Consequently, for F > 1, we exploit an alternative representation of the data model enabling the extraction of the R-D signal tones separately. Therefore, the third contribution of this paper is the derivation of a new algorithm for estimating parameters of R-D damped signals in which the results of the previous contribution apply. The effectiveness of the new algorithm for multiple R-D tones is also analyzed. One very interesting by-product of this approach is that the pairing of R-D parameters is achieved for free, without any further association stage.

The usual way to assess the performances of an estimation method is to compare the variance of the estimates to the CRLB. In [START_REF] Hua | Estimating two-dimensional frequencies by matrix enhancement and matrix pencil[END_REF] Y. Hua derived the CRLB for 2-D frequencies, i.e., undamped 2-D exponentials; no damped signals are considered. Closed-form expressions of the CRLB for the general undamped R-D case are derived in [START_REF] Boyer | Deterministic asymptotic Cramér-Rao bound for the multidimensional harmonic model[END_REF]. CRLB for 2-D damped signals are derived in [START_REF] Clark | Two-dimensional modal anlysis based on maximum likelihood[END_REF]. Therefore, to the best of our knowledge, no compact expressions of the CRLB's are available for the general R-D damped model. Thus, another contribution of the paper is the derivation of the CRLB's for the frequency, damping factor, amplitude and phase of this model. The remainder of this paper is organized as follows. In section 2, we introduce notation and present the R-D modal retrieval problem. In section 3, we formulate the R-D modal estimation problem as R simultaneous sparse estimation problems, show how to construct a modal dictionary on a uniform grid and then describe the new fast multigrid strategy. In section 4, we give sufficient conditions for convergence of the multigrid dictionary refinement in the case of single tone R-D signals. In light of these new results, we propose in section 5 a new efficient algorithm for multiple tones R-D signals. In section 6, we derive the expressions of the CRLB's for the parameters of R-D damped exponentials in Gaussian white noise. We then give the CRLB in the cases of single damped and undamped R-D cisoids. The effectiveness of the proposed method is demonstrated using simulation signals in section 7. Finally, conclusions are drawn in section 8.

Notation and Problem Statement

Notation

Scalars are denoted as lower-case letters (a, b, α), column vectors as lower-case bold-face letters (a, b), matrices as boldface capitals (A, B), and tensors as calligraphic bold-face letters (A, B). Notations (•) T , (•) H and (•) † stand for the transpose, the Hermitian transpose and the pseudo-inverse, respectively. The symbols " " and " " will denote the Khatri-Rao product (column-wise Kronecker) and the Kronecker product, respectively. Both words "mode" and "tone" are used to refer to a component of the multidimensional signal. The tensor operations used here are consistent with [START_REF] Comon | Tensors: a brief introduction[END_REF]:

• the outer product of two tensors A ∈ C M 1 ו••×M R and B ∈ C K 1 ו••×K N is given by: C = A ⊗ B ∈ C M 1 ו••×M R ×K 1 ו••×K N , c(m 1 , . . . , m R , k 1 , . . . , k N ) = a(m 1 , . . . , m R )b(k 1 , . . . , k N ) (1)
• the contraction product acting on the r-th index of a tensor

A ∈ C M 1 ו••×M R and the second index of a matrix U ∈ C K×M r is: B = A • r U ∈ C M 1 ו••×M r-1 ×K×M r+1 ו••×M R , b(m 1 , m 2 , . . . , m r-1 , k r , m r+1 , . . . , m R ) = M r m r =1 a(m 1 , m 2 , . . . , m R )u(k r , m r ) (2) • the matrix A (r) ∈ C M r ×(M 1 •••M r-1 M r+1 •••M R )
represents the unfolding (dimension-r matricization) of the tensor A and corresponds to the arrangement of the dimension-r fibers of A in the columns of the resulting matrix.

• A denotes the the Frobenius norm for tensors.

• The

concatenation of two tensors A 1 ∈ C M 1 ו••×M r-1 ×K 1 ×M r+1 ו••×M R and A 2 ∈ C M 1 ו••×M r-1 ×K 2 ×M r+1 ו••×M R
along the rth dimension is denoted by A 1 r A 2 and obtained by stacking A 1 and A 2 along the rth dimension.

Finally, throughout this paper, the tilde symbol ( ˜) denotes a noisy signal; e.g. ỹ(•) = y(•) + e(•).

Problem Formulation

An R-D modal signal is modeled as the superposition of F multidimensional damped complex sinusoids:

ỹ(m 1 , . . . , m R ) = F f =1 c f R r=1 a m r -1 f,r + e(m 1 , . . . , m R ) (3)
where m r = 1, . . . , M r for r = 1, . . . , R. M r denotes the sample support of the r-th dimension,

a f,r = exp (α f,r + jω f,r ) ∈ C is the f -th mode in the r-th dimension, {α f,r } F,R f =1,r=1 , α f,r ∈ R -, are the damping factors, {ω f,r = 2πν f,r } F,R
f =1,r=1 are the angular frequencies, and c f = λ f exp( jφ f ) is the complex amplitude of the f -th mode where λ f = |c f | denotes the magnitude and φ f the phase. e(m 1 , m 2 , . . . , m R ) is a zero-mean complex Gaussian white noise with variance σ 2 and mutually independent in all dimensions.

In a tensor form, the R-D signal in (3) may be written as

Y = Y + E (4)
where

{ Y, Y, E} ∈ C M 1 ×M 2 ו••×M R .
The problem consists in estimating the set of parameters {a f,r } F,R f =1,r=1 and {c f } F f =1 from the R-D signal samples.

Simultaneous Sparse Approximation for R-D Modal Signals

Tensor Formulation of the Data Model

The noise-free data tensor Y in (4) can be written in the following form:

Y = F f =1 c f a f,1 ⊗ a f,2 ⊗ • • • ⊗ a f,R (5)
where a f,r = [1, a f,r , . . . , a M r -1 f,r ] T , r = 1, . . . , R. Equation ( 5) is called the Canonical Polyadic (CP) decomposition form, or the Candecomp/Parafac decompostion of the tensor Y [START_REF] Comon | Tensors: a brief introduction[END_REF][START_REF] Kolda | Tensor decompositions and applications[END_REF]. The CP model ( 5) can be concisely denoted by Y = c; A 1 , A 2 , . . . , A R where A r = [a 1,r , a 2,r , . . . , a F,r ], r = 1, . . . , R, and c = [c 1 , c 2 , . . . , c F ] T is the vector of complex amplitudes. Using these definitions, the matricized form of Y along the r-th dimension is given by

Y (r) = A r ∆ c (A R • • • A r+1 A r-1 • • • A 1 ) T (6)
where ∆ c = diag(c). Then, we can write

Y (r) = A r H r + E (r) (7)
where

H r ∈ C F×M r is H r def = ∆ c (A R • • • A r+1 A r-1 • • • A 1 ) T (8)
and

M r = R k=1 k r M k . Therefore Y (r) def = [y (r),1 , . . . , y (r),M r ] =         F f =1 h r ( f, 1)a f,r , . . . , F f =1 h r ( f, M r )a f,r         (9)
where h r ( f, m r ) is the ( f, m r ) entry of the matrix H r , for f = 1, . . . , F and m r = 1, . . . , M r .

Simultaneous Sparse Approximation

Assuming 1 that M r > F, ∀r, it is easy to see from ( 9) that the mode coordinates {a f,r } F r f =1 (F r ≤ F) in the r-th dimension are identifiable from any column of Y (r) . This process may be repeated on each dimension r = 1, . . . , R to get all the modes coordinates. In practice, we have to replace the matrix Y (r) by its noisy counterpart Y (r) accounting for the additive white noise. In this case, (9) holds only approximately. Consequently, for each column ỹ(r),m r , m r = 1, . . . , M r , the modal estimation problem can be formulated as a sparse approximation problem corresponding to the following constrained optimization:

x m r = arg min x x 0 subject to ỹ(r),m r -Q r x 2 2 ≤ (10)
where Q r ∈ C M r ×N , N M r , is a (known) modal dictionary, x ∈ C N is a (sparse) vector containing the coefficients of the activated columns in Q r , and is a small reconstruction error related to the noise variance. The pseudo-norm x 0 counts the number of nonzero elements in x. The design of Q r is discussed in section 3.3. The fact that each vector ỹ(r),m r corresponds to a 1-D signal generated by the same modes implies that the position of nonzero entries in x m r should be the same for m r = 1, 2, . . . , M R . Let X be the matrix defined by

X = [x 1 , x 2 , . . . , x M r ], (11)
then the sparsity of X may be measured by computing the Euclidian norms of the rows: those providing a nonzero norm define the rows of active atoms (which are estimations of modes a f,r in the dimension r) in the dictionary Q r . Therefore, we are facing a simultaneous sparse approximation problem:

X r = arg min X X 0,2 subject to Y (r) -Q r X 2 F ≤ (12)
where X 0,2 is the mixed 0 / 2 -norm of X (the number of rows with nonzero 2 -norm). The simultaneous sparse representation models, called also "multiple measurement vectors" (MMV), have been studied from several angles of view, and different approaches have been proposed (see [START_REF] Rakotomamonjy | Surveying and comparing simultaneous sparse approximation (or group-Lasso) algorithms[END_REF] and references therein). As the goal of the present paper is to develop a fast approach well adapted to large signals, we restrict our attention to the SOMP algorithm [START_REF] Tropp | Algorithms for simultaneous sparse approximation[END_REF]. However, it is worth mentioning that, in more intricate cases and/or small size signals, much more efficient simultaneous sparse algorithms may be used at the price of an increased computational burden. A straightforward way to get the R-tuples {(a f,1 , . . . , a f,R )} F f =1 consists in estimating the modes a f,r in the R dimensions using matrices Y (r) , r = 1, . . . , R, which requires a further pairing step to form the R-D modes. To get accurate estimates using the described scheme, two conditions have to be satisfied, 1) the dictionary should contain all modes in the signal, 2) the sparse approximation method should have sufficient guarantees for selecting the true atoms from the dictionary: this is known as "exact recovery guarantees". These problems are discussed in the following sections and an alternative representation of the data is used to avoid the pairing stage in the multiple tones case.

Modal Dictionary Design and Multigrid Strategy 3.3.1. Uniform Modal Dictionary

The dictionary Q r ∈ C M r ×N can be defined from a discretization of the (ν, α) plane. Each point of the grid corresponds to a hypothetic mode. Let N µ be the number of points of a uniform grid covering the frequency interval [0, 1). Similarly, let N β be the number of points of a uniform grid covering the damping factor interval (β min , 0], where β min is a lower bound on {α f,r } F f =1 . Then Q r is given by

Q r = [q r (0, 0), . . . , q r ((N µ -1)δ µ , 0), q r (0, δ β ), . . . , q r ((N µ -1)δ µ , δ β ), . . . , q r ((N µ -1)δ µ , (N β -1)δ β )] ( 13
)
where

q r (µ, β) = a r (µ,β) ||a r (µ,β)|| 2 , a r (µ, β) =
[1, e (β+ j2πµ) , . . . , e (β+ j2πµ)(M r -1) ] T , δ β = β min /N β , and δ µ = 1/N µ . The total number of columns in Q r is N = N µ N β F, each of them is called atom. In the aim of reducing the computational complexity, we propose to estimate frequencies and then damping factors by calling twice the sparse approximation method. At the first step, the frequencies are estimated using a harmonic dictionary. In the second step, the damping factors are estimated using a modal dictionary formed by the already estimated frequencies and a damping factor grid. These two steps are explained in section 4.

Multi-Grid Dictionary Refinement

To achieve a high-resolution modal estimation, a possible way is to define uniform grids as before and selecting very small values for δ µ and δ β . However, the resulting dictionaries will lead to prohibitive calculation cost and memory capacities requested. Rather, we propose to start with a coarse grid (N µ and N β low) and to adaptively refine it through a multigrid scheme as sketched on Figure 1. Let be the current grid level ( = 0, . . . , L -1). At level , we first restore the signal X r ( ) related to the dictionary Q r ( ) by applying the SOMP method. Then we refine the dictionary by inserting atoms inbetween pairs of Q r ( ), in the neighborhood of each activated atom, and we apply again the SOMP method to restore X r ( +1) with respect to the refined dictionary Q r ( + 1). This process is

for i = 1 : numel(Ω) do d i,1 = linspace (d(Ω(i) -1), d(Ω(i)), η) d i,2 = linspace(d(Ω(i)), d(Ω(i) + 1), η) d i = [d T i,1 , d T i,2 (2 : η)] T end d updated = union(d 1 , . . . , d numel(Ω) ) return d updated
repeated until the desired level of resolution is reached. This procedure is applied for both frequencies and damping factors. Algorithm 1 presents the one-step dictionary refinement (DICREF), from level to + 1, where, for a and b reals, linspace(a, b, η) generates a set of η equispaced points in the interval [a, b]. The difference between the present framework and that in [START_REF] Sahnoun | Sparse multidimensional modal analysis using a multigrid dictionary refinement[END_REF] is the following. In [START_REF] Sahnoun | Sparse multidimensional modal analysis using a multigrid dictionary refinement[END_REF] the multigrid algorithm refines jointly R 2-D grids, which leads to expensive computations when R ≥ 2, without convergence guarantees. The present mutigrid scheme refines linear grids, which leads to low computational complexity with convergence guarantees as will be shown in the next section.

Finding the convergence conditions of the new multigrid strategy in the general case (multiple tones) is not easy and depends on the selected sparse approximation algorithm. By contrast, it is possible to show that, under mild conditions, the convergence may be guaranteed in the single tone case. This issue is discussed in the next section. In section 5, we make use of an alternative representation of the data model in the case of multiple tones and we propose a method allowing one to retrieve the signal tones separately.

Single R-D Mode Estimation

In the previous section, we have shown how the R-D modal retrieval problem may be tackled using a sparse approximation algorithm by estimating the set of parameters in each dimension r = 1, . . . , R. Here, we give the sufficient conditions for convergence of the multigrid dictionary refinement for F = 1.

Without loss of generality, we set R = 1. For notation simplicity, we omit reference to the dimension index r.

According to (3), the 1-D modal signal containing a single mode can be written as follows:

y(m) = c 1 a m-1 1 = c 1 e (α 1 + j2πν 1 )(m-1) , m = 1, . . . , M. (14)
Let Q be a normalized modal dictionary Q = [q 1 , . . . , q N ], with

q n = 1 m |q n | 2m [1, q n , . . . , q M-1 n ] T , (15)
q n = exp(β n + j2πµ n ), µ n ∈ [0, 1), β n ∈ (β min , 0], for n = 1, . . . , N.
The single tone sparse approximation of y with respect to Q is the solution of the criterion:

min x J(x) = ||y -Qx|| 2 s.t. ||x|| 0 = 1. ( 16
)
The optimal solution is given by

x * n = q H n y, x * {1,...,N}\n = 0, J(x * ) = ||y|| 2 -y H q n q H n y ( 17
)
where n is the selected column number in Q. Finally, the minimum J(x * ) is reached for an atom q n that maximizes

J (q n ) = y H q n q H n y = |q H n y| 2 , n = 1, . . . , N.
Remark 1. At this point we notice that, for 1-D single-tone harmonic signals, maximizing J (q n ) leads to the well-known beamforming method (or more precisely to the periodogram in the single snapshot case). Hence, an estimate of the frequency ν 1 may be obtained from the peak of J (q n ). Here, as in [START_REF] Duarte | Spectral compressive sensing[END_REF], we use the sparse approximation framework. In our case, it allows to estimate the frequency and then the damping factor from a unified point of view and using a unique algorithm based on SOMP and DICREF.

Estimating the Frequency: The Harmonic Dictionary

First, we estimate frequency ν 1 using a harmonic dictionary (i.e. assuming β n = 0, ∀n). In this case, we have:

J (µ n ) = |c 1 | 2 M 1 -e α 1 M+ j2π(ν 1 -µ n )M 1 -e α 1 + j2π(ν 1 -µ n ) 2 . ( 18
)
The following theorem gives a sufficient condition for the multigrid dictionary refinement scheme to converge to the global maximum of J .

Theorem 1. Let y(m) be a single tone (F = 1) noiseless signal of length M and Q( = 0) = [q 1 q 2 . . . q N(0) ] T be the initial harmonic dictionary in which the columns are sorted in increasing order of µ n (0), n = 1, 2, . . . , N(0) and covering the frequency interval [0,1): µ 1 (0) = 0 and µ N(0) (0) = 1 -1/M. Then the refinement scheme is convergent (i.e. ∃n ∈ {1, . . . , N( )} s.t. lim →∞ µ n ( ) = ν 1 ) if the following condition is satisfied:

max n∈{1,...,N(0)-1} |µ n+1 (0) -µ n (0)| < 2ζ M ( 19
)
where ζ M is a constant depending only on M. Proof It is easy to check that the global maximum of J (µ n ) is reached for µ n = ν 1 , ∀α 1 . Figure 2 shows the shape of J (µ n ) for β n = 0. For α 1 = 0, J (µ n ) reduces to a Fejér kernel of order M which has exactly one local maximum in the interval [ν 1 + k/M, ν 1 + (k + 1)/M], k 0. Let J 1 be the amplitude of the first sidelobe and

ν 1 -1 M ν 1 ν 1 + 1 M ν 1 + 2 M 0 1 4 1 2 3 4 1 ζ M µ n J (µ n )/||y|| 2 α 1 = 0.0 α 1 = -0.3
ν 1 + ζ M be the value of µ n such that J (µ n = ν 1 + ζ M ) = J 1 in the interval [ν 1 , ν 1 +1/M] (we assume 2 that M > 2).
For the dictionary refinement strategy to converge to the global maximum, it is sufficient to the sparse approximation algorithm to select, at a given level , an atom whose frequency satisfies

|µ n * ( ) -ν 1 | < ζ M < 1/M, where µ n * ( ) = arg max n J (µ n ). Indeed, if µ n * ( ) ∈ (ν 1 -ζ M , ν 1 + ζ M )
then adding two atoms whose frequencies are located on both sides of µ n * ( ) will lead to the selection, at level + 1, of an atom that satisfies

|µ n * ( + 1) -ν 1 | ≤ |µ n * ( ) -ν 1 |:
the distance between the selected atom and the true frequency is a monotonically decreasing sequence. Finally, the convergence is guaranteed if the initial dictionary contains an atom n such that |µ n (0

)-ν 1 | < ζ M , which is satisfied if max n∈{1,...,N(0)-1} |µ n+1 (0) -µ n (0)| < 2ζ M . ( 20
)
given the fact that the sequence {µ n (0)} covers the interval [0, 1). For α 1 < 0, the main lobe of J (µ n ) becomes broader and ζ M larger than for α 1 = 0. Consequently, condition (20) is also sufficient for α 1 < 0.

Corollary 1. In the single tone case, the harmonic dictionary refinement is convergent if the initial frequency grid ( = 0) is the Fourier grid.

Proof Fourier bins are obtained for N = M and µ n (0

) = (n - 1)/M. Since ζ M > 1/2M, the proof is straightforward.
It is important to note that condition (20) is sufficient but not necessary. Moreover, this condition is established when adding a single atom on both sides of the selected one (i.e. η = 1 in Algorithm 1). When η 1, the condition may be relaxed and the rate of convergence is expected to be higher.

Estimating the Damping Factor: The Modal Dictionary

Assume that the previous sparse approximation using a harmonic dictionary has converged to select an atom with µ n = ν 1 . Now, we have to estimate the damping factor α 1 . We form a modal dictionary using the damping factor grid and the frequency ν 1 , i.e. q n = exp(β n + j2πν 1 ). Consequently,

J (β n ) = |c 1 | 2 (1 -e 2β n ) 1 -e 2β n M 1 -e (α 1 +β n )M 1 -e (α 1 +β n ) 2 . ( 21
)
Theorem 2. Let y(m) be a single tone (F = 1) noiseless signal of length M and Q(0) = [q 1 q 2 . . . q N(0) ] T be the initial modal dictionary formed using the frequency ν 1 , i.e., q n = exp(β n (0) + j2πν 1 ), where ν 1 is the frequency of signal y. The columns are sorted in increasing order of β n (0), n = 1, 2, . . . , N and covering the damping factor interval (β min , 0]. Then the refinement scheme is convergent (i.e. ∃n s.t.

lim →∞ β n ( ) = α 1 ) if α 1 ∈ (β min , 0].
Proof Let g(β n ) be the derivative of J (β n ) in ( 21) with respect to β n . It is easy to check that g(β n ) > 0 for β n < α 1 , g(β n ) < 0 for β n > α 1 , and g(β n ) = 0 when β n = α 1 . In other words, J (β n ) is monotonically increasing before the maximum reached at α 1 and monotonically decreasing after α 1 . Therefore, the multigrid algorithm converges to α 1 if

β min < α 1 .
As a consequence of Theorem 2, the initial modal dictionary can be formed using only two points in the damping factor grid:

β 1 (0) = β min and β 2 (0) = 0.
We can now state that the multigrid algorithm based on two sparse approximations (for frequency and then damping factor) converges in the single tone case under some conditions. Note that in the noisy case when the SNR is sufficiently high, the convergence analysis is still valid as in the noiseless case, and the proposed multigrid sparse scheme for single tone converges to the global maximum of the Fejér kernel. The extension to the single tone R-D modal retrieval problem is straightforward and can be performed according to the formulation presented in Section 3.2. The details of this approach (STSM: Single Tone Sparse Method) are presented in Algorithm 2. The algorithm takes as input a noisy single tone R-D signal, and a couple of integers η ν and η α that correspond respectively to the number of frequency and damping factor atoms to be added on both sides of the corresponding selected ones. Next, for each dimension r = 1, . . . , R, we run two tasks to estimate the frequency and then the damping factor: in each step we apply SOMP combined to DICREF using corresponding dictionaries and taking into account the convergence conditions discussed previously. Then parameters of a r , i.e., ν r and α r , are given by the corresponding selected atoms.

Multiple R-D Modes Estimation

In the multiple tones case, sparse approximation algorithms yield suboptimal solutions when the coherence of the dictionary is high [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF]. This is a crucial point because the refinement procedure will increase the coherence with increasing , which may prevent convergence even in the noiseless case. In the following, we present a low complexity algorithm that is accurate and robust in the presence of noise. The idea is to begin by an Algorithm 2: Single tone sparse method (STSM) based on a multigrid refinement

input : A tensor Y ∈ C M 1 ו••×M R , (η ν , η α ) ∈ N × N output: Parameters of the single R-D mode: a 1 , . . . , a R initialization: (k ν , k α ) = (0, 0) initialize d (0)
ν and d (0) α using ζ for r = 1 : R do while halting criterion false do

k ν = k ν + 1 Ω (kν) ν = SOMP(Q(d (kν) ν , 0), Y (r) , Iter = 1) d (kν+1) ν = DICREF(d (kν) ν , Ω (k) ν , η ν ) end while halting criterion false do k α = k α + 1 Ω (kα) α = SOMP Q d (kν) ν (Ω (kν) ν ), d (kα) α , Y (r) d (kα+1) α = DICREF(d (kα) α , Ω (k) α , η α ) end a r = exp(d (kα) α (Ω (kα) α ) + 2π d (kν) ν (Ω (kν) ν )) end return a 1 , . . . , a R
initialization step where F single tone modal signals of order R -1 are extracted from the multiple tones R-D signal. Then an iterative technique is proposed to improve this decomposition and estimate more accurately the underlying parameters.

It is assumed that the frequencies are distinct in at least one dimension with M r > F. Then dimensions are permuted such that the dimension with distinct frequencies becomes the first one (r = 1).

From Multiple Tones to Multiple Single-Tone Signals

According to [START_REF] Sacchini | Two-dimensional Prony modeling and parameter estimation[END_REF], Y can be written as [START_REF] Comon | Tensors: a brief introduction[END_REF][START_REF] Kolda | Tensor decompositions and applications[END_REF]

Y = I • 1 A 1 • 2 • • • • R A R • R+1 c T (22) = I • 2 A 2 • 3 • • • • R A R • R+1 c T • 1 A 1 = S • 1 A 1 ( 23
)
where I ∈ C F×Fו••×F is a tensor of order R + 1 containing ones on its diagonal and zeros otherwise, and

S = I • 2 A 2 • 3 • • • • R A R • R+1 c T ( 24
)
is a complex tensor of order R and size

F × M 2 ו • • M R .
Similar expressions are evoked in, among others, [START_REF] Sun | Accurate and computationally efficient tensor-based subspace approach for multi-dimensional harmonic retrieval[END_REF]. The matricization of S along the first dimension is given by [30]

S (1) = I (1) (c T A R • • • A 2 ) T (25)
= I F (c T A R • • • A 2 ) T (26)
where I (1) and I F denote respectively the matricization of I along the first dimension and the identity matrix of size F × F.

Then,

S (1) = (c T A R • • • A 2 ) T (27) =                  c 1 (a 1,R • • • a 1,2 ) T c 2 (a 2,R • • • a 2,2 ) T . . . c F (a F,R • • • a F,2 ) T                  (28)
We can see that each row f = 1, . . . , F of S (1) represents the matricization along the first dimension of a sub-tensor

S f of size 1 × M 2 × • • • × M R that contains a single (R -1)-D tone,
where:

S f = c f a f,2 ⊗ a f,3 ⊗ • • • ⊗ a f,R . (29)
The tensor S can then be written as the concatenation of the F sub-tensors S f along the first dimension

S = S 1 1 S 2 1 • • • 1 S F (30)
This is the key property to transform the multiple tone R-D signal into F single tone signals, thus allowing to use the singletone estimation of Section 4. To do so, we need firstly to estimate A 1 . Now, we show how to estimate the modes of the first dimension, i.e. A 1 , using the matricized form of Y along the first dimension Y (1) (equation ( 6)). The singular value decomposition (SVD) of Y (1) yields

Y (1) = UΣV H (31)
where matrices U ∈ C M 1 ×L and V ∈ C M 1 ×L are orthonormal and contain respectively the left and right singular vectors of Y (1) , with L = min{M 1 , M 1 }. Σ is a diagonal matrix containing the singular values σ i , i = 1, . . . , L, sorted in a decreasing order. As the number of components in Y is equal to F, we can decompose the SVD of Y (1) in [START_REF] Rakotomamonjy | Surveying and comparing simultaneous sparse approximation (or group-Lasso) algorithms[END_REF] as follows

Y (1) = U F Σ F V H F + U n Σ n V H n ( 32
)
where U F (resp. V F ) stands for the matrix formed with the first F columns of U (resp. V) and Σ F contains nonzero singular values on its diagonal Σ F = diag(σ 1 , . . . , σ F ). U n (resp. V n ) is formed by the remaining columns associated with zero singular values Σ n = 0. It can be established from ( 7) and ( 32) that A 1 and U F span the same subspace, and thus there exists an unknown nonsingular matrix T that satisfies

A 1 = U F T. (33)
Denote by M (resp. M) the matrix obtained from M by deleting the first (resp. last) row. By harnessing the Vandermonde structure of A 1 , there exists a diagonal matrix

D = diag(a 1,1 , . . . , a F,1 ) such that A 1 = A 1 D. Since A 1 = U F T and A 1 = U F T, then U F T = U F
TD, which proves that matrix T can be estimated by the eigenvectors of U † F U F . Therefore, using this estimate of T, exponentials of A 1 can be estimated by [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF] where U F is obtained from [START_REF] Duarte | Spectral compressive sensing[END_REF].

In the presence of noise, the SVD of Y (1) is given by

Y (1) = ŨF ΣF ṼH F + Ũn Σn ṼH n . (34)
Due to the noise, all the quantities in the right-hand side (RHS) of ( 34) may be perturbed versions of those in the RHS of [START_REF] Duarte | Spectral compressive sensing[END_REF].

In this case we may express ŨF by ŨF = U F + ∆U F . Then, neglecting Σn , an approximation of Y (1) , denoted by Ŷ(1) , can be obtained using the first F principal components of the SVD of Y (1) :

Ŷ(1) = ŨF ΣF ṼH F . (35)
Thereby S can be estimated from the noisy data and Â1 using equation [START_REF] Adalbjörnsson | High resolution sparse estimation of exponentially decaying two-dimensional signals[END_REF] as follows

Ŝ = Y • 1 Â † 1 , (36)
then Ŝ f , f = 1, . . . , F are extracted from Ŝ according to [START_REF] Kolda | Tensor decompositions and applications[END_REF].

Each Y f = c f a f,1 ⊗ • • • ⊗ a f,R can be estimated by Ȳ(0) f = Ŝ f • 1 â f,1 .
The sparse multigrid algorithm for single tone (STSM) can be applied on each Ȳ(0) f , f = 1, . . . , F to estimate the parameters. However, we propose in the following to improve the separated components using an iterative technique.

Improving the Estimation Accuracy

It is clear from [START_REF] Kay | Fundamentals of statistical signal processing. Estimation theory[END_REF] that, in the noisy case, the error in estimating S (due to the estimation of A 1 ) will propagate when estimating the parameters a f,2 , . . . , a f,R . Hence, we propose to improve iteratively the mode estimates. The following procedure is used to update estimates at each iteration i = 0, . . . , K

apply STSM to estimate a

f,2 , . . . , a f,R , f = 1, . . . , F {â f,2 , . . . , â f,R } = STSM( Ȳ(i) f , η ν , η α , r = 2, . . . , R) (37)
2. estimate c f a f,1 , f = 1, . . . , F by least squares using the already estimated a f,2 , . . . ,

a f,R , f = 1, . . . , F c f a f,1 = Ȳ(i) f (1) (â f,R • • • â f,2 ) T † (38) 3. compute Ŷ(i) f Ŷ(i) f = c f a f,1 ⊗ â f,2 ⊗ • • • ⊗ â f,R (39)
where

Ȳ(i) f = Ŷ(i-1) f + R (i) f -1 , R (i) ( f ) = R (i) f -1 + Ŷ(i-1) f - Ŷ(i) f , f = 1, . . . , F, R (i) 0 def = R (i-1) F , and R (0) F = Ỹ -F f =1
Ŷ(0) f . This iterative scheme will be analyzed in the next section.

Finally, the algorithm we propose (MTM: Multiple Tones Method) is summarized in Algorithm 3. Note that no association step of R-D modes is required. The initialization step consists in initializing: i) Â1 and Ŝ using [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF], [START_REF] Kay | Fundamentals of statistical signal processing. Estimation theory[END_REF] and [START_REF] Kolda | Tensor decompositions and applications[END_REF], ii) the estimated single tones Ȳ(0) f , f = 1, . . . , F.

Algorithm 3: Multiple tones method (MTM)

input : A tensor Y ∈ C M 1 ו••×M R , (η ν , η α ) ∈ N × N output: Parameters of the multiple R-D modes : {a f,r } F,R f =1,r=1
initialization:

1. Compute Â1 and Ŝ f , f = 1, . . . , F using ( 33), ( 36) and ( 30)

2. Ȳ(0) f = Ŝ f • 1 â f,1 , f = 1, . . . , F
For f = 1, . . . , F, compute Ŷ(0) f using ( 37), ( 38) and ( 39) 37), ( 38) and ( 39)

R (0) F def = R (1) 0 = Ỹ -F f =1 Ŷ(0) f for i = 1 : K do for f = 1 : F do Ȳ(i) f = Ŷ(i-1) f + R (i) f -1 compute Ŷ(i) f using (
R (i) f = Ȳ(i) f - Ŷ(i) f , if f = F, then R (i+1) 0 def = R (i) F end end For f = 1, . . . , F, extract a f,1 using a f,1 = STSM( Ŷ(K) f + R (K) F , η ν , η α , r = 1) return {â f,r } F,R f =1,r=1
Note that the columns of Â1 are iteratively updated without extracting the related modes, whereas the modes of the other dimensions are extracted at each iteration using [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. Solely after the last iteration (i = K), the parameters of the first dimension are extracted using STSM algorithm. K denotes the maximum number of iterations, which is fixed to 2 in the simulations since no improvement was observed for K > 2.

Analysis of the Algorithm

Following the separation step described in ( 31)-( 36), we can state that the algorithm yields the expected solution when the SNR is sufficiently high. We want now to prove that the second stage (next iterations), in addition to estimating the parameters from the single tones, is also improving the estimation accuracy. The general idea is inspired from greedy forward/backward sparse approximation, where the solution is refined by adding/removing atoms to/from the set of activated atoms. The improvement of the estimates is stated by the following theorem.

Theorem 3. Assuming that the noise E is sufficiently small such that the ordering of the singular values in Σ in [START_REF] Rakotomamonjy | Surveying and comparing simultaneous sparse approximation (or group-Lasso) algorithms[END_REF] is the same as the ordering of the corresponding singular values when E = 0. Using the procedure expressed by ( 37), ( 38) and [START_REF] Djermoune | Perturbation analysis of subspace-based methods in estimating a damped complex exponential[END_REF] to estimate Y f at iteration i = 0, . . . , K

Ŷ(i) f = arg min X∈H Ȳ(i) f -X ( 40
)
where

H = {X ∈ C M 1 ו••×M R |X = b 1 ⊗ b 2 ⊗ • • • ⊗ b R , b r ∈ P for r 1} with P = {v ∈ C M r |v = [1, v, . . . , v M r -1 ] T , v = exp(β + jω), β ∈ R -, ω ∈ [0, 2π)}.
Then, at each iteration i, the residual is decreased:

Ỹ - Ŷ(i) ≤ Ỹ - Ŷ(i-1) (41)
where

Ŷ(i) = F f =1 Ŷ(i) f . Proof See Appendix A.

Identifiability

Based on the assumptions under which Algorithm 3 is operating, the identifiability condition can be stated as F < M 1 and min{M 2 , ...., M R } ≥ 2. In [START_REF] Jiang | Almost-sure identifiability of multidimensional harmonic retrieval[END_REF], the condition is M r ≥ 4, r = 1, . . . , R, and

F ≤ M 1 2 R r=1 M r 2 .
We note that, when M r ≥ 4, r = 1, . . . , R, the number of identifiable modes is slightly smaller than in [START_REF] Jiang | Almost-sure identifiability of multidimensional harmonic retrieval[END_REF], but the proposed algorithm is able to outperform the conventional methods in terms of computational complexity and accuracy. In addition, another advantage of the proposed algorithm is clear when the number of samples in one or more dimensions is less than 4 (i.e. M r < 4), where identifiability in [START_REF] Jiang | Almost-sure identifiability of multidimensional harmonic retrieval[END_REF] is not satisfied. This latter case (i.e. ∃r, M r < 4) can be encountered in signal processing applications when the size of one or more diversities (dimensions in our formulation problem) is less than 4.

Cramér-Rao Lower Bounds for R-D Cisoids in Noise

In this section, we derive the expressions of the CRLB for the parameters of R-D damped exponentials in Gaussian white noise. We then give the CRLB in the cases of single damped and undamped R-D cisoids. We consider the R-D sinusoidal model given in [START_REF] Nion | Tensor algebra and multidimensional retrieval in signal processing for MIMO radar[END_REF]. Let

θ = [ω 1,1 . . . ω 1,R ω 2,1 . . . ω F,R α 1,1 . . . α 1,R α 2,1 . . . α F,R λ 1 . . . λ F φ 1 . . . φ F ] T
be the unknown parameter vector. The aim here is to derive the CRLB of the parameters in θ.

The joint probability density function (pdf) of ỹ is

p(ỹ; θ) = 1 (σ 2 π) M exp - 1 σ 2 (ỹ -µ(θ)) H (ỹ -µ(θ)) (42)
where µ(θ) is the noise-free part of y and ỹ = [ỹ(1, . . . , 1, 1), . . . , ỹ(1, . . . , 1, M R ), ỹ(1, . . . , 2, 1),

. . . , ỹ(1, . . . , 2, M R ), . . . , ỹ(M 1 , . . . , M R )] T . ( 43
)
The i-th entry of µ(θ) can be written as:

µ(θ) i = F f =1 c f R r=1 a t i,r f,r , (44)
for i = 1, . . . , M, where

t i,r =          i -1 R =r+1 M          mod M r , ( 45
)
and • is the floor function. In the following, we derive the expressions of the CRLB in the general case (F > 1) and then we deduce the result corresponding to a single R-D modal signal (F = 1).

Derivation of the CRLB

Given the joint pdf in (42), the (k, l) entry of the Fisher information matrix is [START_REF] Yao | Cramér-Rao lower bounds for a damped sinusoidal process[END_REF][START_REF] Kay | Fundamentals of statistical signal processing. Estimation theory[END_REF]:

[F(θ)] kl = 2 σ 2 Re        ∂µ(θ) ∂θ k H ∂µ(θ) ∂θ l        . ( 46
)
After some lengthy calculations, the M × (2RF + 2F) matrix ∂µ(θ)/∂θ can be expressed as

∂µ(θ) ∂θ = [ jZ Φ Z Φ Zφ jZφ] V • blkdiag(Λ, Λ, I F , λ) S ( 47
)
where

Z = [Z 1 , . . . , Z F ] ∈ C M×RF , with Z f (i, l) = t i,l R r=1 a t i,r f,r , (48)
Λ = blkdiag(λ 1 I R , . . . , λ F I R ) ∈ R RF×RF , (49)
Φ = blkdiag(e jφ 1 I R , . . . , e jφ F I R ) ∈ C RF×RF , ( 50
) Z = [z 1 , . . . , z F ] ∈ C M×F , with z f (i) = R r=1 a t i,r f,r , (51)
λ = diag([λ 1 , . . . , λ F ]) ∈ R F×F , (52)
φ = diag([e jφ 1 , . . . , e jφ F ]) ∈ C F×F . (53)
Finally, the inverse of the Fisher information matrix is

F -1 (θ) = σ 2 2 S -1 Re{V H V} -1 S -1 = σ 2 2 S -1 WS -1 (54)
where Re{•} stands for the real part. The CRLB of θ k is given by [F -1 (θ)] kk . More explicitly, for f = 1, . . . , F and r = 1, . . . , R:

CRLB(ω f,r ) = 2σ 2 W R( f -1)+r,R( f -1)+r λ 2 f (55) CRLB(α f,r ) = 2σ 2 W RF+R( f -1)+r,RF+R( f -1)+r λ 2 f (56) CRLB(λ f ) = 2σ 2 W 2RF+ f,2RF+ f (57) CRLB(φ f ) = 2σ 2 W 2RF+F+ f,2RF+F+ f λ 2 f ( 58
)
Theorem 4. For the general R-D exponential process, the CRLB's for f = 1, . . . , F and r = 1, . . . , R satisfy

CRLB(ω f,r ) = CRLB(α f,r ) (59) CRLB(λ f ) = λ 2 CRLB(φ f ) ( 60
)
Proof It is based on the special block structure of matrix Re{V H V} (see for instance [START_REF] Yao | Cramér-Rao lower bounds for a damped sinusoidal process[END_REF]).

Single Mode Case

In this section, the CRLB's will be simplified in the case of a single R-D modal signal (F = 1) to obtain more precise details on their parameter dependency. For the sake of simplicity, the subscripts denoting the mode f = 1 will be omitted. First, assume that |a r | = exp(α r ) < 1. We shall express the products Z H Z , Z H Z and Z H Z. After some calculations, we get:

[Z H Z ] nk = R r=1 r n,k 1 -|a r | 2M r 1 -|a r | 2 ×                    M n -1 m=0 m|a n | 2m M k -1 m=0 m|a k | 2m , if n k M n -1 m=0 m 2 |a n | 2m , if n = k (61)
Z H Z = R r=1 1 -|a r | 2M r 1 -|a r | 2 (62) [Z H Z] n = R r=1 r n 1 -|a r | 2M r 1 -|a r | 2 × M n -1 m=0 m|a n | 2m . (63)
Denoting

M (α) = R r=1 (1 -|a r | 2M r )/(1 -|a r | 2 ), q 1 (n) = M n -1 m=0 m|a n | 2m / M n -1 m=0 |a n | 2m and q 2 (n) = M n -1 m=0 m 2 |a n | 2m / M n -1 m=0 |a n | 2m
, we then obtain:

[P] nk = M (α) ×        q 1 (n)q 1 (k), if n k q 2 (n), if n = k (64)
G = M (α) (65) [Q] n = M (α) q 1 (n), (66)
and

Re{V H V} =                P 0 0 Q 0 P Q 0 0 Q T G 0 Q T 0 0 G                . ( 67
)
The inversion of Re{V H V} yields the following expressions of the CRLB's:

CRLB(ω r ) = CRLB(α r ) = σ 2 2λ 2 M (α) × (1 -|a r | 2 ) 2 (1 -|a r | 2M r ) 2 -M 2 r |a r | 2M r (1 -|a r | 2 ) 2 + |a r | 2 (1 -|a r | 2M r ) 2 , (68) CRLB(λ) λ 2 = CRLB(φ) = σ 2 2λ 2 M (α)        1 + R r=1 q 2 1 (r) q 2 (r) -q 2 1 (r)        .
(69) Finally, for a single R-D purely harmonic signal (α r = 0, ∀r), we have M (α) = R r=1 M r = M and taking the limit of the CRLB's when α r → 0 leads to:

lim α r →0 CRLB(ω r ) = 6σ 2 λ 2 M(M 2 r -1) (70) lim α r →0 CRLB(λ) λ 2 = σ 2 2λ 2 M        1 + 3 R r=1 M r -1 M r + 1        . (71)
Hence, for the undamped case, our result in (70) is consistent with [START_REF] Sun | Accurate and computationally efficient tensor-based subspace approach for multi-dimensional harmonic retrieval[END_REF].

Simulation Results

Numerical simulations have been carried out to assess the performances of the proposed method for 2-D and 3-D modal signals in the presence of white Gaussian noise. The performances are measured by the total root-mean square error (RMSE) on estimated parameters and the computational time. The total RMSE is defined as 2 where ξ f,r is an estimate of ξ f,r , and E p is the average on p Monte-Carlo trials. In our simulations, ξ f,r can be either a frequency or a damping factor.

RMSE total = 1 RF E p R r=1 F f =1 (ξ f,r -ξ f,r )

RMSE for 2-D and 3-D Signals

Experiment 1. To show the effectiveness of the multigrid scheme, this experiment presents the results obtained on Signal #1 (see Table 1) with different multigrid levels and different initial grids. Signal #1 is a single tone 2-D modal signal of size 10 × 10. The number of multigrid levels is fixed to L = 2, i.e., = 0, 1, 2. Then the results are presented as a function of the number of atoms in the initial dictionaries N(0) and the number of atoms η ν or η α added at each level . The results we obtain for the first step, i.e., for the harmonic estimation, are presented in Figure 3. We can observe that the frequency RMSE obtained with the R-D sparse algorithm can reach the CRLB using a uniform initial harmonic dictionary of 10 atoms and η ν = 31 (Figure 3.a). Figure 3.b shows that the frequency RMSE is improved at low SNR if the initial dictionary contains 20 atoms, and reaches the optimal estimates with η ν = 21. Figure 4 shows the damping factor RMSE obtained by R-D sparse algorithm using different settings of the initial damping factor dictionary and η α . We can observe that the damping factor RMSE depends on the number of atoms in the dictionary, the more atoms the better. At low SNR, the RMSE also depends β min . Therefore, it is better to choose β min with small absolute value if we have a prior knowledge of the interval of damping factors in the signal. In general, the estimation error is of order 1 N(0)η 2 . For instance, in the frequency step estimation, we recommend to chose N(0) to be greater than or equal to 3 2 M r if we want a good accuracy at lower SNR levels. Otherwise, we can set N(0) = M r . Once N(0) is set, η can be chosen with respect to the desired accuracy. Let ε be the desired estimation error, then ε = 1 N(0)η 2 and we can set η = 1 √ εN(0) .

In the rest of this section, the proposed algorithms are compared with 2-D ESPRIT [START_REF] Rouquette | Estimation of frequencies and damping factors by two-dimensional ESPRIT type methods[END_REF], Tensor-ESPRIT [START_REF] Haardt | Higher-order SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems[END_REF], PUMA [START_REF] So | An efficient approach for twodimensional parameter estimation of a single-tone[END_REF] and TPUMA [START_REF] Sun | Accurate and computationally efficient tensor-based subspace approach for multi-dimensional harmonic retrieval[END_REF]. If the R-D signal contains one tone then Algorithm 2 (STSM) is used, otherwise Algorithm 3 ( MTM) is used. Thus, to facilitate notation, both proposed algorithms, Algorithm 2 and Algorithm 3, will be called R-D sparse. For the proposed method, the initial grid used to build the harmonic dictionary is the same for all dimensions; it contains 50 frequency points uniformly distributed over the interval [0, 1) and 10 damping factors β ∈ [-0.05, 0]. To simulate a random dictionary, at each run, the frequency grid is perturbed by a small random quantity. As a consequence of experiment 1, we use the Since the proposed method is applied directly on data without using spatial smoothing, i.e., it does not require the construction of a large matrix or an augmented order tensor, then a relevant comparison will be with algorithms that do not use spatial smoothing. Thereby, in the next experiments, the proposed algorithm is compared to PUMA [START_REF] So | An efficient approach for twodimensional parameter estimation of a single-tone[END_REF] and TPUMA [START_REF] Sun | Accurate and computationally efficient tensor-based subspace approach for multi-dimensional harmonic retrieval[END_REF], which are algorithms that do not require spatial smoothing. We also report comparisons with 2-D ESPRIT [START_REF] Rouquette | Estimation of frequencies and damping factors by two-dimensional ESPRIT type methods[END_REF] and Tensor-ESPRIT [START_REF] Haardt | Higher-order SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems[END_REF], which need spatial smoothing.

Single tone R-D modal signal

Experiment 2. This experiment tends to show the efficiency of the proposed algorithm in estimating parameters of single tone R-D modal signals. We use the same signal as before (Signal #1). Our R-D sparse algorithm is compared to 2-D ESPRIT [START_REF] Rouquette | Estimation of frequencies and damping factors by two-dimensional ESPRIT type methods[END_REF] and PUMA [START_REF] So | An efficient approach for twodimensional parameter estimation of a single-tone[END_REF]. For each level of noise, 1000 Monte-Carlo trials are performed. Figure 5 shows the obtained results. We can observe that: i) the proposed algorithm and PUMA reach the CRLB and outperform 2-D ESPRIT, ii) R-D sparse outperform PUMA in SNR less than 3 dB.

Multiple tones R-D modal signals

Several configurations are studied in the case of multiple tones to compare the proposed algorithm with Tensor-ESPRIT [START_REF] Haardt | Higher-order SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems[END_REF] and TPUMA [START_REF] Sun | Accurate and computationally efficient tensor-based subspace approach for multi-dimensional harmonic retrieval[END_REF]. These configurations (Experiments 3-4) are summarized in Table 2, in which the number of modes and the distance between frequencies in different dimensions are varied. ∆ Fr denotes the Rayleigh frequency resolution limit, which has the same value in all dimensions because M 1 = M 2 = M 3 . In Experiment 5 we examine the case when the size of only one dimension is larger than 4, i.e., the identifiability condition of [START_REF] Jiang | Almost-sure identifiability of multidimensional harmonic retrieval[END_REF] is not satisfied. The parameters of the used signals are given in Table 1. Experiment 3. In this experiment, we simulate a 3-D signal (Signal #2) of size 8×8×8 and containing two modes whose frequencies in each dimension are well separated. Parameters of the signal are given in Table 1. Figure 6 shows the obtained results. Here, the proposed method performs as TPUMA. Tensor-ESPRIT yields a slightly higher RMSE. 8. We observe that the proposed method outperforms TPUMA mainly at low SNR levels.

Numerical Complexity

It is known that in the case of 1-D signals of size M, OMP costs O(NF M) in terms of multiplications [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]; F is the sparsity (number of components) and N is the number of atoms in the dictionary. For a M-measurements R-D signal, the complexity of the STSM algorithm over a set of L multigrid levels is O(MNLR), assuming that the number of dictionary atoms is maintained constant (equal to N) over all levels. Regarding the approach proposed in Algorithm 3, the main operations are the call of STSM and the update

of c f a f,1 = Ȳ(i) f (1) (â f,R • • • â f,2 ) T † which has a complex- ity of O(M) since (â f,R • • • â f,2 ) T † is a row vector of length R r=2 M r and Ȳ(i) f (1
) is a matrix of size M 1 × R r=2 M r . Therefore, the whole complexity of the proposed algorithm is O ((NL(F(R -1)K + 1) + FK)M), which is linear in the number of measurements M. The complexity of the Tensor-ESPRIT algorithm with spatial smoothing is mainly related to that of the SVD which is at least O(k t F(R + 1)PM) where k t is a constant depending on the implementation of the SVD algorithm. Here P = R r=1 P r where {P r } R r=1 are design parameters used to get smoothed measurements (see [START_REF] Haardt | Higher-order SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems[END_REF]). The accuracy of the estimates provided by ESPRIT depends on these parameters. Since the optimal value for P r is a fraction of M r (e.g. [START_REF] Lemma | Analysis of joint anglefrequency estimation using ESPRIT[END_REF][START_REF] Djermoune | Perturbation analysis of subspace-based methods in estimating a damped complex exponential[END_REF][START_REF] Hua | Matrix pencil method for estimating parameters of ex-ponentially damped/undamped sinusoids in noise[END_REF]), the complexity of the SVD step is, in fact, O(M 2 ). The complexities of PUMA and TPUMA algorithms are O(M 3 ) and O(k t M(R + F -1)) + R r=1 O(K(F + 1)M 3 r ), respectively. Compared to PUMA and TPUMA, the proposed algorithm has an attractive computational complexity for large size signals. We observe that the proposed method involves a low computational complexity compared to TPUMA and Tensor-ESPRIT when M 1 is large.

Conclusion

We presented an efficient sparse estimation approach for the analysis of multidimensional (R-D) damped or undamped modal signals. The idea consists in exploiting the simultaneous sparse approximation principle to separate this joint estimation problem into R multiple measurements problems. To be able to handle large size signals and yield accurate estimates, a multigrid dictionary refinement scheme is associated with the simultaneous orthogonal matching pursuit (SOMP) algorithm. We gave the convergence proof of the the refinement procedure in the single tone case. Then, for the general multiple tones R-D case, the signal tensor model is decomposed in order to handle each tone separately in an iterative scheme so that the pairing of the R-D parameters is automatically achieved. Also, the CRLB of the R-D modal signal parameters were derived. The tests performed on simulated signals showed that the proposed algorithm attains the CRLB and outperforms state-of-the-art subspace algorithms. We also have shown that the complexity of the algorithm is linear with respect to the number of measurements, which allows the processing of large size signals. Finally, it is worth mentioning that this approach can be straightforwardly applied to other multidimensional array processing problems.

Appendix A. Proof of Theorem 3

We begin the proof by introducing the following lemma.

Lemma 1. Consider Y = Y + ∆Y, where Y is the perturbed version of the data tensor Y and ∆Y is the perturbation. Assuming that ∆Y is sufficiently small such that the ordering of the F singular values in Σ in [START_REF] Rakotomamonjy | Surveying and comparing simultaneous sparse approximation (or group-Lasso) algorithms[END_REF] is the same as the ordering of the corresponding singular values when ∆Y = 0. Then the perturbation ∆Y f contains a linear combination of all Y f , f = 1, . . . , F:

∆Y f = D f + F i=1 v f,i Y i where v T f = [v f,1 , . . . , v f,F ] = ∆A † 1 ( f, :)A 1 and D f = ∆Y • 1 a f,1 A † 1 ( f, :) + Y s, f • 1 ∆a f,1 .
Proof From (36) S = Y •

∆Y f = F p=1 v f,p Y p + S f • 1 ∆a f,1 + ∆Y • 1 a f,1 A † 1 ( f, :). 12
Using the previous lemma

Ȳ(0) f = S f • 1 (a f,1 + v f, f a f,1 + ∆a f )+ F p=1,p f v f,p Y p + ∆Y • 1 a f,1 A † 1 ( f, :)
Therefore, a f,2 , . . . , a f,R , f = 1, . . . , F can be estimated using STSM algorithm since Ȳ(0)

f (r) = c f a f,r (a f,R • • • a f,r+1 a f,r-1 • • • (a f,1 + v f, f a f,1 + ∆a f,1 )) +         F p=1,p f v f,p Y p + ∆Y • 1 a f,1 A † 1 ( f, :)         (r)
.

Since Ȳ(0) f (1) has the following form Ȳ(0)

f (1) = c f (a f,1 + v f, f a f,1 + ∆a f,1 )(a f,R • • • a f,2 ) +         F p=1,p f v f,p Y p + ∆Y • 1 a f,1 A † 1 ( f, :)         (1)
, we estimate c f a f,1 by least squares once a f,2 , . . . , a f,R are estimated using STSM

c f a f,1 = min a Ȳ(0) f -a ⊗ â f,2 ⊗ • • • ⊗ â f,R = Ȳ(0) f (1) (â f,R • • • â f,2 ) T †
So, we put Ŷ(0)

f = c f a f,1 ⊗ â f,2 ⊗ • • • ⊗ â f,R and R f = Ȳ(0) f - Ŷ(0)
f . Therefore, the procedure to estimate Y f at iteration i = 0, . . . , K can be summarized in [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF], [START_REF] Lemma | Analysis of joint anglefrequency estimation using ESPRIT[END_REF] and [START_REF] Djermoune | Perturbation analysis of subspace-based methods in estimating a damped complex exponential[END_REF]. Note that this procedure is optimal because STSM and the least squares are optimal when they are used to estimate a f,2 , . . . , a f,R , f = 1, . . . , F and c f a f,1 , f = 1, . . . , F, respectively. Now we present the technique for improving the estimation of Y f . Let R (0)

F = R (1) 0 = Ỹ -F f =1 Ŷ(0) f and Ŷ(1) f = arg min X∈H Ŷ(0) f + R (1) f -1 -X (A.1)
where R (1) f = Ŷ(0) f + R (1) f -1 -Ŷ(1) f , f = 1, . . . , F, and Ŷ(i) f is an improved estimate of Y f . We follow the same procedure as described in equations ( 37), [START_REF] Lemma | Analysis of joint anglefrequency estimation using ESPRIT[END_REF] and [START_REF] Djermoune | Perturbation analysis of subspace-based methods in estimating a damped complex exponential[END_REF] to calculate Ŷ(1) f . We can state that there is improvement in the estimation of

Y f if Ỹ - F f =1 Ŷ(1) f ≤ Ỹ - F f =1 Ŷ(0) f (A.2)
We have R (1)

0 = Ȳ(0) 1 - Ŷ(0) 1 + F f =2 R f +V where V = Ỹ - Ȳ and Ȳ = F f =1
Ȳ(0) f . It can be verified that

R (1) f =         Ȳ(0) f + f -1 p=1 ( Ȳ(0) p - Ŷ(1) p ) + F p= f +1 R p + V         - Ŷ(1) f R (1) f -1 =         Ȳ(0) f + f -1 p=1 ( Ȳ(0) p - Ŷ(1) p ) + F p= f +1 R p + V         - Ŷ(0)
f However, from equation (A.1), Ŷ(1) f is the minimizer with respect to X ∈ H of

        Ȳ(0) f + f -1 p=1 ( Ȳ(0) p - Ŷ(1) p ) + F p= f +1 R p         -X Therefore, R (1)
f ≤ R (1) f -1 , f = 1, . . . , F. As consequence, R

F ≤ R (0) (1)
F , which we are seeking in expression (A.2). Similarly, we can prove that R (i) F ≤ R (i-1) F , i > 1, using the general forms of R (i) f and R (i)

f -1 R (i) f =         Ȳ(0) f + f -1 p=1 ( Ȳ(0) p - Ŷ(i) p ) + f -1 p= f +1 ( Ȳ(0) p - Ŷ(i-1)

Level µ or β µ or β µ or β 0 + 1 Figure 1 : 1 :

 0111 Figure 1: The multigrid dictionary refinement procedure with η = 1. ( ) atoms in the dictionary; ( • ) activated atoms; ( ) new atoms

Figure 2 :

 2 Figure 2: J (µ n ) in the single mode case with β n = 0

Figure 3 :

 3 Figure 3: Frequency RMSE using R-D sparse algorithm with different η ν . 2-D signal containing a single tone (Signal #1). (M 1 , M 2 ) = (10, 10). 1000 Monte-Carlo trials. (a) The initial harmonic dictionary contains N(0) = 10 atoms, (b) N(0) = 20 atoms.

Experiment 4 .Experiment 5 .

 45 3-D signal of size 10 × 10 × 10 containing three 3-D modes ( Signal #3). Note that there exists identical modes in two dimensions and frequencies in the first dimension are separated by less than 1/M 1 . The results are shown on Figure 7. In this experiment, the proposed R-D sparse approach performs better than TPUMA and Tensor-ESPRIT. Observe also that Tensor-ESPRIT outperforms TPUMA in this configuration (close frequencies and identical modes in dimensions 2 and 3). Results on Signal #4 of size 10 × 3 × 3 containing 4 modes are given in Figure

Figure 4 :Figure 5 :

 45 Figure 4: Damping factor RMSE using R-D sparse algorithm with different η α , β min and N(0). 2-D signal containing a single tone (Signal #1). (M 1 , M 2 ) = (10, 10). 1000 Monte-Carlo trials.

  Figure 9 shows the CPU time results of the proposed, Tensor-ESPRIT and TPUMA algorithms versus M 1 for a 3-D damped signal containing two modes with M 2 = M 3 = 4.

Figure 6 :Figure 7 :

 67 Figure 6: Frequency total root-mean square error for a 3-D signal containing two 3-D modes (Signal #2). (M 1 , M 2 , M 3 ) = (8, 8, 8). 1000 Monte-Carlo.

Figure 8 :

 8 Figure 8: Frequency total root-mean square error for a 3-D signal containing 4 modes with (M 1 , M 2 , M 3 ) = (10, 3, 3) (Signal #4). 200 Monte-Carlo.

Figure 9 :

 9 Figure 9: Average CPU time for a single run under M 2 = M 3 = 4 and F = 2.

1 A † 1 (

 11 get ∆S f = ∆Y • f, :) + F p=1 v f,p S p . Y f is estimated using Y f = S f • 1 a f,1, we differentiate and obtain

Table 1 :

 1 2-D and 3-D parameters of Signal #1 through #4

	Signal	ν f,1	α f,1	ν f,2	α f,2	ν f,3	α f,3	c f
	#1	0.22 -0.011	0.34 -0.015	-	-	1
	#2	0.40 -0.01	0.1	-0.01	0.1	-0.01	1
		0.20 -0.01	0.3	-0.15	0.25 -0.01	1
	#3	0.28 -0.01	0.31 -0.01	0.22 -0.01	1
		0.12 -0.01	0.45 -0.015	0.11 -0.01	1
		0.20 -0.01	0.31 -0.01	0.11 -0.01	1
	#4	0.30 -0.01	0.1	-0.01	0.1	-0.01	1
		0.13 -0.01	0.45 -0.015	0.4	-0.01	1
		0.20 -0.01	0.31 -0.01	0.1	-0.01	1
		0.42 -0.012	0.22 -0.01	0.32 -0.01	1

Table 2 :

 2 Different configurations for experiments 3 and 4 Fr ∃ identical modes ∃ identical modes following settings (L, η µ , η β ) =[START_REF] Gershman | Space-time processing for MIMO communications[END_REF][START_REF] Sahnoun | Sparse modal estimation of 2-D NMR signals[END_REF][START_REF] Haardt | Higher-order SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems[END_REF]. The number of iterations in Algorithm 3 is set to K = 2 because no improvement was observed for K > 2.

		F	Dim. 1	Dim. 2	Dim. 3
	Exp. 3	2 ∆ν > ∆ Fr	∆ν > ∆ Fr	∆ν > ∆ Fr
	Exp. 4	3 ∆ν < ∆

Note that this assumption is considered only in this section.

The case of M ≤ 2 in not of practical interest but the theorem is still valid by setting ζ M = 1 2 because J (µ n ) is a monotonically decreasing function in the interval [ν 1 , min{[START_REF] Stoica | Introduction to Spectral Analysis[END_REF] 2 , 1 2 + ν 1 }].

Acknowledgment

The authors would like to thank PhD Weize Sun for providing them with the TPUMA algorithm code.

$ This work is funded by the European Research Council under the European Community's Seventh Framework Programme FP7/2007-2013 Grant Agreement no. 320594, DECODA project.