Graph Signal Processing for Narrowband Direction of Arrival Estimation

For direction of arrival (DOA) estimation based on graph signal processing (GSP), it has been assumed that there is a phase shift between adjacent snapshots of the received signals. However, this assumption does not hold for narrowband signals and thus affects the performance of the corresponding algorithms. To improve the performance, a new GSP-based DOA estimation method is proposed. By building a periodic directed graph based on a graph shift operator and computing the spectrum using the Kronecker product, the relationship between the input narrowband signals and the graph adjacency matrix of different direction coefficients is constructed. Simulation results show that this method performs better than existing algorithms based on GSP.


I. INTRODUCTION
Since the introduction of graph signal processing [1], a large number of studies on various applications of graph signal processing (GSP) have emerged, especially for massive data with irregular sensor structures [2].For example, the filtering operation in GSP can be applied to smoothing [3], denoising [4] and classification [5] of irregular signals and the clustering of graphs can be used to identify weak sources in dense sensor arrays [6].The result from the classical signal time-frequency uncertainty principle was extended to graph signal processing and the use of eigenvectors of the graph Laplacian matrix is justified as the base for the graph Fourier transform in [7].
Direction of arrival (DOA) estimation is a classic array signal processing problem and has been solved by various methods, such as MUSIC, ESPRIT and sparsity based ones [8]- [10].Recently, based on the graph Fourier transform, a directed graph was built to represent characteristics in the spatial domain, and the relationship between eigenvectors of the adjacency matrix and the input angle for DOA estimation is derived in [11].On top of this, by adding an adjacency matrix representing the time domain, the estimation accuracy is improved by combining signals from different snapshots [12].As the Fourier transform can also be used to perform distance estimation, it is possible to extend the GSP-based DOA algorithm to two-dimensional estimation by combining DOA and distance information [13].In [14], two graph adjacency matrices of the signal were formed and spectral peak searching was then performed according to the orthogonality property of its subspaces so that the incoming wave direction of the target can be obtained with a reduced computational cost.
Currently, existing works about DOA estimation based on GSP assume that there is a single-frequency based phase shift between signals received by the same sensor at different snapshots and derive the targeted angle from the graph adjacency matrix based on this assumption, such as algorithms from [12] and [14].However, this assumption does not usually hold for narrowband signals as we normally work at the baseband and the delay between adjacent baseband snapshots of the signal is much higher than the considered delay based on the carrier frequency, which leads to model mismatch issues and affects the estimation performance of the derived algorithms.In [15], the graph shift operator is proposed to represent a periodic time series signal.Similar to this representation, we use the Kronecker product to combine the time adjacency matrix based on this graph shift operator and the space adjacency matrix based on phase shift between sensors in the same array.With eigenvalue decomposition (EVD) of the above combined matrix, we then derive a relationship between its eigenvectors and received narrowband signals.Simulation results show that this method improves the accuracy of DOA estimation in a range of SNR values and different numbers of snapshots, compared to existing solutions based on GSP.
The rest of the paper is structured as follows.The basics of DOA estimation and GSP and their relationship are presented in Sec.II.In Sec.III, the new GSP based DOA estimation method is introduced.Simulation results are shown in Sec.VI, and conclusions drawn in Sec.V.

II. PRELIMINARIES
The uniform linear array (ULA) model with M sensors is shown in Fig. 1.Consider a narrowband complex plane-wave signal s(t)e jωt arriving from angle θ, where θ ∈ [− π the broadside, s(t) is a baseband signal, and ω is the carrier frequency, which is much larger than the bandwidth of s(t).

A. Array Model
Assume that the signal received by the first sensor is x 0 (t) = s(t)e jωt ; then, at the mth sensor the signal is x m (t) = s(t − τ m )e jω(t−τm) ≈ s(t)e jω(t−τm) (1) for m = 1, ..., M − 1, where τ m is the propagation delay for the signal from sensor 0 to sensor m and is a function of θ.
The approximation s(t − τ m ) ≈ s(t) is valid since s(t) has a narrow bandwidth and changes little over the time delay τ m .
Thus, there is a phase shift between signals from different sensors at the same snapshot (e.g., x 2 (t) = x 1 (t)e −j(2πd sin θ)/λ ).However, signals at different snapshots are independent of this phase shift, since the delay between adjacent baseband snapshots of the signal is much larger than the considered delay based on the carrier frequency.This discrepancy leads to a model mismatch problem in existing GSP-based DOA estimation algorithms and affects their estimation performance.We will propose a new solution based on the narrowband model in Sec.III.

B. Graph Signal Processing
An unweighted graph G = (V, B) is defined as a set of nodes V connected by edges B. If each node of a directed graph is connected to only one preceding node and one succeeding node, then such a graph is said to be a directed circular graph [16], as shown in Fig. 2(a).For the element A ij of the graph adjacency matrix, assume that A ij = 0 if nodes i and j are not connected by edges, and A ij ∈ (0, 1] if nodes i and j are connected by edges.Then the adjacency matrix of the graph is shown in Fig. 2  The Kronecker product of two graphs The adjacency matrix A ⊗ of the new graph is given by

C. Graph Model for Array Signal Processing
Signals received at the same snapshot from different sensors are related by e −jσ with σ = 2πd sin θ/λ (Fig. 3).
The zeros on the diagonal ensure that its rank is M , which is the number of array sensors.Thus, the adjacency matrix has different eigenvalues and each eigenvector is orthogonal to the other eigenvectors.The coefficient 1 2 is used to form the following relationship for the case of one source where the effect of noise has been ignored.Therefore, the input signal x is the eigenvector corresponding to the adjacency matrix A s when the eigenvalue is 1.

A. Time Adjacency Matrix
For signals received from the same sensor at different snapshots, the periodic time series relationship can be represented by a unit directed cyclic graph [15] as shown in Fig. 4. Its corresponding time adjacency matrix is given by where A t is an N × N matrix, and N is the number of snapshots.Moreover, A t can be seen as a graph shift operator, and since this is a time series signal of period N , it will return to the original position after shifting N times.Thus, it will be a unit matrix after shifting N − 1 times,

B. Time-space Adjacency Matrix
The Kronecker product of two disjointed graphs is defined in Eq. (7).Therefore, the adjacency matrix of the time-space graph, which is the Kronecker product of the time and space graphs, is given by Based on Eqs. ( 8), ( 10) and ( 12), the time-space adjacency matrix is given by where A ⊗ is an M N × M N matrix.Based on Eqs. ( 11) and ( 13), the time-space adjacency matrix after shifting N − 1 times is given by Assume that the signals received by array sensors at all snapshots are given by According to Eq. ( 9), the product of the input signal and shifted time-space adjacency matrix is given by Accordingly, the input signal x is also the eigenvector of the shifted time-space adjacency matrix A N ⊗ corresponding to the unit eigenvalue.

C. Spectral Decomposition
Since both A s and A t are diagonalizable matrices [17], through eigendecomposition, they can be decomposed into the product of matrices consisting of their eigenvalues and eigenvectors, as shown in Eq. (17): where V is the matrix whose columns represent the eigenvectors of adjacency matrix, and Λ is the diagonal matrix consisting of corresponding eigenvalues.Based on Eqs. ( 12) and ( 17), the spectral decomposition of the time-space adjacency matrix A ⊗ is given by As both A s and A t are normal matrices, A N ⊗ is also a normal matrix.As a result, for matrix A N ⊗ , eigenvectors corresponding to different eigenvalues are orthogonal to each other.In addition, both A s and A t are of full rank, so A N ⊗ has full rank as well.Thus, the eigenvalues of A N ⊗ vary and each eigenvector is orthogonal to the others.Furthermore, since the input signal x is the eigenvector of A N ⊗ , corresponding to a unit eigenvalue, x will be orthogonal to other eigenvectors with non-unit eigenvalues.

D. GSP Applied to DOA Estimation
Similar to the definition of the Discrete Fourier Transform (DFT), a graph Fourier transform (GFT) is defined as where each column of V N is a different eigenvector of A N ⊗ .Each element of x f is equal to the dot product of each eigenvector and x.According to Section III-C, x is orthogonal to the eigenvectors except for itself, which means that all elements in x f are 0 except for the element corresponding to the unit eigenvalue.
Since eigenvectors of the matrix A N ⊗ can be changed by changing the azimuth angle θ, x f of different angles can be calculated, and when x f satisfies the condition that only one element is non-zero and the others are all zeros, the corresponding angle is the desired result.
The absolute values |x f i | of the GFT coefficients obtained from a 5-sensor ULA with two snapshots is shown in Fig. 5.As shown in Fig. 5(a), there is only one peak with a large =30.3° (Correct angle) GFT coefficient, while there are several peaks with smaller GFT coefficients in Fig. 5(b).For eigenvector matrix V N at the correct angle, the energy will be concentrated on the element corresponding to unit eigenvalue.For V N at other angles, the energy will be spread over different elements.
Thus, a cost function can be constructed as where i uni represents the eigenvalue index corresponding to the unit eigenvalue.The cost function f (θ) achieves its maximum only if the eigenvector matrix V N corresponds to the correct angle.Note that although the above cost function is derived based on the one source assumption, it is applicable to multiple sources, as the GFT matrix is independent of source number and the received data vector x can be decomposed into the sum of sub-data vectors corresponding to each source (for the correct DOA angle corresponding to one source, the other subdata vectors will simply be treated as noise).

IV. SIMULATION RESULTS
In the following simulations, the number of sensors M = 5, the array element distance d = λ/2, the real angle θ r = −30.3• for a single source signal, and θ r = −30.3The GSP method in [12] uses phase shifts to build up the time adjacency matrix A t , while the improved GSP method proposed in this work uses a unit directed cyclic graph instead.As shown in Fig. 6, the performance of our proposed GSP method is better than that of the method in [12].With the white Gaussian noise, some disturbance can be seen in both methods, but the positions of dominant peaks produced by the proposed GSP method are a significantly better match to the true angles of arrival than those for the method in [12].
As shown in Fig. 7 (a) and (b), when the signal to noise ratio (SNR) is 10 dB, with the increasing number of snapshots N , the root-mean-square error (RMSE) of the estimated angle in degrees has been reduced for all three methods considered.The proposed GSP method performs better than the method in [12], but still falls short compared to the MUSIC algorithm.The gap between GSP and MUSIC decreases as N increases.For N = 40, with increasing SNR, the RMSE decreases for all three methods.The proposed method shows more accurate DOA estimation than the method in [12].MUSIC has an advantage over the two GSP-based methods for SNRs less than 5, but when the SNR is greater than 5, the results of the three methods are almost the same.
As shown in Fig. 7 (c) and (d), for two source signals, all three methods perform better as the number of snapshots increases, while the proposed method shows a significant improvement compared to the method in [12].With a larger N , the proposed method narrows the gap between itself and MUSIC.In Fig. 7 (d), a higher SNR leads to less disturbance to the array signals, resulting in higher accuracy for all three methods.With increasing SNR, the proposed method performs better than the method in [12], while it also has better noise immunity compared to MUSIC when the SNR equals 0 dB.This result demonstrates the potential of the proposed method for applications in noisy and complex environments.

V. CONCLUSION
A new DOA estimation method based on GSP has been proposed.Using the phase shift caused by the time delay in propagation, the space adjacency matrix is constructed to describe the relationship between array sensors; as signals are not phase shifted between different snapshots, we use a unit directed cyclic graph to represent this time series relationship.Then, by calculating the Kronecker product of the graphs above, we obtain the overall graph A ⊗ for array signals.Exploiting the relationship that the input signal x is an eigenvector of A N ⊗ corresponding to the unit eigenvalue, we establish a mapping of the incident angle θ to the GFT coefficients to estimate DOA.As demonstrated by simulations, the proposed method outperforms an existing GSP-based DOA estimation algorithm in various environments for both single and multiple sources, and it seems to be more robust against noise than MUSIC when SNR is very low.

Fig. 3 .
Fig. 3. Graph representation of signals at different sensors.Its corresponding space adjacency matrix is given as
• , 0 • for two source signals; the source signals are random following the standard normal distribution, and noise n(t) is additive white Gaussian.