Limiting spectral distribution of the sample covariance matrix of the windowed array data

In this article, we investigate the limiting spectral distribution of the sample covariance matrix (SCM) of weighted/windowed complex data. We use recent advances in random matrix theory and describe the distribution of eigenvalues of the doubly correlated Wishart matrices. We obtain an approximation for the spectral distribution of the SCM obtained from windowed data. We also determine a condition on the coefficients of the window, under which the fragmentation of the support of noise eigenvalues can be avoided, in the noise-only data case. For the commonly used exponential window, we derive an explicit expression for the l.s.d of the noise-only data. In addition, we present a method to identify the support of eigenvalues in the general case of signal-plus-noise. Simulations are performed to support our theoretical claims. The results of this article can be directly employed in many applications working with windowed array data such as source enumeration and subspace tracking algorithms.

The distribution of the eigenvalues of the sample covariance matrix (SCM) of data has important impact on the performance of signal processing algorithms. Over the last decade, the properties of complex Wishart matrices are used in the analysis and design of many signal processing algorithms such as in array processing. Our knowledge about the distribution of eigenvalues, eigenvectors and determinants of complex Wishart matrices and their limiting behavior is emerging as a key tool in a number of applications, e.g., in data compression and analysis of wireless MIMO channels [1,2], array processing, source enumeration and identification [3-5], adaptive algorithms [6,7]. The densities of the singular values of random matrices and their asymptotic behavior (as the matrix size tends to infinity) has been employed in some applications [8-10]. The eigenvalues of the SCM are often used to describe many signal processing problems. For example in [8], they are used as sufficient statistics for array source enumeration.

Let X₁, …, X_N be N independent zero mean Gaussian random vectors with covariance matrix of A, i.e., N M ( 0 , A ) , where A is a nonnegative M × M Hermitian matrix. The SCM R_N is defined as R N = 1 N ∑ i = 1 N X i X i H = 1 N XX H , where X = [X₁,…,X_N] contains N snapshots of the received data. In this article, we refer to this SCM as the SCM with rectangular window (SCM-R) as all data samples have equal weights, i.e., a rectangular window is used. In this case R_N has a Wishart distribution [11] and for more than four decades, it has been known that the joint probability density function (PDF) of its eigenvalues, can be expressed in terms of hyper-geometric functions [12]. More recently, a simpler form of this joint PDF was derived in terms of the product of two determinants [13]. However, this form is applicable if the array is small and the eigenvalues of the covariance matrix of the observed data are distinct. Several articles have investigated the behavior of the eigenvalues of R_N when M, N → ∞ assuming M N → c > 0 [14,15]. This is a more realistic assumption than assuming M is finite and N is infinite, because in most practical applications the covariance matrix A slowly varies, hence, the effective window length could not be arbitrary long. For instance, the eigenvalue estimators that are consistent in this asymptotic regime are more robust to finite sample size than other estimators which are only guaranteed to converge for fixed M and N→∞[9]. There are many works on the distribution of eigenvalues in this asymptotic regime, such as information-plus-noise [16] and spiked models where all eigenvalues are equal excluding a small number of fixed eigenvalues (spikes) [17]. Specifically, the distribution of the largest noise eigenvalue is widely studied [18,19].

Some signal processing algorithms process a batch of data together and deal with the SCM-R. In addition, the existing results in literature about the behavior of the eigenvalues mainly consider the rectangular window. However in a number of practical signal processing algorithms, the SCM is estimated by applying a window as follows

where {w_i ≥ 0,i = 1,…,N} is a non-negative sequence. Hereafter, we refer to R_N as the SCM. The SCM-R is obtained using a rectangular window, i.e., where w_i is non-zero and constant for i = 1,…,N. These weights allow to flexibly emphasize or deemphasize some of the observations. For example smaller weights for old data samples allows to improve the agility of the algorithms. For instance in cognitive radio, it is important to detect the activities of users and the idle channels as fast as possible, thereby reducing the detection time and improving the agility of the system [20,21]. Among all windows, the exponential window, w_i = w₀pⁱ, is commonly used. Two reasons for this popularity are (1) this window allows to develop fast recursive algorithms which are considerably less expensive in terms of computational complexity, thereby facilitate the real-time implementation of these algorithms (e.g. see [22,23]) and (2) allows to forget the old data, thereby improving the tracking ability in non-stationary environments. For instance exponentially windowed data is used in most of the existing subspace tracking algorithms[24,25]. That is because only a rank-one update is required for each new data vector to update the underlying SCM, which leads to simple low cost subspace tracking algorithms.

In this article, we study the effects of windowing on the distribution of the eigenvalues of the SCM. In this case, the SCM in (1) has a doubly correlated Wishart distribution [26-30]. We must note that, there are numerous research results for the case of Wishart matrices, however, the spectral properties in the doubly correlated case has not been sufficiently studied.

Manipulating the joint PDF of the eigenvalues which is a very complex function is not practical, particularly for large matrices. An alternative approach used in the literature, is to employ the following empirical spectral distribution (e.s.d.) of a square matrix A ∈ C M × M

where λ₁, λ₂, … ,λ_M are eigenvalues of A and #{.} denotes the cardinality of a set. Note that, in this definition all eigenvalues of A are assumed to be real. Although this formulation is less explicit than the joint PDF of eigenvalues, it describes the statistical behavior of the eigenvalues. In many practical cases A is a random matrix and the e.s.d. F^A(x) is a random function which converges almost surly to a deterministic cumulative distribution function as the dimension of the system grows. In such cases, lim M → ∞F^A(x) is referred to as the limiting spectral distribution (l.s.d.) of A.

In recent years, some results have been obtained on the limiting behavior of the e.s.d. of correlated Wishart matrices. In this article, for the white noise case, we study the behavior of eigenvalues of the SCM. In particular for the exponential window, we extend the results previously demonstrated in [31] and give more details along with the proofs of the required theorems. We then consider the case of signal plus noise and present a method to determine the support of eigenvalues. The main contributions in this article are

• A method is proposed to approximate the spectral distribution of the SCM using arbitrary windows with that of an equivalent Wishart Distribution. For the especial case of white noise (noise only), this approximation is the Marchenko–Pastur (M–P) distribution, which is the known distribution for the case of a rectangular window.

• In Theorem 2, we derive an accurate and explicit equation for the l.s.d. of the SCM of noise-only data for the exponential window. Many simulations are performed to show the accuracy of this l.s.d.

• In Theorem 3 we present a systematic method to compute the support of eigenvalues in the signal plus noise data case using an exponentially weighted window. In addition to the results, we follow up a different and novel approach in proving this theorem compared with the existing proof for the rectangular window case where the Stieltjes transform m(z) has the explicit inverse [15]. This approach can be easily utilized for other window types where the Stieltjes transform is expressed explicitly or implicitly as a function of z.

The demonstrated results provide a key step toward characterization of the distribution of eigenvalues in the general Covariance matrix of windowed data. The results of this work are useful in the design and implementation of robust algorithms using windowed snapshots. Our derivations in Theorems 2 and 3 can be directly used to design unbiased eigenvalue and eigenvector estimators. These estimators are important especially because the exponential window is used in numerous applications. They can be used as a basis to improve the performance and accuracy of many existing algorithms which are based on exponentially windowed data, in many fields such as subspace tracking, DOA estimation and source enumeration.

The remainder of this article is organized as follows: Section 2 introduces the system model and some important mathematical tools. We derive an approximation for the Stieltjes transform of l.s.d. of eigenvalues of weighted windowed array data in Section 3. Asymptotic spectrum of the eigenvalues in noise-only data case is analyzed in Section 4. The signal plus noise case is studied in Section 5. Section 6 provides simulation results. Finally, we conclude this work and suggest future works in Section 7.

We assume that X i ∈ C M in (1) is a circularly symmetrical independent Gaussian random vector process with zero mean and covariance matrix of A ∈ C M × M , i.e., X i ∼ N M ( 0 , A ) . In this case, we can rewrite (1), supposing that SCM is estimated using a window of size N with positive coefficients w₁,…,w_N,

where U = [U₁,…,U_N] is an M × N matrix contains i.i.d. zero-mean unit-variance complex Gaussian entries and W N = Δ diag ( w 1 , … , w N ) . The matrix R_N has a doubly correlated Wishart distribution. In practice, it is very complex to directly characterize the e.s.d. of R_N thus, we use the Stieltjes transform of this distribution and indirectly characterize the behavior of the eigenvalues. Then, in the asymptotic regime as M,N→∞ given M N → c > 0 , the inverse transform of the limit gives the l.s.d of SCM.

Definition 1

[15]Stieltjes transformm(z), z ∈ C + ≡ { z ∈ C : Im ( z ) > 0 } of a distribution functionF^R(x) is defined as

The inverse Stieltjes transform formula is as follows:

Hence, in order to characterize the asymptotic distribution of the sample eigenvalues, we alternatively characterize the asymptotic behavior of the corresponding Stieltjes transform, and then use the Stieltjes inversion formula in (5) to obtain l.s.d. of SCM f^R(x). We use the following theorem which gives the Stieltjes transform of the correlated Wishart matrix [29] and is the basis for derivations in this article.

Theorem 1

For a finite length window with length ofN, consider the matrix defined by R N = 1 N A N 1 2 U W N U H A N 1 2 . Assume that all elements of U ∈ C M × M are i.i.d. random variables with zero-mean, unit variance and finite E {|U_ij|⁴}. In addition, suppose that A N ∈ C M × M is a Hermitian nonnegative definite matrix,W_N = diag(w₁, …, w_N), F A N → D F A , F W N → D F W whenM,N → ∞ with M N → c > 0 . In this case, the empirical distribution F R N , with probability 1, converges weakly to a probability distribution functionF^Rwhose Stieltjes transform m(z), for z ∈ C + , is given by

where e(z) is the unique solution of the following equation in C +

Proof 1

See[29]for proof. Similar results are also demonstrated in[26], [28]with some differences in the assumptions on correlation matrices.□

We emphasize that (6) and (7) give the exact distribution in the asymptotic regime as M,N → ∞ with M N → c > 0 . Since in practice, the array dimension and/or sample size are usually finite numbers, this method gives a deterministic approximation for the actual sample eigenvalue distribution.

To show how this method works, we now consider the simplest case (where the distribution is well known) using a rectangular window and white Gaussian noise, i.e., W = I_N × N and A = σ²I_M × M. In this case, we have dF^W(w) = δ(w − 1)dw and dF^A(x) = δ(x − σ²)dx, where δ(x) is the Dirac delta function. Thus with straightforward manipulations of (6) and (7), the Stieltjes transform is found to be the solution of

In this case, as expected the e.s.d. of the SCM-R, F R N ( x ) , converges to the M–P distribution [14] as follows,

where a ± = σ 2 1 ± c 2 and π a , b ( x ) = 1 a ≤ x ≤ b 0 otherwise .

Now, let us consider an arbitrary window and white noise A = σ²I_M×M. In this case from (6), (7) and dF^A(x) = δ(x − σ²)dx, we obtain m ( z ) = 1 σ 2 e ( z ) . Thus, from (6), we obtain

In this section, we define the effective length of a window which allows to approximate the distribution of the eigenvalues of windowed SCM with that of a rectangular window with an equivalent length, assuming that the covariance matrix of data A satisfies the assumptions of Theorem 1. In several existing articles some intuitive equivalent length are defined simply to extend the previously existing results for the rectangular case in order to analyze the behavior of the eigenvalues in the weighted window cases [22], [23].

Consider a window w_i>0 of length N and denote W N = diag ( w 1 , … , w N ) ∈ R N × N with a converging distribution, i.e., lim N → ∞ F W N = F W . We assume that the sample size N is much larger than the array dimension M, i.e., c is small. It is known that m(z) is bounded for z ∈ C⁺[15]. Thus for 0 < c ≪ 1 we have 0 < cσ²|m| sup |w| < β < 1 where β is some constant number. It is easy to show that^a, we have σ 2 w 1 + cw σ 2 m − 1 cm ∑ i = 1 I ( − cw σ 2 m ) i ≤ | cm | I | w σ 2 | I + 1 1 − β . This yields

where E { . } = ∫ ( . ) d F W ( w ) and | O | ≤ | m | I | σ 2 | I + 1 1 − β E { w I + 1 } .

Since 0 < c ≪ 1, for I = 2 and defining c e = c E { w 2 } E 2 { w } and w_e=E{w} as the effective parameters, we can rewrite (11) as

where using E{w³} < sup{w²}E{w} and E{w²} < sup{w}E{w} it is easy to show that the approximation error is bounded by σ 2 E { w } 2 β 2 1 − β .

Definition 2

The expression (12) represents the M–P distribution as in (8) for a rectangular window of length

with all coefficients equal tow_e. The average weightw_eis a scale parameter for the eigenvalues of covariance matrix of the received data. Although we have derived the effective length for the noise only data, our results reveal that this effective window length gives accurate results for the signal plus noise case.

For the white noise data, the l.s.d. of SCM can be approximated by the M–P distribution defined in (9) by substituting c and σ², with c_e and w_eσ², respectively. Note that the effective window length is always smaller than the number of samples N. This approximation can be intuitively interpreted as a Wishart approximation where the effect of “windowing” is approximated with a rectangular window with an effective number of samples of N_e and the covariance matrix of the received data is scaled to A_e=w_eA”.

Now, we compute the effective length of the triangular and the exponential windows. A triangular window is defined by w i = 2 ( 1 − i − 1 N − 1 ) for i = 1,…,N ≫ 1 and has the average weight of 1 N ∑ i = 1 N w i = 1 . Using (13), the effective length of the triangular window is

The exponential window is very popular in signal processing applications due to its simple implementation and is defined by w_i=w₀pⁱ for i=1,2,…, where p∈(0,1) and w₀ is a normalization constant. We note that the exponential window is inherently an infinite length window. Interestingly, in Theorem 1 the window length and the array dimension jointly tend to infinity where lim M , N → ∞ M N = c > 0 . Here for a finite array dimension M, we first approximate the exponential window with N coefficients, which is only accurate if N is large enough such that the omitted coefficients are negligible. Asymptotically as M,N jointly tend to ∞, the results from this truncated window become accurate for describing the underlying distributions using the exponential window. In this case, with some calculations we obtain N e = 1 − p 2 ( 1 − p ) 2 ( 1 − p N − 1 ) 2 1 − p 2 N − 2 . Thus the effective length of the exponential window (for N≫1) becomes

which is not a function of N. As expected the effective length of the window increases as the forgetting factor p approaches one.

In this section, for the windowed data case, the l.s.d. of the SCM is characterized more accurately. In practice, the array dimension and the effective window length are both finite. However, we are interested in the impact of the weights of the window f^W(w), on the limiting distribution of the eigenvalues as M,N→∞ employing Theorem 1. We use two approaches to model f^W(w), Discrete and Continuous. The former considers f^W(w) as a finite sum of discrete masses at the coefficients of the window. The discontinuous distribution function modeling is useful to analyze the support of eigenvalues and its connectivity. The latter approach, approximates f^W(w) as a continuous function allowing to derive some explicit equations for the Stieltjes transform.

Let S_F denote the support of the function F^R(x) and S F c shows its complement. From (5), we see that S_F consists of points on the real axis where the imaginary part of m(z) is positive, i.e., the support is the union of some subintervals. Thus to find the support of the distribution of eigenvalues, we must determine such intervals. In [15], it is shown that limy→0⁺m_F(x+iy) exists for all x≠0, and therefore we can define

The following lemma is the key to determine these intervals on real axis [15].

Lemma 1 ([32], Lemma 6.1)

For any c.d.f.F, letS_Fdenote its support and S F c be the complement ofS_F. For x ∈ S F c ,m=m_F(x) is the only real solution ofx = z (m) which satisfies

wherez(m) is the inverse function ofm(z). Also conversely, for any realmin the domain ofz(m) if dz ( m ) dm > 0 thenx = z (m) is outside the support ofF.

This simply means that the support S_F, is the union of intervals on the vertical axis where z(m) is increasing for real values of m. According to (10), for noise only data z(m) can be written as follows

4.1 Discrete distribution function approach

Suppose the window consists of N_d distinct weights w_i,i=1,…,N_d, each with multiplicity n = N N d . Therefore as (M,N→∞ and M N → c > 0 ), we can evaluate (18) in terms of the weights

Figure 1, represent a typical case of the function on the right-hand side of (19). Lemma 1 states that the support of the distribution of eigenvalues is the complement of the set of all values x ∈ R⁺ for which x = z(m) is increasing for real values of m, i.e., ( dz ( m ) dm > 0 ) . The function z(m) has poles at m=0, − 1 c w 1 σ 2 , … , − 1 c w N d σ 2 . In addition, z(m) is an analytic function and we have

For c < 1, i.e., where the length of the window is more than the array dimension, lim m → ± ∞ z ( m ) = 0 ± , thus as Figure 1 shows, dz ( m ) dm = 0 has at least two solutions which we denote them m u ∈ − 1 c σ 2 max ( w i ) , 0 and m_l∈(0,∞). For c > 1, as Figure 2 shows typically, from lim_m→±∞z(m) = 0^∓ we conclude that m_l must be in − ∞ , − 1 c σ 2 min ( w i ) . We must note that, for c > 1, the SCM has M − N zero eigenvalues expressed with a probability mass of 1 − 1 c in the l.s.d. of SCM, that is not counted as a cluster in these derivations, i.e. for c >, the PDF of the distribution includes a term of 1 − 1 c δ ( x ) . If the weights are widely separated, the support of eigenvalues may become fragmented into union of a number of disjoint intervals.

In many signal processing applications the white noise subspace is separated from the signal subspace based on the eigenvalues of the SCM. Such a fragmentation of the support of noise eigenvalues misleads the subspace based algorithms and leads to noise eigenvalues to be mistaken as signal ones.

In fact, it is desirable that the support of eigenvalues be as compact as possible. To avoid such an undesirable fragmentation, the equation dz ( m ) dm = 0 should not have real solution for m ∈ − 1 c σ 2 min ( w i ) , − 1 c σ 2 max ( w i ) , i.e.,

Under this connectivity condition, the support of eigenvalues is the interval [x_l = z(m_l),x_u=z(m_u)], which can be calculated, numerically. Our simulations show that this condition is satisfied for popular window types especially for N_d ≫ 1 used in practice. Figure 2 shows a typical case for c > 1 where dz ( m ) dm = 0 has an even number of real-valued solutions (counting multiplicities) which we denote them by m 1 − ≤ m 1 + < ⋯ < m q − ≤ m q + (in addition to m_l,m_u). Each pair of these solutions determines a sub-interval for the support of eigenvalues, i.e., we have S F = [ x l , x u ] − { [ x 1 − , x 1 + ] ∪ ⋯ ∪ [ x q − , x q + ] } , where x i − = z ( m i − ) , x i + = z ( m i + ) . Reducing c or reducing the gap between weight values {w_i} makes the support more compact at the expense of using more temporal samples.

4.2 Continuous function approach

The goal of this approach is to find closed form expressions of Stieltjes integrals of the l.s.d. This approach could be used for any window shapes. However, we start with the triangular window and then consider the exponential window which are more popular. Here, we model the function f^W(w) with a continuous distribution and evaluate (18) to found the Stieltjes transform.

For a triangular window w i = 2 1 − i − 1 N − 1 , i = 1 , … , N , we have F W N ( w ) = 1 N ∑ i = 1 N U ( w − w i ) where U(w) is the unit step function. In this case it is easy to show that F W N ( w ) converges to a uniform distribution as N increases, i.e.,

Substituting f^W(w) in (18), we get

for m ∈ ( − 1 2 c σ 2 , ∞ ) and m≠0.

Again, we first use Lemma 1 and determine the support of eigenvalues (by plotting z(m) for real m and finding the intervals on the vertical axis where z(m) is not increasing). Figure 3 plots the lower and upper boundaries of support of eigenvalues for a triangular window for different values of c. It can be seen that the discrete distribution F W N ( w ) (assuming N_d=50) and the continuous approach result in almost the same boundaries. We also observe that these boundaries are close to those obtained by the Wishart approximation assuming the effective window length in (14). In this figure using the rectangular window with same length as the triangular window, the distribution is referred to as the M–P distribution. Also we observe that the eigenvalues tend to more concentrate around their real value σ²=1 as the window length increases. In addition from this figure, we conclude that the support using the triangular window is looser than than that of the rectangular window for a given value of c, because the effective length of the triangular window is less than that of a rectangular window.

For the exponential window, first we introduce the new parameter γ as the ratio of smallest to largest weights of the truncated exponential window. The coefficients of the window can be redefined as as a function of γ as w i = w 0 γ i N , i = 1 , … , N . Therefore from f W N ( w ) = ∑ i = 1 N 1 N δ ( w − w i ) , from i = N ln γ ln w i w 0 , it is easy to show that

where ⌊.⌋ is the floor function. This increasing staircase function takes values on 0 , 1 N , 2 N , … , 1 . To satisfy the constraints of Theorem 1 for the exponential window, we assume that the ratio of smallest to largest weights of the window, γ=p^N > 0, is an arbitrary small real constant. In other words, the forgetting factor of the window p = γ 1 N ∈ ( 0 , 1 ) approaches to 1, as M,N→∞. The smaller γ, the better this truncated exponential model fits the exponential window with the forgetting factor p. From, lim N → ∞ w 0 = ln γ γ − 1 , we conclude that lim N → ∞ F W N = F W ( w ) where

is a continuous function, independent of window size N and satisfies the assumptions of Theorem 1. Thus, this theorem is applicable to the exponential window truncated at some large integer N.

Substituting f^W(w) in (18), in the asymptotic regime of Theorem 1 as γ→0, such that M n 0 → c 0 , z(m) satisfies

for all m ∈ − 1 c 0 σ 2 , ∞ ∖ { 0 } where n 0 = − 1 ln ( p ) .

One can use the same method as in the discrete distribution function approach and identify the support of the distribution S_F. However, the function z(m) in (26) is simple and the following theorem gives the explicit distribution.

Theorem 2

For the exponentially weighted window, the l.s.d. of SCM,f^R(x), is given by

and upper and lower boundaries of the support are

respectively, whereω_k(x) is the branch of Lambert W function^b[33]withk=−1 andk=0.

Proof 2

According to the Lemma 1, boundaries of the support of eigenvalues are the real solutions ofz^′(m) = 0, i.e., with some simple calculations, are the solutions of

Denotingy= ln(1+c₀σ²m)−c₀−1, we obtain

This equation has two real solutionsm₋ and m₊ expressed using Lambert W function as

Using (26), the boundariesz(m₋) andz(m₊) are obtained as in (28a) and (28b) which determine the support of eigenvalues as the interval [ z ( m − ) , z ( m + ) ] ⊂ R .

To obtain the l.s.d. of SCM, we should findm(z) with positive imaginary part for allz∈[z(m₋),z(m₊)]. In (26), we denote

and obtain

Therefore, the solutions are

According to (16), (33), for the values ofzin the interval of the obtained support on the real axis, due to the properties of the complex logarithm function, the imaginary part ofvis in [−Π,Π], thus only the branches withk = 0 andk = − 1 are acceptable solutions. It is easy to see that forz ∈ [z(m₋),z(m₊)], the expression on the right-hand side of (34) belongs to − e c 0 − 1 , − e − 1 . From (33) and properties of Lambert W function, we also deduce that Im{m} and sin(−Im{v}) have the same signs, and for x ∈ − e c 0 − 1 , − e − 1 the function sin(−Im{ω_k(x)}) is positive fork = − 1 and is negative fork = 0. Therefore, the Stieltjes transform of the l.s.d. of SCM is obtained from (35) and (33) as

Using the inverse formula in (18), the l.s.d of SCM is

Dropping the real terms inside the brackets and applying some simplifications, we obtain (27). □

We can define a second effective window length by employing and comparing the boundaries of the support of eigenvalues in (28a) and in (9) for a rectangular window which is only in terms of σ² and c. Equating the length of the supports in (9) (28a), i.e., a₊−a₋ = x₊−x₋, we can find a rectangular window to match the support as same as that of the exponential window and define the length of this rectangular window as another effective length for the exponential window. In some array signal processing applications, the effective length of the exponential window has been considered to be N e = 1 1 − p [22], [23]. Figure 4 compares these effective lengthes in terms of the forgetting factor p and reveals that the effective window length defined in (13) gives an accurate approximation for the exponential window. We also see a large gap between the traditional approximation for the effective length in [22] and what is obtained in this article using random matrix theory.

Remark 1

In the economic literature, other methods have been proposed to approximate the spectral density function of exponentially weighted financial covariance matrices for Portfolio Optimization[34], [35]. These methods that are used in other articles (e.g., in[36], [37]) are based on numerical calculations rather than developing some closed form expressions. Pafka et al.[34]supposed that the density of the eigenvalues is aproximated by ρ ( u ) = Qν Π where Q = 1 M ( 1 − p ) andνis the root of

In contrast to these methods for the exponential window, we derive an accurate explicit closed form expression which can be easily employed in many applications such as in signal processing and economy.

In this section, we consider the case of white noise plus some signal sources, i.e., where the eigenvalues of A are not equal. In the general case, let λ_q > ⋯ > λ₁ >0 denote the set of q distinct eigenvalues of the covariance matrix and the multiplicity of λ_ℓ is denoted by k_ℓ (we must have M = ∑ ℓ = 1 q k ℓ ). For example suppose a real phased array communication system with q − 1 independent signals impinging on it simultaneously on the same frequency band from different directions where q < M. The smallest eigenvalue λ₁ can be interpreted as the noise eigenvalue and other q − 1 larger eigenvalues are referred to as signal eigenvalues. In the asymptotic regime, when N,M are growing large, we assume that k ℓ M → α ℓ > 0 , where α_ℓ,ℓ= 1, …,q are multiplicity ratios of eigenvalues. In this case the spectral distribution of the matrix A in Theorem 1 can be expressed as sum of Dirac delta functions, i.e. d F A ( a ) = ∑ i = 1 q α i δ ( a − λ i ) da .

In what follows, we present an approach to determine the support of eigenvalues and also the l.s.d. of exponentially weighted SCM of signal plus noise data in the asymptotic regime. The first in determining the distribution of the eigenvalues is to determine its support on the real positive axis.

The definition of the Stieltjes transform in (4) implies that for any distribution F and real x outside the support of F, m(x) is well defined and its derivative, m ′ ( x ) = ∫ dF ( y ) ( y − x ) 2 , is obviously real and positive. Thus, m(z) is increasing on intervals on real line outside the support of its distribution function F[15]. Therefore, the inverse function theorem proves that its inverse exists on these intervals and shall also be increasing. For the one sided correlated Wishart matrices, where the inverse of m(z) has an explicit expression, Lemma 1 shows that the converse of the above statements are also true [15], i.e. for any real m in the domain of z(m), if dz ( m ) dm > 0 then x = z (m) is outside the support of the distribution. Therefore, the support of eigenvalues is a Borel subset of R + for which z(m) is increasing which can be determined by simply plotting the inverse function z(m) for real m. Paul and Silverstein ([29], page 2) suggested the same method for doubly correlated Wishart matrices if there exists an explicit inverse z = z(m) for the limiting Stieltjes transform m(z). Unfortunately, for non-rectangular windows, the inverse of m(z) in general has no explicit expression [29]. Fortunately, by introducing two auxiliary variables u and h in what follows for the exponential window, we found z(h) which implicitly expresses z as a function of m. Then, we prove that the same method can be extended for the exponential window case, while the main difference here is that we are able to use the implicit expressions to determine this Borel set. Although the exponential window case is studied in this article, the same approach may be used for some other window types, to determine the support of eigenvalues.

From (6) and (7) we obtain

Substituting (25) in (39), we get

Substituting dF^A(a) in (7) changes the integral to a summation and we obtain e(z) as

According to Theorem 1, for any z ∈ C + , there is a unique solution e=e(z) for (41) in C + . In this case the Stieltjes transform m(z) is calculated from (40) as

This expression gives the implicit relation between m and z, which cannot be sorted to express m as an explicit function of z or conversely, z as a function of m. Defining the auxiliary variable/function u = c₀ (1+zm(z)) which provides a bijective relation between e and m for all z≠0, we have

This equation reveals that the imaginary parts of u and e have the same signs. In addition since γ <1, c₀ is real and using the properties of the complex logarithm function in (43), we deduce that u always lies in a strip of the positive complex plane where its imaginary part is less than Π, i.e., the domain of u is defined as D_u = {u|0 < Im{u}<Π}. Equation (43) also provides a bijective relation between e and u, therefore according to Theorem 1 for any z ∈ C + , there is a unique u∈D_u, satisfying

Defining the second auxiliary variable/function as

and define D_h as its range for all z ∈ C + . Resorting (44), we have

Proposition 1

The auxiliary variableh, as a function ofuandz, has some interesting properties as:

(1) halways lies in the subset D h ⊂ C + for all z ∈ C + .

(2) forh∈D_h, zcan be explicitly expressed as a function ofh

(1) For any z ∈ C + , a uniquehsatisfying (47) exists in

Proof 3

The first property can is simply implied from (46) as the imaginary part ofh and uhave the same sign. Using (45) and (46), we can easily find (47). The third property is proved as follows. The constraint in (48) is obtained from Im{u} ∈ (0,Π) and (46). According to Theorem 1, for any z ∈ C + , there is a uniqueu∈D_u, satisfying (44). The unique pair (z,u) gives anhin C + according to (45). In order to prove the uniqueness ofh, suppose thath₁andh₂in C + satisfy (47) and (48). Thus, (46) yieldsu₁,u₂∈D_usatisfying (44). In addition, we must haveu₁ = u₂since for any z ∈ C + , there exists a uniqueu₁∈D_u. Thus forzandu₁ = u₂, (45) yields thath₁ = h₂. □

Although z(h) in (47) is defined only for h∈D_h, it is an analytic function for all h ∈ C ∖ 0 , − 1 λ 1 , … , − 1 λ q . In addition note that z ( h ) = − 1 h at the roots of ∑ i = 1 q α i λ i h + 1 − 1 = 0 .

Also using (46) and (47) we express m as a function of h as follows

for h∈D_h. Similar to z(h), the complex function m_h(h) is an analytic function for all h ∈ C except at the set of real values 0 , − 1 λ 1 , … , − 1 λ q and the points where z(h) = 0.

The inverse Stieltjes transform in (5) reveals that the l.s.d. depends on the behavior of m(z) in the vicinity of the real axis, i.e. for z → x 0 ∈ R . Proposition 1 shows that the z(h) in (47) is injective over h ∈ D h ⊂ C + and allows us to treat h(z) as its inverse for z ∈ C + . To determine the range of h(z) denoted by D_h we can evaluate z(h) for all h in (48) and take only those values of h for which z is in C + . As an example Figure 5 shows this region for c₀ = 0.2 and a covariance matrix with two distinct eigenvalues 2 and 1 with the multiplicity ratios α 1 = 1 3 and α 2 = 2 3 . The white regions in Figure 5 shows the values of h for which z(h) has negative imaginary part, and the blue parts are the values where Im {z(h)} > 0. We observe that some parts of the positive complex plane are not in the domain of z(h) as we restrict the range to z ( h ) ∈ C + .

5.1 Support of eigenvalues

Theorem 3

For the exponentially weighted window defined in Theorem 2, under the assumptions of Theorem 1, the complement of support of eigenvalues, is the set of values ofx = z(h) on the vertical axis where dz ( h ) dh > 0 for some h ∈ R , wherez(h) is defined in (47).

Proof 4

LetS_Fdenotes the support of the functionF^R(x) and S F c shows its complement. To prove Theorem 3, first we show that for any x ∈ S F c , there exist a h ∈ R where dz ( h ) dh > 0 . Then, we prove the converse, i.e. if dz ( h ) dh > 0 for some h ∈ R , thenx=z(h) is a real number outside the support of eigenvalues.

From (5), we see thatS_Fconsists of points on the real axis where Im {m(x+iy)} tends to a positive number wheny → 0⁺. Thus to findS_F, we must determine such subintervals on the real axis, or equivalently we can determine S F c by finding the intervals on the real axis where limy → 0⁺{m(x+iy)} is real. Consider any x₁ and x₂ such that ( x 1 , x 2 ) ⊂ S F c ⊂ R + . According to the definition of Stieltjes transform in (4), m(z) and u ( z ) = c 0 ( 1 + zm ( z ) ) = c 0 ∫ λdF ( λ ) λ − z are both real and well defined for any z ∈ (x₁, x₂). In addition du dz = c 0 ∫ λdF ( λ ) ( λ − z ) 2 is nonnegative on this interval. Thus u(z) is a real invertible function on (x₁,x₂), and its inversez(u) is also real and increasing on the interval ( u ( x 1 ) , u ( x 2 ) ) ∈ R , i.e. dz du > 0 .

Lemma 2

For any given z ∈ R + , the function h(u,z) in (45) is monotonically increasing versus u ∈ R .

Proof

Defining h ( 0 , z ) = γ − 1 z = lim u → 0 h ( u , z ) , the function h(u,z) is continuous for all u ∈ R and all z ∈ R + , and for all z ∈ R + we have □

Since dz dh = dz du du dh and du dh = c 0 M ∑ i = 1 q k i λ i ( λ i h + 1 ) 2 are positive for all h ∈ R , Lemma 2 implies that the signs of dz du and dz dh are identical. Thus ifzis an increasing function of u ∈ R , it is also an increasing function of h ∈ R as well, and vice versa, i.e., the intervals for which z is increasing versus u is equal to the intervals for whichzis increasing versus h. This proves the direct part of the theorem.

To prove the converse part, consider that Theorem 3 implies that dz ( h 0 ) d h 0 is real and non-negative for some h 0 ∈ R . Sincez(h) and m_h(h) are both real at pointh = h₀, it is sufficient to show that the pointh₀ belongs to the boundary ofD_h. In this case, as the functionm(h) is continuous in the complex plane (excluding few points as stated after (49)), we conclude that lim_{y → 0⁺} Im {m(h₀+iy)} = Im {m(h₀)} = 0. To show that h₀ is on the boundary of D_h, we prove that the points in the vicinity of h₀ in the positive complex plane, belong to D_h. Let {h_n} be any complex sequence with positive imaginary part converging toh₀ as N → ∞. Sincez(h) is continuous, the sequence {z_n} = {z(h_n)} exists and converges toz(h₀).

Lemma 3

Let z(h) be an analytic function ofh over an open set G, and h(t) ∈ Gbe a differentiable curve at t. Then if dz ( h ) dh is a positive real number, we have arg d dt z ( h ) = arg { h ′ ( t ) } .

Proof 5

This lemma is obtained from the Chain rule; since z(h(t)) is differentiable at t and d dt z ( h ( t ) ) = h ′ ( t ) z ′ ( h ) . Thus for positive real dz ( h ) dh , the argument of d dt z ( h ) and h^′(t) are the same.□

We use Lemma 3 which implies that if dz ( h ) dh is positive and real at the pointh = h₀ then arg d dt z ( h ) = arg { h ′ ( t ) } for any differentiable curve h(t) at t = t₀ where h(t₀) = h₀, i.e. the slope of the curveh(t) in the complex plane is the same as the slope of z(h(t)) ath(t) = h₀. Now for sufficiently large n, consider the line L_n = h(t) = (1 − t)h₀ + th_n, 0 ≤ t ≤ 1 in C + which originates from h₀and ends at h_n. The transformation of L_n, z(L_n), is also a line in the positive imaginary part of complex plane with the same slope as L_n, as we have supposed that dz ( h ) dh ≈ z ′ ( h 0 ) for the points on L_n. Thus the point z_n also lies in the positive complex plane. In the other words, for sufficiently large n, the sequence {z_n} lies in C + ; hence the sequence {h_n} is inD_h. Finally, we conclude that the Stieltjes transform is defined on any such sequences and the sequence of Stieltjes transform {m(z_n)} = {m_h(h_n)} is also in C + for those values of n and converges to m_h (h₀) which is a real number. Thus z(h₀) is outside the support of eigenvalues.□

Remark 2

We must note that we use a different approach in proving Theorem 3 comparing with proof exists for the rectangular window case where the Stieltjes transformm(z) has the explicit inverse[15]. This approach is very simple and can be used in other cases where the Stieltjes transform is expressed explicitly or implicitly as a function of z.

Theorem 3 states that in order to find the support of eigenvalues, we could first find the intervals on the real line where z(h) is increasing. In a sufficiently small vicinity of these intervals on the positive imaginary part of the complex plane, it is discussed in the proof that the imaginary part of z(h) is also positive for all h in this vicinity, therefore this vicinity lies in D_h. Having a closer look at Figure 5, we find that D h ⊂ C + approaches real axis only for some values of h which can be easily studied that these are the intervals for which z^′(h)>0. Thus according to this theorem the support of eigenvalues consists of three disjoint intervals for the setting of Figure 5.

Employing Theorem 3 and plotting z(h) for h < 0 one can determine the support of eigenvalues of the SCM in the asymptotic regime. The function z(h) has asymptotes at − 1 λ 1 , … , − 1 λ q with the following one-sided limits

Figure 6 shows a typical representation of the support of eigenvalues in the signal plus noise case when c₀ = 0.1 and the covariance matrix has four distinct eigenvalues 5, 3, 2, 1 with multiplicities α₁ = α₂ = α₃ = 0.1 and α₄ = 0.7. It can be studied that in general, z(h)→ + ∞ as h→0⁻ and z(h)→0⁺ as h→−∞ and also analogous with the rectangular window case [38] the number of extrema of z(h) (counting the multiplicities) is even and are the solutions of dz dh = 0 . Generally, in order to determine the support of eigenvalues, we identify all intervals on the vertical axis where z(h) is increasing and in general case denote them by S F , b c , b ∈ { 1 , … , s } . Removing these intervals from R , what is left is S_F and according to the proof of Theorem 3. these intervals will not overlap each other. To see this, we note for each x ∈ S F c , there is a unique h∈D_h, such that x=z(h). Assume that I H , b c , b ∈ {1,…,s} are the subintervals in the h domain where z(h) is increasing. Therefore, I H , b c uniquely determines S F , b c , which is an interval in S F c . The complement of these intervals are the points determine the support of eigenvalues. It can be seen in Figure 6 that the support of the distribution is the union of four clusters where each of them represents the support of the distribution of only one of the eigenvalues. This is analogous with the results proven in the literature for rectangular windows [38], i.e., in this case all eigenvalues are separable on the vertical axis.

Figure 7 illustrates the same curves for c₀ = 0.4, i.e., the forgetting factor p is reduced compared with Figure 6. We observe that the smaller the forgetting factor of the exponential window the larger the width of the subintervals associated to distinct eigenvalues. In some cases, some of adjacent subintervals may overlap, e.g. in Figure 7, the support associated to λ₄ = 5 and λ₃ = 3 have overlap whereas the two smaller ones are separable. Figure 8 shows D h ⊂ C + , domain of h in the complex plane, using the same setting as in Figure 7. It has been shown that D_h approaches real axis only for the values of h for which z^′(h)>0 in Figure 7 which identifies the regions on the real axis outside of the support of eigenvalues. We observe that D_h has no intersection with the real axis between h = − 1 , − 1 2 , − 1 3 which reveals that the subintervals of support associated with three smallest eigenvalues λ=1,2,3, are not disjoint.

Figure 9, demonstrates the support of l.s.d. of SCM identified using Theorem 3 for c₀ ∈ {0.1,0.3} and λ₂ ∈ [1,4] with multiplicity of α₂ = 0.1 and λ₁ = 1 with multiplicity of α₁ = 0.9. We observe that for large values of λ₂, the support associated with two eigenvalues are disjoint intervals. However, these two disjoint intervals become connected as the distance between λ₂ and λ₁ reduces. In practice, the value of c₀ determines the window shape and has an important impact on the width of these intervals and on the location of the breakpoint. The location of breakpoint determines the capability of the window to identify two distinct eigenvalues. Figure 9 illustrates that the larger the value of c₀, the smaller the breakpoint of the support, i.e., by increasing p, we may be able to separate closer eigenvalues.

5.2 Limiting spectral distribution

In the noise only case, we find an explicit equation for the l.s.d. of the exponentially weighted SCM employing Lambert W function. However in the signal plus noise case, the l.s.d. can not be obtained explicitly and should be calculated numerically using (5) and (47). It is the same as the rectangular window case where the l.s.d of noise only data has M–P distribution, however there is no explicit equation for the signal plus noise case.

To find the imaginary part of the Stieltjes transform, one could alternatively find the complex roots with positive imaginary part of the inverse function z(m) for all z in the support of the eigenvalues, i.e., z ∈ S_F. Since the imaginary parts of m(z) and h(z) have the same sign and there is no explicit expression for z(m), we find the complex roots of z(h) using (47) and (48) for any real x_h = z(h)∈S_F, where Re {h}∈(h_b−,h_b+), b∈{1,…,s}. This can be done by finding ν = Im{h} for which Im {z(h)} = 0. By inserting the calculated h in (49), we obtain the Stieltjes transform for x_h∈S_F. Finally F^R(x) is obtained using (5). According to Proposition 1, for any z ∈ C + there exists a unique h satisfying (47) and (48), thus the above procedure results in the desired value of h and m.