Member, IEEE Øyvind Ryan

Senior Member, IEEE Mérouane Debbah
email: mer-ouane.debbah@supelec.fr

Asymptotic Behaviour of Random Vandermonde Matrices with Entries on the Unit Circle

Keywords: Vandermonde matrices, Random Matrices, deconvolution, limiting eigenvalue distribution, MIMO

Analytical methods for finding moments of random Vandermonde matrices with entries on the unit circle are developed. Vandermonde Matrices play an important role in signal processing and wireless applications such as direction of arrival estimation, precoding, and sparse sampling theory, just to name a few. Within this framework, we extend classical freeness results on random matrices with independent, identically distributed (i.i.d.) entries and show that Vandermonde structured matrices can be treated in the same vein with different tools. We focus on various types of matrices, such as Vandermonde matrices with and without uniform phase distributions, as well as generalized Vandermonde matrices. In each case, we provide explicit expressions of the moments of the associated Gram matrix, as well as more advanced models involving the Vandermonde matrix. Comparisons with classical i.i.d. random matrix theory are provided, and deconvolution results are discussed. We review some applications of the results to the fields of signal processing and wireless communications.

I. INTRODUCTION

Vandermonde matrices have for a long time had a central position in signal processing due to their connections with important tools in the field such as the FFT [START_REF] Burrus | DFT/FFT and Convolution Algorithms[END_REF] or Hadamard [START_REF] Golub | Matrix Computations[END_REF] transforms, to name a few. Vandermonde matrices occur frequently in many applications, such as finance [START_REF] Norberg | On the Vandermonde matrix and its application in mathematical finance[END_REF], signal array processing [START_REF] Schmidt | Multiple emitter localization and signal parameter estimation[END_REF], [START_REF] Wax | Detection of signals by information theoretic criteria[END_REF], [START_REF] Johnson | Array Signal processing: Concepts and Techniques[END_REF], [START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF], [START_REF] Porat | Analysis of the asymptotic relative efficiency of the MUSIC algorithm[END_REF], ARMA processes [START_REF] Klein | On Stein's equation, Vandermonde matrices and Fisher's information matrix of time series processes. part I: The autoregressive moving average process[END_REF], cognitive radio [START_REF] Sampaio | Vandermonde frequency division multiplexing[END_REF], security [START_REF]Vandermonde matrices for security applications[END_REF], wireless communications [START_REF] Wang | Vandermonde-Lagrange mutually orthogonal flexible transceivers for blind CDMA in unknown multipath[END_REF], and biology [START_REF] Waterfall | Sloppymodel universality class and the Vandermonde matrix[END_REF], and have been much studied. The applied research has been somewhat tempered by the fact that very few theoretical results have been available.

A Vandermonde matrix with entries on the unit circle has the following form:

V = 1 √ N      1 • • • 1 e -jω 1 • • • e -jω L . . . . . . . . . e -j(N -1)ω1 • • • e -j(N -1)ω L      (1)
We will consider the case where ω 1 ,...,ω L are i.i.d., taking values in [0, 2π). Throughout the paper, the ω i will be called phase distributions. V will be used only to denote Vandermonde matrices with a given phase distribution, and the dimensions of the Vandermonde matrices will always be N × L.

Known results on Vandermonde matrices are related to the distribution of the determinant [START_REF] Girko | Theory of Random Determinants[END_REF]. The large majority of known results on the eigenvalues of the associated Gram matrix concern Gaussian matrices [START_REF] Mehta | Random Matrices[END_REF] or matrices with independent entries. Very few results are available in the literature on matrices whose structure is strongly related to the Vandermonde case [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF], [START_REF] Bordenave | Eigenvalues of Euclidean random matrices[END_REF]. Known results depend heavily on the distribution of the entries, and do not give any hint on the asymptotic behaviour as the matrices become large. In the realm of wireless channel modeling, [START_REF] Muller | A random matrix model of communication via antenna arrays[END_REF] has provided some insight on the behaviour of the eigenvalues of random Vandermonde matrices for a specific case, without any formal proof.

In many applications, N and L are quite large, and we may be interested in studying the case where both go to ∞ at a given ratio, L N → c. Results in the literature say very little on the asymptotic behaviour of [START_REF] Burrus | DFT/FFT and Convolution Algorithms[END_REF] under this growth condition. The results, however, are well known for other models. The factor 1 √ N , as well as the assumption that the Vandermonde entries e -jω i lie on the unit circle, are included in [START_REF] Burrus | DFT/FFT and Convolution Algorithms[END_REF] to ensure that the analysis will give limiting asymptotic behaviour. Without this assumption, the problem at hand is more involved, since the rows of the Vandermonde matrix with the highest powers would dominate in the calculations of the moments for large matrices, and also grow faster to infinity than the 1 √ N factor in [START_REF] Burrus | DFT/FFT and Convolution Algorithms[END_REF], making asymptotic analysis difficult. In general, often the moments, not the moments of the determinants, are the quantities we seek. Results in the literature say also very little on the moments of Vandermonde matrices (however, see [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF]), and also on the mixed moments of Vandermonde matrices and matrices independent from them. This is in contrast to Gaussian matrices, where exact expressions [START_REF] Thorbjørnsen | Mixed moments of Voiculescu's Gaussian random matrices[END_REF] and their asymptotic behaviour [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF] are known through the concept of freeness [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF], which is central for describing the mixed moments.

The framework and results presented in this paper are reminiscent of similar results concerning i.i.d. random matrices [START_REF] Marchenko | Distribution of Eigenvalues for Some Sets of Random Matrices[END_REF] which have shed light on the design of many important wireless communication problems such as CDMA [START_REF] Shamai | The impact of frequency-flat fading on the spectral efficiency of CDMA[END_REF], MIMO [START_REF] Telatar | Capacity of multi-antenna gaussian channels[END_REF], or OFDM [START_REF] Debbah | MMSE Analysis of Certain Large Isometric Random Precoded Systems[END_REF]. This contribution aims to do the same. We will show that, asymptotically, the moments of the Vandermonde matrices depend only on the ratio c and the phase distribution, and have explicit expressions. The expressions are more involved than what was claimed in [START_REF] Muller | A random matrix model of communication via antenna arrays[END_REF]. Moments are useful for performing deconvolution. Deconvolution for our purposes will mean retrieving "moments" tr L ((D(N )) i ), ..., tr L ((D(N )) i ) (where D(N ) are unknown matrices), from "mixed moments" of D(N ) and matrices on the form [START_REF] Burrus | DFT/FFT and Convolution Algorithms[END_REF], with L and N large. We are only able to perform such deconvolution when the D(N ) are square diagonal matrices independent from V. We will see that such deconvolution can be very useful in many applications, since the retrieved moments can give useful information about the system under study. Deconvolution has previously been handled in cases where V is replaced with a Gaussian matrix [START_REF] Anderson | Asymptotic theory for principal component analysis[END_REF], [START_REF] Abed-Meraim | A subspace algorithm for certain blind identification problems[END_REF], [START_REF] Thorbjørnsen | Mixed moments of Voiculescu's Gaussian random matrices[END_REF], [START_REF] Ryan | Free deconvolution for signal processing applications[END_REF]. As will be seen, the way the phase distribution influences these moments can be split into several cases. Uniform phase distribution plays a central role in that it minimizes the moments. When the phase distribution has a "nice" density (for our purposes this means that the density of the phase distribution is continuous), a nice connection with the moments for uniform phase distribution can be given. When the density of the phase distribution has singularities, for instance when it has point masses, it turns out that the asymptotics of the moments change drastically.

We will also extend our results to generalized Vandermonde matrices, i.e. matrices where the columns do not consist of uniformly distributed powers. Such matrices are important for applications to finance [START_REF] Norberg | On the Vandermonde matrix and its application in mathematical finance[END_REF]. The tools used for standard Vandermonde matrices in this paper will allow us to find the asymptotic behaviour of many generalized Vandermonde matrices as well.

While we provide the full computation of lower order moments, we also describe how the higher order moments can be computed. Tedious evaluation of many integrals is needed for this. It turns out that the first three limit moments coincide with those of the Marchenko Pastur law [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF], [START_REF] Tulino | Random Matrix Theory and Wireless Communications[END_REF]. For higher order moments this is not the case, although we state an interesting inequality involving the Vandermonde limit moments and the moments of the classical Poisson distribution and the Marchenko Pastur law.

The paper is organized as follows. Section II provides background essentials on random matrix theory needed to state the main results. Section III states the main results of the paper. It starts with a general result for the mixed moments of Vandermonde matrices and matrices independent from them. Results for the uniform phase distribution are stated next, both for the asymptotic moments, and the lower order moments. After this, the nice connection between uniform phase distribution and other phase distributions is stated. The case where the density of ω has singularities is then handled. The section ends with results on generalized Vandermonde matrices, and mixed moments of (more than one) independent Vandermonde matrices. Section IV discusses our results and puts them in a general deconvolution perspective, comparing with other deconvolution results, such as those for Gaussian deconvolution. Section V presents some simulations and useful applications showing the implications of the presented results in various applied fields, and discusses the validity of the asymptotic claims in the finite regime. First we apply the presented Vandermonde deconvolution framework to wireless systems, where we estimate the number of paths, the transmissions powers of the users, the number of sources, and what is commonly referred to as wavelength. Finally we apply the results on Vandermonde matrices to the very active field of sparse signal reconstruction.

II. RANDOM MATRIX BACKGROUND ESSENTIALS

In the following, upper (lower boldface) symbols will be used for matrices (column vectors), whereas lower symbols will represent scalar values, (.) T will denote transpose operator, (.) conjugation, and (.) H = (.) T hermitian transpose. I L will represent the L × L identity matrix. We let T r be the (non-normalized) trace for square matrices, defined by,

T r(A) = L i=1 a ii ,
where a ii are the diagonal elements of the L × L matrix A. We also let tr L be the normalized trace, defined by tr L (A) =

L T r(A).

Results in random matrix theory often refer to the empirical eigenvalue distribution of matrices:

Definition 1: With the empirical eigenvalue distribution of an L × L hermitian random matrix T we mean the (random) function

F L T (λ) = #{i|λ i ≤ λ} L , (2)
where λ i are the (random) eigenvalues of T.

In the following, D r (N ), 1 ≤ r ≤ n will denote nonrandom diagonal L × L matrices, where we implicitly assume that L N → c. We will assume that the D r (N ) have a joint limit distribution as N → ∞ in the following sense:

Definition 2: We will say that the {D r (N )} 1≤r≤n have a joint limit distribution as N → ∞ if the limit D i1,...,is = lim

N →∞ tr L (D i1 (N ) • • • D is (N )) (3)
exists for all choices of i 1 , ..., i s ∈ {1, .., n}.

The matrices D i (N ) are assumed to be non-random throughout the paper. However, all presented formulas extend naturally to the case when D i (N ) are random and independent from the Vandermonde matrices. The difference when the D i (N ) are random is that expectations of products of traces also come into play, in the sense that, similarly to expressions on the form (3), expressions of the form

lim N →∞ E[ tr L (D i1 (N ) • • • D is (N )) • • • tr L (D j 1 (N ) • • • D j r (N ))]
also enter the picture. Our framework can also be extended naturally to compute the covariance of traces, defined in the following way: Definition 3: By the covariance C i,j of two traces tr L (A i ) and tr L (A j ) of an L × L random matrix A, we mean the quantity

C i,j (A) = E tr L A i tr L A j -E tr L A i E tr L A j . ( 4
)
When A is replaced with an ensemble of matrices, A L , the limits lim L→∞ LC i,j (A L ) are also called second order moments.

The normalizing factor L is included in order to obtain a limit. It will be explained later why this is the correct normalizing factor for the matrices we consider. The term second order moment is taken from [START_REF] Collins | Second order freeness and fluctuations of random matrices: III. higher order freeness and free cumulants[END_REF], where different matrix ensembles were considered. For these matrices, the second order moments were instead defined as lim L→∞ L 2 C i,j (A L ) (i.e. a higher order normalizing factor was used), since these matrices displayed faster convergence to a limit. We will present expressions for the second order moments lim

N →∞ LC i,j (D(N )V H V).
Most theorems in this paper will present expressions for various mixed moments, defined in the following way:

Definition 4: By a mixed moment we mean the limit

M n = lim N →∞ E[tr L ( D 1 (N )V H VD 2 (N )V H V • • • × D n (N )V H V)], (5)
whenever this exists.

A joint limit distribution of {D r (N )} 1≤r≤n is always assumed in the presented results on mixed moments. Note that when

D 1 (N ) = • • • = D n (N ) = I L , the M n compute
to the asymptotic moments of the Vandermonde matrices themselves, defined by

V n = lim N →∞ E tr L V H V n = lim N →∞ E λ n dF L V H V (λ) . Similarly, when D 1 (N ) = • • • = D n (N ) = D(N ), we will also write D n = lim N →∞ tr L (D(N ) n ). (6)
Note that this is in conflict with the notation D i1,...,is , but the name of the index will resolve such conflicts.

To prove the results of this paper, the random matrix concepts presented up to now need to be extended using concepts from partition theory. We denote by P(n) the set of all partitions of {1, ..., n}, and use ρ as notation for a partition in P(n). Also, we will write ρ = {W 1 , ..., W k }, where W j will be used repeatedly to denote the blocks of ρ, |ρ| = k will denote the number of blocks in ρ, and |W j | will denote the number of elements in a given block. Definition 2 can now be extended as follows.

Definition 5:

For ρ = {W 1 , ..., W k }, with W i = {w i1 , ..., w i|W i | }, we define D W i = D i w i1 ,...,i w i|W i | (7) D ρ = k i=1 D W i . ( 8
)
To better understand the presented expressions for mixed moments, the notion of free cumulants will be helpful. They are defined in terms of noncrossing partitions [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF].

Definition 6: A partition ρ is called noncrossing if, whenever we have i < j < k < l with i ∼ k, j ∼ l (∼ meaning belonging to the same block), we also have i ∼ j ∼ k ∼ l (i.e. i, j, k, l are all in the same block). The set of noncrossing partitions of {1, , , ., n} is denoted N C(n).

The noncrossing partitions have already shown their usefulness in expressing what is called the freeness relation in a particularly nice way [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF].

Definition 7: Assume that A 1 , ..., A n are L × L-random matrices. By the free cumulants of A 1 , ..., A n we mean the unique set of multilinear functionals κ r (r ≥ 1) which satisfy

E [tr L (A i 1 • • • A i n )] = ρ∈N C(n) κ ρ [A i 1 , ..., A i n ] (9)
for all choices of i 1 , ..., i n , where

κ ρ [A i 1 , ..., A i n ] = k j=1 κ W j [A i 1 , ..., A i n ] κ W i [A i 1 , ..., A i n ] = κ |Wi| [A i w i1 , ..., A i w i|W i | ],
where ρ = {W 1 , ..., W k }, with W i = {w i1 , ..., w i|W i | }. By the classical cumulants of A 1 , ..., A n we mean the unique set of multilinear functionals which satisfy [START_REF] Klein | On Stein's equation, Vandermonde matrices and Fisher's information matrix of time series processes. part I: The autoregressive moving average process[END_REF] with N C(n) replaced by the set of all partitions P(n).

We have restricted our definition of cumulants to random matrices, although their general definition is in terms of more general probability spaces (Lecture 11 of [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF]). ( 9) is also called the (free or classical) moment-cumulant formula. The importance of the free moment-cumulant formula comes from the fact that, had we replaced Vandermonde matrices with Gaussian matrices, it could help us perform deconvolution. For this, the cumulants of the Gaussian matrices are needed, which asymptotically have a very nice form. For Vandermonde matrices, it is not known what a useful definition of cumulants would be. However, from the calculations in Appendix A, it will turn out that the following quantities are helpful.

Definition 8: For ρ ∈ P(n), define

K ρ,ω,N = 1 N n+1-|ρ| × (0,2π) |ρ| n k=1 1-e jN (ω b(k-1) -ω b(k) ) 1-e j(ω b(k-1) -ω b(k) ) dω 1 • • • dω |ρ| , (10)
where ω W 1 , ..., ω W |ρ| are i.i.d. (indexed by the blocks of ρ), all with the same distribution as ω, and where b(k) is the block of ρ which contains k (notation is cyclic, i.e. b(0

) = b(n)). If the limit K ρ,ω = lim N →∞ K ρ,ω,N
exists, then it is called a Vandermonde mixed moment expansion coefficient. These quantities do not behave exactly as cumulants, but rather as weights which tell us how a partition in the moment formula we present should be weighted. In this respect our formulas for the moments are different from classical or free moment-cumulant formulas, since these do not perform this weighting. The limits K ρ,ω may not always exist, and necessary and sufficient conditions for their existence seem to be hard to find. However, it is easy to prove from their definition that they do not exist if the density of ω has singularities (for instance when the density has point masses). On the other hand, Theorem 3 will show that they exist when the same density is continuous. P(n) is equipped with the refinement order ≤ [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF], i.e. ρ 1 ≤ ρ 2 if and only if any block of ρ 1 is contained within a block of ρ 2 . The partition with n blocks, denoted 0 n , is the smallest partition within this order, while the partition with 1 block, denoted 1 n , is the largest partition within this order. In the following sections, we will encounter the complementation map of Kreweras (p. 147 of [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF]), which is an order-reversing isomorphism of N C(n) onto itself. To define this we need the circular representation of a partition: We mark n equidistant points 1, ..., n (numbered clockwise) on the circle, and form the convex hull of points lying in the same block of the partition. This gives us a number of convex sets H i , equally many as there are blocks in the partition, which do not intersect if and only if the partition is noncrossing. Put names 1, ..., n on the midpoints of the 1, ..., n (so that ī is the midpoint of the segment from i to i + 1). The complement of the set ∪ i H i is again a union of disjoint convex sets Hi .

Definition 9: The Kreweras complement of ρ, denoted K(ρ), is defined as the partition on { 1, ..., n} determined by i ∼ j in K(ρ) ⇐⇒ ī, j belong to the same convex set Hk .

An important property of the Kreweras complement is that (p. 148 of [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF])

|ρ| + |K(ρ)| = n + 1. (11)

III. STATEMENT OF MAIN RESULTS

We first state the main result of the paper, which applies to Vandermonde matrices with any phase distribution. It restricts to the case when the expansion coefficients K ρ,ω exist. Different versions of it adapted to different Vandermonde matrices will be stated in succeeding sections.

Theorem 1: Assume that the {D r (N )} 1≤r≤n have a joint limit distribution as N → ∞. Assume also that all Vandermonde mixed moment expansion coefficients K ρ,ω exist. Then the limit

M n = lim N →∞ E[tr L ( D 1 (N )V H VD 2 (N )V H V • • • × D n (N )V H V)] (12
) also exists when L N → c, and equals

ρ∈P(n) K ρ,ω c |ρ|-1 D ρ . ( 13
)
The proof of Theorem 1 can be found in Appendix A. Theorem 1 explains how "convolution" with Vandermonde matrices can be performed, and also provides us with an extension of the concept of free convolution to Vandermonde matrices. It also gives us means for performing deconvolution. Indeed, suppose

D 1 (N ) = • • • = D n (N ) = D(N )
, and that one knows all the moments M n . One can then infer on the moments D n by inspecting [START_REF] Waterfall | Sloppymodel universality class and the Vandermonde matrix[END_REF] for increasing values of n. For instance, the first two equations can also be written

D 1 1 = M 1 K 11,ω D 12 = M 2 -cK 0 2 ,ω D 0 2 K 1 2 ,ω ,
where we have used [START_REF] Porat | Analysis of the asymptotic relative efficiency of the MUSIC algorithm[END_REF], that the one-block partition 0 1 = 1 1 is the only partition of length 1, and that the two-block partition 0 2 and the one-block partition 1 2 are the only partitions of length 2. This gives us the first moments D 1 and D 2 defined by [START_REF] Johnson | Array Signal processing: Concepts and Techniques[END_REF], since D 11 = D 1 , D 02 = D 2 1 , and D 12 = D 2 .

A. Uniformly distributed ω

For the case of Vandermonde matrices with uniform phase distribution, it turns out that the noncrossing partitions play a central role. The role is somewhat different than the relation for freeness. Let u denote the uniform distribution on [0, 2π).

Proposition 1: The Vandermonde mixed moment expansion coefficient

K ρ,u = lim N →∞ K ρ,u,N
exists for all ρ. Moreover, 0 < K ρ,u ≤ 1, the K ρ,u are rational numbers for all ρ, and K ρ,u = 1 if and only if ρ is noncrossing.

The proof of Proposition 1 can be found in Appendix B. The same result is proved in [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF], where the K ρ,u are given an equivalent description. The proof in the appendix only translates the result in [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF] to the current notation. Due to Proposition 1, Theorem 1 guarantees that the mixed moments (12) exist in the limit for the uniform phase distribution, and are given by [START_REF] Waterfall | Sloppymodel universality class and the Vandermonde matrix[END_REF]. The K ρ,u are in general hard to compute for higher order ρ with crossings. It turns out that the following computations suffice to obtain the 7 first moments.

Proposition 2:

The following holds:

K {{1,3},{2,4}},u = 2 3 K {{1,4},{2,5},{3,6}},u = 1 2 K {{1,4},{2,6},{3,5}},u = 1 2 K {{1,3,5},{2,4,6}},u = 11 20
K {{1,5},{3,7},{2,4,6}},u = 9 20
K {{1,6},{2,4},{3,5,7}},u = 9 20 .
The proof of Proposition 2 is given in Appendix C. Combining Proposition 1 and Proposition 2 one can prove the following:

Proposition 3: Assume D 1 (N ) = • • • = D n (N ) = D(N ),
and that the limits

m n = (cM ) n = c lim N →∞ E tr L D(N )V H V n ( 14
)
d n = (cD) n = c lim N →∞ tr L (D n (N )) . ( 15
)
exist. When ω = u, we have that

m 1 = d 1 m 2 = d 2 + d 2 1 m 3 = d 3 + 3d 2 d 1 + d 3 1 m 4 = d 4 + 4d 3 d 1 + 8 3 d 2 2 +
d 3 d 2 2 + 175d 3 d 2 d 2 1 + 35d 3 d 4 1 + 77d 3 2 d 1 + 280 3 d 2 2 d 3 1 + 21d 2 d 5 1 + d 7 1 .
Proposition 3 is proved in Appendix D. Several of the following theorems will also be stated in terms of the scaled moments ( 14)- [START_REF] Mehta | Random Matrices[END_REF], rather than M n , D n . The reason for this is that the dependency on the matrix aspect ratio c can be absorbed in m n , d n , so that the result itself can be expressed independently of c, as in the equations of Proposition 3. The same usage of scaled moments has been applied for large Wishart matrices [START_REF] Ryan | Free deconvolution for signal processing applications[END_REF]. Similar computations to those in the proof of Proposition 3 are performed in [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF], although the computations there do not go up as high as the first seven mixed moments. To compute higher order moments, K ρ,u must be computed for partitions of higher order also. The computations performed in Appendix C and D should convince the reader that this can be done, but that it is very tedious.

Following the proof of Proposition 1, we can also obtain formulas for the second order moments of Vandermonde matrices. Since it is easily seen that

C 1,n (D(N )V H V) = C n,1 (D(N )V H V) = 0,
the first nontrivial second order moment is the following:

Proposition 4: Assume that V has uniform phase distribution, let d n be as in [START_REF] Mehta | Random Matrices[END_REF], and define

m i,j = c lim L→∞ LC i,j D(N )V H V) . ( 16
)
Then we have that

m 2,2 = 4 3 d 2 2 . ( 17
)
Proposition 4 is proved in Appendix E, and relies on the same type of calculations as those in Appendix C. Following the proof of Proposition 1 again, we can also obtain exact expressions for moments of lower order random Vandermonde matrices with uniform phase distribution, not only the limit. We state these only for the first four moments.

Theorem 2: Assume D 1 (N ) = D 2 (N ) = • • • = D n (N ), set c = L
N , and define

m (N,L) n = cE tr L D(N )V H V n (18) d (N,L) n = ctr L (D n (N )) . ( 19
)
When ω = u we have that

m (N,L) 1 = d (N,L) 1 m (N,L) 2 = 1 -N -1 d (N,L) 2 + (d (N,L) 1 ) 2 m (N,L) 3 = 1 -3N -1 + 2N -2 d (N,L) 3 +3 1 -N -1 d (N,L) 1 d (N,L) 2 + (d (N,L) 1 ) 3 m (N,L) 4 = 1 - 20 3 N -1 + 12N -2 - 19 3 N -3 d (N,L) 4 + 4 -12N -1 + 8N -2 d (N,L) 3 d (N,L) 1 + 8 3 -6N -1 + 10 3 N -2 (d (N,L) 2 ) 2 +6 1 -N -1 d (N,L) 2 (d (N,L) 1 ) 2 + (d (N,L) 1 ) 4 .
Theorem 2 is proved in Appendix F. Exact formulas for the higher order moments also exist, but they become increasingly complex, as higher order terms N -k also enter the picture. These formulas are also harder to prove for higher order moments. In many cases, exact expressions are not what we need: first order approximations (i.e. expressions where only the N -1 -terms are included) can suffice for many purposes. In Appendix F, we explain how the simpler case of these first order approximations can be computed. It seems much harder to prove a similar result when the phase distribution is not uniform.

An important result building on the results we present is the following, which provides a major difference from the limit eigenvalue distributions of Gaussian matrices.

Proposition 5: The asymptotic mean eigenvalue distribution of a Vandermonde matrix with uniform phase distribution has unbounded support.

Proposition 5 is proved in Appendix G.

B. ω with continuous density

The following result tells us that the limit K ρ,ω exists for many ω, and also gives a useful expression for them in terms of K ρ,u and the density of ω.

Theorem 3: The Vandermonde mixed moment expansion coefficients K ρ,ω = lim N →∞ K ρ,ω,N exist whenever the density p ω of ω is continuous on [0, 2π). If this is fulfilled, then

K ρ,ω = K ρ,u (2π) |ρ|-1 2π 0 p ω (x) |ρ| dx . ( 20
)
The proof is given in Appendix H. Although the proof assumes a continuous density, we remark that it can be generalized to cases where the density contains a finite set of jump discontinuities also. In Section V, several examples are provided where the integrals (20) are computed. An important consequence of Theorem 3 is the following, which gives the uniform phase distribution an important role.

Proposition 6: Let V ω denote a Vandermonde matrix with phase distribution ω, and set

V ω,n = lim N →∞ E tr L V H ω V ω n .
Then we have that

V u,n ≤ V ω,n .
The proof is given in Appendix I. An immediate consequence of this and Proposition 5 is that all phase distributions, not only uniform phase distribution, give Vandermonde matrices with unbounded mean eigenvalue distributions in the limit. Besides providing us with a deconvolution method for finding the mixed moments of the {D r (N )} 1≤r≤n , Theorem 3 also provides us with a way of inspecting the phase distribution ω, by first finding the moments of the density, i.e. 2π 0 p ω (x) k dx. However, note that we can not expect to find the density of ω itself, only the density of the density of ω. This follows immediately by noting that 2π 0 p ω (x) k dx remains unchanged when the phase distribution ω is cyclically shifted.

C. ω with density singularities

The asymptotics of Vandermonde matrices are different when the density of ω has singularities, and depends on the density growth rates near the singular points. It will be clear from the following results that one can not perform deconvolution for such ω to obtain the higher order moments of the {D r (N )} 1≤r≤n , as only their first moment can be obtained. The asymptotics are first described for ω with atomic density singularities, as this is the simplest case to prove. After this, densities with polynomic growth rates near the singularities are addressed.

Theorem 4: Assume that p ω = r i=1 p i δ α i is atomic (where δ αi (x) is dirac measure (point mass) at α i ), and denote by p

(n) = r i=1 p n i . Then lim N →∞ E[T r( D 1 (N ) 1 N V H VD 2 (N ) 1 N V H V • • • × D n (N ) 1 N V H V)] = c n-1 p (n) lim N →∞ n i=1 tr L (D i (N )) .
Note here that the non-normalized trace is used. The proof can be found in Appendix J. In particular, Theorem 4 states that the asymptotic moments of 1 N V H V can be computed from p (n) . The theorem is of great importance for the estimation of the point masses p i . In blind seismic and telecommunication applications, one would like to detect the locations α i . Unfortunately, Theorem 4 tells us that this is impossible with our deconvolution framework, since the p (n) , which are the quantities we can find through deconvolution, have no dependency to them. This parallels Theorem 3, since also there we could not recover the density p ω itself. Having found the p (n) through deconvolution, one can find the point masses p i , by solving for

p 1 , p 2 , ... in the Vandermonde equation    p 1 p 2 • • • p r p 2 1 p 2 2 • • • p 2 r . . . . . . . . . . . .       1 1 . . .    =    p (1)
p (2) . . .

   .
The case when the density has non-atomic singularities is more complicated. We provide only the following result, which addresses the case when the density has polynomic growth rate near the singularities.

Theorem 5: Assume that

lim x→α i |x -α i | s p ω (x) = p i for some 0 < s < 1
for a set of points α 1 , ..., α r , with p ω continuous for ω = α 1 , ..., α r . Then

lim N →∞ E[T r( D 1 (N ) 1 N s V H VD 2 (N ) 1 N s V H V • • • × D n (N ) 1 N s V H V)] = c n-1 q (n) lim N →∞ n i=1 tr L (D i (N ))
where

q (n) = 2Γ(1 -s) cos (1-s)π 2 n p (n) × [0,1] n n k=1 1 |x k-1 -x k | 1-s dx 1 • • • dx n , (21)
and p (n) = i p n i . Note here that the non-normalized trace is used.

The proof can be found in Appendix K. Also in this case it is only the point masses p i which can be found through deconvolution, not the locations α i . Note that the integral in [START_REF] Marchenko | Distribution of Eigenvalues for Some Sets of Random Matrices[END_REF] can also be written as an m-fold convolution. Similarly, the definition of K ρ,ω,N given by [START_REF] Sampaio | Vandermonde frequency division multiplexing[END_REF] can also be viewed as a 2-fold convolution when ρ has two blocks, and as a 3-fold convolution when ρ has three blocks (but not for ρ with more than 3 blocks).

A useful application of Theorem 5 we will return to is the case when ω = k sin(θ) for some constant k (see [START_REF] Cardoso | Cooperative Spectrum Sensing Using Random Matrix Theory[END_REF]), with θ uniformly distributed on some interval. This case is simulated in Section V-A. It is apparent from [START_REF] Haykin | Cognitive Radio: Brain Empowered Wireless Communications[END_REF] that the density goes to infinity near ω = ±k, with rate x -1/2 . Theorem 5 thus applies with s = 1/2. For this case, however, the "edges" at ±π/2 are never reached in practice. Indeed, in array processing [START_REF] Krim | Two decades of array signal processing research[END_REF], the antenna array is a sector antenna which scans an angle interval which never includes the edges. We can therefore restrict ω in our analysis to clusters of intervals [α i , β i ] not containing ±1, for which the results of Section III-B suffice. In this way, we also avoid the computation of the cumbersome integral [START_REF] Marchenko | Distribution of Eigenvalues for Some Sets of Random Matrices[END_REF].

D. Generalized Vandermonde matrices

We will consider generalized Vandermonde matrices on the form

V = 1 √ N      e -j N f (0) ω 1 • • • e -j N f (0) ω L e -j N f ( 1 N ) ω 1 • • • e -j N f ( 1 N ) ω L . . . . . . . . . e -j N f ( N -1 N ) ω1 • • • e -j N f ( N -1 N ) ω L      , ( 22
)
where f is called the power distribution, and is a function from [0, 1) to [0, 1). We will also consider the more general case when f is replaced with a random variable λ, i.e.

V = 1 √ N      e -jN λ 1 ω 1 • • • e -jN λ 1 ω L e -jN λ2ω1 • • • e -jN λ2ω L . . . . . . . . . e -jN λ N ω 1 • • • e -jN λ N ω L      , ( 23
)
with the λ i i.i.d. and distributed as λ, defined and taking values in [0, 1), and also independent from the ω j . We will define mixed moment expansion coefficients for generalized Vandermonde matrices also. The difference is that, while we in Definition 8 simplified using the geometric sum formula, we can not do this now since we do not assume uniform power distribution anymore. To define expansion coefficients for generalized Vandermonde matrices of the form [START_REF] Shamai | The impact of frequency-flat fading on the spectral efficiency of CDMA[END_REF], define first integer functions

f N from [0, N -1] to [0, N -1] by f N (r) = N f r N
. Let p f N be the corresponding density for f N . The procedure is similar for matrices of the form [START_REF] Telatar | Capacity of multi-antenna gaussian channels[END_REF]. The following definition captures both cases:

Definition 10: For ( 22) and ( 23), define

K ρ,ω,f,N = 1 N 1-|ρ| × (0,2π) |ρ| n k=1 N -1 r=0 p f N (r)e jr(ω b(k-1) -ω b(k) ) dω 1 • • • dω |ρ| K ρ,ω,λ,N = 1 N 1-|ρ| × (0,2π) |ρ| n k=1 1 0 N e jN λ(ω b(k-1) -ω b(k) ) dλ dω 1 • • • dω |ρ| , ( 24
)
where ω W1 , ..., ω W |ρ| are as in Definition 8. If the limits

K ρ,ω,f = lim N →∞ K ρ,ω,f,N K ρ,ω,λ = lim N →∞ K ρ,ω,λ,N ,
exist, then they are called Vandermonde mixed moment expansion coefficients.

Note that (1) corresponds to [START_REF] Shamai | The impact of frequency-flat fading on the spectral efficiency of CDMA[END_REF] with f (x) = x. The following result holds:

Theorem 6: Theorem 1 holds also with Vandermonde matrices (1) replaced with generalized Vandermonde matrices on either form [START_REF] Shamai | The impact of frequency-flat fading on the spectral efficiency of CDMA[END_REF] or [START_REF] Telatar | Capacity of multi-antenna gaussian channels[END_REF], and with K ρ,ω replaced with either

K ρ,ω,f or K ρ,ω,λ .
The proof follows the same lines as those in Appendix A, and is therefore only explained briefly at the end of that appendix. As for matrices of the form (1), it is the case of uniform phase distribution which is most easily described how to compute for generalized Vandermonde matrices also. Appendix B shows how the computation of K ρ,u boils down to computing certain integrals. The same comments are valid for matrices of the form [START_REF] Shamai | The impact of frequency-flat fading on the spectral efficiency of CDMA[END_REF] or [START_REF] Telatar | Capacity of multi-antenna gaussian channels[END_REF] in order to compute K ρ,ω,f and K ρ,ω,λ . This is further commented at the end of that appendix.

We will not consider generalized Vandermonde matrices with density singularities.

E. The joint distribution of independent Vandermonde matrices

When many independent random Vandermonde matrices are involved, the following holds:

Theorem 7: Assume that the {D r (N )} 1≤r≤n have a joint limit distribution as N → ∞. Assume also that V 1 , V 2 , ... are independent Vandermonde matrices with the same phase distribution ω, and that the density of ω is continuous. Then the limit

lim N →∞ E[tr L ( D 1 (N )V H i 1 V i 2 D 2 (N )V H i 2 V i 3 • • • × D n (N )V H i n V i 1 )
] also exists when L N → c. The limit is 0 when n is odd, and equals

ρ≤σ∈P(n) K ρ,ω c |ρ|-1 D ρ , ( 25
)
where σ = {σ 1 , σ 2 } = {{1, 3, 5, ..., }, {2, 4, 6, ...}} is the partition where the two blocks are the even numbers, and the odd numbers.

The proof of Theorem 7 can be found in Appendix L. That appendix also contains some remarks on the case when the matrices D i (N ) are placed at different positions relative to the Vandermonde matrices. From Theorem 7, the following corollary is immediate:

Corollary 1: The first three mixed moments

V (2) n = lim N →∞ E tr L V H 1 V 2 V H 2 V 1 n of independent Vandermonde matrices V 1 , V 2 are given by V (2) 1 = I 2 V (2) 2 = 2 3 I 2 + 2I 3 + I 4 V (2) 3 = 11 20 I 2 + 4I 3 + 9I 4 + 6I 5 + I 6 ,
where

I k = (2π) k-1 2π 0 p ω (x) k dx .
In particular, when the phase distribution is uniform, the first three mixed moments are given by

V (2) 1 = 1 V (2) 2 = 11 3
V (2) 3 = 411 20
The results here can also be extended to the case with independent Vandermonde matrices with different phase distributions:

Theorem 8: Assume that {V i } 1≤i≤s are independent Vandermonde matrices, where V i has continuous phase distribution ω i . Denote by p ω i the density of ω i . Then Equation ( 25) still holds, with K ρ,ω replaced by

K ρ,u (2π) |ρ|-1 2π 0 s i=1 p ωi (x) |ρi| dx,
where ρ i consists of all numbers k such that i k = i.

The proof is omitted, as it is a straightforward extension of the proofs of Theorem 3 and Theorem 7.

IV. DISCUSSION

In the recent work [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF], the Vandermonde model ( 1) is encountered in reconstruction of multidimensional signals in wireless sensor networks. The authors also recognize a similar expression for the Vandermonde mixed moment expansion coefficient as in Definition 8. They also state that, for the case of uniform phase distribution, closed form expressions for the moments can be found, building on an analysis of partitions and calculation of volumes of convex polytopes described by certain constraints. This is very similar to what is done in this paper.

We will in the following discuss some differences and similarities between Gaussian and Vandermonde matrices.

A. Convergence rates

In [START_REF] Thorbjørnsen | Mixed moments of Voiculescu's Gaussian random matrices[END_REF], almost sure convergence of Gaussian matrices was shown by proving exact formulas for the distribution of lower order Gaussian matrices. These deviated from their limits by terms of order 1/N 2 . In Theorem 2, we see that terms of order 1/N are involved. This slower rate of convergence may not be enough to make a statement on whether we have almost sure convergence for Vandermonde matrices. However, [START_REF] Bryc | Spectral measure of large random Hankel, Markov and Toeplitz matrices[END_REF] shows some almost sure convergence properties for certain Hankel and Toeplitz matrices. These matrices are seen in that paper to have similar combinatorial descriptions for the moments, when compared to Vandermonde matrices in this paper. Therefore, it may be the case that the techniques in [START_REF] Bryc | Spectral measure of large random Hankel, Markov and Toeplitz matrices[END_REF] can be generalized to address almost sure convergence of Vandermonde matrices also. Figure 1 shows the speed of convergence of the moments of Vandermonde matrices (with uniform phase distribution) towards the asymptotic moments as the matrix dimensions grow, and as the number of samples grow. The differences between the asymptotic moments and the exact moments are also shown. To be more precise, the MSE values in Figure 1 are computed as follows:

1) K samples V i are independently generated using (1).

2) The 4 first sample moments vji =

1 L tr n V H i V i j (1 ≤ j ≤ 4)
4 j=1 (E j -A j ) 2 .
Figure 1 is in sharp contrast with Gaussian matrices, as shown in Figure 2. First of all, it is seen that the asymptotic moments can be used just as well instead of the exact moments (for which expressions can be found in [START_REF] Ryan | Channel capacity estimation using free probability theory[END_REF]), due to the O(1/N 2 ) convergence of the moments. Secondly, it is seen that only 5 samples were needed to get a reliable estimate for the moments.

B. Inequalities between moments of Vandermonde matrices and moments of known distributions

We will state an inequality involving the moments of Vandermonde matrices, and the moments of known distributions. The classical Poisson distribution with rate λ and jump size α is defined as the limit of

1 - λ n δ 0 + λ n δ α * n
as n → ∞ [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF], where * denotes classical (additive) convolution, and * n denotes n-fold convolution with itself. For our analysis, we will only need the classical Poisson distribution with rate c and jump size 1, denoted ν c . The free Poisson distribution with rate λ and jump size α is defined similarly N XX H with X a complex standard Gaussian matrix. The MSE of the first 4 exact moments from the asymptotic moments is also shown.

as the limit of

1 - λ n δ 0 + λ n δ α n
as n → ∞, where is the free probability counterpart of * [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF], [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF], and where n denotes n-fold free convolution with itself. For our analysis, we will only need the free Poisson distribution with rate 1 c and jump size c, denoted µ c . µ c is the same as the better known Marchenko Pastur law, i.e. it has the density [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF]

f µc (x) = (1 - 1 c ) + δ 0 (x) + (x -a) + (b -x) + 2πcx , ( 26
)
where

(z) + = max(0, z), a = (1- √ c) 2 , b = (1+ √ c) 2 .
c , then φ(a n 1 ) = ρ∈N C(n) c n-|ρ| = ρ∈N C(n) c |K(ρ)|-1 = ρ∈N C(n) c |ρ|-1 . ( 27
)
Here we have used the Kreweras complementation map and (11), with φ denoting the expectation in a non-commutative probability space [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF]. Also, if a 2 has the distribution ν c , then

E(a n 2 ) = ρ∈P(n) c |ρ| . ( 28
)
We immediately recognize the c |ρ|-1 -entry of Theorem 1 in ( 27) and ( 28) (with an additional power of c in ( 28)). Combining Proposition 1 with [START_REF] Ryan | Free deconvolution for signal processing applications[END_REF], and (28), we thus get the following corollary to Proposition 1: Corollary 2: Assume that V has uniform phase distribution. Then the limit moment

D 1 (N ) = • • • = D n (N ) = I L ,
V n = lim N →∞ E tr L V H V n satisfies the inequality φ(a n 1 ) ≤ V n ≤ 1 c E(a n 2 ),
where a 1 has the distribution µ c of the Marchenko Pastur law, and a 2 has the Poisson distribution ν c . In particular, equality occurs for m = 1, 2, 3 and c = 1 (since all partitions are noncrossing for m = 1, 2, 3). Corollary 2 thus states that the moments of Vandermonde matrices with uniform phase distribution are bounded above and below by the moments of the classical and free Poisson distributions, respectively. The left part of the inequality in Corollary 2 was also observed in Section VI in [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF]. The different Poisson distributions enter here because their (free and classical) cumulants resemble the c |ρ|-1 -entry in Theorem 1, where we also can use that K ρ,u = 1 if and only if ρ is noncrossing to get a connection with the Marchenko Pastur law. To see how close the asymptotic Vandermonde moments are to these upper and lower bounds, the following corollary to Proposition 3 contains the first moments:

Corollary 3: When c = 1, the limit moments

V n = lim N →∞ E tr L V H V n ,
the moments f p n of the Marchenko Pastur law µ 1 , and the moments p n of the Poisson distribution ν 1 satisfy

f p 4 = 14 ≤ V 4 = 44 3 ≈ 14.67 ≤ p 4 = 15 f p 5 = 42 ≤ V 5 = 146 3 ≈ 48.67 ≤ p 5 = 52 f p 6 = 132 ≤ V 6 = 3571 20 ≈ 178.55 ≤ p 6 = 203 f p 7 = 429 ≤ V 7 = 2141
3 ≈ 713.67 ≤ p 7 = 877. The first three moments coincide for the three distributions, and are 1, 2, and 5, respectively.

The numbers f p n and p n are simply the number of partitions in N C(n) and P(n), respectively. The number of partitions in N C(n) equals the Catalan number

C n = 1 n+1
2n n [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF], and are easily computed. The number of partitions of P(n) are also known as the Bell numbers B n [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF]. They can easily be computed from the recurrence relation

B n+1 = n k=0 B k n k .
In Figure 3, the mean eigenvalue distribution of 640 samples of a 1600 × 1200 (i.e. c = 0.75) Vandermonde matrix with uniform phase distribution is shown. While the Poisson distribution ν 1 is purely atomic and has masses at 0, 1, 2, and 3 which are e -1 , e -1 , e -1 /2, and e -1 /6 (the atoms consist of all integer multiples), the Vandermonde histogram shows a more continuous eigenvalue distribution, with the peaks which the Poisson distribution has at integer multiples clearly visible, although not as sharp. We remark that the support of V H V for a fixed N goes all the way up to N , but lies within [0, N ]. It is unknown whether the peaks at integer multiples in the Vandermonde histogram grow to infinity as we let N → ∞. From the histogram, only the peak at 0 seems to be of atomic nature. The effect of decreasing c amounts to stretching the eigenvalue density vertically, and compressing it horizontally, just as the case for the different Marchenko Pastur laws. An eigenvalue histogram for Gaussian matrices which in the limit give the corresponding (in the sense of Corollary 2) Marchenko Pastur law for Figure 3 (i.e. µ 0.75 ) is shown in Figure 4. Figure 5 shows an eigenvalue histogram in the case of a non-uniform phase distribution. Here we have taken 640 samples of a 1600×1200 Vandermonde matrix with phase distribution with density [START_REF] Haykin | Cognitive Radio: Brain Empowered Wireless Communications[END_REF], with λ = 2d, α = π 4 . This density, also shown in Figure 6, is used in the applications of Section V-A. Experiments show that the eigenvalue histogram tends to flatten when the phase distribution becomes "less uniform", with a higher concentration of larger eigenvalues.

It is unknown whether the inequalities for the moments can be extended to inequalities for the associated capacity. If X is an N ×N standard, complex, Gaussian matrix, then an explicit expression for the asymptotic capacity exists [START_REF] Tulino | Random Matrix Theory and Wireless Communications[END_REF]:

lim N →∞ 1 N log 2 det I N + ρ 1 N XX H = 2 log 2 1 + ρ -1 4 √ 4ρ + 1 -1 2 -log 2 e 4ρ √ 4ρ + 1 -1 2 . ( 29
)
In Figure 7(a), several realizations of the capacity are computed for Gaussian matrix samples of size 36 × 36. The asymptotic capacity ( 29) is also shown. In Figure 7(b), several realizations of the capacity are computed for Vandermonde matrix samples of the same size, for the case of uniform phase distribution. It is seen that the variance of the Vandermonde capacities is higher than for the Gaussian counterparts. This should come as no surprise, due to the slower convergence to the asymptotic limits for Vandermonde matrices. Although the capacities of Vandermonde matrices with uniform phase N log 2 det I N + ρ 1 N XX H when X is standard, complex, Gaussian. The asymptotic capacity ( 29) is also shown. distribution and Gaussian matrices seem to be close, we have no proof that the capacities of Vandermonde matrices are even finite due to the unboundedness of its support.

C. Deconvolution

Deconvolution with Vandermonde matrices (as stated in [START_REF] Waterfall | Sloppymodel universality class and the Vandermonde matrix[END_REF] in Theorem 1) differs from the Gaussian deconvolution counterpart [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF] in the sense that there is no multiplicative [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF] structure involved, since K ρ,ω is not multiplicative in ρ. The Gaussian equivalent of Proposition 3 (i.e. V H V replaced with

1 N XX H , with X an L × N complex, standard, Gaussian matrix) is m 1 = d 1 (30) m 2 = d 2 + d 2 1 (31) m 3 = d 3 + 3d 2 d 1 + d 3 1 (32) m 4 = d 4 + 4d 3 d 1 + 2d 2 2 + 6d 2 d 2 1 + d 4 1 ( 33
)
m 5 = d 5 + 5d 4 d 1 + 5d 3 d 2 + 10d 3 d 2 1 + 10d 2 2 d 1 + 10d 2 d 3 1 + d 5 1 (34) m 6 = d 6 + 6d 5 d 1 + 6d 4 d 2 + 15d 4 d 2 1 + 3d 2 3 + 30d 3 d 2 d 1 + 20d 3 d 3 1 + 5d 3 2 + 10d 2 2 d 2 1 + 15d 2 d 4 1 + d 6 1 (35) m 7 = d 7 + 7d 6 d 1 + 7d 5 d 2 + 21d 5 d 2 1 + 7d 4 d 3 + 42d 4 d 2 d 1 + 35d 4 d 3 1 + 21d 2 3 d 1 + 21d 3 d 2 2 + 105d 3 d 2 d 2 1 + 35d 3 d 4 1 + 35d 3 2 d 1 + 70d 2 2 d 3 1 + 21d 2 d 5 1 + d 7 1 , ( 36
)
where the m i and the d i are computed as in ( 14)-( 15). This follows immediately from asymptotic freeness [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF], and from the fact that 1 N XX H converges to the Marchenko Pastur law µ c . In particular, when all D i (N ) = I L and c = 1, we obtain the limit moments 1, 2, 5, 14, 42, 132, 429, which also were listed in Corollary 3. One can also write down Gaussian equivalents to the second order moments of Vandermonde matrices (17) using techniques from [START_REF] Collins | Second order freeness and fluctuations of random matrices: III. higher order freeness and free cumulants[END_REF]. However the formulas look quite different, and the asymptotic behaviour is different. We have for instance

lim L→∞ L 2 D 1,1 D(N ) 1 N XX H = cd 2 , ( 37
)
where it is not needed that the matrices D(N ) are diagonal.

Similarly, one can write down an equivalent to Theorem 2 for the exact moments. For the first three moments (the fourth moment is dropped, since this is more involved), these are

m 1 = d 1 m 2 = d 2 + d 2 1 m 3 = 1 + N -2 d 3 + 3d 1 d 2 + d 3 1 .
This follows from a careful count of all possibilities after the matrices have been multiplied together (see also [START_REF] Ryan | Channel capacity estimation using free probability theory[END_REF], where one can see that the restriction that the matrices D i (N ) are diagonal can be dropped in the Gaussian case). It is seen, contrary to Theorem 2 for Vandermonde matrices, that the second exact moment equals the second asymptotic moment [START_REF] Krim | Two decades of array signal processing research[END_REF], and also that the convergence is faster (i.e. O(N -2 )) for the third moment (this will also be the case for higher moments).

The two types of (de)convolution also differ in how they can be computed in practice. In [START_REF] Ryan | Free deconvolution for signal processing applications[END_REF], an algorithm for free convolution with the Marchenko Pastur law was sketched. A similar algorithm may not exist for Vandermonde convolution. However, Vandermonde convolution can be subject to numerical approximation: To see this, note first that Theorem 3 splits the numerics into two parts: The approximation of the integrals p ω (x) |ρ| dx, and the approximation of the K ρ,u . A strategy for obtaining the latter quantities could be to randomly generate many numbers between 0 and 1 and estimate the volume as the ratio of the solutions which satisfy (63) in Appendix B. Implementations of the various Vandermonde convolution variants given in this paper can be found in [START_REF] Ryan | Tools for convolution with Vandermonde matrices[END_REF].

In practice, one often has a random matrix model where independent Gaussian and Vandermonde matrices are both present. In such cases, it is possible to combine the individual results for both of them. In Section V, examples on how this can be done are presented.

V. APPLICATIONS

The applications presented here all use the deconvolution framework for Vandermonde matrices. Since additive, white, Gaussian noise also is taken into account, Vandermonde deconvolution is combined with Gaussian deconvolution. Matlab code for running the different simulations can be found in [START_REF] Ryan | Tools for convolution with Vandermonde matrices[END_REF].

In the eigenvalue histograms for Vandermonde matrices shown in figures 3 and 5, large matrices were used in order to obtain something close to the asymptotic limit. In practical scenarios, and in the applications we present, N and L are much smaller than what was used in these figures, which partially explains the uncertainty in some of the simulations. In particular, the uncertainty for non-uniform phase distributions such as those in Section V-A is high, since exact expressions for the lower order moments are not known, contrary to the case of uniform phase distribution. In all the following, d is the distance between the antennas whereas λ is the wavelength. The ratio d λ is a figure of the resolution with which the system will be able to separate (and therefore estimate the position of) users in space.

A. Detection of the number of sources

Let us consider a basestation equipped with N receiving antennas, and with L mobiles (each with a single antenna) in the cell. The received signal at the base station is given by

r i = VP 1 2 s i + n i . ( 38
)
Here r i is the N × 1 received vector, s i is the L × 1 transmit vector by the L users which is assumed to satisfy

E s i s H i = I L , n i is N × 1 additive, white, Gaussian noise of variance σ √
N (all components in s i and n i are assumed independent). In the case of a line of sight between the users and the base station, and considering a Uniform Linear Array (ULA), the matrix V has the following form:

V = 1 √ N      1 • • • 1 e -j2π d λ sin(θ1) • • • e -j2π d λ sin(θ L ) . . . . . . . . . e -j2π(N -1) d λ sin(θ1) • • • e -j2π d λ sin(θ L )      (39
) Here, θ i is the angle of the user in the cell and is supposed to be uniformly distributed over [-α, α]. P 1 2 is an L × L diagonal power matrix due to the different distances from which the users emit. In other words, we assume that the phase distribution has the form 2π d λ sin(θ) with θ uniformly distributed on [-α, α]. The fact that the phase has the form 2π d λ sin(θ) is a well known result in array processing [START_REF] Krim | Two decades of array signal processing research[END_REF]. The user's distribution can be known (in the case of these simulations, the uniform distribution has been accounted for without loss of generality) through measurements in wireless systems up to some parameters (here, α typically). This is usually done to have a better understanding of the user's behaviour. It is easily seen, by taking inverse functions, that the density is, when 2d sin α λ < 1,

p ω (x) = 1 2α 4π 2 d 2 λ 2 -x 2 (40)
on [-2πd sin α λ , 2πd sin α λ ], and 0 elsewhere (see Figure 6). Throughout the paper we will assume, as in Figure 5, that λ = 2d, α = π 4 when model ( 39) is used. With this assumption, 2d sin α λ < 1 is always fulfilled. The goal is to detect the number of sources L and their respective power based on the sample covariance matrix supposing that we have K observations, of the same order as N . When the number of observation is quite higher than N (and the noise variance is known), classical subspace methods [START_REF] Comon | Independent Component Analysis, a new concept?[END_REF] provide tools to detect the number of sources. Indeed, let R be the true covariance matrix given by

VPV H + σ 2 I N ,
where σ 2 is the noise variance. This matrix has N -L eigenvalues equal to σ 2 and L eigenvalues strictly superior to σ 2 . One can therefore determine the number of source by counting the number of eigenvalues different from σ 2 . However, in practice, one has only access to the sample covariance matrix given by

W = 1 K YY H , with Y = [r 1 , ...r K ] = VP 1 2 [s 1 , ..., s K ] + [n 1 , ..., n K ]. (41)
If one has only the sample covariance matrix W, we have three independent parts which must be dealt with in order to get an estimate of P: the Gaussian matrices S = [s 1 , ..., s K ] and N = [n 1 , ..., n K ], and the Vandermonde matrix V. It should thus be possible to combine Gaussian deconvolution [START_REF] Ryan | Channel capacity estimation using free probability theory[END_REF] and Vandermonde deconvolution by performing the following steps:

1) Estimate the moments of 1 K VP 1 2 SS H P 1 2 V H using multiplicative free convolution as described in [START_REF] Ryan | Free deconvolution for signal processing applications[END_REF]. This is the denoising part.

2) Estimate the moments of PV H V, again using multiplicative free deconvolution. 3) Estimate the moments of P using Vandermonde deconvolution as described in this paper.

Putting these steps together, we will prove the following: Proposition 7: Define

I n = (2π) n-1 2π 0 p ω (x) n dx, ( 42
)
and denote the moments of P and the sample covariance matrix, respectively, by

P i = tr L (P i ) W i = tr N (W i ).
Then the equations

W 1 = c 2 P 1 + σ 2 W 2 = c 2 P 2 + (c 2 2 I 2 + c 2 c 3 )(P 1 ) 2 +2σ 2 (c 2 + c 3 )P 1 + σ 4 (1 + c 1 ) W 3 = c 2 P 3 + (3c 2 2 I 2 + 3c 2 c 3 )P 1 P 2 + c 3 2 I 3 + 3c 2 2 c 3 I 2 + c 2 c 2 3 (P 1 ) 3 +3σ 2 (1 + c 1 )c 2 P 2 +3σ 2 ((1 + c 1 )c 2 2 I 2 + c 3 (c 3 + 2c 2 ))(P 1 ) 2 +3σ 4 (c 2 1 + 3c 1 + 1)c 2 P 1 +σ 6 (c 2 1 + 3c 1 + 1
) provide an asymptotically unbiased estimator for the moments P i from the moments of W i (or vice versa) when

lim N →∞ N K = c 1 , lim N →∞ L N = c 2 , lim N →∞ L K = c 3 .
The proof of this can be found in Appendix M. Note that c 3 = c 1 c 2 , so that the definition of c 3 is really not necessary. We still include it however, since c 1 , c 2 and c 3 are matrix aspect ratios which represent different deconvolution stages, so that they all are used when these stages are implemented and combined serially. In the simulations, Proposition 7 is put to the test when P has three sets of powers, 0.5, 1, and 1.5, with equal probability, with phase distribution given by [START_REF] Cardoso | Cooperative Spectrum Sensing Using Random Matrix Theory[END_REF]. Both the number of sources and the powers are estimated. For the phase distribution [START_REF] Cardoso | Cooperative Spectrum Sensing Using Random Matrix Theory[END_REF], the integrals I 2 and I 3 can be computed exactly (for general phase distributions they are computed numerically), and are [START_REF] Rottmann | Mathematische Formelsammlung. B.I. Wissenschaftsverlag[END_REF]

I 2 = λ 4dα 2 ln 1 + sin α 1 -sin α I 3 = λ 2 tan α 4d 2 α 3 .
Under the assumptions λ = 2d, α = π 4 used throughout this paper, the integrals above take the values

I 2 = 40 π 2 ln 2 + √ 2 2 - √ 2 I 3 = 1600 π 3 .
For estimation of the powers, knowing that we have only three sets of powers with equal probability, it suffices to estimate the three lowest moments in order to get an estimate of the powers (which are the three distinct eigenvalues of P). Therefore, in the following simulations, Proposition 7 is first used to get an estimate of the moments of P. Then these are used to obtain an estimate of the three distinct eigenvalues of P using the Newton-Girard formulas [START_REF] Seroul | Programming for Mathematicians[END_REF]. These should then lie close to the three powers of P. Power estimation for the model ( 39) is shown in the first plot of Figure 10. In the plot, K = L = N = 144, and σ = √ 0.1. Experiments show that when the phase distribution becomes "less" uniform, larger matrix sizes are needed in order for accurate power estimation using this method. This will also be seen when we perform power estimation using uniform phase distribution in the next section.

For estimation of the number of users L, we assume that the power distribution of P is known, but not L itself. Since L is unknown, in the simulations we enter different candidate values of it into the following procedure:

1) Computing the moments P i = tr L (P i ) of P.

2) The moments tr L (P i ) are fed into the formulas of Proposition 7, and we thus obtain candidate moments W i of the sample covariance matrix W. 3) Compute the sum of the square errors between these candidate moments, and the moments Ŵi of the observed sample covariance matrix Ŵ, i.e. compute

The estimate L for the number of users is chosen as the one which gives the minimum value for the sum of square errors after these steps.

In Figure 8, we have set σ = √ 0.1, N = 100, and L = 36.

We tried the procedure described above for 1 all the way up to 100 observations. It is seen that only a small number of observations are needed in order to get an accurate estimate of L. When K = 1, it is seen that more observations are needed to get an accurate estimate of L, when compared to K = 10.

B. Estimation of the number of paths

In many channel modeling applications, one needs to determine the number of paths in the channel [START_REF] Rappaport | Wireless Communications, Principles and Practice[END_REF]. For this purpose, consider a multi-path channel of the form:

h(τ ) = L i=1 s i δ(τ -τ i )
Here, s i are i.d. Gaussian random variables with power P i and τ i are uniformly distributed delays over [0, T ]. The s i represent the attenuation factors due to the different reflections. L is the total number of paths. In the frequency domain, the channel is given by

H(f ) = L i=1 s i G(f )e -j2πf τ i .
Sampling the continuous frequency signal at f i = i W N where W is the bandwidth, the model becomes (for a given channel realization)

H = VP 1 2 s where V = 1 √ N       1 • • • 1 e -j2π W τ 1 N • • • e -j2π W τ L N . . . . . . . . . e -j2π(N -1) W τ 1 N • • • e -j2π(N -1) W τ L N       , (43)
We will here set W = T = 1, which means that the ω i of (1) are uniformly distributed over [0, 2π). The corresponding eigenvalue histogram was shown in Figure 3. When additive noise (n) again is taken into consideration, our model again becomes that of [START_REF] Rappaport | Wireless Communications, Principles and Practice[END_REF], the only difference being that the phase distribution of the Vandermonde matrix now is uniform. L now is the number of paths, N the number of frequency samples, and P is the unknown L×L diagonal power matrix. Taking K observations we arrive at the same form as in [START_REF] Ganesan | Coping with Irregular Spatio-temporal Sampling in Sensor Networks[END_REF]. In this case with uniform phase distribution, we can do even better than Proposition 7, in that one can write down estimators for the moments which are unbiased for any number of observations and frequency samples:

Proposition 8: Assume that V has uniform phase distribution, and let P i be the moments of P, and W i = tr N (W i ) the moments of the sample covariance matrix. Define also

c 1 = N K , c 2 = L N , and c 3 = L K . Then E [W 1 ] = c 2 P 1 + σ 2 E [W 2 ] = c 2 1 - 1 N P 2 + c 2 (c 2 + c 3 )(P 1 ) 2 +2σ 2 (c 2 + c 3 )P 1 + σ 4 (1 + c 1 ) E [W 3 ] = c 2 1 + 1 K 2 1 - 3 N + 2 N 2 P 3 + 1 - 1 N 3c 2 2 1 + 1 K 2 + 3c 2 c 3 P 1 P 2 + c 3 2 1 + 1 K 2 + 3c 2 2 c 3 + c 2 c 2 3 (P 1 ) 3 +3σ 2 (1 + c 1 )c 2 + c 1 c 2 2 KL 1 - 1 N P 2 +3σ 2 c 1 c 3 2 KL + c 2 2 + c 2 3 + 3c 2 c 3 (P 1 ) 2 +3σ 4 c 2 1 + 3c 1 + 1 + 1 K 2 c 2 P 1 +σ 6 c 2 1 + 3c 1 + 1 + 1 K 2
Just as Proposition 7, this is proved in Appendix M. In the following, this result is used in order to determine the number of paths as well as the power of each path. The different convergence rates of the approximations are clearly seen in the plots.

In Figure 9, the number of paths is estimated based on the procedure sketched above. We have set σ = √ 0.1, N = 100, and L = 36. The procedure is tried for 1 all the way up to 100 observations. The plot is very similar to Figure 8, in that only a small number of observations are needed in order to get an accurate estimate of L. When K = 1, it is seen that more observations are needed to get an accurate estimate of L, when compared to K = 10.

For the estimation of powers simulation, we have set K = N = L = 144, and σ = √ 0.1, following the procedure also described above, up to 1000 observations. The second plot in Figure 10 shows the results which confirms the usefulness of the approach.

C. Estimation of wavelength

In the field of MIMO cognitive sensing [START_REF] Cardoso | Cooperative Spectrum Sensing Using Random Matrix Theory[END_REF], [START_REF] Haykin | Cognitive Radio: Brain Empowered Wireless Communications[END_REF], terminals must decide on the band on which to transmit and in particular sense which band is occupied. One way of doing so is to find the wavelength λ in [START_REF] Cardoso | Cooperative Spectrum Sensing Using Random Matrix Theory[END_REF], based on some realizations of the sample covariance matrix. In our simulation, we have set K = 10, L = 36, N = 100, and σ = √ 0.1, in addition to

λ = 2, d = 1, α = π 4 .
We have tried values between 0 and 5 as candidate wavelengths (to be more precise, the values 0.05, 0.1, 0.15, ..., 5 are tried), and chosen the one which gives the smallest deviation (in the same sense as above, i.e. the sum of the squared errors of the first three moments are taken) from a different number of realizations of sample covariance matrices. The resulting plot is shown in Figure 11, and shows that the Vandermonde deconvolution method can also be used for wavelength estimation.

D. Signal reconstruction and estimation of the sampling distribution

For signal reconstruction, one can provide a general framework where only the sampling distribution matters asymptotically. The sampling distribution can be estimated with the help of the presented results. Several works have investigated how irregular sampling affects the performance of signal reconstruction in presence of noise in different fields, namely sensor networks [START_REF] Ganesan | Coping with Irregular Spatio-temporal Sampling in Sensor Networks[END_REF], [START_REF] Nordio | Performance of Field Reconstruction Techniques with Noise and Uncertain Sensor Locations[END_REF], image processing [START_REF] Strohmer | How to Recover Smooth Object Boundaries in Noisy Medical Images[END_REF], [START_REF] Early | Image Reconstruction and Enhanced Resolution Imaging from Irregular Samples[END_REF], geophysics [START_REF] Rauth | A frequency Domain Approach to the Reconvery of Geophysical Potenials[END_REF], and compressive sampling [START_REF] Candès | Sparsity and Incoherence in Compressive Sampling[END_REF]. The usual Nyquist Theorem states that for a signal with maximum frequency f max , one needs to sample the signal at a rate which is at least twice this number. However, in many cases, this can not be performed, or one has an observation of a signal at only a subset of the frequencies. Moreover, one feels that if the signal has a sparse spectrum, one can take fewer samples and still have the same information on the original signal. One of the central motivations of sparse sampling is 39) and ( 43), for various number of observations. K = N = L = 144, and σ = √ 0.1. The actual powers were 0.5, 1, and 1.5, with equal probability. exactly to understand under which condition one can still have less samples and recover the original signal up to an error of [START_REF] Feichtinger | Efficient Numerical Methods in Non-Uniform Sampling Theory[END_REF]. Let us consider the signal of interest as a superposition of its frequency components (this is also the case for a unidimensional bandlimited physical signal), i.e.

r(t) = 1 √ N N -1 k=0 s k e -j2πkt N
and suppose that the signal is sampled at various instants [t 1 , ..., t L ] with t i ∈ [0, 1]. This can be identically written as

r(ω) = 1 √ N N -1 k=0 s k e -jkω ,
or r = V T s. In the presence of noise, one can write where r = [r(ω 1 ), ...r(ω L )] T , s and n are as in [START_REF] Rappaport | Wireless Communications, Principles and Practice[END_REF], and with V on the form (1). A similar analysis for such cases can be found in [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF].

r = V T s + n, (44)
In the following, we suppose that one has K observations of the received sampled vector r:

Y = [r 1 , ...r K ] = V T [s 1 , ..., s K ] + [n 1 , ..., n K ] ( 45
)
The vector r is the discrete output of the sampled continuous signal r(w) for which the distribution is unknown (however, c is known). This case happens when one has an observation without the knowledge of the sampling rate for example. The difference in [START_REF] Rauth | A frequency Domain Approach to the Reconvery of Geophysical Potenials[END_REF] from the model [START_REF] Ganesan | Coping with Irregular Spatio-temporal Sampling in Sensor Networks[END_REF] lies in that the adjoint of a Vandermonde matrix is used, and in that there is no additional diagonal matrix P included. The following result can now be stated and proved similarly to Proposition 7 and 8: Proposition 9:

E [tr n (W)] = 1 + σ 2 (46) E tr n W 2 = c 2 I 2 + (1 + c 3 )(1 + σ 2 ) 2 (47) E tr n W 3 = 1 + 3c 2 (1 + c 3 )I 2 3c 3 + c 2 3 + c 2 2 I 3 3σ 2 (1 + 3c 3 + c 2 3 + c 2 (1 + c 3 )I 2 ) 3σ 4 c 2 (c 2 3 + 3c 3 + 1) σ 6 (c 2 3 + 3c 3 + 1), (48)
where lim N →∞

N K = c 1 , lim N →∞ L N = c 2 , lim N →∞ L K = c 3 , I n is
defined as in Proposition 7, and W = 1 K YY H . The proof of Proposition 9 is commented in Appendix M. We have tested ( 46)-( 48) by taking a phase distribution ω which is uniform on [0, α], and 0 elsewhere. The density is thus 2π

α on [0, α], and 0 elsewhere. In this case we can compute that

I 2 = 2π α I 3 = 2π α 2 .
The first of these equations, combined with ( 46)-( 48), enables us to estimate α. This is tested in Figure 12 for various number of observations. In Figure 13 we have also tested estimation of I 2 , I 3 from the observations using the same equations. When one has a distribution which is not uniform, the integrals I 3 , I 4 , ... would also be needed in finding the characteristics of the underlying phase distribution. Figure 13 shows that the estimation of I 2 requires far fewer observations than the estimation of I 3 . In both figures, the values K = 10, L = 36, N = 100, and σ = √ 0.1 were used and α was π 4 . It is seen that the of I 3 is a bit off even for higher number of observations. This is to be expected, since an asymptotic result is applied.

VI. CONCLUSION AND FURTHER DIRECTIONS

We have shown how asymptotic moments of random Vandermonde matrices with entries on the unit circle can be computed analytically, and treated many different cases. Vandermonde matrices with uniform phase distribution proved to be the easiest case, and it was shown how the case with more general phases could be expressed in terms of this. The case where the phase distribution has singularities was also handled, as this case displayed different asymptotic behaviour. Also, mixed moments of independent Vandermonde matrices were investigated, as well as the moments of generalized Vandermonde matrices. In addition to the general asymptotic expressions stated, exact expressions for the first moments of Vandermonde matrices with uniform phase distribution were also stated. We have also provided some useful applications of random Vandermonde matrices. The applications concentrated on deconvolution and signal sampling analysis. As shown, many useful system models use independent Vandermonde matrices and Gaussian matrices combined in some way. The presented examples show how random Vandermonde matrices in such systems can be handled in practice to obtain estimates on quantities such as the number of paths in channel modeling, the transmission powers of the users in wireless transmission, or the sampling distribution for signal recovery. The paper has only touched upon a limited number of applications, but the results already provide benchmark figures in the nonasymptotic regime.

From a theoretical perspective, it would also be interesting to find methods for obtaining the generalized expansion coefficients K ρ,ω,λ from K ρ,u,u , similar to how we found the expansion coefficients K ρ,ω from K ρ,u . This could also shed some light on whether uniform phase-and power distribution also minimizes moments of generalized Vandermonde matrices, similarly to how we showed that it minimizes moments in the non-generalized case.

Throughout the paper, we assumed that only diagonal matrices were involved in mixed moments of Vandermonde matrices. The case of non-diagonal matrices is harder, and should be addressed in future research. The analysis of the maximum and minimum eigenvalue is also of importance. The methods presented in this paper can not be used directly to obtain explicit expressions for the p.d.f. of the asymptotic mean eigenvalue distribution, so this is also a case for future research. A way of attacking this problem could be to develop for Vandermonde matrices analytic counterparts to what one has in free probability, such as the R-, S-, and the Stieltjes transform [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF]. Interestingly, certain matrices similar to Vandermonde matrices, have analytical expressions for the moments: in [START_REF] Bordenave | Eigenvalues of Euclidean random matrices[END_REF], analytical expressions for the moments of matrices with entries of the form A i,j = F (ω i -ω j ) are found. This is interesting for the Vandermonde matrices we consider,

since 1 N V H V i,j = sin N 2 (ω i -ω j ) N sin 1 2 (ω i -ω j )
.

Unfortunately, the function F N (x) = sin( N 2 x) N sin( 1 2 x) depends on the matrix dimension N , so that we can not find a function F which fits the result from [START_REF] Bordenave | Eigenvalues of Euclidean random matrices[END_REF].

Finally, another case for future research is the asymptotic behaviour of Vandermonde matrices when the matrix entries lie outside the unit circle.

APPENDIX A THE PROOF OF THEOREM 1

We can write

E tr L D 1 (N )V H VD 2 (N )V H V • • • D n (N )V H V (49) as L -1 i 1 ,...,i n j 1 ,...,j n E( D 1 (N )(j 1 , j 1 )V H (j 1 , i 2 )V(i 2 , j 2 ) D 2 (N )(j 2 , j 2 )V H (j 2 , i 3 )V(i 3 , j 3 ) . . . D n (N )(j n , j n )V H (j n , i 1 )V(i 1 , j 1 )) (50)
The (j 1 , ..., j n ) uniquely identifies a partition ρ of {1, ..., n}, where each block W j of ρ consists of the positions of the indices which equal j, i.e.

W j = {k|j k = j}.
We will also say that (j 1 , ..., j n ) give rise to ρ. Write

W j = {w j1 , w j2 , ..., w j|W j | }.
When (j 1 , ..., j n ) give rise to ρ, we see that since

j w j1 = j w j2 = • • • = j w j|W j | ,
we also have that

ω jw j1 = ω jw j2 = • • • = ω jw j|W j | ,
and we will denote their common value by ω W j as in Definition 8. With this in mind, it is straightforward to verify that (50) can be written as

ρ∈P(n) (i 1 ,...,i n ) (j 1 , ..., j n ) giving rise to ρ N -n L -1 × |ρ| k=1 E e j k∈W j i k-1 -k∈W j i k ω W k ×D 1 (N )(j 1 , j 1 ) × • • • × D n (N )(j n j n ), (51)
where i 1 , ..., i n takes values between 0 and N -1. We will in the following switch between the form (51) the form

ρ∈P(n) (j 1 , ..., j n ) giving rise to ρ (i 1 ,...,i n ) N |ρ|-n-1 c |ρ|-1 L -|ρ| ×E n k=1 e j(ω b(k-1) -ω b(k) )i k ×D 1 (N )(j 1 , j 1 ) × • • • × D n (N )(j n , j n ), (52)
where we also have reorganized the powers of N and L in (51), and changed the order of summation (i.e. summed over the different i 1 , ..., i n first). Noting that (i1,...,in)

N |ρ|-n-1 E n k=1 e j(ω b(k-1) -ω b(k) )i k (53) = N |ρ|-n-1 E   (i1,...,in) n k=1 e j(ω b(k-1) -ω b(k) )i k   (54) = N |ρ|-n-1 E n k=1 N -1 i k =0 e j(ω b(k-1) -ω b(k) )i k (55) = N |ρ|-n-1 E n k=1 1 -e jN (ω b(k-1) -ω b(k) ) 1 -e j(ω b(k-1) -ω b(k) ) (56)
= N |ρ|-n-1 × (0,2π) |ρ| n k=1 1 -e jN (ω b(k-1) -ω b(k) ) 1 -e j(ω b(k-1) -ω b(k) ) dω 1 • • • dω |ρ| (57) = K ρ,ω,N , (58)
Definition 8 of the Vandermonde mixed moment expansion coefficients comes into play, so that (52) can also be written

ρ∈P(n) (j 1 , ..., j n ) giving rise to ρ c |ρ|-1 L -|ρ| K ρ,ω,N ×D 1 (N )(j 1 , j 1 ) • • • × ×D n (N )(j n , n ). (59)
The notation for a joint limit distribution simplifies (52). Indeed, add to (52) for each ρ the terms

ρ ∈P(n),ρ >ρ (j 1 , ..., j n ) giving rise to ρ c |ρ|-1 L -|ρ| K ρ,ω,N ×D 1 (N )(j 1 , j 1 ) • • • × D n (N )(j n , j n ). (60)
These go to 0 as N → ∞, since they are bounded by

c |ρ|-1 L -|ρ| K ρ,ω,N L |ρ | = K ρ,ω,N c |ρ|-1 L |ρ |-|ρ| = O(L -1 ).
After this addition, the limit of (59) can be written

ρ∈P(n) c |ρ|-1 K ρ,ω D ρ , ( 61
)
which is what we had to show. We also need to comment on the statement of Theorem 6, where generalized Vandermonde matrices are considered. In this case, the derivations after (52) are different since the power distribution is not uniform. For the case of ( 22), we can in (55) replace

n i k =1 e j(ω b(k-1) -ω b(k) )i k with N -1 r=0 N p f N (r)e jr(ω b(k-1) -ω b(k) )
, since the number of occurrences of the power e jr(ω b(k-1) -ω b(k) ) is N p f N (r). The rest of the proof of Theorem 6 follows by canceling n powers of N after this replacement. The details are similar for the case [START_REF] Telatar | Capacity of multi-antenna gaussian channels[END_REF], where the law of large numbers is applied to arrive at the second formula in [START_REF] Debbah | MMSE Analysis of Certain Large Isometric Random Precoded Systems[END_REF].

APPENDIX B THE PROOF OF PROPOSITION 1

Note that for each block W j , E e

j k∈W j i k-1 -k∈W j i k ω W j = 0 when k∈W j i k-1 = k∈W j i k ,
and 1 if

k∈W j i k-1 = k∈W j i k . ( 62
)
If we denote by S ρ,N the set of all n-tuples 62), and define |S ρ,N | to be the cardinality of S ρ,N , it is clear that

(i 1 , ..., i n ) (0 ≤ i k ≤ N -1, 1 ≤ k ≤ n) which solve (
K ρ,u = lim N →∞ K ρ,u,N = lim N →∞ 1 N n+1-|ρ| |S ρ,N |.
It is straightforward to show that the solution set of (62) has n + 1 -|ρ| free variables. After dividing the equations (62) by N and letting N go to infinity, K ρ,u can thus alternatively be expressed as the volume in R n+1-|ρ| of the solution set of

k∈W j x k-1 = k∈W j x k , ( 63
) with 0 ≤ x k ≤ 1.
It is clear that the volume of this solution set computes to a rational number. It is the form (63) which will be used in the other appendices to compute K ρ,u for certain lower order ρ. Appendix D of [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF] states the same equations for finding quantities equivalent to Vandermonde mixed moment expansion coefficients for the uniform phase distribution. The fact that K ρ,u ≤ 1 follows directly from Appendix D of [START_REF] Nordio | Reconstruction of multidimensional signals from irregular noisy samples[END_REF].

The same applies for the fact that K ρ,u = 1 if and only if ρ is noncrossing.

For any ρ, we can define a partition of {1, ..., n} into n+1-|ρ| blocks, where two elements are defined to be in the same block if and only if the corresponding variables in solutions to (63) are linearly dependent. When ρ is noncrossing, it is straightforward to show that two such variables are dependent if and only if they are equal, and also that this partition is the Kreweras complement K(ρ) of ρ. This fact is used elsewhere in this paper.

We will also briefly explain why the computations in this appendix are useful for generalized Vandermonde matrices with uniform phase distribution. For [START_REF] Shamai | The impact of frequency-flat fading on the spectral efficiency of CDMA[END_REF], the number of solutions i 1 , ..., i k to (62) needs to be multiplied by

N p f N 1 ) • • • N p f N (i k ),
since each i j now may occur N p f N (i j ) times. This means that K ρ,ω,f can be computed as the integrals in this appendix, but that we also need to multiply with the density p f for each variable. The computations of these new integrals become rather involved when f is not uniform, and are therefore dropped.

APPENDIX C THE PROOF FOR PROPOSITION 2

We will in the following compute the volume of the solution set of (63), as a volume in [0, 1] n+1-|ρ| ⊂ R n+1-|ρ| , as explained in the proof of Proposition 1. These integrals are very tedious to compute, and many of the details are skipped. The formula r!s! (r s + 1)! = 1 0

x r (1 -x) s dx can be used to simplify some of the calculations for higher values of n.

A. Computation of K {{1,3},{2,4}},u

This is equivalent to finding the volume of the solution set of

x 1 + x 3 = x 2 + x 4
in R 3 . Since this means that x 4 = x 1 + x 3 -x 2 lies between 0 and 1, we can set up the following integral bounds: When x 1 + x 3 ≤ 1, we must have that 0 ≤ x 2 ≤ x 1 + x 3 , so that we get the contribution

1 0 1-x 1 0 x 1 +x 3 0 dx 2 dx 3 dx 1 ,
which computes to 1 3 . When 1 ≤ x 1 + x 3 , we must have that x 1 + x 3 -1 ≤ x 2 ≤ 1, so that we get the contribution

1 0 1 1-x1 1 x1+x3-1 dx 2 dx 3 dx 1 ,
which also computes to 1 3 . Adding the contributions together we get 2 3 , which is the stated value for K {{1,3},{2,4}},u . It turns out that when the blocks of ρ are cyclic shifts of each other, the computation of K ρ,u can be simplified. Examples of such ρ are {{1, 3}, {2, 4}} (for which we just computed K ρ,u ), {{1, 3, 5}, {2, 4, 6}}, and {{1, 4}, {2, 5}, {3, 6}}. We will in the following describe this simplified computation. Let a (m) l (x) be the polynomial which gives the volume in R m-1 of the solutions set to

x 1 + • • • + x m = x (constrained to 0 ≤ x i ≤ 1) for l ≤ x ≤ l + 1. It is clear that these satisfy the integral equations a (m+1) l (x) = l x-1 a (m) l-1 (t)dt + x l a (m) l (t)dt, (64)
which can be used to compute the a m l (x) recursively. Note first that a For m = 3 we have

a (3) 0 (x) = x 0 a (2) 0 (t)dt = 1 2 x 2 a (3) 1 (x) = 1 x-1 a (2) 0 (t)dt + x 1 a (2) 1 (t)dt = 1 - 1 2 (x -1) 2 - 1 2 (2 -x) 2 a (3) 2 (x) = 2 x-1 a (2) 1 (t)dt = 1 2 (3 -x) 2 .
B. Computation of K {{1,3,5},{2,4,6}},u For m = 3, integration gives 1 0 (a

(3) 0 ) 2 (t)dt + 2 1 (a (3) 1 ) 2 (t)dt + 3 2 (a (3)
2 ) 2 (t)dt, which computes to 11 20 . This is the stated expression for K {{1,3,5},{2,4,6}},u .

C. Computation of K {{1,4},{2,5},{3,6}},u

This is equivalent to finding the volume of the solution set of

x 1 + x 4 = x 2 + x 5 = x 3 + x 6 in R 4 , which is computed as 1 0 (a (2) 0 ) 3 (t)dt + 2 1 (a (2)
1 ) 3 (t)dt, which computes to 1 2 . This is the stated expression for K {{1,4},{2,5},{3,6}},u .

D. Computation of K {{1,4},{2,6},{3,5}},u

This is equivalent to finding the volume of the solution set of

x 1 + x 4 = x 2 + x 5 x 2 + x 6 = x 3 + x 1 in R 4 .
Since this means that x 5 = x 1 -x 2 + x 4 lies between 0 and 1, x 6 = x 1 -x 2 + x 3 lies between 0 and 1, we can set up the following integral bounds: For x 2 ≥ x 1 we must have x 2 -x 1 ≤ x 3 , x 4 ≤ 1, so that we get the contribution

1 0 1 x 1 1 x 2 -x 1 1 x 2 -x 1 dx 4 dx 3 dx 2 dx 1 ,
which computes to 1 4 . It is clear that for x 1 ≥ x 2 we get the same result by symmetry, so that the total contribution is This is equivalent to finding the volume of the solution set of

x 1 + x 5 = x 2 x 6 x 3 + x 7 = x 4 + x 1 in R 5 , or
x 6 = x 5 + x 1 -x 2 lies between 0 and 1, x 7 = x 4 + x 1 -x 3 lies between 0 and 1 .

(65)

This can be split into the following volumes: 1)

x 1 ≤ x 2 ≤ x 3 , 2) x 1 ≤ x 3 ≤ x 2 , 3) x 3 ≤ x 2 ≤ x 1 , 4) x 2 ≤ x 3 ≤ x 1 , 5) x 2 ≤ x 1 ≤ x 3 , 6) x 3 ≤ x 1 ≤ x 2 .
Each of these volumes can be computed by setting up an integral with corresponding bounds. Computing these we get the values This is equivalent to finding the volume of the solution set of

x 1 + x 6 = x 2 + x 7 x 2 + x 4 = x 3 + x 5
in R 5 , or x 6 = x 7 + x 2 -x 1 lies between 0 and 1, x 5 = x 4 + x 2 -x 3 lies between 0 and 1, .

This can be obtained from (65) by a permutation of the variables, so the contribution from K {{1,6},{2,4},{3,5,7}},u must also be 9 20 , which proves the claim. PROPOSITION 3 Note first that multiplying both of sides of (13) with c gives

APPENDIX D THE PROOF FOR

cM n = ρ∈P(n) K ρ,ω (cD) ρ , ( 66
)
where we now can substitute the scaled moments ( 14)- [START_REF] Mehta | Random Matrices[END_REF].

With D 1 (N ) = D 2 (N ) = • • • = D n (N ) = D(N ), D ρ
as defined in Definition 2 does only depend on the block cardinalities |W j |, so that we can group together the K ρ,ω for ρ with equal block cardinalities. If we group the blocks of ρ so that their cardinalities are in descending order, and set

P(n) r 1 ,r 2 ,...,r k = {ρ = {W 1 , ..., W k } ∈ P(n)||W i | = r i ∀i},
where r 1 ≥ r 2 ≥ • • • ≥ r k , and also write

K r 1 ,r 2 ,...,r k = ρ∈P(n) r 1 ,r 2 ,...,r k K ρ,ω , (67)
(66) can be written

m n = r 1 ,...,r k r 1 +•••+r k =n K r 1 ,r 2 ,...,r k k j=1 d r j . ( 68
)
For the first 5 moments this becomes

m 1 = K 1 d 1 (69) m 2 = K 2 d 2 + K 1,1 d 2 1 (70) m 3 = K 3 d 3 + K 2,1 d 2 d 2 1 + K 1,1,1 d 3 1 ( 71
)
m 4 = K 4 d 4 + K 3,1 d 3 d 1 + K 2,2 d 2 2 + K 2,1,1 d 2 d 2 1 + K 1,1,1,1 d 4 1 (72) m 5 = K 5 d 5 + K 4,1 d 4 d 1 + +K 3,2 d 3 d 2 + K 3,1,1 d 3 d 2 1 + K 2,2,1 d 2 2 d 1 + K 2,1,1,1 d 2 d 3 1 + K 1,1,1,1,1 d 5 1 . ( 73
)
Thus, to prove Proposition 3, we have to compute the K r1,r2,...,r k by going through all partitions. We will have use for the following result, taken from [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF]: Lemma 1: The number of noncrossing partitions in N C(n) with r 1 blocks of length 1, r 2 blocks of length 2 and so on (so that r 1 + 2r

2 + 3r 3 + • • • nr n = n) is n! r 1 !r 2 ! • • • r n !(n + 1 -r 1 -r 2 • • • r n )! .
Using this and a similar formula for the number of partitions with prescribed block sizes, we obtain cardinalities for noncrossing partitions and the set of all partitions with a given block structure. These numbers are the used in the following calculations. For the proof of Proposition 3, we need to compute (67) for all possible block cardinalities (r 1 , ..., r k ), and insert these in (69)-(73). The formulas for the three first moments are obvious, since all partitions of length ≤ 3 are noncrossing. For the remaining computations, the following two observations save a lot of work:

• If ρ 1 ∈ P(n 1 ), ρ 2 ∈ P(n 2 )
with n 1 < n 2 , and ρ 1 can be obtained from ρ 2 by omitting elements k in {1, ..., n 2 } such that k and k + 1 are in the same block, then we must have that K ρ1,u = K ρ2,u . This is straightforward to prove since it follows from the proof of Proposition 1 that i k+1 can be chosen arbitrarily between 0 and N -1 in such a case.

• K ρ 1 ,u = K ρ 2 ,

A. The moment of fourth order

The result is here obvious except for the case for the three partitions with block cardinalities (2, 2) (for all block cardinalities, all partitions are noncrossing, so that K r 1 ,r 2 ,...,r k is simply the number of noncrossing partitions with block cardinalities (r 1 , ..., r k ). this number can be computed from Lemma 1). Two of the partitions with blocks of cardinality (2, 2) are noncrossing, the third one is not. We see from Proposition 2 that the total contribution is

K 2,2 = 2 + K {{1,3},{2,4}},u = 2 + 2 3 = 8 3 .
The formula for the fourth moment follows.

B. The moment of fifth order

Here two cases require extra attention:

1) ρ = {W 1 , W 2 } with |W 1 | = 3, |W 2 | = 2:
There are 10 such partitions, and 5 of them have crossings and contribute with K {{1,3},{2,4}},u . The total contribution is therefore

5 + 5 × K {{1,3},{2,4}},u = 5 + 5 × 2 3 = 25 3 . 2) ρ = {W 1 , W 2 , W 3 } with |W 1 | = |W 2 | = 2, |W 3 | = 1:
There are 15 such partitions, of which 5 have crossings. The total contribution is therefore 10 + 5 × K {{1,3},{2,4}},u = 10 + 5 × 2 3 = 40 3 . The computations for the sixth and seventh order moments are similar, but the details are skipped. These are more tedious in the sense that one has to count the number of partitions with a given block structure, and identify each partition with one of the coefficients listed in Proposition 2.

APPENDIX E THE PROOF OF PROPOSITION 4

C i,j (D(N )V H V) is computed as in Appendix A. Since some terms in E tr L A i tr L A j cancel those in E tr L A i E tr L A j , we can restrict to summing over partitions of 1, 2, ..., i + j where at least one block contains elements from both [1, ..., i] and [i+1, ..., i+j]. We denote this set by P(i, j), and set n = i + j. In our new calculations,(52) now instead takes the form

L ρ∈P(i,j) (j 1 , ..., j n ) giving rise to ρ (i 1 ,...,i n ) N |ρ|-i-j-1 L -1 c |ρ|-1 L -|ρ| ×E n k=1 e j(ω b(k-1) -ω b(k) )i k ×D 1 (N )(j 1 , j 1 ) × • • • × D n (N )(j n , j n ),( 74
)
where the normalizing factor L from Definition 4 has been included. Simplifying this as in Appendix A, and restricting to uniform phase distribution, we obtain

lim L→∞ LC i,j (D(N )V H V) = ρ∈P(i,j) c |ρ|-1 K 2,ρ,u D ρ ,
where K 2,ρ,u is the volume of the solution set of

k∈W j x σ -1 (k-1) = k∈W j x k , ( 75
)
where σ is the permutation which shifts [1, i] and [i+1, ..., i+ j] to the right cyclically so that the result is contained within the same interval. Thus, when the normalizing factor L is included, we see that the second order moments exist. 17) is computed by noting that K 2,{{1,3},{2,4}},u and K 2,{{1,4},{2,3}},u both equal 2 3 , and that there are 9 other partitions in P(2, 2), and K 2,π,u = 1 for all these π (all these values are computed as in Appendix C). By adding up for the different block cardinalities we get that

C 2,2 (D(N )V H V) in (
c lim L→∞ LC 2,2 (D(N )V H V) = 4 3 d 2 2 ,
and using the substitution ( 16) we arrive at the desired result.

APPENDIX F THE PROOF OF THEOREM 2

In order to get the exact expressions in Theorem 2, we now need to keep track of the K ρ,u,N defined by [START_REF] Sampaio | Vandermonde frequency division multiplexing[END_REF], not only the limits K ρ,u (if we had not assumed ω = u, the calculations for K ρ,ω,N would be much more cumbersome). When ρ is a partition of {1, ..., n} and n ≤ 4, we have that K ρ,u,N = K ρ,u = 1 when ρ = {{1, 3}, {2, 4}}. We also have that

K {{1,3},{2,4}},u,N = 2 3 + 1 3N 2 , ( 76
)
where we have used that

N i=1 i 2 = N 3 (N + 1)(N + 1
2 ) [START_REF] Rottmann | Mathematische Formelsammlung. B.I. Wissenschaftsverlag[END_REF]. We also need the exact expression for the quantity

T ρ = (j 1 ,...,jn) giving rise to ρ L -|ρ| D 1 (N )(j 1 , j 1 )ו • •×D n (N )(j n , j n ) from (
T ρ = D (N,L) ρ - ρ >ρ L |ρ |-|ρ| T ρ , ( 77
)
which can be used recursively to express the T ρ in terms of the D (N,L) ρ

. We obtain the following formulas for n = 4:

T {{1,2,3,4}} = D (N,L) 4 (78) T {{1,2,3},{4}} = (N,L) 3 D (N,L) 1 -L -1 D (N,L) 4 (79) T {{1,2},{3,4}} = (D (N,L) 2 ) 2 -L -1 D (N,L) 4 (80) T {{1,2},{3},{4}} = D (N,L) 2 (D (N,L) 1 ) 2 -2L -1 (D (N,L) 3 D (N,L) 1 -L -1 D (N,L) 4 ) -L -1 (D (N,L) 2 ) 2 -L -1 D (N,L) 4 -L -2 D (N,L) 4 = D (N,L) 2 (D (N,L) 1 ) 2 -L -1 (D (N,L) 2 ) 2 -2L -1 D (N,L) 3 D (N,L) 1 +2L -2 D (N,L) 4
(81)

T {{1},{2},{3},{4}} = (D (N,L) 1 ) 4 -6L -1 (D (N,L) 2 (D (N,L) 1 ) 2 -L -1 (D (N,L) 2 ) 2 -2L -1 D (N,L) 3 D (N,L) 1 +2L -2 D (N,L) 4 ) -3L -2 (D (N,L) 2 ) 2 + 3L -3 D (N,L) 4 -4L -2 D (N,L) 3 D (N,L) 1 +4L -3 D (N,L) 4 -L -3 D (N,L) 4 = -6L -3 D (N,L) 4 +L -2 (8D (N,L) 3 D (N,L) 1 +3(D (N,L) 2 ) 2 ) -6L -1 D (N,L) 2 (D (N,L) 1 ) 2 + (D (N,L) 1 ) 4 . ( 82
)
For n = 3 and n = 2 the formulas are

T {{1,2,3}} = D (N,L) 3 (83) T {{1,2},{3}} = D (N,L) 1 D (N,L) 2 -L -1 D (N,L) 3 (84) T {{1},{2},{3}} = (D (N,L) 1 ) 3 -3L -1 D (N,L) 1 D (N,L) 2 +2L -2 D (N,L) 3 (85) T {{1,2}} = D (N,L) 2 (86) T {{1},{2}} = (D (N,L) 1 ) 2 -L -1 D (N,L) 2 . ( 87
)
It is clear that (78)-( 82) and ( 83)-(87) cover all possibilities when it comes to partition block sizes. Using ( 14)-( 15), and putting (76), ( 78)-(82), and (83)-( 87) into (59) we get the expressions in Theorem 2 after some calculations.

If we are only interested in first order approximations rather than exact expressions, (77) gives us

T ρ ≈ D ρ - ρ >ρ |ρ|-|ρ |=1 L -1 D ρ ,
which is easier to compute. Also, we need only first order approximations to K ρ,u,N , which is much easier to compute than the exact expression. For (76), K {{1,3},{2,4}},u,N ≈ 2 3 is already a first order approximation. Inserting the approximations in (59) gives a first order approximation of the moments.

where we have converted to Lebesgue measure, and where we have also written

dx = dx 1 • • • dx |ρ| . Consider first the set U = {ω||x b(k-1) -x b(k) | ≤ π∀k}. When N ≤ |ω b(k-1) -ω b(k) | ≤ π, it is clear that sin N (x b(k-1) -x b(k) )/2 sin (x b(k-1) -x b(k) )/2 ≤ 4 x b(k-1) -x b(k) , ( 90
) since sin N (x b(k-1) -x b(k) )/2 ≤ 1, and since | sin(x)| ≥ | x 2 | when |x| ≤ π 2 . When |x b(k-1) -x b(k) | ≤ 2π N we have that sin N (x b(k-1) -x b(k) )/2 sin (x b(k-1) -x b(k) )/2 ≤ N. ( 91
)
Let k 1 , ..., k |ρ| ∈ Z, and assume that k |ρ| = 0. By using the triangle inequality, it is clear that on the set

D k 1 ,...,k |ρ|-1 = {ω| x i - 2k i π N ≤ π N ∀1 ≤ i ≤ |ρ|}, when |k r -k s | ≥
N n+1-|ρ| U |F (x)|dx ≤ 1 N n+1-|ρ| N n × 0≤k 1 ,...,k|ρ|-1<N all k i different n r=1 D |k b(r-1) -k b(r) | 2π(2π) |ρ|-1 N 1-|ρ| = (2π) |ρ| D n 0≤k 1 ,...,k|ρ|-1<N all k i different n r=1 1 |k b(r-1) -k b(r) | , ( 92
) where we have integrated w.r.t. x |ρ| also (i.e. k |ρ| is kept constant in (92)). A similar analysis as for U applies for the complement set

V = {ω|π ≤ |x b(k-1) -x b(k) | ≤ 2π for some k}, so that we can find a constant C such that 1 N n+1-|ρ| [0,2π) |ρ| |F (x)|dx ≤ C 0≤k 1 ,...,k|ρ|-1<N all k i different n r=1 1 |k b(r-1) -k b(r) | , (93)
It is clear this sum converges: First of all, this is only needed to prove for ρ = 0 n , since the summands for ρ = 0 n is only a subset of the summands for ρ = 0 n . Secondly, for ρ = 0 n , (93) can be bounded by considering convolutions of the following function with itself:

f (x) = 1 |x| for |x| > 1 0 for |x| ≤ 1 (94)
The assumption that f (x) = 0 in a neighbourhood of zero is due to the fact that the k i are all different. Note that |f (x)| ≤ 1 |x| 1-for any 0 < < 1. Also, the n -2-fold convolution (we wait with the n -1'th convolution till the end) of with itself exist outside 0 whenever 0 < (n -2) < 1, and is on the form r

1 |x| 1-(n-2)
for some constant r [START_REF] Rottmann | Mathematische Formelsammlung. B.I. Wissenschaftsverlag[END_REF]. Therefore, (93) is bounded by

|x|>1 r 1 |x| 1-(n-2) 1 |x| dx = |x|>1 r 1 |x| 2-(n-2) dx = 2r (n -2) -1 .
This proves that the entire sum (93) is bounded, and thus also the statement on the existence of the limit K(ρ, ω) in Theorem 3 when the density is continuous.

For the rest of the proof of Theorem 3 , we first record the following result:

Lemma 2: For any > 0, lim

N →∞ 1 N n+1-|ρ| B ,r F (ω)dω = 0, (95)
where This means that in the integral for K ρ,ω,N , we need only integrate over the ω which are arbitrarily close to the diagonal, (where ω 1 = • • • = ω |ρ| ). We thus have

B ,r = {(ω 1 , ..., ω |ρ| )||ω b(r-1) -ω b(r) | > }.
K ρ,ω = lim N →∞ 1 N n+1-|ρ| [0,2π) |ρ| F (x) |ρ| r=1 p ω (x r )dx = lim N →∞ 1 N n+1-|ρ| [0,2π) |ρ| F (x)p ω (x |ρ| ) |ρ| dx = lim N →∞ 1 N n+1-|ρ| 2π 0 p ω (x |ρ| ) |ρ| [0,2π) |ρ|-1 F (x)dx 1 • • • dx |ρ|-1 dx |ρ| .
We used here that the density is continuous. Using that

lim N →∞ 1 N n+1-|ρ| [0,2π) |ρ|-1 F (x)dx 1 • • • dx |ρ|-1 = (2π) |ρ|-1 K ρ,u (96)
when x |ρ| is kept fixed at an arbitrary value (this is straightforward by using the methods from the proof of Proposition 1 and (89)), we get that the above equals

K ρ,u (2π) |ρ|-1 2π 0 p ω (x |ρ| ) |ρ| dx |ρ| ,
which is what we had to show.

This has it's highest order when |ρ| = n, so that we can restrict to looking at 0 n . Note also that we may just as well assume that p ω (x) is identical to p i -ω i | -s at an interval around ω i , since lim x→α i |x -α i | s p ω (x) = p i implies that

p ω (x) = p i |x -ω i | -s + k(x)|x -ω i | -s (100)
where lim x→ωi k(x) = 0. It is straightforward to see that the contribution of the second part in (100) to (99) vanishes as N → ∞, so that we may just as well assume that p ω (x) is identical to p i |x -ω i | -s at an interval around ω i , as claimed. Also, since lim n→∞ |x|>

x -s e jnx dx = 0 for all > 0, and since the contributions from large n dominate in (101) below (since n |n| -s diverges), it is clear that we can restrict to an interval around ω i when computing the limit also (since p ω is continuous outside the singularity points, this follows from Theorem 3, and due to the additional 1 N s -factor added to (1)). After restricting to 0 n , multiplying both sides with N , summing over all singularity points, and using (98), we obtain the approximation It therefore suffices to prove that lim N →∞ q (n,N ) = q (n) . To see this, write

N -s |i k-1 -i k | 1-s = 1 N 1 1 N 1-s |i k-1 -i k | 1-s = 1 N 1 i k-1 N -i k N 1-s .
Summing over all 1 ≤ i 1 , ..., i n ≤ N , it is clear from this that q (n,N ) can be viewed as a Riemann sum which converges to q (n) as N → ∞.

APPENDIX L THE PROOF OF THEOREM 7 AND COROLLARY 1

Proof of Theorem 7: we define S j to be the blocks of σ, i.e. S j = {k|i k = j}.

Note that Theorem 3 guarantees that the limit K ρ,ω = lim N →∞ K ρ,ω,N exists. The partition ρ simply is a grouping of random variables into independent groups. It is therefore impossible for a block in ρ to contain elements from both S 1 and S 2 , so that any block is contained in either S 1 or S 2 . As a consequence, ρ ≤ σ.

Until now, we have not treated mixed moments of the form

D 1 (N )V i 2 V H i2 D 2 (N )V i 3 V H i3 • • • × D n (N )V i 1 V H i1 ,
which are the same as the mixed moments of Theorem 7 except for the position of the D i (N ). We will not go into depths on this, but only remark that this case can be treated in the same vein as generalized Vandermonde matrices by replacing the density p f (or p λ in case of continuous generalized Vandermonde matrices) with functions p D i (x) defined by p Di (x) = D i (N )( Lx , Lx ) for 0 ≤ x ≤ 1. This also covers the case of mixed moments of independent, generalized Vandermonde matrices (and, in fact, there are no restrictions on the horizontal and vertical phase densities p ω i and p λ j for each matrix. They may all be different). The proof for this is straightforward.

Proof of Corollary 1: this follows in the same way as Proposition 3 is proved from Proposition 2, by only considering ρ which are less than σ, and also by using Theorem 3. σ are for the listed moments {{1}, {2}}, {{1, 3}, {2, 4}}, and {{1, 3, 5}, {2, 4, 6}}, respectively.

APPENDIX M THE PROOFS OF PROPOSITION 7 AND 8

The moments E tr n W i will be related to the moments P i through three convolution stages: 1) relating the moments of W with the moments of

Γ = VP 1 2 1 K SS H P 1 2 V H , ( 102
)
from which we easily get the moments of

S = 1 K SS H P 1 2 V H VP 1 2 , ( 103
)
2) relating the moments of S with the moments of

T = PV H V, (104)
3) relating the moments of T with the moments of P.

+σ 6 c 2 1 + 3c 1 + 1 + 1 K 2 , , (107)
which are obtained by replacing R in [START_REF] Ryan | Channel capacity estimation using free probability theory[END_REF] by VP

+σ 6 c 2 1 + 3c 1 + 1 + 1 K 2 . ( 116
)
Up to now, all formulas have provided exact expressions for the expectations. For the next step, exact expressions for the expectations are only known when the phase distributions are uniform, in which case the formulas are given by Theorem 2: where the approximation is O(N -1 ), and where I k is defined by [START_REF] Nordio | Performance of Field Reconstruction Techniques with Noise and Uncertain Sensor Locations[END_REF]. Proposition 8 is proved by combining (114)-( 116) with (117)-(122), while Proposition 7 is proved by combining (114)-( 116) with ( 123)-(128). Proposition 9 is proved by first observing that the roles of L and N are interchanged, since the Vandermonde matrix is replaced by its transpose. This means that we obtain the formulas (114)-(116), with c 1 and c 3 interchanged, and c 2 replaced with c2 . The matrix T is now instead VV H , and these can be scaled to obtain the moments of V H V. Finally the integrals I n or the angle α can be estimated from these moments, using (123)-(128) with the moments of P replaced with 1 (since no additional power matrix is included in the model).

Matlab code for implementing the steps (105)-( 107), ( 110)-(113), and (117)-(122) can be found in [START_REF] Ryan | Tools for convolution with Vandermonde matrices[END_REF].

4 j=1 Vj -E j 2 . 7 )

 427 are computed from the samples. 3) The 4 first estimated moments Vj are computed as the mean of the sample moments, i.e. Vj = 1 K K i=1 mji . 4) The 4 first exact moments E j are computed using Theorem 2. 5) The 4 first asymptotic moments A j are computed using Proposition 3. 6) The mean squared error (MSE) of the first 4 estimated moments from the exact moments is computed as The MSE of the first 4 exact moments from the asymptotic moments is computed as

  Fig.1. MSE of the first 4 estimated moments from the exact moments for 80 and 320 samples for varying matrix sizes, with N = L. Matrices are on the form V H V with V a Vandermonde matrix with uniform phase distribution. The MSE of the first 4 exact moments from the asymptotic moments is also shown.

  Fig. 2. MSE of the first 4 estimated moments from the exact moments for 5 samples for varying matrix sizes, with N = L. Matrices are on the form 1

Fig. 3 .Fig. 4 . 1 N

 341 Fig. 3. Histogram of the mean eigenvalue distribution of 640 samples of V H V, with V a 1600 × 1200 Vandermonde matrix with uniform phase distribution.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Histogram of the mean eigenvalue distribution of 640 samples of V H V, with V a 1600 × 1200 Vandermonde matrix with phase distribution pω defined in (40) with λ = 2d, α = π 4 .

Fig. 7 .

 7 Fig. 7. Realizations of the capacity for Gaussian and Vandermonde matrices of size 36 × 36.

10 Fig. 8 .

 108 Fig. 8. Estimate for the number of users. Actual value of L is 36. Also, σ = √ 0.1, N = 100. The powers were 0.5, 1, and 1.5, with equal probability.

10 Fig. 9 .

 109 Fig. 9. Estimate for the number of paths. Actual value of L is 36. Also, σ = √ 0.1, N = 100.

Fig. 10 .

 10 Fig.10. Estimation of powers for the two models[START_REF] Cardoso | Cooperative Spectrum Sensing Using Random Matrix Theory[END_REF] and[START_REF] Strohmer | How to Recover Smooth Object Boundaries in Noisy Medical Images[END_REF], for various number of observations. K = N = L = 144, and σ = √ 0.1. The actual powers were 0.5, 1, and 1.5, with equal probability.

Fig. 11 .

 11 Fig. 11. Estimation of wavelength. Deconvolution was performed for varying number of observations, assuming different wavelengths, In the true model (39), λ = 2, d = 1, α = π 4 , K = 10, L = 36, N = 100, and σ = √ 0.1.

Fig. 12 . 3 Fig. 13 .

 12313 Fig. 12.Estimated values of α using (46)-(48), for various number of observations, and for K = 10, L = 36, N = 100, σ = √ 0.1. The actual value of α was π 4 .

  )dt = 2 -x.

  59) (i.e. we can not add (60) to obtain the approximation (61) here). Setting D (N,L) n = tr L (D n (N )), and D (N,L) ρ = k i=1 D (N,L) W i , we see that

Proof:

  The set B ,r corresponds to those k 1 , ..., k |ρ| in (93) for which |k b(r-1) -k b(r) | > N 2π . Thus, for large N , we sum over k 1 , ..., k |ρ| in (93) for which |k b(r-1) -k b(r) | is arbitrarily large. By the convergence of the Fourier integral of 1 |x| , it is clear that this converges to zero. Define B = {(ω 1 , ..., ω |ρ| )||ω i -ω j | > for some i, j}. If ω ∈ B , there must exist an r so that |ω b(r-1) -ω b(r) | > 2 n , so that ω ∈ B r,2 /n . This means that B ⊂ ∪ r B r,2 /n , so that by Lemma 2 also lim N →∞ 1 N n+1-|ρ| B F (ω)dω = 0.

( 1 × 1 - 2 n(

 112 i1,...,in) a N -ns c n-2p a Γ(1 -s) cos (i k-1 -i k )αa |i k-1 -i k | 1-s ×tr L (D 1 (N )) × • • • × tr L (D n (N )) (101) to (52). Since n k=1 e j(i k-1 -i k )α a = 1, we recognize q (n,N ) = 2Γ(1 -s) cos (1-s)π a p n a ) × (i 1 ,...,i n ) N -ns n k=1 1 |i k-1 -i k | 1-s ,as a factor in (101) such that the limit of (101) as N → ∞ can be writtenc n-1 lim N →∞ q (n,N ) lim N →∞ n i=1tr L (D i (N )) .

For 2 +2σ 2 ( 1 + 3 +3σ 2 ( 1 + c 1 )E tr N Γ 2 +3σ 2 c 1 E 2 +3σ 4 c 2 1 + 3c 1 + 1 + 1 K 2 E

 2213211212112 the first stage, the moments of Ŵ and Γ relate through the formulasE [tr n (W)] = E [tr N (Γ)] + σ 2 (105) E tr n W 2 = E tr N Γ c 1 )E N (Γ)] +σ 4 (1 + c 1 ) (106) E tr n W 3 = E tr N Γ (tr N (Γ)) [tr N (Γ)]

1 2c 3 2 =c 3 2 +2σ 2 (c 2 +E tr n W 3 = c 2 1 + 1 K 2 2 +c 2 c 2 3 E 3 +3σ 2 ( 1 + c 1 )c 2 + c 1 c 2 2 KL E tr L T 2 +3σ 2 c 3 (c 3 + 2 +3σ 4 c 2 1 + 3c 1 + 1 + 1 K 2 c 2 E

 1323222312233211222332112 S, with c = c 1 = N K . For the second part of the first stage, note thatE tr N Γ k = c 2 E tr L Sk (108) E (tr N (Γ)) k = c k 2 E tr L S k , (109)where c 2 = L N . We can now apply Theorem 2 to obtainc 3 E tr L S = c 3 E [tr L (T)] E tr L S3 = 1 + K -2 c 3 E tr L T 3 +3c 2 3 E (tr L T) tr L T 2 +c 3 3 E (tr L (T)) 3(112)E tr L S E (tr L (T)) = L K , and T = PV H V. (105)-(107), (108)-(109), and (110)-(113) can be combined toE [tr n (W)] = c 2 E [tr L (T)] + σ 2 (114) E tr n W 2 = c 2 E tr L T 2 + c 2 c 3 E (tr L (T)) c 3 )E [tr L (T)] + σ 4 (1 + c 1 )(115) E tr L T 3 +3c 2 c 3 E (tr L (T)) tr L T (tr L (T)) 2c 2 )E (tr L (T)) [tr L (T)]

c 2 E 2 = 3 = 2 =E

 2232 [tr L (T)] = c 2 tr L (P) (117) c 2 E tr L T 2 = 1 -N -1 c 2 tr L (P 2 ) +c 2 2 (tr L (P)) 2 (118) c 2 E tr L T 3 = 1 -3N -1 + 2N -2 c 2 tr L (P 3 ) +3 1 -N -1 c 2 2 tr L (P)tr L (P 2 ) +c 3 2 (tr L (P)) 3 (119) E (tr L (T)) tr L (P) 2 (120) E (tr L (T)) tr L (P) 3 (121) E (tr L (T)) tr L T 2 = 1 -N -1 tr L (P )tr L (P 2 ) + c 2 (tr L (P)) 3 . (122)If the phase distribution ω is not uniform, Theorem 1 and Theorem 3 gives the following approximation:c 2 E [tr L (T)] = c 2 tr L (P) (123) c 2 E tr L T 2 ≈ c 2 tr L (P 2 ) + c 2 2 I 2 (tr L (P)) 2 (124) c 2 E tr L T 3 ≈ c 2 tr L (P 3 ) + 3c2 2 I 2 tr L (P)tr L (P 2 ) +c 3 2 I 3 (tr L (P)) 3 (125) E (tr L (T)) (tr L (T)) tr L T 2 ≈ tr L (P)tr L (P 2 ) + c 2 I 2 (tr L (P)) 3 , (128)

  u if the set of equations (63) for ρ 1 can be obtained by a permutation of the variables in the set of equations for ρ 2 . Since the rank of the matrix for (63) equals the number of equations -1, we actually need only have that |ρ 1 | -1 of the |ρ 1 | equations can be obtained from permutation of |ρ 2 | -1 equations of the |ρ 2 | equations in the equation system for ρ 2 .

i=1 |W i -Ŵi | 2 .

ACKNOWLEDGMENT

The authors would like to thank the anonymous reviewers for their insightful and valuable comments, which have helped improve the quality of the paper. They would also like to thank the Associate Editor Prof. G. Taricco for a very professional processing of the manuscript.

This project is partially sponsored by the project BIONET (INRIA) This work was supported by Alcatel-Lucent within the Alcatel-Lucent Chair on flexible radio at SUPELEC as well as the ANR project SESAME This paper was presented in part at the 1st Workshop on Physics-Inspired

APPENDIX G THE PROOF OF PROPOSITION 5 We only state the proof for the case c = 1. In [START_REF] Bryc | Spectral measure of large random Hankel, Markov and Toeplitz matrices[END_REF] it is stated that the asymptotic 2n-moment (m 2n ) of certain Hankel and Toeplitz matrices can be expressed in terms of volumes of solution sets of equations on the form (63), with ρ restricted to partitions with all blocks of length 2. Rephrased in our language of Vandermonde mixed moment expansion coefficients, this means that m 2n = ρ ∈ P(2n) ρ has two elements in each block

In the language of [START_REF] Bryc | Spectral measure of large random Hankel, Markov and Toeplitz matrices[END_REF], the formula is not stated exactly like this, but rather in terms of volumes of solution sets of equations of the form (63). This translates to (88), since we in Appendix B interpreted K ρ,u as such volumes. In Proposition A.1 in [START_REF] Bryc | Spectral measure of large random Hankel, Markov and Toeplitz matrices[END_REF], unbounded support was proved by showing that (m 2n ) 1/n → ∞. Again denoting the asymptotic moments of Vandermonde matrices with uniform phase distribution by V n , we have that m 2n ≤ V 2n , since we sum over a greater class of partitions than in (88) when computing the Vandermonde moments. This means that (V 2n ) 1/n → ∞ also, so that the asymptotic mean eigenvalue distribution of the Vandermonde matrices have unbounded support also.

APPENDIX H THE PROOF OF THEOREM 3

We will use the fact that

where integration is w.r.t. Lebesgue measure.

For ρ = 1 n Theorem 3 is trivial. We will thus assume that ρ = 1 n in the following. We first prove that lim N →∞ K ρ,ω,N exists whenever p ω is continuous. To simplify notation, define

and set ω = (ω 1 , ..., ω |ρ| ) and dω = dω 1 • • • dω |ρ| . Since ω is continuous, there exists a p max such that p ω (ω i ) ≤ p max for all ω i . Then we have that PROPOSITION 6 Proposition 6 will follow directly if we can prove the following result:

n ] and define ω λ1,...,λn (0

Proof: This follows immediately by noting that

, where we have used that

APPENDIX J THE PROOF OF THEOREM 4

The contribution in the integral K ρ,ω,N comes only from when the ω i coincide with the atoms of p. Actually, we evaluate 1-e jN ω

1-e jω in points on the form ω = α i -α j . This evaluates to N n p n i when all ω are chosen equal to the same atom α j . Since lim N →∞ 1-e jN ω N (1-e jω ) = 0 for any fixed ω = 0, lim N →∞ K ρ,ω,N N -n = 0 when ω is chosen from nonequal atoms. (52) (with additional 1/N -factors) thus becomes

where lim N →∞ a ρ,N = 0. Multiplying both sides with N and letting N go to infinity gives

It is clear that this converges to 0 when ρ = 0 n (since |ρ| < n in this case), so that the limit is

which proves the claim APPENDIX K THE PROOF OF THEOREM 5 We need the following identity [START_REF] Rottmann | Mathematische Formelsammlung. B.I. Wissenschaftsverlag[END_REF]:

where sgn(x) = 1 if x > 0, sgn(x) = -1 if x < 0, and 0 otherwise. From this it follows that

Note that the measure with density p, has the same asymptotics near α i as the measure with density p i |x -

As in the proof in Appendix J, the integral for the expansion coefficients is dominated by the behaviour near the points (α i , ..., α i ). To see this, note that the behaviour near the singular points on the diagonal is O (s(|ρ| -n) -1) when polynomic growth of order s of the density near the singular points is assumed. This is very much related to (93) in Appendix H, since K ρ,ω here in a similar way can be bounded by (taking into account new powers of N )