DOA estimation for periodically modulated sources

This paper considers the problem of direction-of-arrival estimation for periodically modulated signals using one uniform linear array of sensors. By means of modulating the sources with periodic modulation sequences, we can form a series of linear equations relating the autocorrelation matrices of the received data and the outer products of the scaled steering vectors. Solving these linear equations yields a group of Hermitian matrices formed from the outer products of the scaled steering vectors. Then taking the eigendecomposition of these Hermitian matrices, we can obtain all the scaled steering vectors. By utilizing a special structure of the scaled steering vectors, we can find the directions of signals impinging on the array. We also examine the relation of the modulation sequences and the estimation performance, and a design of the modulation sequences to resist the effect of spatial noise is proposed. One merit of the proposed method is that it can be used in the scenarios of more sources than sensors. The simulation result also shows that it has the capacity to distinguish the closely spaced sources.


Introduction
The subject of array signal processing is concerned with the extraction of information from signals collected using an array (or arrays) of sensors [1,2]. One important information is the direction of arrival (DOA) of the incident signals. Take wireless communications for example, the information of DOA can be used for mobile localization for directional transmission in the downlink [3]. Hence, the DOA estimation of sources is one important research topic, and various algorithms in this field over the past decades have been proposed [1][2][3][4][5][6][7]. One of the famous algorithms is the multiple signal classification (MUSIC) algorithm proposed in [4]. The merit of the algorithm is that the accuracy of estimation can be obtained for large data samples or at high signal-to-noise ratio (SNR) scenarios. Another famous algorithm is the deterministic method developed by Van Der Veen [5]. Instead of requiring the statistical data and the search procedure in the angular spectrum inherently in the MUSIC algorithm, the deterministic method estimates the DOA directly in terms of eigenvalues of a certain matrix obtained from the received data. Due to the limitation of space, we cannot Correspondence: yischen@fcu.edu.tw Department of Communications Engineering, Feng Chia University, Taichung 40724, Taiwan introduce all algorithms of DOA estimation, and we refer the interested readers to [1,2,7] and the references therein for a detailed review.
In digital communications, although the actual transmitted symbol stream is unknown to the receiver, some a priori information of the transmitted signals, for example, the modulation scheme, is available to the receiver. The receiver can take advantage of this extra information with the received data to carry out some tasks including source separation and channel estimation [8][9][10][11]. Particularly in [9,10], the signal sources transmitted are first multiplied at a symbol rate by known amplitude-variation sequences, called modulation sequences, to aid symbol recovery or channel estimation at the receiver. However, there is little research for DOA estimation using modulation sequences. This motivates our research interest in developing a new DOA estimation method for one uniform linear array (ULA) based on modulation sequences.
Our idea and method are shown as follows: By means of modulating the sources with periodic modulation sequences, we can form a set of autocorrelation matrices of the received data. Then the set of autocorrelation matrices allows us to formulate a series of linear equations relating the outer products of the scaled steering vectors. Solving the set of linear equations produces a group of http://asp.eurasipjournals.com/content/2013/1/110 Hermitian matrices, which are the outer products of the scaled steering vectors. Then taking the eigendecomposition of these Hermitian matrices, we can obtain all the scaled steering vectors. By utilizing a special structure of the scaled steering vectors, we can find the directions of signals impinging on the array. We also examine the relation of the modulation sequences and the estimation performance, and a design of the modulation sequences to resist the effect of spatial noise is proposed. The merit of the proposed method is that it can be used in the case of more sources than sensors. In addition, the method has the capacity to distinguish the closely spaced sources.
It is worth to mention that since the periodically modulated signals are artificial, the proposed method is suitable for communication signals. One possible application of the proposed method is mobile localization for directional transmission in the downlink since the modulation formats of the mobile units are available to the base station in the uplink [3,11].
This paper is organized as follows: Section 2 briefly reviews the system model and provides basic assumptions. In Section 3, we derive the estimation method and discuss some properties of the proposed algorithm. Simulation results are given in Section 4. Section 5 concludes this paper.

Notations
(·) * , (·) T , and (·) H denote the complex conjugate, transpose, and conjugate transpose operations, respectively. The notation · 2 is the 2-norm. The symbols R and C represent the set of real numbers and the set of complex numbers, respectively. I M is the identity matrix of dimension M×M. A•B is the Hadamard product of matrices A ∈ C m×n and B ∈ C m×n ( [12], p. 190), and A(:, k) and A(l, :) are the kth column vector and the lth row vector of A, respectively. For a vector b ∈ C m , b(r 1 : r 2 ) is the subvector formed from the r 1 th element to the r 2 th element of b. In addition, for any m × m matrix G = [ g k,l ] 0≤k,l≤m−1 , we define the operation j (G) = [ g 0,j g 1,j+1 · · · g m−1−j,m−1 ] T , for 0 ≤ j ≤ m − 1, i.e., j (G) is the column vector formed from the jth superdiagonal of G.

Problem statement
Consider a uniform linear array of M sensors where adjacent sensors are separated with equal distance d. The sensor array receives N narrowband sources from farfield whose directions of arrival are θ 1 , θ 2 , · · · , and θ N , respectively. Before transmission, each source s i (k), ∀i = 1, 2, · · · , N, is multiplied by a real and periodic modulation sequence c i (k). Then the standard system model is shown as follows: y(k) = As p (k) + w(k), (2.1) where y(k) ∈ C M is the received vector at time k. A = [a(θ 1 ) a(θ 2 ) · · · a(θ N )] ∈ C M×N is the array response matrix with the ith column being the steering vector and λ ≥ 2d is the signal wavelength. w(k) ∈ C M is the spatial noise vector at time k. s p (k) ∈ C N is the transmitted vector defined as follows: The purpose of this paper is to develop a method of estimating θ i for i = 1, 2, · · · , N, using the second-order statistics of the received data based on the following assumptions: (i) The source vector s(k) is a zero-mean, temporally and spatially uncorrelated, and wide-sense stationary vector with E[ |s i (k)| 2 ] = d 2 i . The noise vector w(k) is zero-mean, wide-sense stationary, and is the Kronecker delta function. In addition, the source signal is uncorrelated with the noise, i.e.,

DOA estimation
In this section, we first derive the estimation method when noise is absent. The design of the modulation sequences when noise is present is given in Section 3.2. Some further discussions about the proposed method are given in Section 3.3.

The proposed method
Using (2.2), the system model (2.1) for the noiseless case can be expressed as Taking the expectation of y(k)y(k) H for k = 1, 2, · · · , P and using assumption (i), we obtain the following P autocorrelation matrices: Since D and C(k) are real and diagonal matrices, (3.2) can be further expressed as (3.3) http://asp.eurasipjournals.com/content/2013/1/110 Here we let AD Since C(k) is a diagonal matrix formed from the periodic modulation sequences c 1 (k), c 2 (k), · · · , c N (k) with period P by assumption (iii), we know that C(k) 2 is also periodic with period P, i.e., C(k) 2 = C(k + P) 2 , which implies that R k+P = R k , for example, R 1+P = R 1 . In addition, for the purpose of DOA estimation, we need the following proposition to aid our derivation of the proposed method. Proposition 1. For k = 1, 2, · · · , P, the ith upper diagonal of R k can be expressed as follows: Proof. Please see Appendix.
With the aid of Proposition 1, we obtain the vectors (3.5), it is obvious that each column of V i can be written as a linear equation shown as follows: Since P > N by assumption (iii), we can appropriately design the modulation sequences {c n (1), c n (2), · · · , c n (P)}, n = 1, 2, · · · , N, to make the matrix W be full column rank. Then for each i, i = 0, 1, · · · , M−1, the least squares solutions of (3.6) are shown as follows: Then each column of X i , solved by (3.7), is used to form the matrix X i , ∀i = 0, 1, · · · , M − 1. Taking the transpose In addition, from (3.4), we know that By writing down the elements in Z i with the aid of the Hadamard product, we observe that (3.10) Since for each n, h n h H n is a Hermitian matrix, the vectors in (3.10) obtained from (3.7) allow us to form N rankone Hermitian matrices Q n = h n h H n , which is the outer product matrix of h n , n = 1, 2, · · · , N. Then each column vector h n of H is estimated up to a scalar ambiguity α n by computing the unit-norm eigenvector associated with the maximal eigenvalue of the matrix Q n , i.e., h n = h n α n , n = 1, 2, · · · , N. (3.11) Since h n = a(θ n )d n , we have (3.12) Then we divide each scaled steering vector h n into two subvectors, namely h n1 = h n (1 : M − 1) and h n2 = h n (2 : M). It is clear that h n1 = h n2 e jφ(θ n ) and e jφ(θ n ) can be obtained from the least squares solution. Then the DOA θ n is thus obtained from the angle of e jφ(θ n ) .

Remark 1.
If the array is only composed of two sensors, i.e., M = 2, then h n is a 2 × 1 vector, and it is obvious that h n (1, :) = h n (2, :)e jφ(θ n ) . In this case, we can divide the first entry of h n by the second one to obtain the DOA for n = 1, 2, · · · , N. In other words, the proposed method can carry out the DOA estimation using only two sensors in theory. http://asp.eurasipjournals.com/content/2013/1/110

The design of the modulation sequences
We have derived the estimation method in Section 3.1. We now discuss how to design the modulation sequences to combat the effect of noise on DOA estimation.
When noise is present, the system model (3.1) becomes y(k) = AC(k)s(k) + w(k) (3.13) and the autocorrelation matrices in (3.3) become (3.14) Since the noise variance only appears on the main diagonal of R k , the groups of vectors i (R k ), i = 1, 2, · · · , M − 1 in (3.4) remain unchanged, and only the group of vectors 0 (R k ) needs to be changed as Then the corresponding M linear equations from (3.16) become V 0 (:, k) = WX 0 (:, k) + σ 2 w 1 P , k = 1, 2, · · · , M. (3.17) From (3.7), we know that the solutions of (3.17) become where q = (W T W) −1 W T 1 P . From (3.18), we know that X 0 (:, k) is the actual X 0 (:, k) plus a perturbation term σ 2 w q due to noise. Since q is formed from the modulation sequences, we need to design the modulation sequences to minimize q 2 and the effect of the resulting perturbation term. However, the high nonlinearity of the modulation sequences contained in q makes it difficult to design. Hence, we adopt another reasonable design criterion which is also used in [10] to tackle this problem.
From (3.18), we know the X 0 (:, k) = X 0 (:, k) if and only if 1 P is orthogonal to every column of W, i.e., where w i = c i (1) 2 c i (2) 2 · · · c i (P) 2 T ∈ R P is the ith column of W, ∀i = 1, 2, · · · , N. If the modulation sequences can be selected to meet the orthogonality condition (3.19), the effect of noise is completely eliminated, but this is impossible since the entries of w 1 , w 2 , · · · , w N and 1 P are positive. Therefore, we seek to choose the modulation sequences such that 1 P is as close to being orthogonal to w i as possible, for i = 1, 2, · · · , N. To this end, we define the following correlation coefficients: and try to choose the modulation sequences to minimize the correlation coefficient γ i subject to the following two constraints: for i = 1, 2 · · · , N. Roughly, constraint (3.21) normalizes the power gain of the modulation sequence of each source to 1, and constraint (3.22) requires that at each instant, the power gain is no less than τ . Note that the optimization problem is identical to the case considered in [10], and the resulting optimal sequences are given by, for any fixed 1 ≤ m i ≤ P, for i = 1, 2, · · · , N. In the study of blind channel estimation using periodic modulation [10], the two-level sequence in (3.23) is also shown to be optimal for mitigating the channel noise effect. In addition, with the optimal solution in (3.23), the corresponding γ i is , ∀i = 1, 2, · · · , N. Note that γ opt decreases as τ decreases, and thus, the noise effect imposed on V 0 is reduced and hence estimation performance improves. We will give a simulation example to illustrate this property. From (3.23), we know that each of the N modulation sequences is a two-level sequence with a single peak in one period. However, to make the matrix W be full column rank such that the least squares solutions (3.7) can be computed, the peak locations of the N modulation sequences in one period need to be distinct with one another, i.e., m k = m l for all k = l. Without loss of generality, we can let m i = i for i = 1, 2, · · · , N. Remark 2. It is worth noting that even for the colored noise case, the design of modulation sequences is the same as the case of white noise. To see this, we know that if the zero-mean, wide-sense stationary noise vector w(k) in assumption (i) now becomes colored, then the M × M autocorrelation matrix R w = E[ w(k)w(k) H ] becomes a Hermitian and Toeplitz matrix with i (R w ) = σ 2 w (i)1 M−i , ∀i = 0, 1, · · · , M − 1. In addition, the autocorrelation matrix of the received vector in (3.14) becomes R k = http://asp.eurasipjournals.com/content/2013/1/110

E[ y(k)y(k) H ] = HC(k) 2 H H + R w , which means for
and The least squares solution of (3.25) is From (3.26), it is obvious that if W T 1 P = w 1 w 2 · · · w N T 1 P = 0, then the effect of noise can be eliminated. From here, the discussion and derivation are the same as the content from (3.19) to (3.23). Hence we know that the optimal sequences given in (3.23) can work well for the colored noise case. We will give a simulation in Section 4 to demonstrate this feature.

Discussions
We now discuss some notable features of the proposed method. First, from the result at the end of Section 3.1, we know that after h n is estimated, the DOA angle θ n can be obtained from the linear equation h n1 = h n2 e jφ(θ n ) . Since h n is an M × 1 vector, the DOA angle θ n can be acquired provided that M ≥ 2, where M is the number of sensors for the ULA. Hence, the proposed method can carry out DOA etimation not only for the case of more sensors (M ≥ N), but also for the case of less sensors (M < N), as long as the number of sensors M ≥ 2. Second, the DOA θ n ∈[ − π 2 , π 2 ], n = 1, 2, · · · , N, may not be distinct with each other since from Section 3.1, we know that the estimates of θ 1 , θ 2 , · · · , and θ N are independently obtained from the corresponding scaled steering vectors h 1 , h 2 , · · · , and h N , respectively. Hence, the proposed method possesses the capacity to distinguish the closely spaced sources. We will give a simulation to demonstrate this feature. The third feature of the proposed method is that it provides a design of the modulation sequences to minimize the effect of noise on DOA estimation and thus improves the accuracy of the solution.
We now summarize the proposed approach as the following algorithm: 1. Collect the received data {y(m)} S m=1 , where S divides P, the period of the modulation sequences.
2. Compute the autocorrelation matrices R 1 , R 2 , · · · R P , via the following time average: 3. Use the autocorrelation matrices in (3.27) to form V 0 , V 1 , · · · , V M−1 , and use the designed modulation sequences (3.23) to form W. 4. Obtain the matrices X 0 , X 1 , · · · , X M−1 using the least squares solutions in (3.7) with the aid of the matrices V 0 , V 1 , · · · , V M−1 , and W obtained from the previous step. 5. Form N rank-one Hermitian matrices Q 1 , Q 2 , · · · , Q N from Z i = X T i , i = 0, 1, · · · , M − 1, with the aid of (3.10). 6. For each Q n = h n h H n , compute the estimate h n as the unit-norm eigenvector associated with the maximal eigenvalue of Q n , ∀n = 1, 2, · · · , N.

Simulation
In this section, we use several simulations to demonstrate the performance of the proposed method. For all simulation examples, the SNR is defined as SNR = E y(k)−w(k) 2 2 E w(k) 2 2 . We use the root-mean-square error (RMSE) of angles as the performance measure, which is defined as RMSE = E 1 N N n=1 (θ n − θ n ) 2 , where θ n is the estimate of θ n . The number of Monte Carlo trials is 500. The source symbols are independent and identically distributed binary phase-shift keying signals. The noise is zero-mean and white Gaussian (except for simulation 3).

Simulation 1 -underdetermined DOA estimation
In this simulation, we execute two underdetermined experiments where the number of sensors is less than the number of sources. There are six signal sources in both experiments with the corresponding DOA: For both experiments, we set d = 0.5λ and we adopt the modulation sequence c i (k) in (3.23) for source i with period P = 7 and peak location m i = i, i = 1, 2, · · · , 6. The number of symbol blocks is S = 700. In the first experiment, we set SNR = 10 dB, and the number of sensors is increasing from two to five. Simulation results in Figure 1a show that increasing the number of sensors http://asp.eurasipjournals.com/content/2013/1/110 improves the performance of DOA estimation. This may be due to more information obtained to average out the computational error and noise effect for more sensors.
In the second experiment, the ULA is only composed of four sensors. From Figure 1b, we observe that the RMSE of angles decreases as SNR increases. In addition, both experiments also show that the estimation performs better for smaller τ , which is consistent with our analysis in Section 3.2.

Simulation 2 -comparison with existing methods
In this simulation, we examine the performance of the proposed method with those of the conventional MUSIC algorithm [4] and the deterministic method [5] for two sources. The ULA is composed of five sensors with d = 0.5λ. For the proposed method, the modulation sequence c i (k) for source i is chosen as in (3.23) with period P = 3, τ = 0.2, and peak location m i = i, i = 1, 2. The number of symbol blocks is S = 900. We first consider a case of http://asp.eurasipjournals.com/content/2013/1/110 nonclosely spaced sources: {θ 1 , θ 2 } = {−15 • , 15 • }. Simulation results in Figure 2a show that the proposed method performs better than both methods in [4,5] except for the deterministic method for SNR > 15 dB. Then we consider the case of two closely spaced sources: {θ 1 , θ 2 } = {−15 • , −12 • }. Since the MUSIC algorithm fails to resolve the closely spaced signals under the simulation setting, we only compare the proposed method with the deterministic method. From Figure 2b, we observe that the proposed method performs better than the deterministic algorithm for this case. It also shows that the capacity of the proposed method to distinguish two closely spaced sources is good.

Simulation 3 -underdetermined DOA estimation in the presence of colored noise
In this simulation, we examine the performance of the proposed method with that of the Khatri-Rao (KR) subspace method [6] for the colored noise case. The additive colored noise w(n) is generated by passing a zero-mean with unit variance white sequence w v (n) through a finite impulse response filter c(z) = 1.2+0.6z −1 +0.3z −2 whose output is w(n) = c(z)w v (n). Assume that there are four signal sources with {θ 1 , θ 2 , θ 3 , θ 4 } = {−60 • , −30 • , 30 • , 60 • } and the ULA is composed of three sensors with d = 0.5λ. For the proposed method, the modulation sequence c i (k) for source i is chosen as in (3.23) with period P = 5, τ = 0.2, and peak location m i = i, i = 1, 2, 3, 4. Since the KR subspace method is suitable for the quasi-stationary signals, we randomly choose the standard deviation of signals following a uniform distribution on [ 0, 2 √ 3] such that the variance is 1. In addition, 50 frames with each period being 100 are used to carry out the simulation for the KR subspace method. In other words, the number of symbol blocks is S = 5, 000. From Figure 3, we observe that the KR subspace results in better performance for SNR < 7 dB. However, for SNR > 7 dB, the proposed method achieves smaller RMSE of angles than that of the KR subspace method.

Conclusion
This paper has proposed a new DOA estimation algorithm for one ULA based on periodic modulation. The proposed method has three notable features. First, the proposed algorithm can handle more sources than sensors, which may be few as two. In addition, the great capacity to distinguish the closely spaced sources is the second feature of the proposed method. The final feature of the proposed method is that the performance of the estimation algorithm depends on the choice of the modulation sequences to resist the noise effects. Hence, we can properly choose the modulation sequences to improve the performance of estimation. Simulation results are used to demonstrate the performance of the proposed method and to compare it with some existing methods. proposed method method in [4] method in [5] (a) 5 10 15 20 25 30 10 -2 10 -1 10 0 10 1 10 2 SNR (dB) RMSE of angles (degree) proposed method method in [5] (b)