A Searchfree DOA Estimation Algorithm for Coprime Arrays
Abstract
Recently, coprime arrays have been in the focus of research because of their potential in exploiting redundancy in spanning large apertures with fewer elements than suggested by theory. A coprime array consists of two uniform linear subarrays with interelement spacings and , where and are coprime integers and is the wavelength of the signal. In this paper, we propose a fast searchfree method for directionofarrival (DOA) estimation with coprime arrays. It is based on the use of methods that operate on the uniform linear subarrays of the coprime array and that enjoy many processing advantages. We first estimate the DOAs for each uniform linear subarray separately and then combine the estimates from the subarrays. For combining the estimates, we propose a method that projects the estimated point in the twodimensional plane onto onedimensional line segments that correspond to the entire angular domain. By doing so, we avoid the search step and consequently, we greatly reduce the computational complexity of the method. We demonstrate the performance of the method with computer simulations and compare it with that of the FDroot MUSIC method.
keywords:
DOA estimation, coprime arrays, coprime sampling, uniform linear array, searchfree1 Introduction
Directionofarrival (DOA) estimation using sensor arrays is a problem that is frequently encountered in many engineering areas including radar, sonar, and wireless communication, and it has been studied for several decades. Recently, the notion of coprime array signal processing has emerged as an area of interest where, as the name suggests, the processing is applied to signals acquired by coprime arrays (1); (2); (3); (4). In general, it has been shown that coprime sampling allows us to sample a signal sparsely and estimate parameters of signals at higher resolutions (2). The parameters can be the frequencies of temporal signals or the DOAs of spatial signals. By sparse sampling, one can reduce the rate of analogtodigital converter for sampling in the temporal domain (4) or reduce the number of array sensor elements for sampling in the spatial domain (1). This entails that systems where coprime sampling is applied can have lower cost than systems with Nyquist sampling, yet without degrading the performance of the latter. For example, the cost of ADC usually grows exponentially with the bandwidth or the sampling rate. By using coprime sampling, one can reduce the requirements on the analog frontend by significantly lowering the sampling rate (4).
In general, the idea behind coprime arrays is to extend the concept of minimally redundant sensor arrays and span large apertures by using far fewer elements than prescribed by textbook antenna theory. A coprime array can be constructed with two uniform linear arrays with and elements with interelement spacings and , respectively, and where and are coprime integers, i.e., integers whose greatest common divisor is one. It has been shown that in the presence of spatially stationary fields, one can match the information extracted with such arrays to that of fully populated arrays with elements (2). As already pointed out, a very important payoff for the use of coprime arrays is simplified and reduced cost array designs. Instead of having an array with elements, one can use an array with elements (where one element is shared).
The reduced cost of coprime arrays represents a strong motivation for studying DOA methods that can process signals for such arrays. An obvious approach is to consider methods that are developed to operate on arbitrary geometry arrays. In this class there are several popular DOA estimation methods including beamspan (5), Capon (6), MUSIC (7) and several versions of maximum likelihood estimators (8). For finding the optimal estimate, however, they all require a computationally expensive search step in a nonconvex space. We note that for uniform linear arrays (ULAs), the array response vector has a Vandermonde form, which allows the search step to be replaced by polynomial rooting. These so called searchfree algorithms include IQML (9), MODE (Method of Direction Estimation) (10), rootMUSIC (11) and ESPRIT (12). An interesting alternative to the methods for arbitrary arrays are the approaches that transform the arbitrary arrays to equivalent ULAs so that one can take the advantage of efficient algorithms for ULAs, e.g., array interpolation (13), manifold separation (14) and FourierDomain (FD) rootMUSIC (15). These approaches basically try to approximate the steering vectors by virtual ULAs. Their drawback is that for achieving satisfactory performance, they require a large number of virtual arrays, which entails increased complexity. Furthermore, they suffer fixed error at high signaltonoise ratios (SNRs) and thus, do not converge to the CramérRao bound as the SNR increases.
The DOA estimation with ULAs has been thoroughly studied, but there has not been much work for designing special algorithms that take advantage of the partial uniform linear structure of coprime arrays. In (1), the authors apply the MUSIC algorithm to coprime processing, where they mainly consider the problem from the perspective of degree of freedom, i.e., how many signals the coprime array can identify. To achieve additional degree of freedom with coprime arrays, that is, to enhance the rank of the covariance matrix of the received signal, they propose the use of spatial smoothing. Nevertheless, the proposed method is not searchfree, and its accuracy and computational complexity have not been investigated in detail. In (4), where processing of time series is addressed, the authors use an iterative method to search for the optimal frequency. Although it is searchfree, it does not use the coprime property and it cannot handle the case of multiple signals.
In this paper, we address the problem of DOA estimation with coprime arrays with the emphasis on reduced computational complexity while preserving estimation accuracy. We propose a searchfree method that exploits the uniform linear structure of the subarrays of the coprime arrays. Considering that the DOA estimation using ULAs is fast and accurate, we first estimate the DOA for each subarray separately. Since the interelement distance for each subarray is larger than halfwavelength of the signal, ambiguity is present and as a result, the DOAs cannot be uniquely determined. To resolve the ambiguity, we use a projectionlike method that combines the estimates from different subarrays, where the correctness of the estimate is guaranteed as a consequence of the coprimeness.
The paper is organized as follows. In Section 2, we introduce the model and briefly review the DOA estimation with ULAs. The proposed algorithm is discussed in detail in Section 3 and the numerical experiments are presented in Section 4. Section 5 provides final remarks.
We use the following notation: and denote the transpose and the conjugate transpose, respectively; refers to the complex space; signifies the Kronecker delta; stands for the expectation operator; and is the trace of the matrix .
2 DOA estimation for ULA
The problem of finding the DOAs of narrowband plane waves impinging on a ULA of sensors can be modeled as follows (16):
(1) 
where is the signal received by the sensors at the th time slot; is the number of snapshots; contains the complex envelopes of the emitter signals; is a noise process; and
(2) 
with
(3) 
is the distance from the th sensor to the reference point; is the DOA of our interest. The number of sources, , is assumed to be known. Thus, is the vector of DOAs we wish to estimate. We assume the signal and the noise are independent zeromean complex Gaussian random processes with the following moments:
(4)  
(5) 
with being the identity matrix, , the signal covariance matrix, and , the noise power. The covariance matrix of the received signal can be written as
(6) 
We denote by the sample covariance matrix
(7) 
The eigendecomposition of in (7) can be written as
(8) 
where and are diagonal matrices that contain the eigenvalues of the signal and the noise subspaces, respectively, whereas and are composed of the eigenvectors of the signal and the noise subspaces, respectively.
For ULAs, the steering vector becomes
(9) 
There exists a large number of fast algorithms for estimating DOAs with ULAs, among which MODE is the leading candidate (17); (18). Since we later use MODE as part of our algorithm, we review it here briefly. The MODE algorithm was originally proposed in (10). It estimates the DOAs via the minimization of the following function:
(10) 
where
(11)  
(12)  
(13)  
(14) 
and is a Toeplitz matrix defined by
(15) 
Let , and denote by the complexvalued polynomial
(16)  
(17) 
The angles of the roots of the polynomial are the estimates of the DOAs. It was shown in (10) that the MODE algorithm is an asymptotically efficient estimator of the DOAs. As per (18); (17), convergence of MODE is guaranteed. To minimize (10), both iterative and noniterative steps have been proposed. See (17); (18) for detail.
3 DOA estimation for coprime arrays
We consider a coprime array with two uniform linear subarrays. We assume that the two subarrays have interelement spacing and , respectively, with and being coprime. Fig. 1 shows the case for and . Because the two subarrays share the first sensor, the total number of sensors is equal to . The corresponding steering vector in (2) becomes
(18) 
For such arrays, the efficient methods for ULAs are not directly applicable. We can of course use the algorithms for arbitrary arrays, but they are slow and the uniform linear structure of the subarrays is wasted. Our objective is to develop an algorithm that is applicable to coprime arrays while at the same time enjoys the efficiency brought by the uniform linear structure of the subarrays. The proposed approach is composed of two steps. We first estimate the DOAs from each subarray separately, which is similar to the approach in (4). We then combine the two estimates in an innovative way.
3.1 DOA estimation using a single subarray
For a subarray, the interelement spacing is with . This has the same effect as we undersample a signal by a factor of in the temporal domain. According to the basic sampling theorem, we will have ambiguities, or aliasing, in the angular domain. Specifically, if we use the same algorithm as described above, the polynomial in (17) becomes
(19) 
This polynomial has roots instead of . If is a root of (17), it is still a root of (19). The problem is that aliasing occurs with the period of ; for are also the roots of (19), as shown in Fig. 2. Therefore, we have no way of distinguishing the we want from the remaining roots.
Note that in the first step, it does not matter which algorithm we use. In the experiment, we selected MODE, which was shown to have the best performance for ULAs (18).
3.2 Combination of the estimates from subarrays
The present ambiguities in the estimates from the two subarrays notwithstanding, we show that we are able to disambiguate them. For example, in the case of one source and no noise, in Fig. 2, we can see that the estimates from the two subarrays coincide at the red arrows. Therefore, we can be very confident that the position of the red arrow is the DOA. In the presence of noise, however, coincidences as the one in Fig. 3 are unlikely, and consequently, it becomes difficult to tell the exact DOA. Furthermore, in the case with multiple sources, obtaining the DOAs in this way becomes impossible.
We propose a projection method that can uniquely determine the exact DOA. Recall that the aliasing period in the angular domain for the subarray is . Thus, we only need to take the values between . The ranges of the outputs of the two subarrays define a rectangle in the twodimensional plane as shown in Fig. 4, where and . It is not difficult to see that the 45 degree oblique line segments colored in blue correspond to the entire angular domain. The Chinese remainder theorem guarantees that the map is onetoone and onto as a result of the coprimeness. This is the key point of our algorithm. In Fig. 4, the entire angular domain consist of four line segments; L1 corresponds to ; L2, ; L3, ; and L4, . For general and , there will be oblique line segments. Suppose the outputs of the first and the second subarray are and , respectively. The two estimates specify a point in the plane. We project the point onto the oblique line segments that corresponds to the entire angular domain. Specifically, we seek the point on the oblique line segments that is nearest to that point specified by the two estimates. This ensures that the combination is optimal and the result is valid. For point B in Fig. 4, we can simply draw a line that is vertical to the oblique line and measure the distance from B to the intersections. For point A, however, one intersection falls outside the rectangle area and should be wrapped around. As a result, the intersection is actually on L3 in the figure. After the nearest point on the line segments is found, it is easy to calculate the corresponding DOA by solving a set of modular linear equations. As an alternative to solving the equations, one may have a precalculated lookup table with solutions of the modular equations.
The main difference between the proposed method and the one in (4) is that here we take advantage of the coprimeness and use a projection process in the search for optimal DOAs, while the method in (4) iteratively tries different combinations of the estimates and find the best one among them. Although solving the modular equations in our method requires an iterative process, it can be avoided by using a lookup table.
3.3 Combination of the estimates in the case of multiple sources
In the case of multiple sources, the method gets more complicated. Suppose that the outputs of the first subarray are and and the outputs of the second subarray, and . We have no way of knowing which value corresponds to the first and which to the second target. Therefore, we have to check all the four points. In Fig. 4, , , and define four points on the twodimensional plane, which are marked A, B, C, and D. Since candidate points are obtained, we can readily evaluate the likelihood probability for these points and take the best ones.
In short, the proposed algorithm can be summarized as follows:

Estimate the angles for each subarray.

Plot the corresponding points on the 2dimensional plane.

Project these points to the nearest diagonal line segments.

Map the projected points back onto the entire angular domain.

Evaluate the likelihood probabilities at the points.

Select points with the largest likelihood probabilities.
3.4 Example
In this subsection, we provide a simple example to better illustrate the proposed method. We consider the case for , and the number of targets . Suppose the outputs of the first subarray are and , and those of the second subarray are and . We denote by and the final results of the estimation. Using the estimates of the subarrays, we can locate the following four points on the twodimensional plane: A, B, C and D. Also, we can see by simple calculation that

A is closer to L1 with distance ;

B is closer to L4 with distance ;

C is closer to L3 with distance ;

D is closer to L4 with distance .
We then project these points back onto the diagonal line segments. For example, B is projected onto L4 and the projected point, say , is . Subsequently we construct the modular equations
(20)  
(21) 
and the solutions will be the value for on the entire angular domain. The equations can be solved using the extended Euclidean algorithm. Alternatively, as we stated before, we can also use a lookup table to solve the equations. After we obtain the four values for the four points, we can then evaluate the likelihood probability at each point and take the largest two values as the final estimates.
3.5 Discussion
Asymptotic performance
In this subsection, we briefly discuss the performance of the proposed algorithm. Our method basically consists of two steps. First, we estimate the DOA separately for each subarray, and second, we combine the estimates by projecting the combined point onto a line segment. For the first step, any classical DOA estimation method can be used. In our simulation, we applied the MODE method, which is a type of maximum likelihood estimators (18).
It is well known that maximum likelihood estimators are asymptotically unbiased and normally distributed. The asymptotical property of the proposed algorithm depends on the SNR level. For high SNRs, the point on the plane defined by the two estimates of the subarrays will probably lie very close to the correct point. After projection, it could be even closer. This can be seen on Fig. 5. The estimation errors from the subarrays are indicated by and , which moves the combined point away from the correct one . In the second step, we project the point back onto the line segment and obtain . We can see that the distance between and is shorter than that of and . It is obvious that as long as does not lie outside the blue threshold, which is the line between the two line segments, is always closer to than is to . Moreover, we notice that the projection is linear. Since the estimates from the subarrays are asymptotically unbiased, we can expect that the final estimate is also unbiased. In other words, the error would not propagate for high SNRs.
Under low SNRs, however, is more likely to fall on the other side of the blue threshold. As a result, will be incorrectly projected to a different line segment and therefore, the final estimate will be poor. This explains why a threshold effect can be seen in the simulations. We point out that this effect is seen in many other unwrapping phasebased methods. After the SNR goes above certain level, the mean square error (MSE) decreases rapidly. We defined four areas in Fig. 5. Given as shown in the figure, if the point lies in Area 1, we have ; for Area 2, we have and ; for Area 3, we have ; for Area 4, is projected to the other line segment, and therefore the error is magnified. We have and .
Maximum number of detectable sources
For ULAs, it is well known that an array with elements can identify up to targets. Thus, the two subarrays in the coprime array can identify and targets, respectively. Consequently, there will be intersections on the twodimensional plane, which entails that the maximum number of sources that a coprime array can identify is .
We point out that our algorithm is particularly suitable for the case of a single target. For multiple sources, it may produce poor results because of the resolution problem in the processing step of each subarray. For ULAs, when two sources are very close to each other, the estimates may be poor. For coprime arrays, it is very likely that the true values of the sources overlap on the wrapped angular domain of one of the subarrays even if on the entire angular domain they do not. In that case, the estimates from the subarray are not accurate at all, and consequently the final estimates are unreliable.
4 Simulations
In this section, we use numerical experiments that demonstrate the performance of the proposed method. We show results of simulations that compare the method with that of the the FD rootMUSIC method as well as with the CramérRao bound (CRB). We note that the FD rootMUSIC method uses Fourier series coefficients to approximate the nullspectrum function (15), and we use it because it was shown to have better performance than other searchfree methods (19).
In the simulations, we used two sets of arrays, one with , and one with elements. For each coprime arrays, we test the cases of one source and three sources.
4.1 In the case of one source
First, we simulated one source with . In the first two experiments, we studied the performance of the method on the coprime array . In the first of them, we fixed the number of snapshots to , and we varied the SNR from dB to dB. For the FD rootMUSIC, we used and Fourier series coefficients, respectively. We computed the mean square error (MSE) of the methods by using 2000 independent realizations. The results and the CRB are presented in Fig. 6. We can see that at SNRs below dB, the FD rootMUSIC outperforms the proposed method. For SNRs greater than dB, the performance of the proposed method approaches the CRB while the FD rootMUSIC suffers from fixed errors due to its approximations. As the number of Fourier series coefficients increases, as expected, the performance of the FD rootMUSIC method improves. However, the method cannot reach the CRB.
In the second experiment, we kept the SNR at dB, and varied the number of snapshots from 100 to 1000 in steps of 100. The remaining parameters of the experiment were the same as before. From Fig. 7 we see that when the number of snapshots is 100, the MSE of the proposed method is far from the CRB. The results are shown in Fig. 7. We see that the FD rootMUSIC is better than our method when the number of snapshots or smaller. On the other hand, our method approaches the CRB when whereas the FD rootMUSIC does not improve for .
In the remaining two experiments, we repeated the simulations with a coprime array with elements. Again, first we kept and varied the SNR from dB to dB. The results of the MSE’s and the CRB are shown in Fig. 8. We observe similar patterns as in Fig. 6. However, now the proposed method reaches the CRB at dB. For SNRs greater than dB, the difference in performance between the two methods is by orders of magnitude and in favor of the proposed method.
Finally, in the last experiment, we kept the SNR at dB, and varied the number of snapshots from 100 to 1000, in steps of 100. The MSEs of the methods and the CRB are displayed in Fig. 9. We observe, that our method follows the CRB as the number of snapshots increases, whereas the MSE of FD rootMUSIC remains constant.
4.2 In the case of three sources
In this subsection, we repeated the simulations under the presence of three targets. We set the angles to be . Similarly we show four figures here. The results of the coprime array with are shown in Fig. 10 and Fig. 11. Fig. 10 shows the MSEs for different values of SNR with . We observe similar patterns as in Fig. 6 and Fig. 8. The proposed method reaches the CRB at about dB. For SNRs greater than dB, the difference in performance between the two methods is by orders of magnitude and in favor of the proposed method. Also we notice that the MSEs of the FD rootMUSIC methods with are almost the same. That is, with the increase of the approximation order, the performance does not really improve. Fig. 11 shows the MSEs for different values of with SNR= dB. Our method follows the CRB when the number of snapshots is larger than 400, whereas the MSEs of FD rootMUSIC remains constant. The results of the coprime array with are shown in Fig. 12 and Fig. 13. Similarly, Fig. 12 shows the MSEs for different values of SNR with . The proposed method outperforms the FD rootMUSIC methods under both low and high SNRs. This is because for arrays with large size, 200 Fourier coefficients is far from enough to provide good approximation. Fig. 13 shows the MSEs for different values of with SNR= dB. Our method follows the CRB when the number of snapshots is larger than 900, and its performance is much better than the FD rootMUSIC method with .
Moreover in Fig. 6, Fig. 8, Fig. 10 and Fig. 12, we can observe the threshold effect that we mentioned in Section 3. For example, in Fig. 8 the MSE of the proposed method rapidly decreases around SNR=dB; in Fig. 10 the threshold effect happens between SNR being dB and dB. These justify our discussion in Section 3.
In summary, we conclude that the proposed method compares favorably with the FD rootMUSIC method. It is important, however, that we also make a computational comparison between the methods. A drawback of the FD rootMUSIC method is that it requires the use of a large number of Fourier coefficients to achieve satisfactory performance. In (15) it was suggested that, in general, , the number of Fourier coefficients, should be approximately equal to , where is the signal wavenumber, and is the largest distance between the array sensors and the origin of the coordinate system. The complexity of FD rootMUSIC is (15). Thus, to achieve satisfactory performance, the FD rootMUSIC method is computationally expensive even for arrays with small sizes. On the other hand, using the result of MODE in (17), we can obtain that the complexity of the proposed method is , which is smaller than that of FD rootMUSIC.
5 Concluding remarks
In this paper, a fast searchfree DOA estimation algorithm for coprime arrays is proposed. It exploits the uniform linear structure of the subarrays by first estimating the DOA separately within each subarray and then combining the estimates from different subarrays. Simulation shows that the performance of the proposed algorithm is close to the FD rootMUSIC method at low SNRs and much better at high SNRs. We point out that the proposed method has significantly lower complexity.
Here, we only discussed the case with two subarrays. For a coprime array with more than two subarrays, the state space will be a multidimensional hyperrectangle, but the entire angular domain can still be mapped to line segments and the same idea for estimation can be reproduced. Therefore, the presented results can be extended to coprime arrays with multiple subarrays.
Appendix A Vitae
Zhiyuan Weng was born in Shanghai, China. He is a Ph.D candidate in the Department of Electrical and Computer Engineering at Stony Brook University. His research interests include statistical signal processing and machine learning.
Petar M. Djurić received his B.S. and M.S. degrees in electrical engineering from the University of Belgrade and his Ph.D. degree in electrical engineering from the University of Rhode Island. Currently, he is a professor in the Department of Electrical and Computer Engineering at Stony Brook University. He was an associate editor of several signal processing journals including IEEE Transactions on Signal Processing. In 2007, he received the IEEE Signal Processing Magazine Best Paper Award and, during 2008 2009, he was Distinguished Lecturer of the IEEE Signal Processing Society. In 2012, he received the European Association for Signal Processing Technical Achievement Award. He is a Fellow of the IEEE.
Footnotes
 An earlier version of this paper can be found at http://arxiv.org/abs/1301.4155.
References
 P. Pal, P. P. Vaidyanathan, Coprime sampling and the MUSIC algorithm, in: Digital Signal Processing Workshop and IEEE Signal Processing Education Workshop (DSP/SPE), 2011 IEEE, IEEE, 2011, pp. 289–294.
 P. P. Vaidyanathan, P. Pal, Sparse sensing with coprime samplers and arrays, Signal Processing, IEEE Transactions on 59 (2) (2011) 573–586.
 P. P. Vaidyanathan, P. Pal, Theory of sparse coprime sensing in multiple dimensions, Signal Processing, IEEE Transactions on 59 (8) (2011) 3592–3608.
 S. Negusse, P. Händel, The dual channel sinewave model: Coprime sparse sampling, parameter estimation, and the CramérRao bound, Measurement 45 (9) (2012) 2254–2263.
 J. Mayhan, L. Niro, Spatial spectral estimation using multiple beam antennas, Antennas and Propagation, IEEE Transactions on 35 (8) (1987) 897–906.
 J. Capon, Highresolution frequencywavenumber spectrum analysis, Proceedings of the IEEE 57 (8) (1969) 1408–1418.
 R. Schmidt, Multiple emitter location and signal parameter estimation, Antennas and Propagation, IEEE Transactions on 34 (3) (1986) 276–280.
 P. Stoica, A. Nehorai, Performance study of conditional and unconditional directionofarrival estimation, Acoustics, Speech and Signal Processing, IEEE Transactions on 38 (10) (1990) 1783–1795.
 Y. Bresler, A. Macovski, Exact maximum likelihood parameter estimation of superimposed exponential signals in noise, Acoustics, Speech and Signal Processing, IEEE Transactions on 34 (5) (1986) 1081–1089.
 P. Stoica, K. C. Sharman, Maximum likelihood methods for directionofarrival estimation, Acoustics, Speech and Signal Processing, IEEE Transactions on 38 (7) (1990) 1132–1143.
 B. Rao, K. V. S. Hari, Performance analysis of rootMUSIC, Acoustics, Speech and Signal Processing, IEEE Transactions on 37 (12) (1989) 1939–1949.
 R. Roy, T. Kailath, ESPRITestimation of signal parameters via rotational invariance techniques, Acoustics, Speech and Signal Processing, IEEE Transactions on 37 (7) (1989) 984–995.
 B. Friedlander, A. J. Weiss, Direction finding using spatial smoothing with interpolated arrays, Aerospace and Electronic Systems, IEEE Transactions on 28 (2) (1992) 574–587.
 F. Belloni, A. Richter, V. Koivunen, DoA estimation via manifold separation for arbitrary array structures, Signal Processing, IEEE Transactions on 55 (10) (2007) 4800–4810.
 M. Rubsamen, A. B. Gershman, Directionofarrival estimation for nonuniform sensor arrays: from manifold separation to Fourier domain MUSIC methods, Signal Processing, IEEE Transactions on 57 (2) (2009) 588–599.
 H. Krim, M. Viberg, Two decades of array signal processing research: the parametric approach, Signal Processing Magazine, IEEE 13 (4) (1996) 67–94.
 J. Li, P. Stoica, Z. S. Liu, Comparative study of IQML and MODE directionofarrival estimators, Signal Processing, IEEE Transactions on 46 (1) (1998) 149–160.
 H. L. Van Trees, Optimum array processing, WileyInterscience, 2002.
 A. B. Gershman, M. Rübsamen, M. Pesavento, One and twodimensional directionofarrival estimation: An overview of searchfree techniques, Signal Processing 90 (5) (2010) 1338–1349.