title: Structured Autocorrelation Matrix Estimation for Coprime Arrays
authors: Dimitris G. Chachlakis Panos P. Markopoulos
text: 
  
  
  Structured Autocorrelation Matrix Estimation for Coprime Arrays
       This research is supported in part by the U.S. National Science Foundation, under grant OAC-1808582, and by the U.S. Air Force Office of Scientific Research, under the Young Investigator Program.
  
  This is a preprint of work that has been submitted for publication to a peer-reviewed journal.
  Dimitris G. Chachlakis and Panos P. Markopoulos^*^*Corresponding author.
  
  Department of Electrical and Microelectronic Engineering
   Rochester Institute of Technology
   Rochester, NY 14623
  E-mail:  and 
  
      December 30, 2023
  ==================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
  
  empty
  
  
  
  A coprime array receiver processes a collection of received-signal snapshots to estimate the autocorrelation matrix of a larger (virtual) uniform linear array, known as coarray. By the received-signal model, this matrix has to be (i) Positive-Definite, (ii) Hermitian, (iii) Toeplitz, and (iv) its noise-subspace eigenvalues have to be equal. Existing coarray autocorrelation matrix estimates satisfy a subset of the above conditions. In this work, we propose an optimization framework which offers a novel estimate satisfying all four conditions. Numerical studies illustrate that the proposed estimate outperforms standard counterparts, both in autocorrelation matrix estimation error and Direction-of-Arrival estimation. 	
  
  	
  Index Terms – Sensor array processing, Coprime arrays, Coarray, Autocorrelation estimation.
  
  PREPRINTPREPRINT
  	
  
  
  § INTRODUCTION
  
  
  In Direction-of-Arrival (DoA) estimation, coprime arrays offer increased Degrees-of-Freedom (DoF) and enable the identification of more sources than sensors compared to equal-length uniform linear arrays <cit.>. Coprime arrays have been successfully employed in applications such as beamforming design <cit.> and space-time processing <cit.>, to name a few.
  Other non-uniform arrays with increased DoF and closed-form expressions are the nested and MISC arrays <cit.>. 
  Intelligent processing of the autocorrelations of the physical array's elements enables the estimation of a signal subspace corresponding to a larger (virtual) array known as “coarray" which, in general, is non-uniform. Commonly, only a segment of the coarray is retained wherein the elements are uniformly spaced. Alternatively, some works employ interpolation methods to fill the “gaps” of the full coarray. In this work, we consider only the uniform segment of the coarray for simplicity. Our result is straightforwardly extended to the full coarray after its gaps are filled by existing interpolation methods (e.g., <cit.>).
  
  A coprime array receiver processes the autocorrelations of the physical-array's elements and estimates the autocorrelation matrix of the coarray. By the received-signal model, the nominal autocorrelation matrix of the coarray has a specific structure: it is (i) Positive Definite (PD), (ii) Hermitian, (iii) Toeplitz, and (iv) its noise-subspace eigenvalues are equal. In practice, the autocorrelations of the physical-array's elements are estimated by processing a collection of received-signal snapshots and diverge from the nominal ones. Accordingly, existing approaches offer autocorrelation-matrix estimates which diverge from the nominal one, while at the same time, violate at least one of the above structure-properties. 
  
  In this work, we propose an optimization framework for computing an improved coarray autocorrelation-matrix estimate that satisfies properties (i)-(iv). In practice, we iteratively solve a sequence of optimization problems, obtaining upon convergence, an improved estimate that satisfies properties (i)-(iv).
  Our studies illustrate that the proposed method outperforms standard counterparts, both in estimation error and DoA estimation. 
  	
  
  
  § SIGNAL MODEL AND PROBLEM STATEMENT
  
  
  We consider coprime integers M<N and design coprime array with L2M+N-1 elements at locations ℒ{(i-1)Md}_i=1^N ⋃{iNd}_i=1^2M-1, where d=λ 2 is the reference inter-element spacing and λ is the wavelength <cit.>. Narrow-band signals impinge on the array from K<MN+M sources with DoAs Θ{θ_1,θ_2,…,θ_K} under carrier frequency f_c and propagation speed c. Under far-field conditions, the kth source-signal impinges on the array from direction θ_k ∈ (-π/2, π/2] with respect to the broadside. Defining the element-location vector psort(ℒ), the array-response vector of source k becomes 𝐬(θ_k)[  v(θ_k)^[ p]_1, v(θ_k)^[ p]_2, … , v(θ_k)^[ p]_L]^⊤, where v(θ_k)exp(-  2π f_c/csin(θ_k)). Accordingly, the receiver collects received-signal snapshots of the form
  	
      𝐲_q ∑_k=1^K s(θ_k) x_q,k+𝐧_q ∈ℂ^L,
  
  	where x_q,k∼𝒞𝒩(0, d_k) and 𝐧_q ∼𝒞𝒩(0_L, σ^2 𝐈_L), model the qth symbol transmitted by the kth signal-source (power-scaled and flat-fading-channel processed) and Additive White Gaussian Noise, respectively. Received-symbols are uncorrelated across snapshots and sources. Noise-variables are uncorrelated from received symbols. The receiver's objective is to identify Θ from the collected snapshots. Next, we briefly review standard coprime array processing.
    
  Physical-Array Autocorrelation Matrix: The nominal received-signal autocorrelation matrix of the physical array is given by 𝐑_y 𝔼{𝐲_q 𝐲_q}=𝐒diag(𝐝)𝐒^H+σ^2𝐈_L, where S[ s(θ_1),…, s(θ_K)] and d[d_1,…,d_K]^⊤ denote the array-response matrix and source-power vector, respectively. Since Θ, d, and σ^2 are in practice unknown to the receiver, 𝐑_y can not be directly computed and is estimated based on Q received-signal snapshots by 𝐑_y1/Q∑_q=1^Q y_q y_q^H.
  	
  Autocorrelation Sampling:
  Nominally, the receiver processes R_y and computes the autocorrelation-vector 𝐫vec (𝐑_y)=∑_k=1^K a(θ_k) d_k+σ ^2  𝐢_L, where 𝐚(θ_k)𝐬(θ_k)^* ⊗𝐬(θ_k) and 𝐢_L vec(𝐈_L). 
  For j ∈ [L^2]{1,…, L^2}, it holds [𝐚(θ_k)]_j=v(θ_k)^n, n ∈𝒜sort({n_1-n_2| n_1,n_2 ∈ℒ}). The element-locations of the uniform segment of the coarray are described by
  
      ℬ{ n∈𝒜| 1-L'≤ n ≤ L'-1},
   where L'  MN+M. For every n ∈ℬ, the receiver discards duplicates by selecting any single index j_n ∈ [L^2] such that [𝐚(θ_k)]_j_n=v(θ_k)^n. That is, the receiver forms selection-sampling matrix
  
      𝐄_sel [	𝐞_j_1-L', L^2 , … , 𝐞_j_L'-1, L^2],
   where for any p ≤ P ∈ℕ_+, 𝐞_p,P is the pth column of 𝐈_P, and computes 𝐫_sel𝐄_sel^⊤𝐫 = ∑_k=1^K 𝐚_sel(θ_k) d_k + σ^2 𝐞_L', 2L'-1, where 𝐚_sel(θ_k) 𝐄_sel^⊤𝐚(θ_k) = [v(θ_k)^1-L', …,  v(θ_k)^L'-1] ^⊤.
  In practice, the autocorrelation-vector 𝐫 is estimated by 𝐫vec(𝐑_y) and  r_sel is estimated by r_sel𝐄_sel^⊤𝐫.
  	
  	
  Coarray Autocorrelation Matrix:
  	The receiver applies a rank-enhancement approach on r_sel (or, 𝐫_sel in practice) to form the autocorrelation matrix of the coarray. Commonly, the Augmented Matrix <cit.> and Spatial Smoothing <cit.> approaches are employed.
  	According to the augmented matrix approach, the receiver computes
  	
      R_am F ( I_L'⊗ r_sel) ∈ℂ^L' × L',
  
  	where
  	𝐅 [𝐅_1, 𝐅_2, …, 𝐅_L'] and, for every m ∈ [L'],  𝐅_m  [0_L' × (L'-m), 𝐈_L', 0_L' × (m-1)].
  	R_am has full-rank, is PD, Hermitian, Toeplitz, and coincides with the autocorrelation matrix of the coarray
  	
      R_co= 𝐒_codiag(𝐝) 𝐒_co^H + σ^2 𝐈_L',
  
  	where 
  	[𝐒_co]_m,k v(θ_k)^m-1, for every m ∈ [L'] and k ∈ [K]. According to the spatial-smoothing approach <cit.>, in the case of known statistics, the receiver computes the spatially-smoothed matrix
  	R_ss1/L' R_am R_am^H
  which is not an autocorrelation matrix but an autocorrelation matrix is extracted from it as a scaled version of its principal square root 
      R_psr√(L' )𝐑_ss^1/2.
  
  	We notice that R_am R_am^H= R_am^2=L' R_ss. Moreover, R_ss admits Singular Value Decomposition (SVD) R_sssvd= U Λ𝐕^H which implies that R_psr= U(√(L')Λ^1/2)𝐕^H= R_am. That is, R_psr and R_am both coincide with . Here, we note that in the (ideal) case of known statistics to the receiver, all estimates above coincide with the nominal autocorrelation matrix of the coarray and satisfy (i)-(iv). However, in the practical case of unknown statistics (case of interest) to the receiver, the estimates above diverge from  and satisfy only a subset of (i)-(iv):
  The augmented matrix approach of <cit.> proposed to substitute the sampling matrix 𝐄_sel by the averaging sampling matrix 	
  
      𝐄_avg [	1/|𝒥_1-L'|∑_j ∈𝒥_1-L'𝐞_j,L^2, … , 1/|𝒥_L'-1|∑_j ∈𝒥_L'-1𝐞_j,L^2],
  
  where for every n ∈ℬ, 𝒥_n={j ∈ [L^2]| [ a(θ_k)]_j=v(θ_k)^n}, substituting r_sel by r_avg𝐄_avg^⊤𝐫.[If the nominal statistics are known, r_sel and r_avg=𝐄_avg^⊤ r coincide. The latter does not hold if r_sel and r_avg are estimated by r_sel and r_avg, respectively.] Accordingly, R_co is estimated by R_am𝐅 (𝐈_L'⊗ r_avg).
  	Importantly, it holds that R_am is Hermitian and Toeplitz, however, it's not guaranteed to be PD. That is, R_am can be an indefinite estimate of  <cit.>.
  	Similarly, R_ss and R_psr are estimated by
  	R_ss1/L'𝐑_am𝐑_am^H and 𝐑_psr√(L') R_ss^1/2, respectively, where R_am= F( I_L'⊗ r_sel).[R_am and R_am denote the augmented matrix approach estimates combined with averaging and selection sampling, respectively.] R_am can be an indefinite matrix. In view of the above,  R_psr is by construction a PD and Hermitian matrix estimate of the coarray autocorrelation matrix, however, it violates the Toeplitz structure-property of .
  	It follows that R_psr and R_am no longer coincide, however, their left-hand singular-valued singular vectors span the same signal subspace. For the unknown statistics case, we summarize the above estimates in Table <ref>, where for each estimate we mention the employed autocorrelation sampling approach. Moreover, for each structure property guaranteed to be satisfied, we place a , otherwise, we place a .  
  	Given a coarray autocorrelation matrix estimate R∈{ R_am,  R_am,  R_psr}, a standard DoA estimation approach–e.g., MUltiple SIgnal Classification (MUSIC)–is applied for identifying the DoAs in Θ.
  
  
  	
  
  
  § PROPOSED AUTOCORRELATION-MATRIX ESTIMATE
  
  
  We propose an algorithm which iteratively solves a sequence of optimization problems returning, upon convegrence, an improved coarray autocorrelation matrix estimate. Motivated by <cit.>, where it was formally proven that averaging autocorrelation sampling attains superior autocorrelation estimates compared to selection sampling with respect to the MSE metric, we propose to initialize the proposed algorithm to P_0=√(1/L' R_am R_am^H).
  	At iteration  i≥ 0, the proposed algorithm computes 
  	
      Q_i=P_i,
          R_i=Q_i,
          P_i+1= R_i,
   
  where for any X= X^H∈ℂ^D × D with Eigen-Value Decomposition (EVD) Xevd= U diag(X)  U^H the following hold.
  	
  		X returns the nearest Toeplitz matrix, in the Euclidean norm sense[Otherwise known as Frobenius norm: ·_F^2 returns the sum of the squared entries of its argument.], to X: X T_0 ∈𝒯^D X - T_0_F^2, where
  		𝒯^D{𝐓∈ℂ^D × D| T is Toeplitz}.
  		
  	 
  	
  		X returns the nearest Positive Semidefinite (PSD) matrix to X: X X_0 ∈𝕊^D_+ X - X_0_F^2, where
  		𝕊_+^D{𝐀∈ℂ^D × D| A= A^H≽0}. 
  		
  	 
  	
  		X performs an eigenvalue-correction operation to the D-d smallest eigenvalues of X. For some general d ∈{1,…,D-1},
  		X U diag(X) U^H,
  		where 
      [X]_i=
      		[X]_i,     i≤ D-d+1, 
      1/D-d∑_i=d+1^D[X]_i,     i> D-d+1.
  
  		
  	  
  	In view of the above, the proposed algorithm seeks to optimize the D-d smallest eigenvalues of the autocorrelation matrix estimate at which it is initialized while preserving the PSD, Hermitian, and Toeplitz structure. Next, we conduct formal convergence analysis of the proposed algorithm. 
  	Consider arbitrary 𝐗= X^H  ∈ℂ^D × D and let 𝐝_i(𝐗) denote a diagonal of 𝐗 (see Fig. <ref>) such that 
  	
      [𝐝_i(𝐗)]_j=
      	[𝐗]_j-i,j,    i ≤ 0, 
      
      	[𝐗]_j,j+i,    i > 0,
  
  	for any j ∈{1,2,…,D-|i|}.
  	The following remarks hold.
  		
  	
  		It holds that 
  		𝐝_i(𝐗)=𝐝_i'^*(𝐗),  if  |i|=|i'|.
  	
  	
  		Let 𝐓X=min_𝐓_0 ∈𝒯^D𝐗-𝐓_0_F^2. For any i ∈{1-D,…, D-1}, it holds that 𝐝_i(𝐓)=1/D-|i|1_D-|i|^⊤𝐝_i(𝐗)1_D-|i|.
  		 
  	
  	
  		Let 𝐓X. It holds that 
  		𝐓_F^2 ≤𝐗_F^2.
  	
  	
  	
      𝐓_F^2   =∑_i=1-D^D-1𝐝_i(𝐓)_2^2
      	=∑_i=1-D^D-11/D-|i|(1_D-|i|^⊤𝐝_i(𝐗))1_D-|i|_2^2
      	=∑_i=1-D^D-1(1_D-|i|^⊤𝐝_i(𝐗))^2/(D-|i|)^21_D-|i|_2^2
      	
         =∑_i=1-D^D-1(1_D-|i|^⊤𝐝_i(𝐗))^2/D-|i|≤∑_i=1-D^D-11_D-|i|_2^2 𝐝_i(𝐗)_2^2/D-|i|=𝐗_F^2.
  
  
  	
  		Let 𝐓X admit EVD[A Hermitian matrix 𝐀 can be expressed as 𝐀=𝐔Λ𝐔^H, where Λ is an upper diagonal with the eigenvalues of 𝐀 in its main diagonal. If 𝐀 is normal (i.e., 𝐀𝐀^H=𝐀^H 𝐀), then Λ is diagonal. Every Hermitian matrix is normal.] 𝐓=𝐔diag(λ_T)𝐔^H. It holds 
  		𝐓_F^2=Tr(𝐔diag(λ_T) diag(λ_T)𝐔^H)
  		=λ_T_2^2.
  		
  	
  	
  		Let 𝐓X admit EVD 𝐓=𝐔diag(λ_T)𝐔^H and define 𝐏 such that 𝐏=𝐔diag(λ_P)𝐔^H,
  		where λ_P=λ_T^+,–i.e., ∀ i ∈ [D], [λ_P]_i=max{[λ_T]_i,0}. It holds that 𝐏 is the solution to 𝐏_0 ∈𝕊_+^Dmin.𝐓-𝐏_0_F^2.
  		A proof for Remark <ref> was first offered for real matrices in <cit.>.
  		For completeness purposes, we offer an analogous proof for complex-valued matrices.
  		
  		
  			
  	Consider Hermitian T ∈ℂ^D × D with EVD 𝐓=𝐔diag(λ_T)𝐔^H, U is unitary (i.e., U U^H= U^H U= I_D). Let H= U^H  P_0 U which implies that P_0= U  H  U^H. It holds 
      min_𝐏_0 ∈𝕊_+^D𝐓-𝐏_0_F^2   =min_ H ∈ℂ^D × Ddiag(λ_T)- H_F^2
      	=min_ H ∈ℂ^D × D∑_i,j
      i≠ j [ H]_i,j^2+∑_i=1^D([λ_T]_i-[ H]_i,i)^2
      	
         ≥∑_i,j
      i≠ j [ H]_i,j^2+∑_i=1^D([λ_T]_i-[ H]_i,i)
      	≥∑_i ∈{1,2,…,D | [λ_T]_i<0}[λ_T]_i^2.
  
  	Similar to <cit.>, the lower bound in (<ref>) is attained by matrix H=diag(λ_P) for λ_P such that [λ_P]_i=max{[λ_T]_i,0}.
  	
  			
  	
  		For P=T, it holds that 
  		𝐏_F^2≤𝐓_F^2. Formally,
  		𝐏_F^2=λ_P_2^2=∑_i=1^D[λ_P]_i^2≤∑_i=1^D[λ_T]_i^2=𝐓_F^2.
  	
  	
  		Let 𝐀=P=𝐔diag(λ_A)𝐔^H,
  		where 
      [λ_A]_i=
      		[λ_P]_i,     i≤ D-d+1, 
      1/D-d∑_j=d+1^D[λ_P]_i,     i>D-d+1,
  
  		for d ∈{1,…,D-1}. It holds that
  		𝐀_F^2≤𝐏_F^2.
  		
  		
  		
  				
      𝐀_F^2   =∑_i=1^D[λ_A]_i^2=∑_i=1^D-d[λ_P]_i^2+d(1/d∑_j=D-d+1^D[λ_P]_j)^2
      				≤∑_i=1^D-d[λ_P]_i^2+∑_j=D-d+1^D[λ_P]_j^2=𝐏_F^2.
  
  		
  		In view of Remarks <ref>-<ref>, the following Proposition derives. 
  
  For Q_i,  R_i, and P_i+1 in (<ref>)-(<ref>), it holds Q_i_F^2≥ R_i_F^2≥ P_i+1_F^2≥ Q_i+1_F^2≥…≥ 0 ∀ i≥ 0.
  		
  	 
  		
  Proposition <ref> states that the proposed algorithm is guaranteed to converge. In practice, one can terminate the iterations when P_i+1- P_i≤ϵ, for some ϵ≥0. For sufficiently small ϵ, Proposition 1 implies that, at convergence, P_i+1= R_i= Q_i which, in turn, implies that the algorithm converged to a PD, Hermitian, and Toeplitz matrix the noise-subspace eigenvalues of which  are equal. A visual illustration of Proposition <ref> and a pseudocode of the proposed algorithm are offered in Fig. <ref> and  Fig. <ref>, respectively. Importantly, ∀ i≥ 0, the Algorithm of Fig. <ref> computes Q_i,  R_i,  P_i+1 by closed-form expressions with cost at most the cost of EVD–i.e., 𝒪(D^3). Overall, the cost of the proposed algorithm is 𝒪(TD^3) where T  is the number of iterations required for convergence.
  	
  	
  
  § NUMERICAL STUDIES
  	
  	
  We consider coprime naturals (M,N)=(3,5) and form coprime array with L=10 elements. Source-signals impinge on the array from K=13 DoAs {θ_k}_k=1^13,  θ_k=(-75+(k-1)12)^∘. The noise variance is set to σ^2=0dB. All sources emit signals with equal power d_k=α^2dB. Accordingly, the Signal-to-Noise Ratio SNR= α^2. The receiver collects Q ∈{150,300,450,600} received-signal snapshots. For every Q, we consider 3000 statistically independent realizations of noise; i.e., { y_r,1,…, y_r,Q}_r=1^3000. At every realization r, we compute coarray autocorrelation matrix estimates corresponding to the augmented matrix approach (AM), principal square root of the spatial smoothed matrix (PSR), nearest Hermitian, PSD, and Toeplitz (H-PSD-T) approach of <cit.>[H-PSD-T seeks a Hermitian-PSD-Toeplitz matrix which fills the gaps of the coarray. When the uniform segment of the coarray is considered, H-PSD-T returns R ∈𝕊_+^L' R- R_am_F^2+μ R_*, where μ R_* is a regularization term that moderates overfitting.], and the proposed structured estimate.  For every method and estimate 𝐑_r at realization r, we compute the normalized squared error NSE(𝐑_r)=𝐑_r-_F^2 _F^-2. Then, we compute the Root Mean Normalized Squared Error  RMNSE=√(1/3000∑_r=1^3000NSE(𝐑_r)). In Fig. <ref> and Fig. <ref>, we plot the  RMNSE versus sample support, Q, for SNR=-4dB and SNR=2dB, respectively. Expectedly, we observe that all methods employing averaging-sampling perform similarly well. The proposed estimate attains superior estimation performance across the board. Moreover, we notice the sensitivity of H-PSD-T with respect to the ad-hoc parameter μ; e.g., for SNR=-4dB, H-PSD-T with μ=1.5 exhibits low performance while for  SNR=2dB it exhibits high estimation performance.
  	
  Thereafter, we consider that the nominal coarray autocorrelation matrix admits SVD = Q_coΣ_co V_co^H+σ^2Q̅_coQ̅_co^H, where Q_co and Q̅_co correspond to the signal and noise subspace bases, respectively. Similarly, every coarray autocorrelation matrix estimate 𝐑_r admits SVD 𝐑_r= Q_rΣ_r  V^H_r +Q̅_rΣ̅_rV̅_r^H, where Q_r denotes the signal-subspace-basis of the K dominant left-hand singular valued singular vectors of 𝐑_r. At each realization and for every value of Q, we compute the normalized squared subspace error  NSSE(𝐐_r)=𝐐_r𝐐_r^H-𝐐_co𝐐_co^H_F^2 (2K)^-1. Then, we compute the Root Mean Normalized Squared Subspace Error  RMN-SSE=√(1/3000∑_r=1^3000NSSE(𝐐_r)). In Fig. <ref> and Fig. <ref>, we plot the RMN-SSE versus sample support for SNR=-4dB and 2dB, respectively. We notice the influence of the ad-hoc parameter μ with respect to H-PSD-T and observe that the proposed structured estimate clearly outperforms all counterparts across the board in subspace estimation performance. 
  		
  Next, for every value of sample support and realization r, we conduct DoA estimation by applying MUSIC on the estimate 𝐑_r which returns estimates {θ_k,r}_k=1^13 <cit.>. Then, we measure the Root Mean-Squared-Error RMSE=√(1/30001/13∑_r=1^3000∑_k=1^13 (θ_k-θ_k,r)^2) and illustrate the corresponding RMSE curves versus sample support Q, in Fig. <ref> and Fig. <ref>, for SNR= -4dB and 2dB, respectively. We include the Cramér Rao Lower Bound (CRLB) curves as benchmarks <cit.>. We notice that the the performances of standard counterparts (AM, PSR) deviate significantly from the CRLB. In contrast, the proposed coarray autocorrelation matrix estimate outperforms all counterparts by at least 0.3^∘ and at most 2^∘. In addition, as Q increases, its performance curves approach the CRLB curves.
  
  
  § CONCLUSIONS
  
  
  We proposed an optimization framework which computes a structured coarray autocorrelation matrix estimate. The proposed algorithm is accompanied by convergence analysis and is guaranteed to return a a coarray autocorrelation matrix estimate satisfying all structure properties of the true autocorrelation matrix. Numerical studies illustrate the enhanced performance of the proposed estimate compared to standard counterparts, both in autocorrelation matrix estimation error and DoA estimation. 
  
  IEEEtran
  
  