Nicolas Forcadel

R Égis Monneau

EXISTENCE OF SOLUTIONS FOR A MODEL DESCRIBING THE DYNAMICS OF JUNCTIONS BETWEEN DISLOCATIONS

Keywords: Dislocation dynamics, non-local equations, junctions, parabolic system of equations AMS subject classifications. 35K15, 74K30

We study a dynamical version of a multi-phase field model of Koslowski and Ortiz for planar dislocation networks. We consider a two-dimensional vector field which describes phase transitions between constant phases. Each phase transition corresponds to a dislocation line, and the vectorial field description allows the formation of junctions between dislocations. This vector field is assumed to satisfy a non-local vectorial Hamilton-Jacobi equation with non-zero viscosity. For this model, we prove the existence for all time of a weak solution.

Introduction.

1.1. Physical motivation. Dislocations are line defects in crystal, and their motion is at the origin of plastic properties of metals. In these crystals, we can observe self-organised structures, like Frank networks, i.e. networks of dislocations related by junctions. See for instance page 190 in Hull, Bacon [START_REF] Hull | Introduction to dislocations[END_REF] for such networks in BCC iron, or page 188 for hexagonal networks in FCC crystals. In the present paper, we consider a special case of a network contained in a single slip plane, where the dislocations can move. We are interested in particular in the motion of the junctions between dislocations, which remains a quite open question, both from the modelling point of view, and from the mathematical analysis point of view (see for instance the work of Rodney, Le Bouar, Finel [START_REF] Rodney | Phase field methods and dislocations[END_REF]). Let us mention, for the stationary case, the work of Garroni [11]. The goal of the present paper is to propose and to study a model for the dynamics of junctions of dislocations.

The question of junctions has several other physical applications and there is various literature on this subject. Let us mention for instance the problem of crystal growth or grain growth (see Taylor [START_REF] Taylor | The motion of multiple-phase junctions under prescribed phase-boundary velocities[END_REF][START_REF]A variational approach to crystalline triple-junction motion[END_REF] and Bronsard, Reitich [START_REF] Bronsard | On three-phase boundary motion and the singular limit of a vector-valued Ginzburg-Landau equation[END_REF]). We also refer to Bonnet [START_REF] Bonnet | On the regularity of the edge set of Mumford-Shah minimizers[END_REF] for problems concerning the minimisation of the Mumford-Shah functional.

1.2. A phase field model for the dynamics of junctions. In a phase field model, the dislocation can be represented as the phase transition of a phase parameter ρ(x) = ρ 1 (x)e 1 +ρ 2 (x)e 2 ∈ R 2 defined for x = x 1 e 1 +x 2 e 2 in the plane R 2 with (e 1 , e 2 ) an orthonormal basis. The energy of the dislocations, in the presence of a constant exterior applied stress σ 0 ∈ R 2 , is then given (see the model of Koslowski and Ortiz [START_REF] Koslowski | A multi-phase field model of planar dislocation networks[END_REF]) by

E(ρ) = R 2 - 1 2 (C 0 ⋆ ρ) • ρ -σ 0 • ρ + W (ρ), (1.1)
where the precise meaning of this expression will be given later.

For any phase transition between two states A and B, the difference B -A needs physically to be the Burgers vector of the dislocation, i.e. a vector of the lattice Λ = Za 1 + Za 2 of the crystal we are considering, with general basis (a 1 , a 2 ). This information is encoded in the potential W : R 2 -→ R + which is is assumed to be minimal on Λ and to have the periodicity of the lattice Λ:

W (ρ + a) = W (ρ) for any a ∈ Λ. When the material is submitted to an exterior shear stress, it makes the dislocations move. The dynamics of a given dislocation line is physically given by its normal velocity, which is called the resolved Peach-Koehler force. This force is the sum of the resolved exterior shear stress and the stress created by all the dislocations lines, including the line itself. In the expression giving the energy (1.1), the kernel C 0 (x) is a 2 × 2 symmetric matrix which takes into account the long range elastic interactions between dislocations and

b 3 ρ ∼ 0 ρ ∼ b 1 b 2 b 1 ρ ∼ b 1 + b 2 = -b 3
(C 0 ⋆ ρ) i = j=1,2 C 0 ij ⋆ ρ j , for i = 1, 2
where ⋆ denotes the usual convolution. In (1.1) and throughout the paper, we denote by A • B the scalar product between two vectors A, B ∈ R 2 . The resolved stress σ[ρ] created by all the dislocations is then formally given by the opposite of the gradient of the energy -E ′ (ρ), and can be expressed as the following non-local quantity

σ[ρ] = σ 0 + C 0 ⋆ ρ -W ′ ρ (ρ). (1.3)
The phase parameter ρ(t, x) ∈ R 2 is then assumed to satisfy the following equation

           (ρ k ) t = |∇ρ| -1 i=1,2 j=1,2 (σ[ρ]) i ∇ j ρ i ∇ j ρ k + ε∆ρ k for k = 1, 2 in (0, T ) × T 2 , ρ(0, x) = ρ 0 (x) on T 2 , (1.4)
where σ is given in (1.3), ρ t = ∂ρ ∂t and ∇ j ρ i = ∂ρ i ∂x j for i, j = 1, 2, and

|∇ρ| 2 = i=1,2 j=1,2 |∇ j ρ i | 2 .
Here the parameter 1 > ε > 0 is a small viscosity introduced in the model, in order to regularise the equation, but which has no real physical meaning. Given any time T > 0, we will study this equation, not on the the whole plane, but on the particular torus T 2 = R 2 \Z 2 , in order to simplify the analysis. This kind of periodic conditions is also meaningful physically, if we want to describe periodic networks of dislocations. This means in particular that σ[ρ] is given by (1.3) where the convolution is done on the torus.

Finally, let us mention that our model (1.4) has some similarities with the model of Allen, Cahn [START_REF] Allen | A microscopic theory for the antiphase boundary motion and its application to antiphase domain coarsening[END_REF] on the motion of curved boundaries in which they consider gradient flow associated with a free-energy functional. This led to the study of scalar Ginzburg-Landau type diffusion equation like

u t = ∆u -W ′ (u).
1.3. Main result. We make the following assumption on the kernel C 0 : T 2 -→ R 2×2 sym (A) We assume that there exists a constant m > 0, such that for any k ∈ Z 2 , the Fourier coefficients of the kernel

C 0 (k) = T 2 dx e -2iπk•x C 0 (x) satisfy C 0 (k) = M (k), where for any ξ = (ξ 1 , ξ 2 ) ∈ R 2 and any p = (p 1 , p 2 ) ∈ R 2              M ∈ C ∞ (R 2 \ {0} ; R 2×2 sym ), M (-ξ) = M (ξ), M (ξ) = |ξ|M ξ |ξ| |ξ| |p| 2 m ≥ - i=1,2 j=1,2 p i • M ij (ξ) • p j ≥ m|ξ| |p| 2 with |p| 2 = i=1,2 (p i ) 2 (1.5)
We also make the following assumption on the potential W :

T 2 -→ R + (B)
We assume that W ∈ C 2 and W satisfies (1.2).

Then we have the following result for the model of dynamics of junctions between dislocations: Theorem 1.1 (Existence of a solution) Under assumptions (A) and (B), if ρ 0 ∈ H 1 (T 2 )

2 , then for any constant applied stress σ 0 ∈ R 2 and for any time T > 0, there exists a solution ρ of (1.4) with

ρ ∈ C 0 [0, T ); L 4 3 (T 2 ) 2 .
The uniqueness of the solution is not known, neither the existence of a solution when ε = 0. Let us mention that equation (1.4) is a non-local system of scalar equations, and can be sketched as the following equation

v t = |∇v| 2 + ∆v (1.6)
Indeed, this comes from our assumption (A) that the convolution with the kernel behaves like a first order operator. A lot of work has been done on equations (or systems) like equation (1.6). Let us mention for instance the works of Boccardo, Murat, Puel [START_REF] Boccardo | Existence de solutions faibles pour des équations elliptiques quasi-linéaires à croissance quadratique[END_REF][START_REF]Résultats d'existence pour certains problèmes elliptiques quasilinéaires[END_REF][START_REF]Existence results for some quasilinear parabolic equations[END_REF] in which they study general equations including equation (1.6) and prove existence result. Equation (1.6) is also similar to the Navier-Stokes equations written for the potential A such that the velocity of the fluid is given by u = curl A (see for instance Leray [16]).

1.4. Organisation of the paper. In Section 2, we give some remarks on the modelling. In Section 3, we study an approximate problem of equation (1.4) where the right hand side is approached by some term at most linear in the solution. The main result is proved in Section 4. In a first subsection, we give some a priori estimates for the solution of the approximate problem and then in a second subsection we pass to the limit in the approximate problem.

Notation.

In what follows, we will denote by C a generic constant, which will then satisfy C + C = C, C • C = C, and so on. We also use the following set:

W 2,1;p (Q T ) = u ∈ L p (Q T ); u t ∈ L p (Q T ) and ∂ 2 u ∂x i , ∂x j ∈ L p (Q T ) for i, j = 1, 2
where

Q T = (0, T ) × T 2 .
2. Preliminary remarks on the modelling.

2.1. Dynamics of a single dislocation. In this section, we consider a special case where a dislocation of Burgers vector b 1 ∈ Λ is described by the phase transition of a scalar parameter ρ such that

ρ = ρ b 1
Then the resolved shear stress that makes the dislocation move, is given by σ

[ρ] = b 1 • σ[ρ b 1 ]
and the dislocation line (described by the phase parameter ρ) moves with normal velocity proportional to this shear stress. More precisely, ρ solves the following equation

ρt = |b 1 | -1 σ[ρ]|∇ρ| + ε∆ρ where ρt = ∂ ρ ∂t and ∇ρ = ∂ ρ ∂x 1 e 1 + ∂ ρ ∂x 1 e 2
. Physically, we expect to have the eikonal equation with ε = 0, but in order to simplify the analysis and get an existence result, we only consider the case with additional viscosity ε > 0. Then, we easily check that ρ = ρ b 1 satisfies (1.4).

Explicit expression of C

0 for isotropic materials. Given a particular Burgers vector b 1 , let us consider c0 = (b 1 ) T • C 0 • b 1 .
In the special case of isotropic linear elasticity with constant coefficients, we recall (see for instance a limit case of the Peierls-Nabarro model in Alvarez et al. [START_REF] Alvarez | Dislocation dynamics: short time existence and uniqueness of the solution[END_REF]) that we have for k

= (k 1 , k 2 ) c0 (k) = - µ|b 1 | 2 2 1 1-ν (k • e) 2 + k ⊥ • e 2 |k| with e = b 1 |b 1 |
where k ⊥ = (-k 2 , k 1 ) is the vector obtained by a rotation of k of angle π/2. Here µ > 0 is a Lamé coefficient and ν ∈ (-1, 1/2) is the Poisson ratio. We deduce that

C 0 (k) = - µ 2|k| 1 1 -ν k ⊗ k + k ⊥ ⊗ k ⊥
which satisfies assumption (A).

3. An approximate problem. We first start to approximate the right hand side of equation (1.4) by some term at most linear in the solution. To this end, we introduce a function h n defined by

h n (r) = h 0 (r -n) with h 0 (r) =            1 if r ≤ 0 1 -r if 0 ≤ r ≤ 1 0 if r ≥ 1
We then look at the following approximate problem:

   ρ t -ǫ∆ρ = f n [ρ] on Q T := (0, T ) × T 2 ρ(0, •) = ρ 0 on T 2 (3.1)
where

f n [ρ] = h n (|∇ρ|) |∇ρ| -1 (∇ρ) T • ∇ρ • σ[ρ]
and σ[ρ] is given in (1.3), and is at most linear in ρ.

The natural idea to find a solution to equation (3.1), is to define the map Φ which associates to any function u, the solution

ρ = Φ(u) of    ρ t -ǫ∆ρ = f n [u] on Q T := (0, T ) × T 2 ρ(0, •) = ρ 0 on T 2 (3.2)
and to prove that Φ has a fixed point in a suitable space. This way, we will prove the following result Theorem 3.1 (Existence of a solution for the approximate problem)

If ρ 0 ∈ H 1 (T 2 )
2 , then for any n ≥ 1 and any T > 0, there exists a solution ρ n of (3.1) with

ρ n ∈ L 2 ((0, T ); H 2 (T 2 ) 2 ) ∩ C 0 ([0, T ); L 2 (T 2 ) 2 ).
In this section, we will make the proof of this theorem. In a first subsection, we will collect some preliminary results, and in a second subsection we will prove that Φ has a fixed point.

Preliminary results.

The following estimate on the stress will be important.

Lemma 3.2 (Estimate on σ[ρ])

For any p ∈ (1, +∞), there exists a constant C (depending on p, on the constant σ 0 , on the potential W and on the constant m defined in assumption (A)), such that for any ρ ∈ W 1,p (T 2 )

2 , we have

σ[ρ] (L p (T 2 )) 2 ≤ C 1 + ∇ρ (L p (T 2 )) 2×2 . (3.3) Partial proof of Lemma 3.2
Let us make the proof for p = 2. We have with σ

= σ[ρ] |σ -σ 0 + W ′ (ρ)| 2 (L 2 (T 2 )) 2 = k∈Z 2 | (C 0 ⋆ ρ)(k)| 2 = k∈Z 2 | C 0 (k) • ρ(k)| 2 ≤ 1 m 2 k∈Z 2 |k| 2 | ρ(k)| 2 ≤ 1 (2πm) 2 k∈Z 2 | ∇ρ(k)| 2 = 1 (2πm) 2 |∇ρ| 2 (L 2 (T 2 )) 2×2 Therefore |σ[ρ]| (L 2 (T 2 )) 2 ≤ |σ 0 | + |W ′ | (L ∞ (T 2 )) 2 + 1 2πm |∇ρ| (L 2 (T 2 )) 2×2
which provides the result in the case p = 2.

The proof for the general case p ∈ (1, +∞) is given in Appendix. 2

We will also need the following result.

Lemma 3.3 (Estimate on f n [u]) If u ∈ (H 1 (T 2 )) 2 , then f n [u] ∈ (L 2 (T 2
)) 2 with the following estimate:

f n [u] (L 2 (T 2 )) 2 ≤ C(n + 1) 1 + ∇u L 2 (T 2 ) 2×2
where the constant C depends on σ 0 , on the potential W and on the constant m defined in assumption (A).

Proof of Lemma 3.3 Since supp(h n ) ⊂ [0, n + 1], the following holds |f n [u]| ≤ (n + 1)|σ[u]|, (3.4)
where we have used the fact that

|B T • B • p| ≤ |B| 2 |p| for B ∈ R 2×2 and p ∈ R 2 . Then f n [u] (L 2 (T 2 )) 2 ≤(n + 1) σ[u] (L 2 (T 2 )) 2 ≤C(n + 1) 1 + ∇u L 2 (T 2 ) 2×2
where we have used Lemma 3.2 2

We now recall some classical results. We start with the following parabolic estimates for the following equation

   g t -ǫ∆g = f on Q T := (0, T ) × T 2 g(0, •) = g 0 on T 2 (3.5)
Proposition 3.4 (Parabolic estimates for the heat equation)

Let g 0 ∈ H 1 (T 2 ) and f ∈ L 2 (Q T ).
Then there exists a unique solution g to (3.5) with

g ∈ L 2 ((0, T ); H 2 (T 2 )) ∩ L ∞ ((0, T ); H 1 (T 2 )), g t ∈ L 2 (Q T ). (3.6)
We have the following estimate

sup 0≤t≤T g(t) H 1 (T 2 ) + g L 2 ((0,T );H 2 (T 2 )) + g t L 2 ((0,T );L 2 (T 2 )) ≤C T f L 2 (QT ) + g 0 H 1 (T 2 ) , (3.7)
where the constant C T only depends on T and ε.

Moreover we have To prove (3.8), we simply multiply equation (3.5) by g and integrate over (0, t) in time, taking the supremum for 0 ≤ t ≤ T . We get

sup 0≤t≤T T 2 g 2 (t) + 4ε T 0 T 2 |∇g| 2 ≤ 2 T 2 (g 0 ) 2 + 4T T 0 T 2 f 2 (3.
sup 0≤t≤T T 2 g 2 (t) 2 + ε T 0 T 2 |∇g| 2 ≤ T 2 (g 0 ) 2 2 + T 0 T 2 |g f |
We now use the fact that

T 0 T 2 |g f | ≤ T 0 T 2 g 2 1 2 • T 0 T 2 f 2 1 2 ≤ T sup 0≤t≤T T 2 g 2 (t) 1 2 • T 0 T 2 f 2 1 2 ≤ 1 4 sup 0≤t≤T T 2 g 2 (t) + T T 0 T 2 f 2
which implies the result. 2

We also recall the Theorem 3.5 (Schaefer's fixed point theorem) Let X be a real Banach space. Suppose that Φ : X → X is a continuous and compact mapping. Assume further that the set {u ∈ X, u = λΦ(u) for some λ ∈ [0, 1]} is bounded. Then Φ has a fixed point.

For the proof of this theorem, we refer to Evans [10, Theorem 4 page 504].

Finally, we will need some compactness argument and a weak continuity property contained in the following two classical results:

Proposition 3.6 (Compactness) We recall that W 2,1;2 (Q T ) = g ∈ L 2 ((0, T ); H 2 (T 2 )), g t ∈ L 2 (Q T ) .
Then the injection

W 2,1;2 (Q T ) -→ L 2 ((0, T ); H 1 (T 2 )) is compact.
For the proof of this result, we refer to Lions [17, Theorem 5.1 page 58].

Proposition 3.7 (Continuity)

With the notation of Proposition 3.6, let us consider a sequence (g m ) m such that

g m ⇀ g weakly in W 2,1;2 (Q T )
We assume also that g m |t=0 = ρ 0 . Then

g |t=0 = ρ 0 .
This result is classical but for the reader's convenience we give the proof in Appendix.

3.2. Proof of Theorem 3.1. We are now ready to make the proof of Theorem 3.1. To this end, for any T > 0, we set

X T = L 2 ((0, T ); H 1 (T 2 )).
In all what follows the index n is assumed fixed. We first remark that if u ∈ X 2 T , then 2 , and then we can consider the solution

f n [u] ∈ (L 2 (Q T ))
ρ of    ρ t -ǫ∆ρ = f n [u] on Q T := (0, T ) × T 2 ρ(0, •) = ρ 0 on T 2 (3.9)
which satisfies ρ ∈ X 2 T because of the parabolic estimates Proposition 3.4. Then we set Φ(u) = ρ, and see that Φ maps X 2 T into X 2 T . We will prove that Φ admits a fixed point using Schaefer's fixed point theorem. We do the proof in four steps.

Step 1 : weak continuity of Φ Let us consider sequences (u m ) m , (ρ m ) m such that

   u m ∈ X 2 T , ρ m = Φ(u m ) u m -→ u in X 2 T From Lemma 3.3, we deduce that ||f n [u m ]|| (L 2 (QT )) 2 ≤ C(n + 1) 1 + ||u m || X 2 T (3.10)
From the parabolic estimates (Proposition 3.4), we deduce that ρ m is bounded in W 2,1;2 (Q T ) 2 , i.e. there exists a constant C > 0 such that

||ρ m || (W 2,1;2 (QT )) 2 ≤ C (3.11)
Therefore, up to a subsequence, we have

ρ m ⇀ ρ in W 2,1;2 (Q T ) 2
and from Proposition 3.7, we deduce that

ρ |t=0 = ρ 0 on T 2
We now claim that

f n [u m ] -→ f n [u] in L 1 (Q T ) (3.12)
Indeed, we can write

f n [u] = g n (∇u) • σ[u] with g n (∇u) := h n (|∇u|) |∇u| -1 (∇u) T • ∇u
From the proof of Lemma 3.2, for p = 2, we already deduce that

σ[u m ] -→ σ[u] in L 2 (Q T ) (3.13)
From the convergence of u m to u in X 2 T , we deduce that up to a subsequence we have ∇u m -→ ∇u a.e. in Q T . Now from the fact that g n is continuous and bounded, we deduce in particular that

g n (∇u m ) -→ g n (∇u) in L 2 (Q T ) (3.14)
Then the convergence (3.12) follows from (3.13) and (3.14). Therefore we conclude that ρ solves (3.9). Finally, by uniqueness of the solutions of (3.9), we deduce that the limit ρ does not depend on the choice of the subsequence, and then that the full sequence converges:

ρ m ⇀ ρ weakly in W 2,1;2 (Q T ) 2 , with ρ = Φ(u)
Step 2 : compactness of Φ The compactness (and the usual continuity) of Φ follows from the compactness of the injection W 2,1;2 (Q T ) 2 -→ X 2 T (see Proposition 3.6).

Step 3 : a priori bounds on the solutions of u = λΦ(u) for T small Let us consider a solution u of u = λΦ(u) for some λ ∈ [0, 1] (3.15)

Then from the parabolic estimates (3.8), we have

sup 0≤t≤T T 2 |u(t)| 2 + 4ε T 0 T 2 |∇u| 2 ≤ 2 T 2 |ρ 0 | 2 + 4T T 0 T 2 |λf n [u]| 2 ≤ 2 T 2 |ρ 0 | 2 + 8T C 2 (n + 1) 2 T + T 0 T 2 |∇u| 2
where in the third line we have used Lemma 3.3 and the fact that |λ| ≤ 1. Therefore for

T ≤ T * := 4C 2 (n + 1) 2 -1 ε (3.16) we have sup 0≤t≤T T 2 |u(t)| 2 + 2ε T 0 T 2 |∇u| 2 ≤ 2 T 2 |ρ 0 | 2 + 2εT
which proves that there exists a constant C > 0 such that any solution of (3.15) satisfies ||u|| X 2 T ≤ C We can then apply the Schaefer's fixed point Theorem (Theorem 3.5), to deduce that Φ has a fixed point on X 2 T , and therefore there is a solution ρ of (3.1) on the time interval (0, T ) if T satisfies (3.16), i.e. if T is small enough independently on the initial data ρ 0 .

Step 4 : solution for any time Let us call ρ ρ 0 , t the function ρ(t, •) obtained at Step 3 as a solution of (3.1) on the time interval [0, T * ) with initial data ρ 0 . From the parabolic estimates (Proposition 3.4), we also know that ρ(t, •) ∈ (H 1 (T 2 )) 2 for any t ∈ [0, T * ). Then we can define with τ = T * /2

u(0) = ρ 0 and u(t) = ρ (u(kτ ), t) if kτ ≤ t < (k + 1)τ with k ∈ N.
Using the fact that u t ∈ L 2 loc ((0, +∞); (L 2 (T 2 ) 2 ), and the fact that the problem is invariant by translation in time, we can easily check that u solves (3.1) for any T > 0 and provides the desired solution ρ n = u of Theorem 3.1. This ends the proof of Theorem 3.1.

4.

A priori estimates and proof of Theorem 1.1.

A priori estimates.

We have the following a priori estimates:

Lemma 4.1 (A priori estimates)

There exists a constant C > 0 such that for all T > 0, n ≥ 1 and 0 < ε < 1, any solution ρ n of (3.1) given by Theorem 3.1, satisfies

ρ n 2 L ∞ ((0,T );(H 1 2 (T 2 )) 2 ) ≤ Ce CT ε , (4.1)
ρ n 2 L 2 ((0,T );(H 3 2 (T)) 2 ) ≤ C ε e CT ε , (4.2)
and

h n (|∇ρ n |)|∇ρ n | -1 2 ∇ρ n • σ[ρ n ] 2 (L 2 (QT )) 2 ≤ Ce CT ε . (4.3)
Proof of Lemma 4.1

Step 1: Preliminaries on the energy We first recall the expression of the energy for a general Z 2 -periodic smooth function

ρ(x) = (ρ 1 (x), ρ 2 (x)) E(ρ) = T 2 - 1 2 (C 0 ⋆ ρ) • ρ -σ 0 • ρ + W (ρ)
For future use, we start to evaluate from below the first term in the energy, using Fourier series

T 2 -(C 0 ⋆ ρ) • ρ = Re k∈Z 2 -(C 0 ⋆ ρ)(k) • ρ * (k) = Re k∈Z 2 -( C 0 (k) • ρ(k)) • ρ * (k) ≥ m k∈Z 2 |k| | ρ(k)| 2
where we have used in the first line the fact that ρ and C 0 are real, and in the last line we have used assumption (A). Then we define

||ρ|| 2 ( Ḣ 1 2 (T 2 )) 2 := k∈Z 2 |k| | ρ(k)| 2 .
Similarly, we compute

T 2 -(C 0 ⋆ (∇ρ) T ) : ∇ρ =(2π) 2 Re k∈Z 2 -( C 0 (k) • ρ(k) ⊗ (ik)) : (ik) * ⊗ ρ * (k) =(2π) 2 Re k∈Z 2 -|k| 2 ( C 0 (k) • ρ(k)) : ρ * (k) ≥(2π) 2 m k∈Z 2 |k| 3 | ρ(k)| 2
where we have used assumption (A) in the last line. Then we define

||ρ|| 2 ( Ḣ 3 2 (T 2 )) 2 := k∈Z 2 |k| 3 | ρ(k)| 2 .
Step 2: Estimate on the time-derivative of the energy Let us fix T > 0. We know that any solution ρ n given by Theorem 3.1 belongs to the space W 2,1;2 (Q T ). In particular, using the following general fact (because of assumption (A))

T 2 -(C 0 ⋆ ρ) • ρ = Re k∈Z 2 -|k|( Ĉ0 k |k| • ρ(k)) • ρ * (k)
we deduce that the energy E(ρ n (t)) is well-defined for almost every t ∈ [0, T ), and that for almost every time t ∈ [0, T ), we have

d dt E(ρ n (t)) = T 2 -σ[ρ n ] • ρ n t = T 2 -h n (|∇ρ n |)|∇ρ n | -1 |∇ρ n • σ[ρ n ]| 2 -εσ[ρ n ] • ∆ρ n = T 2 -h n (|∇ρ n |)|∇ρ n | -1 |∇ρ n • σ[ρ n ]| 2 - T 2 ε W ′′ (ρ n ) : ((∇ρ n ) T • ∇ρ n ) -(C 0 ⋆ (∇ρ n ) T ) : ∇ρ n . Therefore d dt E(ρ n (t)) + T 2 h n (|∇ρ n |)|∇ρ n | -1 |∇ρ n • σ[ρ n ]| 2 ≤Cε T 2 |∇ρ n | 2 + C 0 ⋆ (∇ρ n ) T ) : ∇ρ n (4.4)
But now (up to change the constant line to line)

∇ρ n 2 (L 2 (T 2 )) 2×2 ≤ C k∈Z 2 |k| 2 | ρ n (k)| 2 ≤ C k∈Z 2 |k| 3 2 | ρ n (k)| • |k| 1 2 | ρ n (k)| ≤ C k∈Z 2 1 2α |k| 3 | ρ n (k)| 2 + α 2 |k|| ρ n (k)| 2 ≤ C T 2 - 1 α (C 0 ⋆ (∇ρ n ) T ) : ∇ρ n + T 2 -α(C 0 ⋆ ρ n ) • ρ n ,
where ρ n (k) are the Fourier coefficients of ρ n and α is a constant which will be precised later. We then deduce finally that:

d dt E(ρ n (t)) + T 2 h n (|∇ρ n |)|∇ρ n | -1 |∇ρ n • σ[ρ n ]| 2 ≤ -Cε 1 - 1 α T 2 -(C 0 ⋆ (∇ρ n ) T ) : ∇ρ n + Cεα T 2 -(C 0 ⋆ ρ n ) • ρ n ≤ -Cε ρ n (t) 2 ( Ḣ 3 2 (T 2 )) 2 + Cε 1 + E(ρ n (t)) + |σ 0 | T 2 ρ n (t) . (4.5)
for α chosen large enough, with C a suitable positive constant.

Step 3: Estimate on the time-derivative of the mean-value of the solution Integrating in space equation (3.1), we get

d dt T 2 ρ n (t) = T 2 h n (|∇ρ n |)|∇ρ n | -1 (∇ρ n ) T • ∇ρ n • σ[ρ n ]
and then

d dt T 2 ρ n (t) ≤ T 2 (h n (|∇ρ n |)|∇ρ n |) 1 2 • (h n (|∇ρ n |)|∇ρ n | -1 ) 1 2 |∇ρ n • σ[ρ n ]| ≤ T 2 (1 + |σ 0 |)h n (|∇ρ n |)|∇ρ n | + 1 4(1+|σ 0 |) h n (|∇ρ n |)|∇ρ n | -1 |∇ρ n • σ[ρ n ]| 2 (4.6)
Step 4: Estimate on the energy Setting

F n (t) = 1 + E(ρ n (t)) + (1 + |σ 0 |) T 2 ρ n (t) + (1 + |σ 0 |) 4 , (4.7)
we deduce from (4.5) and (4.6) that

d dt F n (t) + 3 4 h n (|∇ρ n |)|∇ρ n | -1 |∇ρ n • σ[ρ n ]| 2 ≤ -Cε ρ n (t) 2 ( Ḣ 3 2 (T 2 )) 2 + Cε 1 + E(ρ n (t)) + |σ 0 | T 2 ρ n (t) + (1 + |σ 0 |) 2 T 2 h n (|∇ρ n |)|∇ρ n |.
Now we remark that

(1 + |σ 0 |) 2 T 2 h n (|∇ρ n |)|∇ρ n | ≤(1 + |σ 0 |) 2 T 2 |∇ρ n | ≤ Cε 2 T 2 |∇ρ n | 2 + 1 + |σ 0 | 4 2Cε
Using the fact that (since the domain is bounded)

T 2 |∇ρ n | 2 ≤ ρ n (t) 2 ( Ḣ 3 2 (T 2 )) 2 ,
we get

d dt F n (t) + 3 4 T 2 h n (|∇ρ n |)|∇ρ n | -1 |∇ρ n • σ[ρ n ]| 2 + Cε 2 ρ n (t) 2 ( Ḣ 3 2 (T 2 )) 2 (4.8) ≤ Cε 1 + E(ρ n (t)) + |σ 0 | T 2 ρ n (t) + (1 + |σ 0 |) 4 2Cε ≤ C ε F n (t).
This implies, using Gronwall Lemma,

F n (t) ≤ F n (0)e C ε t . (4.9)
Step 5: Estimate on ρ n Let us first remark that

E(ρ n (t)) ≥ - 1 2 (C 0 ⋆ ρ n ) • ρ n -|σ 0 | T 2
ρ n (t) . (4.10)

Using (4.9), (4.10) and the definition of F n (t), yields

T 2 - 1 2 (C 0 ⋆ ρ n ) • ρ n + T 2 ρ n (t) ≤ Ce C ε t .
Using Step 1, we then get

ρ n 2 L ∞ ((0,T );( Ḣ 1 2 (T 2 )) 2 ) ≤ Ce C ε T and T 2 ρ n (t) ≤ Ce C ε T .
This implies (4.1). Taking the integral T 0 in (4.8) and using the fact that ∀t ≤ T , F n (t) ≥ 0, we get

h n (|∇ρ n |)|∇ρ n | -1 2 |∇ρ n • σ[ρ n ]| 2 (L 2 (QT )) 2 + εC ρ n 2 L 2 ((0,T );( Ḣ 3 2 (T 2 ))) ≤ Ce C ε T
which implies (4.2) and (4.3). 2 4.2. Proof of Theorem 1.1. We are now able to prove Theorem 1.1. In this section, we denote by C a generic constant which can depend on ρ 0 , ε and T but which do not depend on n.

Proof of Theorem 1.1 Let T > 0. The idea of the proof is to pass to the limit in Equation (3.1). The only difficulty is to prove that the non-linear term f n [ρ n ] converges in a certain sense to |∇ρ| -1 (∇ρ) T • ∇ρ • σ[ρ], where ρ is the limit of ρ n in an appropriate norm. The proof is decomposed into five steps:

Step 1: a priori bound on f n [ρ n ] We have the following estimate on f n [ρ n ]:

f n [ρ n ] " L 4 3 (QT ) " 2 ≤ C. (4.11)
To prove this, let us write

f n [ρ n ] = |∇ρ n | -1 (∇ρ n ) T • |∇ρ n | 1 2 h n (|∇ρ n |)|∇ρ n | -1 2 ∇ρ n • σ[ρ n ] .
Using (4.3), we have that the last term is bounded in (L 2 (Q T )) 2 by C. Moreover, the first term is bounded by 1 in (L ∞ (Q T )) 2×2 , then we just have to bound the term

|∇ρ n | 1 2 in L 4 (Q T ). Using (4.
2), we have

|∇ρ n | 1 2 L 4 (QT ) = QT |∇ρ n | 2 1 4 ≤ ρ n 1 2 L 2 ((0,T );(H 3 2 (T)) 2 ) ≤ C. (4.12)
This implies (4.11).

Step 2: Strong convergence of ∇ρ n in L 2 ((0, T ); (L where for a Banach space B W 1 2 ;p ((0, T ); B) = g ∈ L p ((0, T ); B),

T 0 T 0 g(t) -g(s) p B |t -s| 1 2 p+1
dt ds < ∞ is equiped with the following norm:

g W 1 2
;p ((0,T );B)

:= T 0 T 0 g(t) -g(s) p B |t -s| 1 2 p+1
dt ds

1 p .
Moreover, using (4.2) we get

∇ρ n L 2 " (0,T ); " H 1 2 (T 2 ) " 2×2 « ≤ C. (4.14)
We then use the following lemma:

Lemma 4.2 (Compactness result) Let (g n ) n be a sequence uniformly bounded in

L 2 (0, T ); H 1 2 (T 2 ) ∩ W 1 2 ; 4 3 (0, T ); L 4 3 (T 2 ) ,
then, for a subsequence,

g n → g strongly in L 2 ((0, T ); L 4 3 (T 2 )).
Formally, the proof uses the fact that H Step

3: Weak convergence of σ[ρ n ] in L 2 (0, T ); L 4 (T 2 ) 2 We have H 1 2 (T 2 ) ⊂ L 4 (T 2 ) with continuous injection (see Adams [1, Theorem 7.57 p217]). So L 2 ((0, T ); H 1 2 (T 2 )) ⊂ L 2 ((0, T ); L 4 (T 2
)) with continuous injection. We then deduce from (4.2) that

∇ρ n L 2 ((0,T );(L 4 (T 2 )) 2×2 ) ≤ C. (4.15)
Using Lemma 3.2, we then get

C 0 ⋆ ρ n L 2 ((0,T );(L 4 (T 2 )) 2 ) ≤ C. (4.16)
Using the fact that the application W

2,1, 4 3 x,t (Q T ) → L 4
3 (Q T ) is compact and the converse of Lebesgue Theorem, we deduce that W ′ (ρ n ) → W ′ (ρ) almost everywhere. This implies that σ

[ρ n ] ⇀ σ[ρ] in L 2 ((0, T ); (L 4 (T 2 )) 2 .
Step 4: Passing to the limit Using Step 2 and Step 3 and the fact that |∇ρ n | -1 ∇ρ n is bounded by 1, we deduce that

f n [ρ n ] → |∇ρ| -1 (∇ρ) T • ∇ρ • σ[ρ] in the distributions sense.
By passing to the limit in (3.1), we obtain

ρ t -ǫ∆ρ = |∇ρ| -1 (∇ρ) T • ∇ρ • σ[ρ] in D ′ ((0, T ) × T 2 ).
(4.17)

Step 5: Initial Condition Using the fact that ρ n t are bounded uniformly in L 4 3 (Q T ) (by parabolic estimates for the heat equation and Step 1), we deduce that (uniformly in n)

ρ n (t + h) -ρ n (t) (L 4 3 (T 2 )) 2 ≤ Ch 1 4 ρ n t L 4 3 ((0,T );(L 4 3 (T 2 )) 2 )
and then ρ ∈ C 0 ((0, T ); (L 4 3 (T 2 )) 2 ) and ρ |t=0 = ρ 0 . This achieves the proof of Theorem 1.1. 2 5. Appendix. Full proof of Lemma 3.2 Here we do the proof for any p ∈ (1, +∞). Under assumption (A), there exists a constant C > 0 only depending on p, such that the following result holds for any

ρ ∈ W 1,p (R 2 ) | C0 ⋆ R 2 ρ| (L p (R 2 )) 2 ≤ C m |∇ρ| (L p (R 2 )) 2×2
where the Fourier transform of C0 satisfies C0 = M with M as in (1.5). This result can be found in the scalar case on R n in Stein [START_REF] Stein | Harmonic analysis: real-variable methods, orthogonality, and oscillatory integrals[END_REF], see proposition 5 page 251, or Coifman, Meyer [START_REF] Coifman | Au delà des opérateurs pseudo-différentiels[END_REF], Theorem 9 page 39 and Proposition 2 page 41. See also Calderon-Zygmund inequalities Theorem 2.7.2 in Morrey [START_REF] Jr | Multiple integrals in the calculus of variations[END_REF]. Here the convolution by C0 is a multiplier operator in the class S 1 of pseudo-differential operators. We then get the result in the vectorial case, summing the scalar components. See also the book of Garroni, Menaldi [START_REF] Garroni | Green functions for second order parabolic integrodifferential problems[END_REF] for complements on integro-differential operators. The fact that the result holds on the torus T 2 is then classical. We prove it for the convenience of the reader. To this end, we consider a smooth function ϕ such that

| C0 ⋆ R 2 ρ| (L p (T 2 )) 2 ≤ C m |B 2 | |B 1 | |∇ρ| (L p (T 2 )) 2×2
i.e.

|C 0 ⋆ T 2 ρ| (L p (T 2 )) 2 ≤ 4C m |∇ρ| (L p (T 2 )) 2×2 with C 0 (x) = k∈Z 2 C0 (x -k).
We then get the final result by density of smooth functions in (W 1,p (T 2 )) 2 . 2

Proof of Proposition 3.7

For simplicity of notation, we denote by g(t) the function x → g(t, x). We have

g m (t) -ρ 0 (L 2 (T 2 )) 2 ≤ t 0 ds g m t (s) (L 2 (T 2 )) 2 ≤ √ t g m t (L 2 (QT )) 2 .
Using the fact that g m is bounded uniformly in W 2,1;2 (Q T ) (this is a consequence of the fact that g m ⇀ g in W 2,1;2 (Q T )), we get Using that ϕ ≥ 0 is arbitrary, we deduce that for almost every t, we have

g m (t) -ρ 0 (L 2 (T 2 )) 2 ≤ C √ t. ( 5
g(t) -ρ 0 2 (L 2 (T 2 )) 2 ≤ √ Ct.
This implies the result. 2

(1. 2 )

 2 In this model, junctions of three dislocations of Burgers vectors b 1 , b 2 , b 3 ∈ Λ with b 1 + b 2 + b 3 = 0 is expected, like for instance as the phase transitions between the states 0, b 1 , -b 3 (see Figure 1.1).

Fig. 1 . 1 .

 11 Fig. 1.1. The junction of three dislocations as phase transitions of ρ.

8 )

 8 Proof of Proposition 3.4 For the proof of (3.6)-(3.7), we refer to Evans [10, Theorem 5 page 360].

4 3 (

 3 T 2 )) 2×2 ) Using the parabolic estimates for the heat equation (see [15, ch 4.3 p 80 and ch 4.9 p 341]) and Step 1, we get

1 2 ⊂ L 4 3

 24 with compact injection in space while the compactness in time comes from(4.13). We refer to Simon [20, Corollary 5, p.86] for a more general result and for the proof of this lemma.Using (4.13), (4.14) and Lemma 4.2, we then deduce that, for a subsequence, ∇ρ n → ∇ρ strongly in L 2 (0, T );

1 |

 1 ϕ(x) = 1 on [-1/3, 1/3] 2 , and supp ϕ ⊂ [-2/3, 2/3] 2 , and 0≤ ϕ ≤ 1 such that k∈Z 2 ϕ(x -k) = 1.For any smooth function ρ : R 2 -→ R 2 which is Z 2 -periodic, we then set for K > 0(S 2K ρ)(x) = |k|≤2K, k∈Z 2 ϕ(x -k)ρ(x).Therefore we get for K > 0 large enough|B K | | C0 ⋆ R 2 ρ| L p ((-1/2,1/2) 2 ) + 0(1/K) ≤| C0 ⋆ R 2 ρ| L p (BK ) ≤| C0 ⋆ R 2 (S 2K ρ)| L p (R 2 ) + | C0 ⋆ (ρ -(S 2K ρ))| L p (BK ) ≤ C m |∇(S 2K ρ)| L p (R 2 ) + +| C0 ⋆ (ρ -(S 2K ρ))| L p (BK ) ≤ C m |B 2K | |∇ρ| (L p (-1/2,1/2)) 2×2 + 0(1/K) + |ρ| (L ∞ (R 2 )) 2 |B K | |z|≥K-C0 (z)|.Using the fact that |z|≥K-1 | C0 (z)| = O(1/K), dividing by |B K | and taking the limit as K -→ +∞, we get

. 1 ) 0 g(s) -ρ 0 2 (L 2 (T 2 )) 2 -

 102222 Now let ϕ ∈ C ∞ c ([0, +∞), R) be such that ϕ ≥ 0. Using (5.1), we get thatt 0 ds g m (s) -ρ 0 2 (L 2 (T 2 )) 2 ϕ(s) ≤ C t 0 ds s ϕ(s).Using Fatou's Lemma, we deduce that t Cs ϕ(s) ≤ 0.

Acknowledgement. The authors thank A. El Hajj for fruitful discussion in the preparation of this paper. This work was supported by the ACI JC 1025.


Helena Peić Tukuljac
email: helena.peictukuljac@epfl.ch

Antoine Deleforge
email: antoine.deleforge@inria.fr

Rémi Gribonval
email: remi.gribonval@inria.fr

MULAN: A Blind and Off-Grid Method for Multichannel Echo Retrieval

This paper addresses the general problem of blind echo retrieval, i.e., given M sensors measuring in the discrete-time domain M mixtures of K delayed and attenuated copies of an unknown source signal, can the echo locations and weights be recovered? This problem has broad applications in fields such as sonars, seismology, ultrasounds or room acoustics. It belongs to the broader class of blind channel identification problems, which have been intensively studied in signal processing. Existing methods in the literature proceed in two steps: (i) blind estimation of sparse discrete-time filters and (ii) echo information retrieval by peak-picking on filters. The precision of these methods is fundamentally limited by the rate at which the signals are sampled: estimated echo locations are necessary on-grid, and since true locations never match the sampling grid, the weight estimation precision is impacted. This is the so-called basis-mismatch problem in compressed sensing. We propose a radically different approach to the problem, building on the framework of finite-rate-of-innovation sampling. The approach operates directly in the parameter-space of echo locations and weights, and enables near-exact blind and off-grid echo retrieval from discrete-time measurements. It is shown to outperform conventional methods by several orders of magnitude in precision.

Introduction

When a wave propagates from a point source through a medium and is reflected on surfaces before reaching sensors, the measured signals consist of mixtures of the direct path signal with delayed and attenuated copies of itself. This physical phenomenon is commonly referred to as echoes and has a wide range of applications in different areas of science, from sonars [START_REF] Kleeman | Sonar sensing[END_REF] to seismology [START_REF] Sato | Seismic wave propagation and scattering in the heterogeneous earth[END_REF], from acoustics [START_REF] Dokmanić | Acoustic echoes reveal room shape[END_REF][START_REF] Crocco | Towards fully uncalibrated room reconstruction with sound[END_REF][START_REF] Crocco | Estimation of TDOA for room reflections by iterative weighted L1 constraint[END_REF] to ultrasounds [START_REF] Achim | Compressive sensing for ultrasound RF echoes using a-stable distributions[END_REF]. For instance, in acoustics, it has been shown that precise knowledge of early echo timing enables the estimation of the positions of reflective surfaces in a room [START_REF] Dokmanić | Acoustic echoes reveal room shape[END_REF][START_REF] Crocco | Towards fully uncalibrated room reconstruction with sound[END_REF]. In [START_REF] Dokmanić | Acoustic echoes reveal room shape[END_REF], the approximate 3D geometry of Lausanne cathedral could be retrieved in this way. On the other hand, echoes' attenuation capture information about the acoustic impedance of surfaces, which is notoriously hard to measure or estimate in practice [START_REF] Antonello | Identification of surface acoustic impedances in a reverberant room using the FDTD method[END_REF][START_REF] Bertin | Joint estimation of sound source location and boundary impedance with physics-driven cosparse regularization[END_REF]. In [START_REF] Dokmanić | Raking the cocktail party[END_REF] and [START_REF] Scheibler | Separake: Source separation with a little help from echoes[END_REF], it is shown that knowing the attenuation and timing of early echoes may improve beamforming and source separation performance, respectively. Systems using echoes for beamforming are commonly referred to as rake receivers in the wireless literature [START_REF] Price | A communication technique for multipath channels[END_REF].

Retrieving echo properties when the emitted signal is known is referred to as active echolocation in biology, and is well exemplified by the sensory system of echoing bats. This principle is for instance 32nd Conference on Neural Information Processing Systems (NeurIPS 2018), Montréal, Canada.

at the heart of active sonar technologies. In the signal processing literature, this problem belongs to the category of system or channel identification, i.e., estimating the filters from a known input to the observed output of a linear system. In the case of echoes, these linear filters consist of streams of Diracs in the continuous-time domain and are hence sparse in the discrete-time domain. The more challenging problem of estimating echoes/filters when the emitted signal is unknown is referred to as passive echolocation or blind system identification (BSI) [START_REF] Xu | A least-squares approach to blind channel identification[END_REF][START_REF] Hua | Fast maximum likelihood for blind identification of multiple FIR channels[END_REF][START_REF] Abed-Meraim | A subspace algorithm for certain blind identification problems[END_REF][START_REF] Aissa-El-Bey | Blind SIMO channel identification using a sparsity criterion[END_REF][START_REF] Lin | Blind channel identification for speech dereverberation using l1-norm sparse learning[END_REF][START_REF] Kammoun | Robustness of blind subspace based techniques using Lp quasi-norms[END_REF][START_REF] Kowalczyk | Blind system identification using sparse learning for tdoa estimation of room reflections[END_REF][START_REF] Crocco | Estimation of TDOA for room reflections by iterative weighted L1 constraint[END_REF][START_REF] Li | Multichannel identification and nonnegative equalization for dereverberation and noise reduction based on convolutive transfer function[END_REF][START_REF] Lee | Fast and guaranteed blind multichannel deconvolution under a bilinear system model[END_REF]. BSI is a long-standing and still active research topic in signal processing, notably due to its fundamental illposedness. In the general setting of arbitrary signals and filters, rigorous theoretical ambiguities under which the problem is unsolvable have been identified [START_REF] Xu | A least-squares approach to blind channel identification[END_REF]. A number of methods for multichannel BSI with general signals and filters have been developed some time ago [START_REF] Xu | A least-squares approach to blind channel identification[END_REF][START_REF] Hua | Fast maximum likelihood for blind identification of multiple FIR channels[END_REF][START_REF] Abed-Meraim | A subspace algorithm for certain blind identification problems[END_REF]. Some wellknown limitations of these approaches are their sensitivity to the chosen length of filters, and their intractability when the filters are too large. Following the compressed sensing wave [START_REF] Candès | An introduction to compressive sampling[END_REF], a number of methods extending these BSI methods to the case of sparse [START_REF] Aissa-El-Bey | Blind SIMO channel identification using a sparsity criterion[END_REF][START_REF] Lin | Blind channel identification for speech dereverberation using l1-norm sparse learning[END_REF][START_REF] Kammoun | Robustness of blind subspace based techniques using Lp quasi-norms[END_REF][START_REF] Kowalczyk | Blind system identification using sparse learning for tdoa estimation of room reflections[END_REF][START_REF] Crocco | Estimation of TDOA for room reflections by iterative weighted L1 constraint[END_REF] or structured [START_REF] Lee | Fast and guaranteed blind multichannel deconvolution under a bilinear system model[END_REF] filters have been developed. They generally extend classical methods using regularizers such as the 1 -norm for sparsity or a bilinear constraint as in [START_REF] Lee | Fast and guaranteed blind multichannel deconvolution under a bilinear system model[END_REF]. Similarly to classical filter estimation methods, they require knowledge of the filters' length and they work in the space of discrete-time filters which are typically thousands of samples long. Because they work in the discrete-time domain, the accuracy at which these methods can recover echo locations is fundamentally limited by the signal's frequency of sampling: the recovered echoes are on-grid. Moreover, the sparsity assumption on filters is invalid in practice due to smoothing and sampling effects at sensors. Interestingly, [START_REF] Chi | Guaranteed blind sparse spikes deconvolution via lifting and convex optimization[END_REF] employs a continuous-time spike model for single-channel blind deconvolution but relies on a strong linear prior on the signal.

In this paper, we propose a drastically different approach to blind echo retrieval based on the framework of finite-rate-of-innovation (FRI) sampling [START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF][START_REF] Blu | Sparse sampling of signal innovations[END_REF][START_REF] Ye | Compressive sampling using annihilating filter-based low-rank interpolation[END_REF]. In stark contrast with existing methods, the approach directly operates in the space of continuous-time echoes, and is hence able to blindly recover their locations off-grid. The proposed method is shown to recover echo delays and attenuation with an accuracy far higher than what the sampling rate would normally allow, using noiseless multichannel discrete-time measurements of an unknown simulated speech emitter in a room. The method does not assume that the filters are finite-length and only requires the number of echoes. The remainder of this paper is organized as follows. In section 2 the signal and measurement models and notations are introduced, and conventional methods are briefly reviewed. In section 3 the proposed approach is presented in the non-blind and blind cases. In section 4 the method is compared and evaluated on both synthetic and room acoustic data. Conclusions and future directions are outlined in section 5.

Background

The signal and measurement models

We start by defining the signal model in the continuous-time domain. Suppose a source emits a band-limited signal s(t) which is reflected and attenuated K times before reaching M sensors. The continuous signal impinging at sensor m is

x m (t) = ( h m * s)(t) (1)
where h m is a linear filter from the source to sensor m and * denotes the continuous convolution operator defined by

(x * y)(t) = +∞ -∞ x(u)y(t -u)du. (2)
The filter consists of the following stream of Diracs:

h m (t) = K k=1 c m,k δ(t -τ m,k ), (3)
where δ denotes the Dirac delta function, {τ m,k } K k=1 denote the K propagation times from the source to sensor m in seconds, i.e. the echo delays or Dirac locations and {c m,k } K k=1 denote the echo attenuations or Dirac weights. In practical applications, continuous time-domain signals are not accessible. They are measured by sensors and discretized to be stored in a computer's memory. Let xm ∈ R N denote N consecutive discrete samples collected by sensor m. Most measurement models assume that the impinging signal undergoes an ideal low-pass filter with frequency support [-F s /2, F s /2] before being regularly sampled at the rate F s in Hz. This is expressed by

xm (n) = ( φ * x m )(n/F s ), n = 0, . . . , N -1 (4)
where φ = sin(πt)/πt is the classical sinc sampling Kernel. The continuous-time model (1) can then be approximated in two different ways, described in the next two sub-sections.

Discrete time-domain model

First, model (1) can be approximated in the discrete, finite-time domain. Let ĥm ∈ R L and ŝ ∈ R N +L-1 denote discrete, sampled versions of the filter h m and signal s respectively. We then have

xm (n) ≈ ( ĥm ŝ)(n) (5)
where the discrete finite convolution operator between two vectors u ∈ R L and v ∈ R D (L ≤ D) is defined by

(u v)(n) = L-1 j=0 u(j)v(L -1 + n -j), n = 0, . . . , D -L. (6)
The following convenient matrix notation will be used in the paper:

u v = Toep 0 (u)v = Toep(v)u = (7)        u L . . .
. . . u 1            v 1 v 2 . . . v D     =      v L v L-1 . . . v 1 v L+1 v L . . . v 2 . . . . . . . . . . . . v D v D-1 . . . v D-L+1          u 1 u 2 . . . u L     ,
where Toep 0 (u) ∈ R (D-L+1)×D and Toep(v) ∈ R (D-L+1)×L . The validity of approximation (5) depends on the way h m and s are sampled. In [START_REF] Van Denboomgaard | Gaussian convolutions numerical approximations based on interpolation[END_REF](Proposition 2), it is showed that if s(t) is band-limited with maximum frequency lower than F s /2 and if we let the number of samples N and the filter length L grow to infinity, then model ( 5) is exact for the following sampling schemes:

ŝ(n) = s(n/F s ), n ∈ Z (8) ĥm (n) = ( φ * h m )(n/F s ), n ∈ Z.
(9) Here, it is important to note that contrary to intuition, even in the idealized case where an infinite number of samples are available, the discrete-time filters { ĥm } M m=1 involved in the measurement model are never streams of Diracs, but non-sparse, infinite-length filters consisting of decimated combinations of sinc functions. This is illustrated in Fig. 1. Recovering the original Dirac positions and coefficients from finitely many samples of such filters is a challenging task in itself.

Discrete frequency-domain model

Alternatively, one may approximate model [START_REF] Kleeman | Sonar sensing[END_REF] in the discrete finite-frequency domain. Let x m ∈ C F denote the discrete Fourier transform (DFT) of xm , defined by

x m (f ) = DFT(x m ) = N -1 n=0 xm (n)e -2πif n/Fs ( 10
)
where f belongs to a set of F regularly-spaced frequencies F = {f 1 , . . . , f F } ⊂]0, F s /2] in Hz. We then have the following approximate model:

x m (f ) ≈ h m (f )s(f ) ≈ K k=1 c m,k e -2πif τ m,k s(f ) (11)
where h m ∈ C F and s ∈ C F denote the DFT of ĥm and ŝ, respectively. Two approximations are made in [START_REF] Price | A communication technique for multipath channels[END_REF]. First, the time-domain convolution between ĥm and ŝ has been transformed into a multiplication through the DFT. This would be exact for a circular convolution, but the actual model is a linear convolution between infinite and non periodic signals, resulting in an approximation error. Second, the formula used for h m in the right hand side of ( 11) is the one that would result from the discrete-time Fourier transform (DTFT) of ĥm which would require infinitely many samples N to be calculated exactly, as opposed to the DFT. Note that the smoothing sinc kernel φ(t) does not impact this formula, since only frequencies below F s /2 are considered. Importantly, both approximations in [START_REF] Price | A communication technique for multipath channels[END_REF] become arbitrarily precise as the number of samples N grows to infinity.

While both the discrete-time model ( 5) and the discrete-frequency model ( 11) become increasingly accurate when N becomes large, the latter directly incorporates the variables of interest {c m,k , τ m,k } M,K m,k=1 , as opposed to the former. In the remainder of this paper, it will be assumed that s(t) is bandlimited with maximum frequency less than F s /2 and that N is sufficiently large such that both models hold very well. This is a reasonable assumption in audio applications, where sensors typically acquire tens of thousands of samples per second. Moreover, we focus on situations where sensor noise is negligible. Hence, the approximation signs will be dropped for convenience.

Existing methods in channel identification

To the best of the authors' knowledge, all existing methods in blind channel identification rely on the discrete-time model [START_REF] Crocco | Estimation of TDOA for room reflections by iterative weighted L1 constraint[END_REF] [START_REF] Xu | A least-squares approach to blind channel identification[END_REF][START_REF] Hua | Fast maximum likelihood for blind identification of multiple FIR channels[END_REF][START_REF] Abed-Meraim | A subspace algorithm for certain blind identification problems[END_REF][START_REF] Aissa-El-Bey | Blind SIMO channel identification using a sparsity criterion[END_REF][START_REF] Lin | Blind channel identification for speech dereverberation using l1-norm sparse learning[END_REF][START_REF] Kammoun | Robustness of blind subspace based techniques using Lp quasi-norms[END_REF][START_REF] Kowalczyk | Blind system identification using sparse learning for tdoa estimation of room reflections[END_REF][START_REF] Crocco | Estimation of TDOA for room reflections by iterative weighted L1 constraint[END_REF][START_REF] Li | Multichannel identification and nonnegative equalization for dereverberation and noise reduction based on convolutive transfer function[END_REF][START_REF] Lee | Fast and guaranteed blind multichannel deconvolution under a bilinear system model[END_REF]. The case of general emitted signals and finite filters was studied both methodologically and theoretically in the 90s [START_REF] Xu | A least-squares approach to blind channel identification[END_REF][START_REF] Hua | Fast maximum likelihood for blind identification of multiple FIR channels[END_REF][START_REF] Abed-Meraim | A subspace algorithm for certain blind identification problems[END_REF], where two main categories of methods emerged, which we briefly review here, focusing on the two-channel (M = 2) case for simplicity. First, the so-called subspace methods rely on the estimation of a time-domain M P × M P covariance matrix where P is a time-window length that must be larger than the filters' length L [START_REF] Abed-Meraim | A subspace algorithm for certain blind identification problems[END_REF]. The filters are estimated by spectral decomposition of this matrix. Second, the more common cross-relation (CR) methods rely on the observation that under noiseless conditions we have ĥm xlĥl xm = 0 N -L+1 for l = m ∈ {1, . . . , M }, by associativity of the convolution. A common approach is therefore to solve a minimization problem of the form:

ĥ * 1 , ĥ * 2 = argmin ĥ1(1)=1 Toep(x 2 ) ĥ1 -Toep(x 1 ) ĥ2 2 2 , (12)
which is a simple least-square problem. The constraint ĥ1 (1) = 1 is used to avoid the trivial solution ĥ1 = ĥ2 = 0 L . Alternatively, the normalization ĥ1 In the case of interest where the goal is to retrieve echo information from the filters, both subspace [START_REF] Kammoun | Robustness of blind subspace based techniques using Lp quasi-norms[END_REF] and to a larger extent CR [START_REF] Aissa-El-Bey | Blind SIMO channel identification using a sparsity criterion[END_REF][START_REF] Lin | Blind channel identification for speech dereverberation using l1-norm sparse learning[END_REF][START_REF] Kowalczyk | Blind system identification using sparse learning for tdoa estimation of room reflections[END_REF][START_REF] Crocco | Estimation of TDOA for room reflections by iterative weighted L1 constraint[END_REF] methods have been extended in order to handle sparse filters. This approach requires two independent steps: first estimating sparse filters, second retrieving echo locations and weights from them, typically using a peak-picking technique. Following the compressed sensing idea [START_REF] Candès | An introduction to compressive sampling[END_REF], sparsity is usually promoted using an 1 -norm penalty term on the filters. For instance in [START_REF] Lin | Blind channel identification for speech dereverberation using l1-norm sparse learning[END_REF], the following LASSO-type [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] problem is considered:

ĥ * 1 , ĥ * 2 = argmin ĥ1(1)=1 Toep(x 2 ) ĥ1 -Toep(x 1 ) ĥ2 2 2 + λ( ĥ1 1 + ĥ2 1 ) (13)
and a Bayesian-learning method for the automatic inference of λ is proposed. Several other approaches relying on similar schemes [START_REF] Aissa-El-Bey | Blind SIMO channel identification using a sparsity criterion[END_REF][START_REF] Kowalczyk | Blind system identification using sparse learning for tdoa estimation of room reflections[END_REF][START_REF] Crocco | Estimation of TDOA for room reflections by iterative weighted L1 constraint[END_REF] have been proposed.

Four important bottlenecks of discrete-time methods for echo retrieval can be identified:

• Although they rely on sparsity-enforcing regularizers, the filters are strictly-speaking nonsparse in practice, due to the sinc kernel (Fig. 1). This general bottle-neck in compressed sensing has been referred to as basis mismatch and was notably studied in [START_REF] Chi | Sensitivity to basis mismatch in compressed sensing[END_REF]. In particular, the true peaks of the filters do not correspond to the true echoes (Fig. 1), even for N → ∞. Though, most existing methods rely on peak-picking [START_REF] Kowalczyk | Blind system identification using sparse learning for tdoa estimation of room reflections[END_REF][START_REF] Crocco | Estimation of TDOA for room reflections by iterative weighted L1 constraint[END_REF]. • For the same reason, these methods are fundamentally on-grid, namely, they can only output echo locations which are integer multiple of the sampling period 1/F s . This prevents subsample resolution, which may be important in applications such as room shape reconstruction from audio signals [START_REF] Dokmanić | Acoustic echoes reveal room shape[END_REF]. • These methods strongly rely on the knowledge of the length L of the filters. However, due to the sinc kernel (Sec. 2.1.1), the true filters are always infinite. • The dimension of the search space is M L-1, which is much larger in practice than the actual number 2M K of unknown variables. This makes the methods computationally demanding and sometimes intractable for large filter lengths (typically in the tens of thousands for acoustic applications).

3 Off-grid echo retrieval by multichannel annihilation

In this section, we introduce a novel method for echo recovery that makes use of the discrete-frequency model [START_REF] Price | A communication technique for multipath channels[END_REF] and overcomes a number of shortcomings of existing approaches. Namely, it works directly in the parameter space, it does not rely on the filters' length but on the number of echoes, and it enables exact off-grid recovery of echoes' locations and weights in the noiseless case. The approach relies on the finite rate of innovation (FRI) sampling paradigm introduced in [START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF]. This is the first time this paradigm is applied to blind channel identification, to the best of the authors' knowledge.

The non-blind case

We start by considering the non-blind case where the emitted signal s ∈ C F in the discrete frequency domain is known. We further assume throughout the paper that this signal is nonzero on the considered frequency grid F = {f 1 , . . . , f F }. We can then transform the discrete-frequency model ( 11) by writing:

h m (f ) = x m (f )z(f ) = K k=1 c m,k e -2πif τ m,k (14)
where the Fourier-inverted signal z ∈ C F is defined by z(f ) = s(f ) -1 . Our goal is to estimate {c m,k , τ m,k } K k=1 from h m = x m z m , where denotes the Hadamard product. If we take our frequency indexes F to be in arithmetic progression with step ∆ f , then the exponential sequence {e -2πifiτ m,k } F i=1 is a geometric progression with ratio r m,k = e -2πi∆ f τ m,k for each m, k. Hence, h m is a weighted sum of geometric progressions. This enables us to use the so called annihilating filter technique [START_REF] Stoica | Introduction to spectral analysis[END_REF]. This technique is based on the observation that

[1, -w] [w 0 , w 1 , w 2 , . . . , w F -1 ] = 0 F -1 , (15)
for any w ∈ C and F ∈ N. We deduce that if we define the filter a m = [a m,0 , . . . , a m,K ] ∈ C K+1 as the following discrete convolution 1 of K filters of size 2:

a m = [1, -r m,1 ] [1, -r m,2 ] • • • [1, -r m,K-1 ] [0 K-1 , 1, -r m,K , 0 K-1 ], (16)
then a m is an annihilating filter for h m , i.e., a m h m = 0 F -K . Importantly, the number of echoes K has to be known upfront in order to define a m . Let us now define the polynomial representation of filter a m by:

P am [y] = K k=0 a m,k y k . ( 17
)
Because a m is an annihilating filter for h m , it follows from the classical interpretation of convolution as polynomial multiplication that P am has exactly K roots, which are the ratios {r m,k } K k=1 . Hence, once an annihilating filter a m for h m has been found, the Dirac locations {τ m,k } K k=1 can be deduced by rooting P am . Once the roots are known, reconstructing the weights is a simple linear problem involving a Vandermonde matrix V(r m ) ∈ C F ×K , obtained by writing [START_REF] Abed-Meraim | A subspace algorithm for certain blind identification problems[END_REF] in matrix form:

    h m (f 1 ) h m (f 2 )
. . .

h m (f F )     =      1 1 . . . 1 r 1 m,1 r 1 m,2 . . . r 1 m,K . . . . . . . . . . . . r F -1 m,1 r F -1 m,2 . . . r F -1 m,K      D m     c m,1 c m,2 . . . c m,K     = V(r m )D m c m . ( 18
)
where D m = Diag(e -2πif1τ m ) ∈ C K×K . The least-square solution of this system is given by c m = D -1 m V(r m ) † h m where {•} † denotes the Moore-Penrose pseudo inverse. In practice, since positive weights are sought, the phases of this complex vector are discarded. General FRI theory [START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF] tells us that F ≥ 2K + 1 is enough to uniquely recover the exact K Dirac locations and weights using this method. In other words, the original echo retrieval problem has been reduced to that of finding an annihilating filter for h m = x m z. In practice, this can be done by solving the following minimization problem for m = 1, . . . , M :

a * m = argmin am 2 2 =1 Toep(x m z)a m 2 2 , ( 19
)
where the unit norm constraint is used to avoid the trivial solution a m = 0 K+1 . The solution of this problem is the eigenvector associated to the lowest eigenvalue of Toep(x m z). Assuming that the true z is given, that model ( 11) holds exactly and that F ≥ 2K + 1, this eigenvalue will be unique and equal to 0.

MULAN: an iterative scheme

In the blind echo retrieval problem of interest, the emitted signal s and hence z are unknown. To solve for all unknown variables jointly, we introduce the following non-convex optimization problem:

z * , a * 1 , . . . , a * M = argmin z 2 2 = a1 2 2 =•••= a M 2 2 =1 M m=1 Toep(x m z)a m 2 2 . ( 20
)
Our strategy to tackle this problem is by alternated minimization with respect to each variable. Minimization with respect to each a m is already covered by the previous section. Minimization with respect to z is also a minimum eigenvalue problem, since the cost function C(z, a) can be rewritten:

C(z, a) = M m=1 Toep(x m z)a m 2 2 = M m=1 Toep 0 (a m )Diag(x m )z 2 2 = Qz 2 2 , (21)
where Q = [Toep 0 (a 1 )Diag(x 1 ); . . .

; Toep 0 (a M )Diag(x M )] ∈ C M (K+1)×F (22)
and [•; •] denotes vertical concatenation. If the algorithm succeeds in bringing down the cost function to zero, it means that appropriate annihilating filters have been found for all channels for a given Fourier-inverted signal z, and the locations and weights of all Diracs can be recovered. We call this method MULAN for MULtichannel ANnihilation. Pseudo-code for the algorithm is given in Alg.

1. Since ( 20) is non-convex, the alternate minimization scheme is at best guaranteed to converge to a stationary point of the cost-function C(z, a). To alleviate this issue, we propose to initialize the method multiple times with random values of z and only keep the run with lowest final C(z, a).

Algorithm 1 MULAN (MULtichannel ANnihilation)

Input: Frequency-domain multichannel measurements {x 1:M (f ); f ∈ F} computed via DFT [START_REF] Scheibler | Separake: Source separation with a little help from echoes[END_REF]; max_iter; conv_thresh. Output: Echo locations and weights {τ m,k , c m,k } M,K m,k=1 . 1: iter := 0; z := random(); // i.i.d. standard complex Gaussian in C F 2: repeat 3:

iter := iter + 1;

Identifiability and ambiguities

The identifiability of blind channel identification for general discrete filters and signals has been studied some time ago [START_REF] Xu | A least-squares approach to blind channel identification[END_REF]. It is known that the filters { ĥm } M m=1 cannot be recovered if their polynomial representations admit at least a common root or if the polynomial representation of the emitted signal ŝ has less than 2L + 1 roots. The latter is ruled out if the emitted signal has a rich enough spectral content (enough nonzero frequencies) which is usually the case for natural signals.

The former has at least one consequence in our case: the problem is unidentifiable if the observed signals are scaled and delayed versions of each other, which may happen in practice. While other common roots may appear in the general setting, it is important to note that MULAN restricts the search of filters to those which are linear combinations of geometrical series in the frequency domain.

There is no complete theoretical study on common roots in this case, to the best of the authors' knowledge. The authors of [START_REF] Choudhary | On the properties of the rank-two null space of nonsparse and canonical-sparse blind deconvolution[END_REF] theoretically studied blind deconvolution of sparse signals, but their results do not apply here since our filters are not sparse (see Sec. 2.2). Another well-known ambiguity is that the filters can only be recovered up to a global time-shift and scaling, because a converse shifting and scaling of the emitted signal yields the same observations. We handle this by adopting the convention τ 1,1 = 0 and c 1,1 = 1. Additionally, we assume that all echoes are located in the first half of temporal filters to avoid time-wrapping ambiguities. Finally, the proposed MULAN algorithm has an extra specific ambiguity. It can be easily shown that multiplying the roots of all polynomials {P am } M m=1 by a complex scalar γ while dividing the Fourier-inverted signal z element-wise by a geometric series of ratio γ does not change the cost function C(z, a). This can be handled by rescaling the roots of all annihilating filters to have unit modulus at each iteration. However, since only the complex arguments of the roots are used in the end, this appeared to be unnecessary in our experiments.

Experiments

On-grid vs. off-grid echo retrieval

We first emphasize the specific ability of the proposed method to recover echo locations off-grid by comparing it to conventional on-grid methods on a simulated room-acoustic scenario and on an artificial scenario with truly sparse discrete filters for reference. For the room-acoustic scenario, there is a point source emitting speech from the TIMIT dataset [START_REF] Garofolo | DARPA TIMIT acoustic-phonetic continous speech corpus CD-ROM[END_REF], and M = 2 microphones are randomly placed inside 100 random shoe-box rooms whose sizes vary from 4m × 6m × 8m to 5m × 7m × 9m. Simulations were performed using the pyroomacoustics library [START_REF] Scheibler | Pyroomacoustics: A python package for audio room simulations and array processing algorithms[END_REF]. The absorption coefficient of each surface of the room is set to 0.2. Only first-order reflections on the 6 surfaces and the direct path are simulated, resulting in K = 7 echoes per channel and filters shorter than 50 ms. For each experiment, it was ensured that the minimum separation of echoes was 1ms. The filters are simulated in the continuous-time domain using the image-source method [START_REF] Jont | Image method for efficiently simulating small-room acoustics[END_REF]. They are then smoothed, sampled and convolved with the source signal at F s = 16kHz according to the measurement model described in Sec. 2.1. The ground-truth echo locations and weights are saved in the time-domain before smoothing and are hence off-grid. The M -channel input signals used are 0.25s long, i.e., N = 0.25F s = 4000 samples. On the other hand, for the artificial scenario, the speech source was discretely convolved with sparse filters of similar length with K = 7 nonzero elements each resulting in N = 4000 samples of M -channel observations. The ground-truth echo locations and weights are hence on-grid in this case. All weights take values between 0 and 1.

For MULAN, the DFT (eq. 10) is applied to each input signal using a grid F of F = 401 regularly spaced frequencies between 200 Hz and 2000 Hz. Such a choice of the frequency range avoids low-frequency bands which are often noisy in real scenario, while focusing on a typical spectral range for speech, but it can be easily adapted depending on the application. An odd number of frequencies was chosen, since it has proven to be a good practice [START_REF] Blu | Sparse sampling of signal innovations[END_REF]. We use 20 random initializations as a good compromise between global convergence and computational complexity, max_iter= 1000 and conv_thresh= 0.1%. The two baseline methods chosen are CR [START_REF] Xu | A least-squares approach to blind channel identification[END_REF] as described in [START_REF] Xu | A least-squares approach to blind channel identification[END_REF] and its LASSO-type extension [START_REF] Lin | Blind channel identification for speech dereverberation using l1-norm sparse learning[END_REF] as described in [START_REF] Hua | Fast maximum likelihood for blind identification of multiple FIR channels[END_REF]. The filters lengths L were always set to the true lengths (which never exceed 0.05F s ) and the sparsity parameter λ for LASSO was manually set to λ = 10 -3 , which empirically showed best performance among the choices {10 -6 , 10 -5 , . . . , 10 2 }, although any value below 10 -2 showed similar performance.

We used two distinct metrics to evaluate Dirac location estimation and Dirac weight estimation. For the first one, a test is counted as successful if the root mean squared error (RMSE) of the 7 × 2 = 14 Dirac locations is below 1 sample (1/F s seconds), and the success rate out of 100 tests is provided. This metric only counts fully successful channel recovery and penalizes tests where some Diracs are missed or completely off. For the second one, we provide the weight RMSE of successful tests only. This is to avoid counting weights estimated at wrong Dirac locations. These metrics for 100 onand off-grid tests and all three methods are showed in Table 1. We can see that for the on-grid case, both CR and MULAN perform well, CR even achieving more location recoveries than MULAN. This is not too surprising since CR is based on the on-grid artificial model, while MULAN uses an off-grid model. We observed that LASSO struggles with the proximity of Diracs and did not perform as well. In terms of weight estimation MULAN yields errors which are 2 to 3 orders of magnitudes smaller than the two competing methods, which is very encouraging. In the more realistic off-grid scenario, we observed that localization errors of CR and LASSO drastically degrades with almost no successful channel estimation. Meanwhile, MULAN achieves near-exact full recovery of locations and weights in 70 out of 100 tests.

Influence of K, M , F on recovery rate

We now conduct further experiments to check the influence of parameters K, M and F on the ability of MULAN to fully recover Dirac locations and weights off-grid. We show results with 20 random initializations, F = 201 or F = 401 in the same frequency range as before, M ∈ {2, . . . , 7} and K ∈ {2, . . . , 7}. The following RMSE thresholds were defined for success of recovery: 1 sample for locations as before and 10 -2 for weights. 100 experiments were performed for every parameter set. Results for F = 201 can be seen in Figures 2 and3, and for F = 401 in Figures 4 and5. As can be seen, a higher recovery rate is generally observed when fewer echoes are present and more frequencies are used. On the other hand, the number of sensors does not significantly affect recovery performance. This is expected since O(KM ) parameters are estimated from O(M F ) observations. Increasing the number of random initializations also showed to increase success by alleviating the non-convexity of the problem, at the cost of an increased computational requirement. 1: Ratio of full Dirac location recovery (RMSE < 1 sample = 1/F s seconds) and weight RMSE (successful cases only) for three methods over 100 on-grid and 100 off-grid tests. Weights take values between 0 and 1.

Conclusion

This paper introduced the first method enabling blind and off-grid recovery of echo locations and weights from discrete-time multichannel measurements, to the best of the authors' knowledge. Future work will include alternative initialization schemes and convex relaxations in the spirit of [START_REF] Chi | Guaranteed blind sparse spikes deconvolution via lifting and convex optimization[END_REF] for the proposed cost function, extensions to sparse-spectrum signals and noisy measurements, and applications to dereverberation and audio-based room shape reconstruction. A better theoretical understanding of recovery guarantees as a function of M , K, F and N will also be sought. The code for this submission can be found at: https://github.com/epfl-lts2/mulan.

Figure 1 :

 1 Figure 1: (a) Continuous-time stream of Diracs h(t), (b) sinc kernel φ(t), (c) smoothed stream ( φ * h)(t), (d) original stream h(t) (red) and its smoothed, sampled version ĥ ∈ R L (blue).

2 2 + ĥ2 2 2 = 1

 221 can be used, leading to a minimum eigenvalue problem.

4 :

 4 for m = 1 → M do: a m := min_eig_vec(Toep(x m z)); end for 5: z := min_eig_vec(Q); // See eq. (22) 6: until iter=max_iter or C(z,a) decreased by less than conv_thresh // See eq. (21) 7: for m = 1 → M do 8: r m := roots(P am ); τ m := -arg(r m )/(2π∆ f ); c m := abs(D -1 m V(r m ) † h m ); // Sec. 3.1 9: end for 10: return Shifted and scaled {τ m,k , c m,k } M,K m,k=1 ; // See Sec. 3.3

Figure 2 :

 2 Figure 2: Rate of location retrieval for F = 201.

Figure 3 :

 3 Figure 3: Rate of weight retrieval for F = 201.

Figure 4 :

 4 Figure 4: Rate of location retrieval for F = 401.Figure 5: Rate of weight retrieval for F = 401.

Figure 5 :

 5 Figure 4: Rate of location retrieval for F = 401.Figure 5: Rate of weight retrieval for F = 401.

The chained discrete convolutions in[START_REF] Lin | Blind channel identification for speech dereverberation using l1-norm sparse learning[END_REF] have to be taken from right to left to be compatible with[START_REF] Achim | Compressive sensing for ultrasound RF echoes using a-stable distributions[END_REF].