Marc Nicodème

Charles Dossal

Flavius Turcu

Yannick Berthoumieu

Lipschitz bounds for noise robustness in compressive sensing: two algorithms

The paper deals with numerical estimations of Lipschitz bounds relating locally the reconstruction error to the measurement error in the compressive sensing framework. Most recent theoretical papers in the field parametrize such bounds relatively to certain families of vectors called dual certificates, which are fundamental to several reconstruction criteria. The paper provides two algorithms for computing dual certificates that optimize their related reconstruction error bounds. We give a greedy algorithm that provides a fast approximate solution, and a convex-projection algorithm that computes the exact optimum.

I. INTRODUCTION

An important topic in compressive sensing, that comes alongside with the central problem of the recoverability of sparse signals, is the robustness to measurement noise for signals that are normally recoverable in the absence of noise. Suppose x 0 ∈ R n is an unknown sparse signal that has to be recovered from an ill-posed linear problem

y = Ax 0 + b (1)
in which y is a known vector, called observation or measure data, A ∈ R m×n is a known linear operator with m < n, called sensing or measuring operator, and b ∈ R m is an unknown vector called sensing or measuring noise. Since the sensing matrix always has a non-trivial kernel, "solving" (1) undergoes a regularization, which classically means finding the solution(s) of (1) that minimize a certain functional. In the noiseless case b = 0, a vector x0 that coincides with the unique solution of the regularized problem is called recoverable or identifiable with respect to the chosen functional.

If one is only interested in sparse signals x 0 , the natural choice for such a functional would be the ℓ 0 -pseudo-norm x 0 = #(supp(x)), i.e. the number of non-zero coordinates of x. Unfortunately solving the ℓ 0 -regularized problem is considered untractable, hence a popular choice is rather to consider a convex ℓ 1 regularization. In the noiseless case b = 0 this amounts to solving the basis pursuit [START_REF] Scott Shaobing Chen | Atomic decomposition by basis pursuit[END_REF]:

min x∈R n x 1 s.t. Ax = y (2)
In the noisy case several variants have been developed, such as the popular Lagrangian formulation, but here we will rather focus on the (equivalent) constrained optimization problem

min x∈R n x 1 s.t. Ax -y 2 ≤ δ. (3)
This model has been explored by many papers and applied in a different applications in signal and image processing. Here the identifiability of the original vector is evaluated by sign consistency and noise robustness.

More precisely, x 0 is considered identifiable with sign consistency if the solution of (3) has the same support than x 0 i.e. the positions of the non zero components. The recovery of the support generally needs strong assumptions and might hold only for small enough noises.

On the other side, noise robustness involves the ℓ 2 recovery error regardless of the support, which in our case is x 0 -x δ 2, for some minimizer x δ of the optimization problem [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF]. Noise robustness is usually expressed in terms of some Lipschitz bound C relating the ℓ 2 recovery error to the ℓ 2 norm of the noise, such as

x δ -x 0 2 ≤ Cδ, b 2 ≤ δ (4)
in the case of the problem (3). Classical theoretical estimations for such constants were given in the well-studied context of random matrices A with various prescribed column distributions. Such matrices are shown to satisfy the restricted isometry property (RIP), hence noise-robustness bounds can be derived in terms of RIP constants.

Unfortunately these RIP assumptions can only be applied for random matrices and are unsuitable in general for some given deterministic matrix. Deterministic matrices occur in various applications, as for instance in tomography, which was our initial target application.

Here sensing matrices can be viewed as discrete versions of Radon transforms, whose behavior does not match the RIP context.

Other known results applying to deterministic matrices give errors bounds estimates using a semi implicit formula of the solution of the noiseless ℓ 1 -minimization problem, but are only valid under smallnoise assumptions.

As an alternative, more recent papers focus on the local behavior of the reconstruction error. One such approach theoretically consists in emphasizing an error bound for each member of a certain family of vectors associated to a local noiseless solution. These vectors are called dual certificates and play a fundamental role in identifiability criteria based on the so-called source condition [START_REF] Burger | Convergence rates of convex variational regularization[END_REF], [START_REF] Grasmair | Linear convergence rates for tikhonov regularization with positively homogeneous functionals[END_REF]. As such, the theoretical bound estimates based on dual certificates are more versatile, since they are local -hence potentially more accurate -and need no a priori assumptions on the magnitude of the measurement noise. However, dual certificates are not unique, and finding them is not straightforward and might have a significant computational cost. The common strategy consists in building "candidates" which can be either generic or specific to a certain context. In the context of super-resolution, Candès et al [START_REF] Candes | Super-resolution from noisy data[END_REF] [3] build a trigonometric polynomial which is a certificate under certain conditions. In the same context, Duval et al [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF] constructed candidates by minimizing their ℓ 2 norm and imposing a zero-derivative whenever x 0 has a non-zero component. Also, Vaiter et al [START_REF] Vaiter | Robust sparse analysis regularization[END_REF] compute a generic candidate.

In this context involving dual certificates, theoretical robustness Lipschitz bound were given throughout [START_REF] Burger | Convergence rates of convex variational regularization[END_REF], [START_REF] Grasmair | Sparse regularization with lq penalty term[END_REF] and [START_REF] Dossal | Consistency of l1 recovery from noisy deterministic measurements[END_REF], as

x δ -x 0 2 ≤ C(η)δ, b 2 ≤ δ, (5)
where the constant C(η) depends, besides on the matrix A and the support of x 0 , on some dual certificate η ∈ R m for x 0 .

In the paper we provide two algorithms for computing dual certificates η which minimize their related Lipschitz error bounds C(η) in [START_REF] Currie | OPTI: Lowering the Barrier Between Open Source Optimizers and the Industrial MATLAB User[END_REF], based on two different approaches:

The first is a greedy algorithm that provides a fast approximate solution, by iteratively minimizing the error bound with respect to a current dual certificate candidate. This produces a fast approximation of the optimal bound, since all the numerical computations only rely on scalar products.

The second algorithm is based on a theoretical geometric result that describes the exact optimal solution. More precisely we showed that the problem of finding the dual certificate minimizing the error bound can be reduced to computing the convex projection of a suitable vector over a certain convex set. As such, it gives the best error bound -which is attained for a unique dual certificate ηopt, but has a higher computational cost than the greedy approximation.

The two algorithms are illustrated and numerically compared for compressive sensing matrices issued from a tomography application framework. The relationship between theoretical bounds and numerical observed bounds is also discussed.

Notations: Throughout the paper, (ei)i=1,...,n stands for the canonical base in R n , A ∈ R m×n denotes a linear operator from R n to R m , and ai = Aei, (i = 1, . . . , n) the columns of A.

For x = [x1, . . . , xn] ∈ R n , its support is denoted

supp(x) = {i ∈ {1, . . . , n} | xi = 0},
and its vector sign by sign(x) = [sign(xi)]i=1,...,n. Also, for some subset I ⊂ {1, . . . , n}, we denote #(I) the cardinal of I, we write xI for the restriction of x to its components indexed by I, and AI for the submatrix of A with the columns ai, i ∈ I.

The notation A * stands for the adjoint of A, while A + denotes the Moore-Penrose pseudoinverse of A, i.e. A + = (A * A) -1 A * . The p, q-operator norm of a matrix A is

A p,q = max x =0 Ax q x p and A p,p = A p (6) Note that A 1,2 = max i ai 2 and A + 2 = (A * A)-1 2.
Finally, for any subset M of R n , we write Aff(M ) and Conv(M ) for the affine, respectively convex envelope of M . In particular, B1 = Conv({±ei}i=1,...,n) stands for the closed unit ball of R n under the ℓ 1 -norm.

II. STATE OF THE ART

Let us begin by recalling the notion of dual certificate, its related source (or range) condition from regularization theory [START_REF] Burger | Convergence rates of convex variational regularization[END_REF] and their connection with identifiability and noise robustness in the framework of ℓ 1 -minimization. Fix a linear operator A ∈ R m×n from R n to R m , a vector x 0 in R n and let I = supp(x 0 ). A vector η in R m is called a dual certificate for x 0 if the vector A * η belongs to the subgradient of the ℓ 1 norm at x 0 . This latter set can be viewed as

∂ x 0 1 = {u ∈ R n : uI = sign(xI ) and u ∞ ≤ 1} (7)
If η satisfies only the first condition above, i.e.

(A * η)I = sign(x 0 I ) (8)
then η is usually called a candidate or a pre-certificate. Also, x 0 is said to satisfy the source condition if it admits at least one dual certificate.

Let us first discuss these notions in connection to the identifiability of x 0 in the noiseless case, which means, as mentioned in the introduction, the property of x 0 to be the unique minimizer of the ℓ 1 -regularized problem [START_REF] Candes | Super-resolution from noisy data[END_REF].

One can easily check that the source condition on x 0 is a sufficient condition for x 0 to be a minimizer of ( 2), but it does not guarantee its uniqueness.

One crucial step towards uniqueness is to turn the second condition in [START_REF] Dossal | Consistency of l1 recovery from noisy deterministic measurements[END_REF] into a strict inequality outside the support of x 0 , i.e. requesting that

(A * η)J ∞ < 1, J = {1, . . . , n} \ I (9)
This strengthened condition insures that no minimizer of (2) can be supported anywhere but inside I. Searching for dual certificates with this strengthened condition has been proposed by several authors. It has been first exploited by Fuchs in [START_REF] Fuchs | On sparse representations in arbitrary redundant bases[END_REF], where the inequality ( 9) is tested for the generic pre-certificate

d(x 0 ) = A + * I sign(x 0 ) (10)
in order to derive a sufficient identifiability criterion. This latter condition is is also known as the irrepresentability condition in statistics, see [START_REF] Dossal | Consistency of l1 recovery from noisy deterministic measurements[END_REF] for more details. Rather than testing a particular candidate, a more general sufficient criterion was given by Vaiter et al in [START_REF] Vaiter | Robust sparse analysis regularization[END_REF], where convex optimization is used to solve the problem

IC(x 0 ) = min η∈R m A * J η ∞ s. t. A * I η = sign(x 0 I ) (11)
whose solutions are dual certificates satisfying (9) whenever

IC(x 0 ) < 1.
Finally, even if x 0 admits a dual certificate satisfying [START_REF] Fuchs | On sparse representations in arbitrary redundant bases[END_REF], in order for x 0 to be the unique minimizer of (2) one has to make sure that there are no other minimizers supported inside I. One can see that this happens if and only if AI is injective.

Summing up, this latter injectivity condition, together with ( 8) and ( 9) represent a strengthened version of the source condition that was considered by Grasmair et al in [START_REF] Grasmair | Linear convergence rates for tikhonov regularization with positively homogeneous functionals[END_REF][13] [START_REF] Grasmair | Sparse regularization with lq penalty term[END_REF]to describe the noiseless ℓ 1 -identifiability.

Dual certificates have a deep geometric interpretation in ℓ 1identifiability regarded in the context of polytope theory, as considered by Donoho in [START_REF] David | Neighborly polytopes and sparse solutions of underdetermined linear equations[END_REF]. In fact, from the main result there it follows that x 0 is identifiable if and only if the convex set

F (x 0 ) := Conv({sgn(x 0 i )ai : i ∈ I}) (12)
is a face of affine dimension #(I) -1 of the image polytope

P := AB1 = Conv({±ai : i = 1, . . . , n}) (13)
Now the fact that F (x 0 ) is a face of P means by definition that there is a linear functional on R m , thus a vector η ∈ R m (via Riesz representation theorem), that separates F (x 0 ) from the rest of the polytope P . This translates by

η, sgn(x 0 i )ai = 1 i ∈ I, (14)
η, ±aj < 1 j ∈ J , (15)
which are nothing but equivalent forms of ( 8) and ( 9) respectively. Moreover the injectivity of AI is equivalent to the fact that F (x 0 ) has affine dimension #(I) -1. This means precisely that dual certificates, in the sense of the strengthened source condition, can be regarded as separating functionals for faces of the image polytope P . This geometric interpretation is also very useful in the matter of noise robustness, that will be discussed next.

Let us suppose now that x 0 is noiseless-identifiable, thus satisfies the strengthened source condition for some dual certificate η. Assume that x 0 is measured in (1) with some noise 0 < b 2 ≤ δ and denote x δ a minimizer of the regularized problem [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF]. Estimating the noise robustness of the reconstruction deals with measuring the reconstruction error x 0 -x δ , which unfortunately has no implicit formula in general. It is much instead easier to measure this error along various directions, each direction producing an error estimate, with more or less accuracy. Among all possible directions, the most "fiable" are those given by the dual certificates, in the sense of scalar products of the type A * η, x 0 -x δ . Such quantities characterize the Bregman distance between x 0 and x δ , which was used by Burger et al. [START_REF] Burger | Convergence rates of convex variational regularization[END_REF] in order to measure the recovery error with a very simple Lipschitz constant. However, the Bregman distance, as a pseudodistance, is a rather weak measure, and one would prefer an ℓ 2estimation instead. This was done by Grasmair et al. [START_REF] Grasmair | Necessary and sufficient conditions for linear convergence of l1-regularization[END_REF] by relating the Bregman distance to the ℓ 2 norm. The resulting Lipschitz constant was afterwards reduced by Dossal et al [START_REF] Dossal | Consistency of l1 recovery from noisy deterministic measurements[END_REF], leading to the following result:

Theorem 1 ([7]
). For any minimizer x δ of (3) and any dual certificate η satisfying the strengthened source condition (9) one has

x δ -x 0 ≤ C(η)δ ( 16
)
with

C(η) = 2 A + I 2 + ( A + I 2 AJ 1,2 + 1)Q(η) (17)
Q(η) = η 2 1 -A * J η ∞ (18)
Clearly the optimal Lipschitz constant is obtained for the dual certificate that minimizes the factor Q(η), which solely depends on η in (17). This leads to the optimization problem

min Q(η) : (A * η)I = sign(x 0 I ), (A * η)J ∞ < 1 (19)
We give in the next section the explicit formula for the solution of this problem and the corresponding algorithm that computes the optimal dual certificate and its related Lipschitz bound.

III. A CONVEX-PROJECTION ALGORITHM FOR THE OPTIMAL LIPSCHITZ BOUND

We saw in the last section that a dual certificate η for an identifiable vector x 0 supported on I has a geometric interpretation as a separating functional for the face F (x 0 ) of the image polytope P (see [START_REF] Grasmair | Sparse regularization with lq penalty term[END_REF],( 13)). This means that the dual affine subspace of η

Hη := {u ∈ R m : η, u = 1} (20)
is a support hyperplane for the face F (x 0 ), i.e. Hη ∩ P = F (x 0 ).

The key fact in solving the problem ( 19) is that the quantity Q(η) also has a strong geometric interpretation. More precisely, consider the "reduced" polytope

PJ = Conv({±aj : j ∈ J}) (J = {1, . . . , n} \ I) (21)
One can show that:

Proposition 2. For any dual certificate η of x 0 one has

Q(η) = 1 D(Hη, PJ ) ( 22
)
where D stands for the affine euclidian distance.

Therefore in order to minimize Q(η) one has to maximize the distance from the reduced polytope PJ to a generic affine hyperplane Hη passing through the face F (x 0 ). Obviously the intersection of all these hyperplanes is the affine closure of the face F (x 0 ), denoted

L(x 0 ) := Aff(F (x 0 )) = Aff({sgn(x 0 i )ai : i ∈ I}) (23)
It is more convenient to consider K the orthogonal complement of L(x 0 ) in R m and project the whole problem onto K. Denote PK the orthogonal projector onto K. More precisely, suppose if I = {i1, ..., ip} and ĀI = ( āi)i∈{i 2 ,...,ip} is the matrix with the columns {ai -ai 1 , i = {i2, ..., ip}}, then

PK = IRm -ĀI Ā+ I . (24)
Put ãj = PK aj and consider the projected polytope PJ = PK PJ = Conv({±ãj : j ∈ J}) and Hη = PK Hη. (25)

Now the entire affine subspace L(x 0 ) projects onto a singleton, which can be shown to be precisely the vector

d(x 0 ) + = d(x 0 ) d(x 0 ) 2 2 , (26)
the dual vector of the generic pre-certificate d(x 0 ) given in [START_REF] Elmer | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF]. It is also straightforward that

Q(η) = 1 D(Hη, PJ ) = 1 D( Hη, PJ ) . ( 27
)
Note that, in contrast to (22), Hη above is now a generic affine hyperplane in K which passes just through the point d(x 0 ) + , and whose distance to the projected reduced polytope PJ needs to be maximized. One can use then a couple of standard convex analysis arguments to show the main result:

Theorem 3. Denote w 0 the point in PJ at minimal distance from d(x 0 ) + . Then min η Q(η) = 1 D(L(x0), PJ ) = 1 D(d(x 0 ) + , PJ ) = = 1 d(x 0 ) + -w 0 2 (28)
Moreover, the minimal value is attained for the unique dual certificate

ηopt = d(x 0 ) 2 2 1 -d(x 0 ), w 0 (d(x 0 ) + -w 0 ). (29)
The algorithmic implementation is straightforward from the last theorem:

Algorithm 1 Convex projection algorithm for optimal Lipschitz bound

One can use any algorithm to compute the convex projection in

Step 4., which clearly is the most costly in terms of complexity. In our implementation we used an m-dimensional version of the Gilbert-Johnson-Keerthi (GJK) algorithm [START_REF] Elmer | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF]. For this algorithm for instance, each iteration comprises two steps, one step which performs a support hyperplane computation, and a second step containing a simplex convex projection. The first step is rather fast, since it relies solely on scalar products. On the contrary, the simplex projection is the costly part, since it is itself recursive.

In the next section we show how one can roughly eliminate the type of costly computation above, at the expense of some precision in the bound estimation.

IV. A GREEDY ALGORITHM FOR FAST LIPSCHITZ BOUNDS

ESTIMATES

In this section we construct an accelerated greedy algorithm 2 which only deals with scalar products, thus fast to compute, for producing approximative Lipschitz bound estimations.

The main idea is roughly to iteratively select a sequence of candidates η which decrease the denominator in (18) until a dual certificate is found. This is done by decreasing the scalar product between columns of A indexed by the cosupport J of the signal and the candidate vector η at each step.

More precisely, as for IC criterion [START_REF] Grasmair | Linear convergence rates for tikhonov regularization with positively homogeneous functionals[END_REF], we aim to construct a candidate η with a minimal A * J η ∞. The algorithm is initiated with the same generic candidate d(x0) in ( 10) which also appears in the previous algorithm.

At each step

k, A * J η ∞ is reduced by adding a new component aj k to the vector η k+1 = η k -t k A + * I k sign k such that J k = J k-1 ∪ {j k } and AJ k injective.
The vector sign k is the appropriated sign for effectively reducing the scalar products with the selected columns. The maximum of the scalar products between these selected columns is scaled down together proportionally to t k and, moreover,

| as, η k | = A * J η k ∞ (s ∈ S)
where S denotes the set of selected columns. When the vector η is modified, its scalar products with A * J are also affected and there are some scalar products which can increase. The descent factor t k is computed to have the best trade-off between the reduced scalar products and those which increase.

Empirically, the sign k is stable by

(A * J k AJ k ) -1 i.e sign((A * J k AJ k ) -1 sign k ) = sign k
Indeed, if the sign of the scalar product of an already selected column is inverted, the scalar product scalar will increase instead of decreasing. Consequently, the factor t k will be smaller and the maximum of the scalar products will be limited. This stability can be an indicator of a column which is no more needed.

Finally, an optional elimination step can be used for cleaning the superfluous columns at the end. For enough sparse x 0 , this elimination step is not necessary. while AI k is injective do 6:

j k ← max j k / ∈Ik | aj k , η | 7:
sign k ← sign k ∪ {sign( aj k , η } 8:

I k ← I k ∪ j k 9: t k ← min i / ∈I k ti = a i ,η -a j k ,η (A * A + * I k sign k ) i -(A * A + * I k sign k ) j k such that ti > 0 10: η ← η -t k A + * I k sign k 11:
end while 12:

I k _bad ← component of I k such that A * I k sign k = sign k 13:
Elimination of components I k (I k _bad)

14:

Elimination of components sign k (I k _bad) 15: end while

V. NUMERICAL RESULTS

In this section, we create synthetic sets of sparse image x 0 ∈ R n of size x = 64 × 64. The operator A of the computed tomography will be a classical Radon transform with parallel lines. There is no real gap with a fan-beam like operator. The Radon transform is performed with the code written by Justin K. Romberg. In the noiseless case, the basis pursuit has been computed using linear programming with the Matlab optimization toolbox and OPTI toolbox [START_REF] Currie | OPTI: Lowering the Barrier Between Open Source Optimizers and the Industrial MATLAB User[END_REF] ( to speed up the calculations). We obtain the same results with those toolbox. In the noisy case, (3) was computed with SPGL1 toolbox [START_REF] Van Den | Probing the pareto frontier for basis pursuit solutions[END_REF] with the best precision.

A. Empirical performance of the greedy algorithm for noiseless identifiability

In this subsection we present the numerical results for the greedy algorithm regarding its ability to produce dual certificates, and compare it with other existing criteria.

We used as a reference an empirical phase transition of the ℓ 1 minimization is obtained, showing the degree of sparsity one can expect for within the identifiability zone. Different sets of sparse images were built for each degree of sparsity. A first one, a nonnegative, where non-zero components are set to 1 and the second where non-zero components can take -1 or 1 as value. At first, a hundred signals are selected for a given sparsity and a certain number of angles. The basis pursuit is performed on each signals of a given sparsity and a probability of success is deduced. The operation is repeated for each sparsity. Results are shown as a function of the sparsity in Figure 1.

For the same set of selected sparse images, the pre-certificates are computed via the greedy algorithm and it is checked whether they are dual certificates. As for the basis pursuit, results are shown as a ratio of the number of dual certificates over the pre-certificates function of sparsity. The label ηIC denotes the candidate produced by IC criterion. d is for the Fuchs criterion that only checks the particular candidate d(x 0 ) -see (10) -, ηg is the greedy algorithm without the elimination step and ηg2 is the greedy algorithm with elimination.

The projection algorithm is not represented in the Figure 2 because it finds a dual certificate for every identifiable x 0 . It is a very good way to estimate the success of sparse approximations by ℓ 1 minimization. There is only one issue, the cost of time computation The ratio of pre-certificates d, ηg, ηg2 and ηIC which are dual certificates as a function of x 0 0 can be near or worst than ℓ 1 minimization itself as it is show in the next section. The IC covers almost the whole set of identifiable signals but it gives the worst Q. The pre-certificate ηg2 is a dual certificate for nearly all of identifiable signals. One can note that the enhanced greedy algorithm has a huge gain over the simple greedy approach. The Fuchs criterion is verified only on signals with a low degree of sparsity.

B. Projection vs Greedy

This subsection presents a comparison between the projection algorithm and the greedy algorithm regarding time complexity and quality of the Lipschitz bounds.

The table I and II present the Q obtained for different degrees of sparsity.

When the degree of sparsity approaches the phase transition, Q is higher and the time computation explodes for the greedy algorithm, Table II: Results for Q obtained by the greedy algorithm because of the elimination of bad components. Close to the phase transition, finding a dual certificate is harder and there are many steps of elimination, thus the algorithm struggles in a long loop. Under a certain degree of sparsity, the fast greedy algorithm leads to dual certificates such that the correspondent Q is neared to those of the projection algorithm.

C. Theoretical bounds versus real bounds

We end this section by a discussion on the quality of the theoretical Lipschitz bounds. One general feature is that these bounds always emphasize a worst case, so they can be very pessimistic for an arbitrary noise. A classical approach for testing how realistic they are is to compute the recovery of a signal with a large number of gaussian noises. In the context of computed tomography, the choice of the directions of the ray, i.e the angles between the rays and an axis, is important. All the results are obtained with a constant angular pitch, except in the worst case where we first choose the directions which produces the best results for our purpose. We also construct a "worst case" noise.

Locally, one knows that the solution x λ with J λ = supp(x λ ) is

x λ J λ = x 0 J λ + A + * J λ b -λ(A * J λ A * J λ ) -1 sign(x λ J λ ) (30)
One can expect to maximize x λ J λ -x 0 J λ 2 by building a noise vector b based on the last support J λ during the execution of (3) where the matrix (A * J λ AJ λ ) -1 is not well-conditioned. It could work for a small noise in the sense of ℓ 2 . Indeed, a higher noise can impact the sign. Following this idea, the algorithm 3 is proposed.

Algorithm 3 Building a worst case noise vector

1: Find the last support J λ such that (A * J λ AJ λ ) -1 is not well- conditioned 2: (U, S, V ) ← SV D((A * J λ AJ λ ) -1
) where SV D is the Singular Value Decomposition. The table III presents the best results with the "worst case" noise vector for small noise and it compares to a hundred gaussian noises and the theoretical boundary. As it is shown, the initial gap in the literature has been filled, while there is still a significant difference.

One sees that the theoretical bounds differ from real bounds by some varying factor. In the case of tomography matrices, this can be explained even in the continuous domain. Indeed, if Rf is the Radon Table III: Ratio between theoretical results and real ℓ 2 -error. T B is the theoretical boundary, mean is the mean ℓ 2 -error with gaussian noises, idem with max and WC is the ℓ 2 -error with the "worst" case noise.

x λ J λ -x 0 J λ 2 Mean Max WC T B 1 / x λ J λ -x 0 J λ 2 from
transform operator, Rf * Rf can be seen as deconvolution operator. Duval et al [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF] show in the context of sparse spike deconvolution that there can be small Dirac movements. Thus one cannot expect a too small ℓ 2 -recovery error. Moreover, in the discrete domain this lead to a further increase of the ℓ 2 -error when the resolution of the discrete grid increases. In this context, even for small noises, where one can generally hope a better approximation, the theoretical estimation of the noise robustness still remains at some factor away. One can note that this bound not only hold for a solution x δ but for every solution x δ such that x δ 1 ≤ x 0 1. Indeed, this theoretical result covers the worst case for this set of solutions of (3). It also explains very well the gap with the reality. Even so, the figure 3 shows a correlation between the theoretical results and the reality. There is approximatively a ratio K.

VI. CONCLUSION

The presented algorithms allow the numerical optimization of some of the most recent and versatile theoretical results concerning noiserobustness in compressive sensing. Numerical results show that they are rather complementary with respect to different sparsity regimes. Also, one can choose one or the other according to the potential application priorities, speed versus quality. In addition, due to their inherently and natural geometrical content, part of the ideas behind the proposed algorithms is susceptible to be adapted for future refinements of theoretical noise-robustness estimations.

1 : 2 ; 2 :; 3 :

 1223 Put d(x 0 ) = A + * I sign(x 0 I ) and d(x 0 ) + = d(x 0 )/ d(x 0 ) 2 Compute the projection matrix onto K: PK = IRm -ĀI Ā+ I Put ãj = PK aj, j ∈ J; 4: Compute the convex projection w0 of d(x 0 ) + onto the convex closure of ±ãj, j ∈ J 5: Compute Qopt = d(x 0 ) + -w 0 -1 2

Algorithm 2 2 :I 3 :

 223 Greedy algorithm to compute a pre-certificate with elimination 1: Calculate η0 = A + * I sign(x 0 I ) sign k ← vector zeros of the length of x 0 I k _bad ← (1) 4: while I k _bad is non empty do 5:

Figure 1 :Figure 2 :

 12 Figure 1: Examples of x 0 ∈ R 64×64

3 :

 3 b ← V (1) where V (1) is the vector with the maximum singular value. 4: normalize b.

Figure 3 :

 3 Figure 3: The theoretical boundary in function of the mean of the recovery error of a thousand different gaussian noises.

Table I :

 I Results for Q obtained by the projection algorithm

	sparsity	80			100		150
		Q	time(sec)	Q	time(sec)	Q	time(sec)
	mean	11.2262	11.9678	15.29	11.38	37.44	15.13
	median	11.1787	11.8658	15.27	11.33	36.71	15.01
	min	9.3402	9.9907	12.45	10.47	26.53	13.15
	max	13.4216	16.3981	19.35	12.85	56.54	17.37
	sparsity	80			100		150
		Q	time(sec)	Q	time(sec)	Q	time(sec)
	mean	11.5985	3.7014	15.84	4.58	42.95	9.64
	median	11.5188	3.5650	15.74	4.44	41.25	9.51
	min	9.5915	2.8845	12.83	3.53	29.14	6.71
	max	13.7379	5.7184	20.02	6.36	72.16	14.54

ACKNOWLEDGEMENT

This work was supported by the TOTAL SA. Authors would like to thank P. Ricoux of TOTAL SA. for his advice and useful discussions.


Martin Devaud
email: martin.devaud@univ-paris-diderot.fr

Éric Bringuier

Thierry Hocquet

Acoustics in the Lagrange picture: an application to the Rayleigh radiation pressure

At the undergraduate level, most lectures and textbooks on hydrodynamics make use of the so-called Euler picture, where the pressure, temperature and velocity of the fluid are treated as continuous fields defined by the value they take at each point of the reference frame the fluid moves in. There nevertheless exists another possible description of the movement which consists in labelling the fluid elements themselves, and keeping this labelling in the course of the motion. This so-called Lagrange picture is scarcely taught for it often introducess more complicated mathematics, as soon as a three-dimensional geometry is considered. Yet it is actually more intuitive than the Euler picture. In this paper, we illustrate the point in the example of the Rayleigh acoustic radiation pressure. An improved physical insight ensues, which is of interest to students graduating in acoustics.

I. INTRODUCTION

Fluid mechanics is often presented [START_REF] Landau | Fluid Mechanics[END_REF][START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] and taught, at least at the undergraduate level, in a way quite similar to electrodynamics. To put it briefly, a (usually Galilean) reference frame is chosen. At a given point r of this frame, and at time t, the physical state of the fluid is described by a set of functions: mass density ρ, pressure P , fluid velocity v (with respect to the frame), temperature T , and so on. So one deals with a set of (coupled) continuous fields ρ( r, t), P ( r, t), v( r, t), T ( r, t), etc. For instance, in the absence of any external force (gravity or else), the movement of an inviscid 19 fluid is governed by the well-known Euler equation,

ρ( r, t) ∂ v ∂t + ( v • --→ grad) v = - --→ grad P. (1)
This description of the fluid motion will be referred to in the present paper as the Euler picture. Such a picture has many advantages. First, it is convenient to describe the dynamics of the fluid by means of local equations coupling fields, exactly as in electrodynamics. The Euler equation ( 1) is an example of such a local description. Second, the Euler picture is particularly well suited to situations in which the fluid really flows: when studying the stream of a river passing under some bridge, we are interested in the very behaviour of the water under this bridge at a time t, whatever the origin or the past behaviour of this water. Nevertheless, the Euler picture has a few drawbacks. First, in equation ( 1), the left-hand side is obviously nonlinear in field v, due to the ( v • --→ grad) v term. A second drawback of the Euler picture shows up when free boundary conditions between two fluids must be imposed. Let us consider the example illustrated in figure 1: two different fluids -say 1 and 2 -are separated at rest by the (infinite) plane at x = 0. Consider a plane pressure wave propagating from x = -∞ in fluid 1 towards the boundary. As is well known, this incident wave splits at the interface into two parts: a reflected wave, travelling back to x = -∞ through medium 1, and a transmitted wave, travelling towards x = + ∞ through medium 2. It is a time-honoured undergraduate level exercise to determine the reflection and transmission coefficients at the interface. In principle, the answer is easy: the continuity of pressure (due to the finite acceleration of a fluid element),

P 1 (interface) = P 2 (interface), (2a)
and of velocitiy (due to the absence of a gap),

v 1 (interface) = v 2 (interface), (2b)
provides two equations enabling us to calculate both coefficients. The above reasoning is undoubtedly correct, but raises a non-trivial difficulty: where is the interface? At x = 0? Certainly not, since the interface itself moves back medium 1

x = 0 x = + ∞ medium 2 x = -∞
FIG. 1: Medium 1 and medium 2 are separated at rest by the plane x = 0. When a pressure wave propagates, the interface does not remain at x = 0, but moves back and forth on either side of the plane x = 0. In the Euler picture, locating the interface at x = 0 appears as a zero-order approximation.

and forth, due to the wave motion. As a matter of fact, the x = 0 plane spends half the time in medium 1 and half the time in medium 2. Of course, locating the interface at x = 0 is the best approximation, and it leads to the correct values of the reflection and transmission coefficients, but it should be regarded only as a zero-order approximation.

In fact, the difficulties of the Euler picture mentioned above can (up to a point) be overcome using a different framework, known as the Lagrange picture [START_REF] Philipp | Methods of Theoretical Physics[END_REF][START_REF] Guyon | Hydrodynamique Physique[END_REF][START_REF] Bennett | Lagrangian Fluid Dynamics[END_REF] . [START_REF]Lagrange picture" are the exact transposition to fluid mechanics of the so-called "Schrödinger picture" and "Heisenberg picture" of quantum mechanics. 21[END_REF] It is precisely the aim of section II to sketch the main features of the Lagrange picture which is hardly taught in academic courses and scarcely used when studying acoustic wave propagation in fluids. The reason is that the Lagrange picture brings in some involved mathematics like tensor calculus and differential geometry as soon as one considers a three-dimensional propagation. Consequently, for the sake of simplicity, we shall restrict ourselves to the one-dimensional case and show in subsection II A that, in the Lagrange picture, linearity or nonlinearity is a purely thermodynamic issue. Subsection II B is devoted to the linear thermodynamic response and the attendant picture of sound propagation. Section III deals with the nonlinearity of the thermodynamic response and tackles the cumbersome problem of the Rayleigh acoustic radiation pressure, which can be given a simpler solution in the Lagrange picture with greater physical insight.

II. AN OUTLINE OF THE LAGRANGE PICTURE

A. What is it all about?

Contrary to the Euler picture, which labels the geometric points of the reference frame disregarding the origin of the fluid elements passing through these points at time t, the Lagrange picture labels the fluid elements disregarding the position they occupy at time t. More specifically, consider a fluid at some time t 0 . We denote by r 0 the fluid element that occurs to stand at point r 0 of the reference frame at time t 0 . We shall henceforth keep this label r 0 to denote this fluid element, whatever its later position. Thus, at time t > t 0 , the fluid element r 0 will be found at some point r given by r( r 0 , t) = r 0 + u( r 0 , t),

where u( r 0 , t) is the displacement undergone by the fluid element r 0 between times t 0 and t (see fig. 2). The physical state of the fluid is still described by a set of continuous fields: mass density, pressure, velocity, temperature, etc. The correspondence between both pictures is very simple. With superscripts E and L respectively standing for "Euler" and "Lagrange", and quantity A standing for whichever parameter ρ, P , v, T , etc., we have

A E r( r 0 , t), t = A L ( r 0 , t), (4)
with r( r 0 , t) given by (3). Concretely, the above equation means that A L ( r 0 , t) denotes the actual value of parameter A taken at time t by the fluid element labelled r 0 which is currently at point r 0 + u( r 0 , t) of the reference frame, i.e.

A L ( r 0 , t) mathematically coincides with to A E r( r 0 , t), t .

In appendix A, the Lagrange picture is applied to the calculation of the reflection and transmission coefficients at an interface between two fluids. It is shown that the Lagrange picture offers several advantages from a technical point of view. To begin with, this picture rids us of the formal nonlinearity associated with the ( v • --→ grad) v term on the left-hand side of the Euler equation (1). This enables us to recognize genuine nonlinearities, thereby allowing a perturbative resolution of the field equations. We shall take advantage of such a simplification in section III, when dealing with the Rayleigh acoustic radiation pressure.

We focus on the simplest situation one may have to face: the one-dimensional problem. Let us therefore consider a fluid at rest occupying a cylindrical volume with axis Ox 0 and cross-sectional area S (fig. 3a), at equilibrium pressure

x 0 x y z y 0 z 0 ρ 0 , P 0 time t 0 reference position r 0 volume V 0 time t position r = r 0 + u volume V ρ, P u (r0 ,t ) -→ -→ -→ -→ -→ -→ FIG.
2: An element of fluid with reference position r0 is displaced and deformed in the course of time. In the Lagrange picture, this fluid element is labelled r0 and keeps this label throughout its motion.

(r 0 , P 0 ) (a) P 0 and mass density ρ 0 . Both ends, labelled x 0 = 0 and x 0 = L, are bounded by pistons that are provisionally supposed to be fixed. As displayed in fig. 3, the slice of fluid located between faces x 0 and x 0 + dx 0 has mass ρ 0 Sdx 0 , where ρ 0 is the equilibrium mass density. At time t, its thickness is

x 0 = 0 x 0 x 0 + dx 0 x 0 = L 0 (b) S S x 0 = 0 x 0 x 0 + dx 0 x 0 + u(x 0 ) x 0 + dx 0 + u(x 0 + dx 0 ) x 0 = L 0
x 0 + dx 0 + u(x 0 + dx 0 , t) -x 0 + u(x 0 , t) = 1 + ∂u ∂x 0 dx 0 , (5a)
so that its mass density is just

ρ(x 0 , t) = ρ 0 1 + ∂u ∂x 0 . (5b)
In the Lagrange picture, the pressure forces undergone by the slice of fluid are respectively SP 0 (x 0 , t) (left end) and -SP (x 0 + dx 0 , t) (right end). Applying Newton's Second Law to the slice, we obtain

ρ 0 ∂ L2 u ∂t 2 = - ∂P ∂x 0 , (6)
where superscript L on the left-hand side recalls that the time derivative is understood at constant x 0 (even if the current position at time t of the face labelled "x 0 " is miles away from the point with abscissa x 0 of the reference frame). Equation ( 6) deserves two comments: (i ) it is exact (no approximation was made); (ii ) it is strictly linear in displacement u or in velocity v = ∂ L u ∂t . If we now want to get a closed-form propagation equation, we have to connect the pressure P (x 0 , t) with the expansion factor ∂u ∂x0 or equivalently the mass density ρ. This connection involves thermodynamics. Throughout the present article, we shall assume, for the sake of simplicity, that any transformation undergone by the fluid is isentropic. In the framework of the Lagrange picture, this means that the entropy of any fluid slice [x 0 , x 0 + dx 0 ] is, at any time, equal to its equilibrium value. So, in the course of the motion, the pressure P (x 0 , t) can be expressed as a function of the sole 21 mass density ρ(x 0 , t). In the isentropic relation of state P = P (ρ), let us expand the extra pressure P (x 0 , t) -P 0 in increasing powers of (ρ(x 0 , t) -ρ 0 )/ρ 0 = 1 + ∂u ∂x0 -1 :

P (x 0 , t) -P 0 = -κ 1 ∂u ∂x 0 + 1 2 κ 2 ∂u ∂x 0 2 + • • • , (7)
where the compressibility κ 1 > 0, due to the Second Law of thermodynamics. Combining the mechanical equation ( 6), in which the superscript L for "Lagrange" is henceforth omitted, with the above thermodynamic relation (7), we get the sound propagation equation

ρ 0 ∂ 2 u ∂t 2 = κ 1 ∂ 2 u ∂x 2 0 1 - κ 2 κ 1 ∂u ∂x 0 + • • • . ( 8
)
The above equation is nonlinear in displacement u, its nonlinearity originating exclusively in the κ 2 , κ 3 , etc. terms in the thermodynamic expansion (7).

B. The linear thermodynamic response and the propagation of sound

In this subsection, we linearize equation ( 7), i.e. we take κ 2 = κ 3 = • • • = 0. The propagation equation ( 8) also becomes linear, and reads

1 c 2 ∂ 2 u ∂t 2 = ∂ 2 u ∂x 2 0 , ( 9a
) with c = κ 1 ρ 0 , (9b)
which is the usual d'Alembert wave equation for sound propagation. Let us look for the associated eigenmodes, i.e. the monochromatic solutions of (9a). Owing to the boundary conditions we have chosen, they necessarily read

u n (x 0 , t) = ℜ A n sin(k n x 0 ) e -iωnt , (10a)
with ℜ denoting the real part. In (10a),

ω n = ck n , k n = nπ L 0 (n = 1, 2, • • • ), (10b)
and A n is a complex amplitude. Note that, insofar as the linearization of the thermodynamic relation ( 7) is relevant, the above solution is exact, contrary to the solution generally proposed in the framework of the Euler picture, which also involves neglecting the ( v • --→ grad) v term. 22 Let us now determine the overall acoustic energy associated with the wave, i.e. the variation (with respect to the rest state) of the total energy of all slices [x 0 , x 0 + dx 0 ]. Since there is neither a heat exchange between neighbouring slices nor an external force, we just have to determine the work done by the pressure force to drive each fluid slice from its equilibrium state to its current state at time t. For the [x 0 , x 0 + dx 0 ] slice, the work is exactly

dE = t 0 dt ′ -SP (x 0 + dx 0 , t ′ ) ∂u(x 0 + dx 0 , t ′ ) ∂t ′ + SP (x 0 , t ′ ) ∂u(x 0 , t ′ ) ∂t ′ = -Sdx 0 t 0 dt ′ ∂ ∂x 0 P ∂u ∂t ′ = -Sdx 0 t 0 dt ′ ∂P ∂x 0 ∂u ∂t ′ + P ∂ 2 u ∂x 0 ∂t ′ . ( 11a
)
Owing to (6) and to the linearized version of (7), the above equation becomes

dE = Sdx 0 1 2 ρ 0 ∂u(x 0 , t) ∂t 2 -P 0 ∂u(x 0 , t) ∂x 0 + 1 2 κ 1 ∂u(x 0 , t) ∂x 0 2 . ( 11b
)
Integrating over the whole fluid and accounting the boundary conditions, we finally get the overall acoustic energy E, which is a constant of the movement:

E = 1 2 ρ 0 S L0 0 dx 0 ∂u ∂t 2 + c 2 ∂u ∂x 0 2 . (11c)
Now, since any solution u(x 0 , t) of the wave equation (9a) is a linear combination of eigenmodes of the type (10a), the above energy E may also be written, according to Parseval theorem,

E = 1 4 ρ 0 SL 0 ∞ n=1 |A n | 2 ω 2 n . (11d)
We complete these results with the following thought experiment. Suppose that, while a given eigenmode (say, n) is established in the cylindrical cavity bounded by the two pistons displayed in fig. 3, we slowly move the piston located at the end labelled "x 0 = L 0 " at, say, a constant velocity V . By "slowly", we mean "adiabatically in the Ehrenfest sense". In this connection, let us define the effective number N n of quanta in mode n such that [START_REF]The reduced Planck constant is introduced here artificially[END_REF] We have discussed at some length this issue in a former paper [START_REF] Devaud | The adiabatic invariant of the n-degree-of-freedom harmonic oscillator[END_REF] , and shown that, in the course of such an adiabatic parametric excitation of the system, the number N n of quanta is conserved. Let us recall that, in the Lagrange picture, the label "x 0 = L 0 " of the fluid in contact with the moving piston remains unchanged, although the total length of the fluid column is obviously L(t) = L 0 + V t. In this respect, it is convenient to split the displacement u(x 0 , t) into two parts, and let

N n ω n = 1 4 ρ 0 SL 0 |A n | 2 ω 2 n , so that E = ∞ n=1 N n ω n .
u(x 0 , t) = x 0 L 0 V t + w(x 0 , t). ( 12
)
The first term on the right-hand side is the displacement of the slice labelled x 0 , associated with a quasistatic expansion (or compression, according to the sign of V ) of the fluid. The second term is the extra displacement of the fluid slice due to the acoustic wave. Observe that the boundary conditions for w(x 0 , t) are the same as for u(x 0 , t): w(x 0 = 0, t) = w(x 0 = L 0 , t) = 0. Now, let us rewrite the wave equation (9a) in terms of w instead of u. From ( 12), we get

1 c 2 ∂ 2 w ∂t 2 = ∂ 2 w ∂x 2 0 , (13)
so that w and u are governed by the same equation. Solution (10a) is consequently unchanged: surprising though it may be, the motion of the piston has strictly no influence upon eigenmode n. In particular, no frequency shift 24 occurs: the wave number k n = nπ L0 as well as the angular frequency ω n = ck n keep their initial values even if L(t) happens to become twice (or half) its initial value L 0 . Moreover, from equation (11d) one concludes that, since neither the quanta number N n (Ehrenfest adiabaticity) nor the angular frequency ω n are modified, the acoustic energy is unchanged. To move the piston, the operator has of course to account for the quasistatic variation of the instantaneous equilibrium pressure P eq (t) of the fluid:

P (x 0 , t) -P 0 = -κ 1 ∂u ∂x 0 = -κ 1 V t L 0 + ∂w ∂x 0 P (x 0 , t) -P eq (t) = -κ 1 ∂w ∂x 0 , (14a)
with

P eq (t) = P 0 -κ 1 V t L 0 , (14b)
but he has no extra work to supply, associated with the acoustic wave itself. Before finishing with our thought experiment, we should emphasize one point: although equation ( 13) holds whatever the value of velocity V , the Ehrenfest adiabaticity is required during the initial acceleration of the piston, from V = 0 to its final speed.

Another interesting point is: what happens when our acoustic wave meets an interface, as illustrated in fig. 1? This point has been raised in the introduction (see equations (2) and the attendant text). The Lagrange picture is shown in appendix A to provide a straightforward answer to this question.

The propagation of longitudinal expansion/compression waves through a mass-distributed spring (e.g. such as those designed as decorative objects or toys for children) is well described by the equations of this subsection, but the ideal fluid considered in the calculations does not exist. In the next subsection, we consider a more realistic approximation of the thermodynamic relation (7), better suited to real fluids.

III. TAKING NONLINEARITY INTO ACCOUNT: THE RAYLEIGH RADIATION PRESSURE

Let us now retain the second-order term in (7):

P (x 0 , t) -P 0 = -κ 1 ∂u ∂x 0 + 1 2 κ 2 ∂u ∂x 0 2 . ( 15
)
Accounting for the nonlinearity κ 2 has the following consequence. Suppose that we move the right piston (labelled "x 0 = L 0 ") by an amount δL, and consider the motion of the fluid with respect to this new equilibrium position. Then we have, replacing V t with δL in (12),

u(x 0 , t) = x 0 L 0 δL + w(x 0 , t), (16a)
so that (15) becomes

P (x 0 , t) -P eq = -κ ′ 1 ∂w ∂x 0 + 1 2 κ 2 ∂w ∂x 0 2 , (16b)
with

P eq = P 0 -κ 1 δL L 0 + 1 2 κ 2 δL L 0 2 , ( 16c
)
κ ′ 1 = κ 1 1 - κ 2 κ 1 δL L 0 . ( 16d
)
Expressions (16b) and ( 15) are alike. They differ in the equilibrium pressure, as was already noticed in the linear case ((16c) is just the generalization of (14b) with V t replaced by δL). They differ also in the linear compressibility coefficient κ 1 being changed into κ ′ 1 , due to the nonzero value of κ 2 . This change of the compressibility causes a frequency shift when moving piston L 0 , as well as an acoustic radiation pressure. We shall come back to this nomenclature later and show that the frequency shift and the radiation pressure are inherently entangled and, so to say, consubstantial.

How can we solve for u(x 0 , t) obeying the nonlinear propagation equation ( 8)? Unless an exact solution can be found, a good approach is the perturbation method, provided that the condition | ∂u ∂x0 | ≪ 1 is fulfilled. For the sake of simplicity, let us start from an eigenmode of the linearized wave equation, say (see (10a))

u (1) (x 0 , t) = A sin(kx 0 ) cos(ωt -ϕ), (17a)
where we have dropped the eigenmode index n. While u [START_REF] Landau | Fluid Mechanics[END_REF] (x 0 , t) is not a solution of the non linear wave equation, the exact solution having u [START_REF] Landau | Fluid Mechanics[END_REF] as its linear approximation can be expanded in increasing powers of amplitude A:

u(x 0 , t) = u (1) (x 0 , t) + u (2) (x 0 , t) + • • • , (17b)
where u (i) ∼ A i . Let us look for u [START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] (x 0 , t). Using ( 8) and (9b), we have

1 c 2 ∂ 2 u (2) ∂t 2 - ∂ 2 u (2) ∂x 2 0 = - κ 2 κ 1 ∂ 2 u (1) ∂x 2 0 ∂u (1) ∂x 0 . ( 18a
)
The above equation means that the first-order solution u [START_REF] Landau | Fluid Mechanics[END_REF] acts as a source term for the second-order displacement u [START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] . From (17a), this source term is:

- κ 2 κ 1 ∂ 2 u (1) ∂x 2 0 ∂u (1) ∂x 0 = κ 2 κ 1 A 2 k 3 2 sin(2k 0 x 0 ) 1 + cos(2ωt -2ϕ) 2 . ( 18b
)
As a result, solution u [START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] is the sum of two contributions: one is static and the other is oscillating at the angular frequency 2ω. Let us focus on the former contribution. Accounting for the boundary conditions u [START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] with index "s" standing for "static". Now it is interesting to calculate, up to the second order in amplitude A, the static extra pressure P [2] s -P 0 associated with the acoustic mode. From (15), we get P [START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] (x 0 , t)

(x 0 = 0, t) = u (2) (x 0 = L 0 , t) = 0, it is u (2) s (x 0 ) = κ 2 κ 1 A 2 k 16 sin(2kx 0 ), (19) (a)
-P 0 = -κ 1 ∂u (1) ∂x 0 + ∂u (2) ∂x 0 + 1 2 κ 2 ∂u (1) ∂x 0 2 . ( 20
)
The first-order extra pressure term -κ 1 ∂u [START_REF] Landau | Fluid Mechanics[END_REF] ∂x0 oscillates at the angular frequency ω, and consequently averages to zero over time. The second-order extra pressure term is the sum of a static and a 2ω-oscillating contribution. Focussing on the former contribution, we find, after all calculations have been carried out, P [START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] s -

P 0 = 1 8 κ 2 A 2 k 2 . ( 21
) Quantity P [2]
s -P 0 is homogeneous (P [2] s does not depend on x 0 ), as expected from a static extra pressure (otherwise it would entail a permanent flow). This static extra pressure is known as the Rayleigh radiation pressure. At this juncture, a comparison with electromagnetic waves is of interest.

The existence of radiation pressure exerted by an electromagnetic wave onto some encountered obstacle is well known and can be easily figured out even at the undergraduate level. The electric field E of the wave moves electric charges in the obstacle (free carriers in a conductor, bound charges in a dielectric), resulting in an electric current density . Then, the magnetic field B of the wave exerts the Lorentz force density  × B onto the charges and the latter force is transmitted to the obstacle's lattice. The calculation of the overall force undergone by the obstacle is easy in the one-dimensional case (plane wave encountering a plane obstacle). We owe to Maxwell the first prediction of this effect in 1871 [START_REF] Clerk | A Treatise on Electricity and Magnetism[END_REF] . This phenomenon can also be looked at in a quantum framework: the incident electromagnetic wave is regarded as a flux of photons, each photon with energy hν carrying a mechanical momentum hν/c. When the wave interacts with matter, the latter momentum (or a part of it) is transferred to the lattice and the classical result is recovered [START_REF] Cagnac | Modern Atomic Physics: Fundamental Principles[END_REF] . Such a radiation pressure was experimentally observed in 1899 by P. N. Lebedev 9 using a device analogous to the Nichols radiometer.

In analogy with electrodynamics, J.W. Rayleigh introduced the acoustic radiation pressure as early as 1902 [START_REF] Rayleigh | On the pressure of vibrations[END_REF][START_REF] Rayleigh | On the momentum and pressure of gaseous vibrations, and on the connexion with the virial theorem[END_REF] . It was observed by W. Altberg in 1903 [START_REF] Wilhem | Ueber die Druckkräfte der Schallwellen und die absolute Messung der Schallintensität[END_REF] . Nevertheless, although an acoustic stress tensor can be defined, which is prima facie analogous to the Maxwell tensor, the comparison stops there: on the one hand, no microscopic mechanism providing us with an acoustic force density analogous to the Lorentz force can be put forward; on the other hand, the photon has no real acoustic counterpart: the phonon is but a quasiparticle carrying no mechanical momentum. Moreover, while there exists only one electromagnetic radiation pressure, there are at least two kinds of acoustic radiation pressure, according as the fluid in which the acoustic wave propagates is bounded (the Rayleigh configuration, illustrated in fig. 4a) or is free to skirt around the obstacle (the Langevin configuration, illustrated in fig. 4b). Since the beginning of the twentieth century, a wealth of studies have been devoted to this rather puzzling and little understood issue [START_REF] Brillouin | Les tenseurs en mécanique et en élasticité[END_REF][START_REF] Erna | Experimental studies on acoustic radiation pressure[END_REF][START_REF] Robert | Radiation pressure-the history of a mislabeled tensor[END_REF][START_REF] Boa | Acoustic radiation pressure produced by a beam of sound[END_REF][START_REF] Wang | Acoustic radiation pressure[END_REF] .

It is clear from figure 4 that the Rayleigh configuration can be implemented in a one-dimensional geometry. On the contrary, the Langevin configuration involves a two-or three-dimensional geometry. For the sake of a pedagogic account, we have focussed in this paper on the simpler one-dimensional case. Moreover we argue that the calculation is easier in the framework of the Lagrange picture.

Let us go back to formula (21). It is noteworthy that P [2] s -P 0 is proportional to coefficient κ 2 , and thus owes its very existence to the nonlinearity of the thermodynamic relation (15): this is the reason why we could not find such an extra pressure in the linear framework of subsection II B. A simple relation exists between the Rayleigh radiation pressure and the (Lagrange) acoustic energy density E that can be derived from (11d):

E = E SL 0 = 1 4 ρ 0 A 2 ω 2 = 1 4 κ 1 A 2 k 2 . ( 22a
)
According to (21), the relation is P [START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] s -P 0 = 1 2

κ 2 κ 1 E. (22b)
Lastly, let us show how deeply the radiation pressure and the frequency shift are entangled. Consider again the thought experiment discussed in subsection II B, in which we slowly move the piston at the end labelled "x 0 = L 0 ". Splitting the displacement u according to (12), we are led to modify equations (14a-b) according to (16b-c), i.e.

P (x 0 , t) -P eq (t) = -κ 1 1 - κ 2 κ 1 V t L 0 ∂w ∂x 0 + 1 2 κ 2 ∂w ∂x 0 2 , (23a)
with

P eq (t) = P 0 -κ 1 V t L 0 + 1 2 κ 2 V t L 0 2 . ( 23b
)
As was already mentioned, the important point is that the linear term -κ 1 ∂w ∂x0 is changed into -κ ′ 1 ∂w ∂x0 , due to the variation δL = V t of the length of the cavity:

κ ′ 1 = κ 1 1 - κ 2 κ 1 V t L 0 , (24)
in agreement with (16d). Consequently, the wave equation governing w becomes

ρ 0 ∂ 2 w ∂t 2 = κ ′ 1 ∂ 2 w ∂x 2 0 1 - κ 2 κ 1 ∂w ∂x 0 , (25)
which is the same equation as for u except that κ 1 is replaced by κ ′ 1 (compare for instance with ( 8)). Linearizing equation (25), we find a wave equation with a modified speed of sound c ′ given by

c ′2 = κ ′ 1 ρ 0 = c 2 1 - κ 2 κ 1 V t L 0 . ( 26a
)
This modification of the speed of sound, associated with an unchanged 25 wavevector k n = nπ L0 , shifts the angular frequency:

ω ′2 = ω 2 1 - κ 2 κ 1 V t L 0 . (26b)
Observe in passing that the amplitude A of mode n is changed, too. Nevertheless, since the piston is moved adiabatically (in the Ehrenfest sense), we have (see (11d) and the discussion thereafter)

A ′2 ω ′ = A 2 ω. (26c)
As a consequence, the change in the acoustic energy of the wave is

δE = 1 4 ρ 0 SL 0 A 2 ωδω. (27a)
Using (26b), (9b), (10b) and ( 21), this is tantamount to

δE = - 1 8 κ 2 A 2 k 2 SδL = -(P [2] s -P 0 )δV, (27b)
which shows that δE is the work the operator has to supply in order to vary the volume of the cavity by an amount δV = SδL. The acoustic pressure present is -δE/δV. This is in line with the usual definition of a pressure, given that the transformation is isentropic.

In appendix B, it is shown that the radiation pressure is related to the static relative expansion through the compressibility κ 1 . It is also shown that the present description of the radiation pressure carries over to solid state physics, where the thermal expansion of a solid may be thought of as the outcome of the Rayleigh radiation pressure.

(10a) would hold unchanged, whereas the wavevector quantification relation (10b) would become k n = n + 1 2 π L0 . As regards the energy balance, (11a) and (11b) would be unchanged while energy E should be replaced by the (conserved) quantity E + P 0 Su(L 0 , t) in (11c) and (11d). Due to the new boundary condition, piston L 0 would move, and its position become L 0 + δL(t). The static part δL s = δL(t) can be obtained by setting P 0 = P eq (δL s ) + 1 8 κ 2 A 2 k 2 , as suggested by (21), with P eq (δL s ) given by (16c). Hence an acoustically induced static expansion δLs L0 of the medium ensues, given in our second-order approximation by

δL s L 0 = 1 8 κ 2 κ 1 A 2 k 2 = 1 κ 1
(P [START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] s -P 0 ). (B1)

Now, let us suppose that not only one, but all modes, are excited simultaneously. Averaging out all interference terms, the static relation (22b) becomes P [START_REF] Georges | An Introduction to Fluid Dynamics[END_REF] s -P 0 = 1 2

κ 2 κ 1 n E n , (B2a)
where E n stands for the overall acoustic energy density in mode n (see (11d) for instance). Relation (B1) can be generalized in the same way, yielding the static expansion caused by the acoustic radiation pressure, namely

δL s L 0 = 1 2 κ 2 κ 2 1 n E n . (B2b)
We now come to the Grüneisen approach of the thermal expansion of matter in solid-state physics [START_REF] Ashcroft | Solid State Physics[END_REF] . This approach consists in regarding solids -which are de facto nonlinear compounds -as linear media with volume-dependent vibrational mode frequencies. That is to say, the angular frequency of mode n is ω n (L) and the so-called Grüneisen parameter g n of a one-dimensional solid is defined as

g n = - d ln ω n d ln L . ( B3
)
Thermal expansion originates precisely in g n being nonzero. Now, in the calculation developed in the present paper, we find (up to the second order, see (16d) or (26a))

ω 2 0 = ω 2 n0 1 - κ 2 κ 1 δL L 0 , (B4a)
whence

g n = - δ ln ω n δ ln L = 1 2 κ 2 κ 1 . (B4b)
In other words, the existence of a Rayleigh acoustic radiation pressure, on the one hand, and the thermal expansion under constant P 0 , on the other hand, originate in a non-vanishing Grüneisen parameter g = 1 2 κ2 κ1 , regardless of n. In a pictorial parlance, one may say that the thermal expansion of a solid is just the outcome of the Rayleigh radiation pressure associated with the thermally excited vibrational modes.

FIG. 3 :

 3 FIG. 3: (a). The fluid at rest, with equilibrium mass density ρ0 and pressure P0. (b). The fluid at time t: both ends, labelled x0 = 0 and x0 = L, are made up of pistons that are provisionally supposed fixed.

FIG. 4 :

 4 FIG. 4: Two kinds of acoustic radiation pressure, according to the geometry. (a) Rayleigh configuration: the fluid that propagates the acoustic wave is bounded by the obstacle. (b) Langevin configuration: the wave can skirt around the obstacle.

Acknowledgments

We acknowledge Professor Jean-Claude Bacri for stimulating conversations during the preparation of this work.

IV. CONCLUSION

In this article, we have promoted the idea that introducing the Lagrange picture of fluid dynamics could be useful in the teaching of acoustics at the undergraduate level. On the one hand, the Lagrange picture can complement the Euler picture, as our alternative derivation of the reflection/transmission coefficients in appendix A shows. On the other hand, the Lagrange picture can be superior to the Euler picture, as our treatment of the Rayleigh radiation pressure has shown. The physical meaning of that pressure if much clearer in the Lagrange picture.

Appendix A: Reflection/transmission coefficients in the Lagrange picture

We go back to the question raised in the introduction and show that the Lagrange picture provides an exact calculation of the reflection/transmission coefficients. Let us consider fig. 1 again. At rest, medium 1 and medium 2 are separated by the plane x = 0. Choosing the rest state of the system to implement the Lagrange labelling of the fluid elements, x 0 = 0 means both the right-hand face of the last slice of medium 1, and the left-hand face of the first slice of medium 2. This labelling will "follow" the motion of the system and (provided of course that no mixing occurs between the two fluids) the Lagrange labelling of the interface will remain x 0 = 0 throughout the propagation of the acoustic wave, whatever the amplitude of the latter and without any approximation. Let us consider an acoustic wave coming from x = -∞. This incident wave is described by the displacement field u i (x 0 , t) = f i t -x0 c1 , where c 1 = κ 11 /ρ 01 is the speed of sound in medium 1 and f i is any (regular) function. When the wave reaches the interface (labelled x 0 = 0, whatever its motion), it splits into a reflected wave u r (x 0 , t) = f r t + x0 c1 and a transmitted wave u t (x 0 , t) = f t t -x0 c2 , where c 2 = κ 12 /ρ 02 is the speed of sound in medium 2. In summary, the displacement field reads

while the (extra) pressure field reads, owing to (7):

where Z 1 = κ 11 /c 1 = √ κ 11 ρ 01 and Z 2 = κ 12 /c 2 = √ κ 12 ρ 02 stand for the acoustic impedances of media 1 and 2.

Writing the continuities of u and p at the interface x 0 = 0 at any time t, we obtain the well-known result

In the framework of the linear thermodynamic response, the Lagrange picture provides thus the simplest quantitative description of the reflection/transmission phenomenon of an acoustic wave at an interface.

Appendix B: Connection with the Grüneisen approach in solid-state physics

In section III, we have considered rigid boundary conditions, namely both pistons were fixed, or moved at a velocity imposed by the operator as regards piston L 0 . One may wonder how our results would be modified if, say, piston L 0 were not fixed, but subject to the external pressure P 0 . Then, the boundary condition on the fluid slice labelled L 0 would no longer be u(x 0 = L 0 ) = 0 at any time, but instead P (x 0 = L 0 ) = P 0 . Of course, equations (5a) through


Anya Traille
email: atraille@laas.fr

Sangkil Kim

Antony Coustou

Herve Aubert

Manos M Tentzeris

A Conformal/Rollable Monolithic Miniaturized Ultra-Portable Ground Penetrating Radar Using Additive and Inkjet Printing

Keywords: Monolithic, conformal, FMCW radar, Ground Penetrating Radar, flexible electronics, inkjet printing additive printing

Typical UWB FMCW Ground Penetrating (GPR) Radars operate at low frequencies that require a wide sweep bandwidth thus necessitating complex architectures and bulky broad-band antennas. This poses unique challenges to the system portability especially for manual, wide-area outdoor measurements. In this paper, we present the first design, fabrication and characterization of a complete conformal and miniaturized radar system to be rolled up in a "poster-like" container using additive printing technology. As the lumped or distributed passives, the active devices and the Rx/Tx antennas may share the same flexible substrate, the proposed radar technology is considered to be monolithic. The presented proofof-concept system performs the most fundamental operations of the FMCW radar including signal generation, as well as the amplification and correlation of the LO and RF signals for GPR frequencies. Specifically, this paper outlines an ultra low cost system integration, packaging and experimental verification of a flexible/conformal monolithic radar system with almost identical performance for different degrees of flexing.

I. INTRODUCTION

Radar is a very mature and established technology dating back to the 1920s. FMCW (frequency modulated continuous wave) radar, in particular, combines the advantages of the high-resolution impulse radar (ΔR ~ 1/B) and the low power CW (continuous wave) radar (E=P av τ p ). Other advantages include high SNR, wide dynamic range, low NF, better angular resolution, cheaper/simple low-power components (Power Amplifier and low speed Analog to Digital Converter), and frequency control/versatility, while its main disadvantages include licensing restrictions, and time gating against reflections and multipath. Both Pulsed and FMCW imaging radars have become popular in ground penetrating radar applications with Pulsed radar being more commercially available and the FMCW radar developed and used mainly by research groups.

Advancements in FMCW radar have been mostly due to major hardware developments, signal processing techniques, waveform design, and system level integration of new components into new architectures, such as the Ultra Wideband FMCW Ground Penetrating Radar (GPR) [START_REF] Holt | Sea ice thickness measurements by ultrawideband penetrating radar: First results[END_REF] which is a recently developed technology for high resolution imaging applications such as concealed weapon detection and detailed subsurface imaging. In the case of GPR, the unique challenge of Ultra wideband FMCW radar is that it requires the use of very low frequencies (<300MHz) to obtain penetration depths of at least several meters (to detect ground water for example) and a very wide sweep bandwidth to achieve submeter resolution for detailed subsurface images [START_REF] Holt | Sea ice thickness measurements by ultrawideband penetrating radar: First results[END_REF]. This makes typical GPR systems too bulky for backpack transportation and simple deployment. In addition, GPR, unlike other radar systems must perform manual wide-area measurements commonly using a vehicle, sled or by walking. There is a severe lack of flexible radar platforms providing the mechanical versatility that would facilitate geophysical studies (particularly in rugged environments such as the arctic, desert, uneven mountainous terrains or space).

In this paper we propose the first additive-printing-based packaging process of a conformal, monolithic radar system for a UWB FMCW applications operating in the VHF/UHF range for deep/shallow subsurface imaging. The main challenges of packaging a UWB FMCW radar are due to the facts that: (1) the high percentage bandwidth at low frequencies complicates the pure, linear, and spurious-free signal generation (for example limited availability steep-skirted filters) and the receiver requirements therefore requiring very complex receiver architectures [START_REF] Jankiraman | Pandora multifrequency FMCW/SFCW radar[END_REF] and, (2) at low GPR frequencies, very large broadband antennas are required which are commonly bulky.

Such a challenge can be met only with packaging modifications in both the electronic system and the antenna, where they both get combined into one single flexible, conformal, rollable and robust system using a cost effective approach. To achieve that in this paper we had to address the most current challenges in flexible electronics which are related to fabrication (metal thickness uniformity, flexible interconnects, accurate alignment), performance distortion, robustness (delamination of IC's) and lifetime (power handling). Flexible PCBs made from flexible polymer do exist however they still require photolithography and lamination fabrication processes [3]. Cost effective inkjet printing of individual components such as passives, interconnects and antennas for sensor modules has been demonstrated in [START_REF] Kim | Inkjet-printed Antennas, Sensors and Circuits on Paper Substrate[END_REF], however this is the first reported design and validation of a complete radar system which can be easily extended to multilayer densely interconnected circuits for more complex multi-mode/multi-channel radar architectures for other radar imaging applications and frequencies. The full system is referred to as monolithic as lumpe passives (hybrid coupler, matching baluns, fi active devices (VCO, LNA and Power Am Rx/Tx antennas may share here the same flex the proposed radar technology.

The following sections will show the de and measurement of the Monolithic F Penetrating Radar and access the challenges o uniformity, IC delamination, alignment interconnects. Upon completion, the entire s the antenna is to be presented at the conferenc II. MONOLITHIC FMCW RADAR DE A simple FMCW architecture using operati customized for high resolution GPR was implemented in Polyethylene terephthalate proof-of-concept demonstrator. Based on th and required output SNR (as well as other The test setup consists of a sign Voltage Controlled Oscillator between -2.5 to +9.5dBm but reduc test frequency range), which then p (Local Oscillator) and RF. The Low Input port receives the RF either fro generator, as it happened in this pap appropriately for the desired sweep mixer's required input LO drive po number of LNAs for our particul specifications were available so we reduce mixer noise (typically aro receiver chain.

The circuit footprints for each co the entire system were designed in in Fig. 2, which clearly demons feature.

III. FABRICAT

The passive interconnects were step electroless plating bath process inkjet printing and curing of cataly (ε r = 2.0, tanδ = 0.04, t=0.1 ẽlectroless growth of the metal C thickness = 0.61~ 3.8 μm).

The Inkjet printing of the catalys plating was chosen for the metalli improve the flexibility and robustn electronics including ICs and pass soldered on the circuit footprint util temperature curing. The soldered cracked unless the stresses are app This is one out of two fabrication ch system (S=signal, radar footprint in he concept of a Radar designed to kly unpacked for Fig. 3. Part of the fabricated mono with components mounted.

GPR operating from 200as calculated to be arounde contains the necessary channel transmission and esting), amplification and hitecture and estimated link nal generator that tunes the (specified output power ced due to harmonic outside produces the necessary LO w Noise Amplifier at the RF om a delay cable or a signal per. The VCO was chosen p frequency range and the ower (< 1dBm). A limited ar frequency and package chose the best candidate to ound 7dBm but varies) in omponent and the traces for CadSoft EAGLE as shown strates its rollable/portable TION fabricated using the twos [START_REF] Cook | Inkjet catalyst printing and electroless copper deposition for low-cost patterned microwave passive devices on paper[END_REF], which consists of the ytic ink onto PET substrate ~0.254mm followed by Cu film (σ = 3.2e6 S/m, st-based ink for electroless ization method in order to ess [START_REF] Kim | Inkjet-printed Antennas, Sensors and Circuits on Paper Substrate[END_REF]. The surface mount sive lump components are lizing a solder paste for low d components cannot be plied to the solder directly. hoices, with the other being olithic FMCW radar system the inkjet printing of a silver nanoparticle ink (σ = 1.1e7 S/m) onto paper, PET or LCP (ε r = 3 ~ 3.2, tan δ = 0.02 ~ 0.05, thickness = 0.254 mm) which offers a more expensive but improved electrical properties. The fabricated system is shown in Fig. 3.

IV. MEASUREMENTS Measurements were performed on the fabricated structure to validate its performance. The signal generated from the Crystek VCO (frequency range of 200 ~ 400 MHz, tuning voltage 0 ~ 5V was tested at the RF output port of the flexible system (Fig. 4).

Due to the absence of a PLL circuit to stabilize the output frequency control phase noise, the VCO needed to be characterized extensively. The VCO output was measured at the RF output port following the hybrid coupler. The full radar system was measured, starting from signal generation all the way to the IF output port (VCO, Hybrid Coupler, LNA and Mixer). The measured IF output frequency (Fig. 5) corresponded exactly to the RF-LO frequency difference s of 22 MHz, 28 MHz and 34 MHz. All components were confirmed to work exactly as specified. The only system weakness observed was the high insertion loss from the mixer. This could be attributed to the matching circuit. The output power level can be improved by adding an external LC matching network such as a T-or -π matching networks at the output port, as such topologies provide more design versatility than L matching networks. In addition, adjusting the two input power levels (LO and RFin) of the mixer (using attenuators) to be closer together may improve the mixer output signal gain.Due to the absence of a VGA equipped Automatic Gain Control circuit found in many transmitters [START_REF] Holt | Sea ice thickness measurements by ultrawideband penetrating radar: First results[END_REF], the normalized output power level changes with frequency.

A. Flexibility Test

A flexibility test was performed in order to benchmark the system performance in a rolled/flexed state. Fig 6 shows the IF output when wrapping the system around a foam cylinder of radius 4.5cm both horizontally and vertically demonstrating an almost flex-independent performance of the radar module. It has to be stressed that output power level (between -70 and -75 dBm in our preliminary experiment) was the same for both rolled and unrolled configurations.

The only parameter limiting the minimum radius of the circuit's bending curvature is limited by the largest surface mount ICs. The largest electronic component of the proposed circuit is the VCO and the whole system (100mm x 80mm) can be easily wrapped around the cylinder of radius is 4.5 cm without damaging the circuit due to delamination. The inkjetprinted catalyst-based electroless plating improves the flexibility of interconnects because the adhesion of the electroless plated conductor to the substrate is much stronger than printing of nanoparticle-based ink [START_REF] Cook | Inkjet catalyst printing and electroless copper deposition for low-cost patterned microwave passive devices on paper[END_REF].

V. CONCLUSION

The first conformal and rollable monolithic radar system has been demonstrated featuring an almost identical performance for rolled/unrolled states. The choice of wide band antenna will be the next step of this work. The antenna gain required to achieve a particular penetration depth must be determined based on the estimated receiver sensitivity and power budget calculations for maximum transmit power without exceeding emission limits. Although we are planning to test antennas such as the cone, bowtie, Vivaldi, horn and spiral are recommended for GPR applications [START_REF] Hertl | UWB Antennas for Ground Penetrating Radar Application[END_REF], we will focus on the integration with inkjet-printed origami-based "accordion" antennas that have recently led to 1-2 orders of size miniaturization and full integrability with inkjet-printed electronic platforms. The simplicity as well as the multilayer/multimaterial capability of the proposed additive printing (fully printed passive, diodes, vias) combined with flexible mounting would enable the realization of highly dense interconnects as well as enhanced miniaturization and conformability of practical GPR systems of the future, that could find numerous applications in wearable biomonitoring, smart skins and geodetection systems.

Fig. 1 .

 1 Fig. 1. Link budget estimation of flexible radar N=noise).

Fig. 2 .

 2 Fig.2. The left shows a portion of the flexible orange, the right shows the integration into th portable/conformal FMCW Ground Penetrating R be easily attached onto a backpack and quick measurements.

Fig. 5 .

 5 Fig. 5. Normalized IF Output power of the complete system.

Fig. 4 .

 4 Fig. 4. The measured VCO output (from RF output port) for various tuning voltages Vt Fig. 6. Measured IF Output of system bent horizontal and vertical.

ACKNOWLEDGEMENT

It is wished to acknowledge IRD (Institut de Recherche pour le Développement) for supporting the research and development of a new FMCW radar as well as the support of National Science Foundation (NSF).