Source separation and localisation via tensor decomposition for distributed arrays

: This study focuses on the problem of power spectra separation and localisation of multiple sources using distributed arrays. First, the array structure and signal model are discussed. By cross-correlating the multi-channel received signals in time-domain, a third tensor is constructed. Then, utilising the multi-dimensional characteristic, the tensor is decomposed to separate the array manifold matrix and the power spectra matrix through alternating least square (ALS) method. Finally, the sources are located using the relative x - y plane relationship between the distributed arrays and the direction of arrival (DOA), which can be estimated by spectrum analysis of each column of the array manifold matrix. The effectiveness and superiority of the proposed method is demonstrated by simulation results.


Introduction
Source separation and localisation are widely used to solve the problem of spectrum management in cognitive radio and localise the active emitter in electronic warfare and radar electronic countermeasure [1][2][3][4]. Spectrum sensing is the first step towards situation awareness. The radio frequency environment can be mapped out to highly efficient reuse in space and time by identifying the power spectra and locations of the active emitters, including enemy radars or jammers.
There is a lot literature on spectrum separation. A parallel frequency bin detection method is proposed in [1]. Exploiting the frequency-domain sparsity, compressive sensing is used in [2,3]. Whereas most works address in reconstructing the Fourier spectrum of the received signals, but in cognitive radio and passive localisation applications only the power spectrum (PS) is necessary. Thus, there is no reason to reconstruct the time-domain signals. Utilising that the PS is the Fourier translation of the autocorrelation function, a finite number of autocorrelation lags can be estimated using received array signals. In [4], the separation of multiple received mixed spectra is treated as non-negative matrix factorisation (NMF) problem. However, the NMF is not unique in general; hence, the identification of the mixed spectra cannot be guaranteed. Moreover, the most existence works consider single received station, which only gives the direction of arrival (DOA) of the each source, not the x-y plane location.
In this paper, we configure a distributed network with multiple stations to separate and localise multiple sources by estimating the PS and DOAs (respect to each station). The array structure and signal model are given in Section 2. In Section 3, first, a third tensor is constructed by cross-correlating the multi-channel received signals in time-domain. Then, utilising the multidimensional characteristic, the tensor is decomposed to separate the array manifold matrix and power spectra matrix through alternating least square (ALS) method. Finally, the sources are localised using the relative x-y plane relationship of the distributed arrays and DOAs. Simulations are reported in Section 4. We conclude in Section 5.
The contribution of this paper lies in formulating the problem of joint power spectra separation and localisation of multiple sources as a problem of tensor decomposition, which benefits from crosscorrelating the multi-channel received signals in time-domain. The identification of both spectrum and DOAs of each source can be guaranteed, as the tensor model is unique under fairly mild conditions.

Array configuration and signal model
Consider a distributed network configuration given in Fig. 1. There are P widely separated stations with index p = 1, 2, …, P, and each station has N closed spaced sub-arrays with index n = 1, 2, …, N, and each sub-array is an uniform linear array (ULA) with M received antennas.
Assuming that there are K uncorrelated narrowband sources in far-field, the received signal of the nth sub-array in the pth station is where x (p) (n, t) ∈ C M × 1 , β n, k (p) is the path loss and phase shift from the kth source to the nth sub-array in the pth station, ] T is the received steering vector of the kth source to the nth sub-array in the pth station, θ n, k (p) is the DOA of the kth source to the nth sub-array in the pth station, d denotes the distance between neighbouring received antennas, λ is the wavelength, s k (t) is narrowband transmitted signal of the kth source. z (p) (n, t) is the noise term. The N subarrays in each station are closed spaced, which means θ n, k can be replaced by Then (1) can be written as Our aim is to estimate power spectra and the DOAs of all K sources, denoted by } with k = 1, 2, …, K and p = 1, 2, …, P.

Temporal correlation
To estimate power spectra and the DOAs of all K sources, we propose to formulate the problem in the temporal correlation domain. The K sources are uncorrelated, so the cross-correlation of the m 1 th and m 2 th antennas of the nth sub-array in the pth station is . δ(x) denotes the Kronecker delta function, i.e. δ(x) = 1 when x = 0 and δ(x) = 0 otherwise. Assuming that the noise at each antenna is white Gaussian, both temporally and spatially, with zeros mean and variance σ 2 , i.e. z (p) (n, t) ∼ N(0, σ 2 ) for all p and n. Then, (3) can be expressed as Applying discrete time Fourier transform on (4), we have The frequency axis in (5) can be discretised as ; ⋯; G n (P) ] ∈ C PQ × F and g n = vec(G n ), stacking g n one after another, we have where (1) ; A (2) ; ⋯; A (P) ] ∈ C PQ × K , B = [β 1 , β 2 , …, β K ] ∈ R NP × K , and β k = vec(B k ), B k is the received path loss matrix with the β n, k (p) 2 on the nth row and the pth column.
Z ∈ C FQP × NP is the noise term.

Parameter matrix estimation via tensor decomposition
Definition 1: The n-mode matrix unfolding of an N-order tensor X ∈ C I 1 × I 2 × ⋯ × I N , denoted by X n ∈ C I n × (I n + 1 ⋯I N I 1 ⋯I n − 1 ) , contains the (i 1 , i 2 , …, i N )th element of X at the position (i n , j), where On the contrary, given X n , X is constructed. Each third-order tensor can be recognised as a data cube, and its matrix unfoldings consist of the slices of the cube along different directions. Using definition 1, a third-order tensor H can be formed by H in (7) [5]. The signal tensor H and its matrix unfoldings are given in Fig. 2 The minimise problem in (9) can be solved via ALS method [7], which successively estimates one of the three matrix unfoldings assuming the other two are known in least squares (LS) principle. Given Â , B, update Given Ŝ, B, update Given Â , Ŝ, update In (10)-(12), the pseudo inverse of the Khatri-Rao product is calculated by [5] It can be seen that the estimated matrix unfolding Ŝ contains the power spectra of sources, and Â contains the DOAs.
Remark 1: Using line search to accelerate the convergence of ALS. Mostly, ALS needs a large number of iterations before converging. The slowness in convergence can be due to the large size of the tensor, or to the bad starting values. Line search is an effective solution proposed to cope with the problem of slow convergence. Some line search methods can be used to speed up ALS for searching the global minimum very quickly. Refer to [8,9] for detail.

DOA estimation and source localisation
Let â k = [â 1, k ; â 2, k ; ⋯; â P, k ] be the kth column of Â , where â p, k = [α −(M − 1) , …, 0, …, α (M − 1) ] T ∈ C Q × 1 and α = e j2πdsin(θ p, k )/λ . The DOA of the kth source respected to the pth station can be estimated via spectrum analysis of â p, k . In this paper, root-MUSIC [10], one of the classical spectrum analysis methods, is used. Calculate the noise subspace U of the covariance matrix R = â p, k â p, k H via eigen decomposition. Then the qth coefficient of polynomial equals the sum of the (q − Q)th diagonal of G = UU H with q = 1, 2, …, (2Q − 1). Then we have θ k (p) = asin[ρλ/(2dπ)], where ρ is the root of the polynomial that is nearest and inside the unit circle.
The (x, y) location of the kth the source can be estimated using the relative triangular relationship of the P stations and the DOAs θ k (p) , where p = 1, 2, …, P.

Steps of the proposed algorithm
The proposed algorithm mainly includes three steps: (i) according to the temporal correlation domain, get the signal matrix H; (ii) Utilising the multi-dimension structure, construct the tensor H, and decompose the tensor by ALS; and (iii) Estimate DOAs from manifold matrix and localise the sources. The outline of the proposed algorithm is summarised in Algorithm 1 (Fig. 3).

Simulations
In this section, we use computer simulations to demonstrate the effectiveness of the proposed algorithm.  Fig. 4 shows the estimated power spectra of the sources by the proposed method. The results are obtained by averaging 50 trials, to enable easy visual assessment of estimate variance, including leakage from one source to the others. We see that the proposed method identifies the power spectra of all three sources fairly well.
Case 2: Fig. 5 shows the accuracy of the estimation power spectra and localisation versus SNR using the proposed method under different number of stations P = 2, 3, 5. Under P = 5, the (x, y) locations of the stations are (−20, 0)Km, (−10, 0)Km, (0, 0)Km, (10, 0)Km, and (20, 0)Km. Under P = 3 and P = 2, the first 3 and 2 stations are used, respectively. The other parameters are the same as Case 1. The number of Monte-Carlo trials is 1000. The rootmean-square error (RMSE) of the estimated power spectra denoted by RMSE PS and the RMSE of the estimated (x, y) location denoted by RMSE x-y are adopted as the performance measure, defined as It can be seen in Fig. 5, the proposed method can jointly separate PS and localise the (x, y) position of multiple sources. RMSE PS and RMSE x-y both reduce as the SNR increases. The proposed method yields reasonable RMSE PS of PS even under lower SNR, e.g. <

Conclusion
The problem of joint power spectra separation and localisation of multiple sources for distributed arrays has been considered in this paper. Utilising the temporal correlation domain, this problem is formulated as a third-order tensor decomposition problem. The tensor and its matrix unflodings are discussed. ALS method is used to solve the tensor, which results the manifold of each array and the power spectra matrix. The DOAs of each source respected to all arrays can be estimated by imposing root-MUSIC method on the estimated array manifold. Simulations illustrated the accuracy and efficacy of the proposed techniques.