Joint Estimation of DOA and Frequency of Multiple Sources with Orthogonal Coprime Arrays

A two-stage method is proposed to jointly estimate the direction-of-arrival (DOA) and carrier frequency (CF) of multiple sources, by using two orthogonal coprime arrays (CPAs). The DOAs of CF-known sources are estimated first by applying a spatial smoothing MUSIC algorithm. The contribution of these source signals is then removed from the originally received signal by applying an orthogonal complement projector. Next, a joint-ESPRIT algorithm is applied to estimate the DOAs and CFs of the remaining CF-unknown sources. With two orthogonal CPA(5, 6), the RMSE of DOA and CF of applying the proposed method to 30 sources, 13 of which have unknown CF, is less than 1% at SNR >5 dB.


Introduction
The direction-of-arrival (DOA) of source signals is a prerequisite in many applications like cognitive radar and antenna beamforming. Typical algorithms for estimating the DOAs include subspace-based MUSIC [1], ESPRIT [2] and sparsity-based compressive sensing (CS) approach [3], among others. However, these algorithms are applicable only when the carrier frequency (CF) of signal sources is known [1][2][3][4]. For subspace-based algorithms, MUSIC [1] can generate multiple steering vectors and calculate their projection norm in the null space of receiving signals. On the other hand, ESPRIT [2] can calculate the phase delays if the sensor array possesses rotational invariance, then the DOA can be estimated under a known carrier frequency. As for compressed sensing-based algorithms [4], minimization of 1 -norm is imposed on the sensing matrix to determine the sparse coding. However, due to Doppler shift [5] or in applications of radar tracking [6] and cognitive radios [7,8], the carrier frequency of signal sources can be ambiguous.
Recently, several algorithms were proposed to jointly estimate the DOA and CF of source signals [7][8][9][10]. In [9], a joint angle-frequency estimation (JAFE) algorithm was proposed, in which both the CF and DOA were iteratively updated from their previous values. In [7,8], a compressed carrier and DOA estimation (CaSCADE) algorithm was proposed. A modulated wideband converter (MWC) was used to estimate the CF at sub-Nyquist rate. Then, a CS approach, a parallel factor (PARAFAC) approach and a joint-ESPRIT (JE) were applied to jointly estimate the CF and DOA. In [10], a two-stage method on an uniform linear array (ULA) was proposed. The frequency components of the received signal were determined by applying Fourier transform and spectrum peak searching. Candidate DOAs at each frequency component were estimated with subspace-based methods like MUSIC. After all the CFs and DOAs were estimated, an exhaustive matching algorithm was applied to match each CF to one DOA, at high computational cost. In practice, the received signal is contributed by a combination of CF-known sources and CF-unknown sources. The performance of estimation algorithms like JE degrade as the number of signal sources increases [7,8].
A Khatri-Rao (KR) subspace was proposed to extend the degrees of freedom (DOFs) of an N-element ULA from N − 1 to 2N − 2 so that DOAs of more signal sources can be estimated [11]. In [12,13], KR subspace-based algorithms were proposed to reduce the uncertainty in gain and phase of sensors. Different array configurations such as nested array [14], coprime array (CPA) [15,16] were proposed to further increase the DOFs. In [4], two types of generalized CPA structure were reviewed and their DOFs were estimated. In [17], the DOF of a CPA was extended by scavenging the fourth-order statistics of the received signals. In [18], a triply primed array (TPA), based on fourth-order statistics, was proposed to increase the DOF beyond CPAs. The effects of mutual interference between sensors in an array were also discussed in [19][20][21].
In applications of beamforming [15,16] or cognitive radios [22], a large number of source signals are usually involved. Previous algorithms for joint estimation of DOA and CF [7][8][9][10] were applicable to a limited number of sources. In this work, a two-stage algorithm is proposed to estimate the DOAs of signal sources with known CF first, followed by the JE algorithm to estimate the DOAs and CFs of the remaining signal sources. Two orthogonal CPAs are used to increase the DOFs. This paper is organized as follows. The signal model is presented in Section 2, the proposed two-stage joint estimation algorithm is presented in Section 3, the simulation results are discussed in Section 4 and, finally, some conclusions are drawn in Section 5.

Sub-Nyquist Sampling Scheme
In this work, a sub-Nyquist sampling scheme [7,8] is used, where the sampling frequency is on the same order as the bandwidth of source signals, which is much lower than the carrier frequency. Figure 1 shows two orthogonal CPAs deployed along x and z axes, respectively, to receive signals from M sources, with the CF and DOA of the m-th source (s m (t)) being f m and θ m , respectively. Without loss of generality, assume that all the signals propagate in the xz plane and are uncorrelated to one another. Also assume the first M c sources have the same CF of f 0 . The bandwidth of each source signal is limited to B = 1/T, where T is the period of baseband signals. The other signals with unknown CF are assumed to fall in disjoint bands, namely, min  Figure 2 shows the configuration of a CPA (N 1 , N 2 ), which is composed of two ULAs at spacings of N 1 d and N 2 d, respectively; where N 1 and N 2 are coprime integers and d is the unit spacing. The first ULA has N 2 sensors at spacing of N 1 d and the second ULA has 2N 1 sensors at spacing of N 2 d. Both ULAs share the same first sensor, which is located at the origin, leading to a total of N = 2N 1 + N 2 − 1 physical sensors in the CPA (N 1 , N 2 ) configuration. The locations of the physical sensors in these two ULAs and the whole CPA are specified by Γ 1 , Γ 2 and Γ p , respectively, with where τ x (θ) = d cos θ/c, and χ n τ x (θ) is the time advance of the n-th sensor with respect to the first one at the origin. The Fourier transform of u n (t) is where S m ( f ) is the spectrum of the m-th baseband signal, with | f | ≤ B/2.

Consider a periodic signal
which is Fourier transformed to havẽ By applying an ideal low-pass filter (LPF) with cutoff frequency f c /2 toX n ( f ), we obtain which is slightly different from that in [8].
Next, take the inverse Fourier transform of X n ( f ) to have which is sampled at the rate of f s = 1/T s to obtain where x n [ ] = x n ( T s ) and w m [ ] =s m ( T s ). Note that we may choose f s = f c . Equation (1) can be put in a matrix form as wherex Similarly, the received signals in the array along the z-axis are processed to havē By taking noise into consideration, (2) and (3) are modified as whereĀ x andĀ z are the abbreviations ofĀ x (f ,θ) andĀ z (f ,θ), respectively;n x [ ] andn z [ ] are uncorrelated additive white Gaussian noise vectors at the th sampling time instant. The noise vectors are assumed to have zero mean and covariance matrix σ 2 nĪN , whereĪ N is an N × N identity matrix.

Second-Order Statistics
Assume the source signals are wide-sense quasi-stationary, which are observed over Q nonoverlapped time frames, with each time frame lasting for LT s [4,11,18]. The average power of the m-th source signal in the q-th time frame is represented as Assume P q m 's from different sources are uncorrelated, with their first and second moments determined as where E q means the expectation value is taken over the q-th time frame. The covariance matrices of the received signals in the q-th time frame are determined as wherex q [ ] is the received signal in the q-th time frame andD q = diag P q 1 , P q 2 , · · · , P q M . In practice,R q x andR q z are estimated as By concatenating all the columns ofR q x andR q z , respectively, into N 2 × 1 vectors, (5) is transformed toȳ difference manifold matrices, and means column-wise Kronecker product. By subtracting the expectation value over time-frame index q fromȳ x [q] andȳ z [q], respectively, we havē Next, by taking the average of entries inȳ x [q] andȳ z [q], respectively, that have the same lags (of same or opposite sign), we derive the observation vectors as where Φ contains the index of entries having lags in In other words,B Φ x andB Φ z serve as steering matrices of effective ULAs deployed along x and z axes, respectively, with virtual sensor locations specified by Γ s . The set Φ associated with the x-direction ULA is the same as that with the z-direction ULA because both ULAs have the same physical configuration.

Joint Estimation of DOA and CF
In this Section, a two-stage algorithm is proposed, in which the DOA of CF-known sources are estimated first, then the DOA and CF of the remaining sources are jointly estimated.

DOA Estimation of Signal Sources with Known CF
The second-order statistics of the received signals can be used to derive Q-manifold signals, which are then processed with subspace-based algorithms like ESPRIT or MUSIC to estimate the DOA and CF of the original sources. If the fourth-order statistics of the source signals are available, the fourth-order covariance matrices can be estimated as The covariance matricesR xx andR zz are then vectorized as By taking the average of entries with same lags inȳ xx andȳ zz , respectively, we obtain where Φ contains the index of entries having lags in [17]. The set Φ associated with the ULA in x-direction is the same as that with the ULA in z-direction because both ULAs have the same physical configuration. Note that subspace-based algorithms are applicable to CPA with contiguous lags in both x and z directions. Next, construct spatial smoothing covariance matrices of the fourth-order as [23] Algorithm 1 lists the major steps of the spatial smoothing MUSIC (SS-MUSIC). In steps 5-7, the conventional MUSIC algorithm is applied toR xx andR zz , respectively, to estimate the DOAs to the nearest grid points, and record them inφ 1x andφ 1z , respectively. Next, for every entry inφ 1z , find the best match inφ 1x with DOA difference smaller than 0.5 • , then record the average of these two DOAs inθ 1 .

Input:
• Q snapshots of manifold vectors,ȳ Φ x andȳ Φ z in (7) • Total number of signal sources, M • Known carrier frequency, f 0 Output: • Estimated DOA vector,θ 1 Algorithm: 1. Estimate fourth-order covariance matrices as 2. Vectorize, average and remove non-contiguous lags to formȳ Φ xx andȳ Φ zz , as in (8) and (9) 3. Construct spatial smoothing covariance matricesR xx andR zz , as in (10) 4. Apply singular-value decomposition (SVD) onR xx andR zz to havē Construct null spacesŪ xn andŪ zn out ofR xx andR zz , respectively, from the last N f − M columns ofŪ x andŪ z 7. Compute the DOA spectra Pick M highest peaks from each spectrum and record their associated angles inφ 1x and φ 1z , respectively. 8. For every entry inφ 1z , find a best match inφ 1x with DOA difference smaller than 0.5 • , then record the average of both DOAs inθ 1 .

Joint-ESPRIT in Projected Subspace
The DOAs of the signal sources with known CF ( f 0 ) are estimated and stored inθ 1 . The steering vectors associated with these sources are then removed by applying the orthogonal complement projectors where N s is the largest lag in Γ s of the second-order virtual arrays,B x andB z are the steering matrices of CPAs in x and z directions, respectively, with contiguous lags specified in Γ s . Next, apply the joint-ESPRIT (JE) algorithm to estimate the DOA and CF of the remaining sources embedded in the residual signals [8]. Construct two subarrays out of each second-order virtual array in x and z direction, respectively, as are the observation vectors from the first 2N s and the last 2N s second-order virtual sensors, respectively, in x direction, at the q-th sampling instant;B Φ x1 andB Φ x2 are the first 2N s rows and the last 2N s rows, respectively, ofB Φ are defined in the same way as their counterparts in x-direction. The matricesB x1 ,B x2 ,B z1 andB z2 satisfy rotational invariance property (RIP), namely, which is prerequisite to applying the JE algorithm. Define covariance matrices which can be implemented as The four matrices in (14) are concatenated to form Then, apply SVD onR c to haveR spans the same subspace asB c does, and M c is the number of signals with known CF, estimated in the first stage. By substituting (16) into (17), we havē Next, take the inverse ofŪ s1 and multiply it toŪ s2 ,Ū s3 andŪ s4 , respectively, to derivē By applying eigenvalue decomposition (EVD) on the average of these three matrices, we obtain whereΛ is a diagonal matrix containing all the eigenvalues ofΩ. The first two equations in (19) imply thatΨ Finally, the DOAs and CFs are estimated as where Ψ x,m and Ψ z,m are the (m − M c )-th diagonal entries ofΨ x andΨ z , respectively. The proposed projected Joint-ESPRIT (PJE) algorithm is listed in Algorithm 2.

Simulation Setup
Two CPA (5,6), each composed of 15 physical sensors, are deployed along x and z axes, respectively. The unit spacing is d = 3 cm, which is λ/2 at 5 GHz. A total of M = 30 sources are considered, including M u CF-unknown sources and M c CF-known sources. For convenience, we assume the first Each source signal is assumed to display quasi-stationary property [11]. The time-frame length of each source signal is set to L = 500. A total of 200 Monte Carlo realizations are simulated in each scenario. The signal-to-noise ratio (SNR) is defined as The CF-known sources emit at f 0 = 5 GHz. The CFs of the remaining sources are chosen from [0, f N − B], with f N = 6 GHz and B = 50 MHz, where the spacing between two adjacent channels is wider than B. To evaluate the performance of the proposed algorithm, define root-mean-square error (RMSE) as |f m − f m | 2 whereθ m and θ m are the estimated DOA and actual DOA, respectively, of source m,f m and f m are the estimated CF and actual CF, respectively, of the m-th CF-unknown source.

Cramer-Rao Bound
In [24,25], the Cramer-Rao bound (CRB) of DOA estimation was analyzed, assuming that there are fewer sources than antennas. In [26,27], the CRB of 1D DOA estimation in [24,25] was extended to 2D DOA estimation. In this work, the approach in [24,25] is extended to derive the CRB of joint DOA and frequency estimation. By concatenatingx[ ] andz[ ] in (4), a signal vector is formed as from which its covariance matrix is estimated as The parameters are listed in a vector as where θ m = θ m /180 and The CRB matrix is formulated as [24,25] where σ 2 n is the noise power defined in (4), LQ is the total number of snapshots, Re{} denotes the real part of its complex argument, The SNR defined in (21) can be rewritten as and the CRBs in DOA and CF estimations can be represented as

Maximum Detectable Number of Sources
In the first scenario, we compare the maximum number of detectable signal sources when there is only one source with unknown CF (1.15 GHz). Figure 3 shows the success probability of DOA and CF estimation under different numbers of signal sources, where the success probability is defined as the number of Monte Carlo realizations with RMSE < 5% in both DOA and CF estimation, divided by the total number of realizations. It is observed that the success probability with conventional JE [7,8] is acceptable when M ≤ 5, and drops dramatically when M > 5. On the other hand, the success probability with the proposed PJE remains about 90% even when M > 5. Figure 4 shows the RMSE of DOA and CF estimation, respectively, under different numbers of signal sources. The RMSE with conventional JE algorithm is higher than 5% when M > 2, while that with the proposed PJE algorithm remains about 1% even at M = 8.
In the second scenario, we try to find the maximum detectable number of CF-unknown signal sources, under a fixed number of sources. The success probability and RMSE with the proposed PJE algorithm are shown in Figures 5 and 6, respectively. The RMSE with the conventional JE algorithm is higher than 7%. The proposed two-stage algorithm works reasonably well when M u ≤ 11, and the RMSE increases dramatically when M u ≥ 13, indicating that the maximum number of CF-unknown sources is about 13 in this case.  Next, we try to find the maximum allowable number of signal sources, including 10 CF-unknown sources. Figure 7 shows that the success probability decreases rapidly when M > 36 and drops to 0 at M = 40. Figure 8 shows that the RMSE of DOA and CF estimation increases significantly when M is increased from 30 to 32, and stays at the same level before M reaches 39.  Figure 9 shows the normalized SS-MUSIC spectrum with Algorithm 1. It is observed that the blue peaks are aligned with the green peaks at the actual DOAs of CF-known sources in both cases of M u = 13 and M u = 14. All the CF-known sources and some of the CF-unknown sources are accurately estimated. Figure 10 shows the estimated DOA and CF, with M u = 13 and 14, respectively. In both cases, the DOA and CF of CF-known sources are accurately estimated. Figure 10a shows that all the CF-unknown sources are accurately estimated, while Figure 10b shows that some CF-unknown sources are poorly estimated, which is consistent with the observation in Figure 6 that the maximum number of detectable CF-unknown sources is 13.   Figure 11a,b show the estimated DOA and CF by using the proposed PJE algorithm and the conventional JE algorithm [8], respectively, with M u = 10 and M = 30. It is observed that the proposed PJE algorithm can detect all the sources accurately, while the RMSE of either DOA or CF with the conventional JE algorithm exceeds 20%. It is interesting to observe in Figure 11b that the error of CF-known sources is higher than that of CF-unknown sources. Since the received signals are contributed by several CF-known sources directionally nearby the CF-unknown sources, a high resolution of signal phase is required to resolve these sources. The joint diagonalization in the conventional JE algorithm probably degrade such phase resolution. Figure 11c shows the results by using the proposed PJE algorithm with M u = 13. All the sources are accurately estimated, with slightly larger error than in the case with M u = 10.

Detail of Proposed Two-Stage Algorithm
The SS-MUSIC algorithm in [17] has been extended to deal with multiple sources, some with known CF and others with unknown CF. By matching the results from x-CPA and z-CPA, the DOAs of CF-known sources are estimated first. Then, the CF-known signals are removed from the received signals by using (11)- (13) in the proposed PJE algorithm. Finally, the conventional JE algorithm is applied to jointly estimate the DOA and CF of the remaining CF-unknown signals.  Figure 12 shows the RMSE of DOA and CF estimation with the proposed PJE algorithm, under different SNR. The RMSE is well below 1% at SNR > 0 dB if M u = 10, and is below 1% at SNR > 5 dB if M u = 13.

Robustness of Proposed Two-Stage Algorithm
In the third scenario, we try to find the minimum angular spacing between two signal sources that can be detected by using the proposed two-stage algorithm. Consider the case of M = 3, with the incident angles of −60 • , 5 • and (5 + x) • , respectively, where x is in the range of [0.1, 1]. Figure 13 shows that the proposed PJE algorithm can detect two signals as close as 0.3 • apart, with success probability higher than 95%, as compared to 0.5 • if conventional JE algorithm is used. Figure 14 shows that the RMSEs of both DOA and CF estimation by using the proposed PJE algorithm are smaller than their counterparts with conventional JE algorithm when the angular spacing is smaller than 0.5 • .

Conclusions
A two-stage method is proposed to estimate the DOA of multiple signal sources by processing the received signals contributed by multiple CF-known sources and multiple CF-unknown sources. In the first stage, the intermediate results derived from the received signals of two orthogonal CPAs are matched to estimate the DOA of CF-known sources. In the second stage, the proposed PJE algorithm is applied to the residual signals, after removing the contribution of CF-known sources, to estimate the DOA and CF of the remaining sources. Simulation results show that by applying the proposed two-stage method to received signals from two orthogonal CPA(5, 6), both DOA and CF of 30 signal sources, 13 of which with unknown CF, can be accurately estimated.