Sequential Estimation of Frequency and Direction-of-Arrival Based on the Relaxed Coprime Array for Multiple Targets

Recently, extensive researches focus on the joint estimation of carrier frequency and direction-of-arrival, whereas the heavy computation burden and affordable hardware cost are always required therein. To alleviate the problem, we propose a relaxed coprime array based two-stage estimator, which can sequentially detect the frequencies and direction-of-arrival for multiple targets. Guided by the closed-form robust algorithm, this relaxation yields a higher array sparsity compared with existing sparse arrays. Specifically, the first stage aims to achieve the classification of the frequency and angle remainders arising from temporal-spatial undersampling, which is realized by combining spectrum correction with pattern clustering. By incorporating the closed-form robust Chinese Remainder Theorem into the remainder classification result, we obtain all the frequencies and direction-of-arrival in the second stage. Both performance analysis and numerical results were conducted to verify the effectiveness of the proposed estimator in the multi-target case.


I. INTRODUCTION
Estimation of the carrier frequency and direction-ofarrival (DOA) of the received signal is a fundamental problem of array signal processing, which is applied in numerous applications including radar, sonar, communication and so on [1]- [3]. Meanwhile, some applications develop in wider band and higher frequency, such as the wideband signal processing in cognitive radio (CR), milimeter wave in 5G wireless communication systems. Both the realizable sampling rate and the affordable power consumption of analogto-digital converters (ADCs) become bottlenecks due to the fact that at least two samples should be collected during a signal period according to the temporal-domain Nyquist sampling theorem [4]. On the other hand, as [5] pointed out, for the sake of reducing mutual coupling between elements, it is attractive to develop approaches to achieve DOA estimation when adjacent elements are spaced farther apart The associate editor coordinating the review of this manuscript and approving it for publication was Gerard-Andre Capolino.
(by N λ/2 or by M λ/2, where N , M ∈ Z + .). As a result, it is urgent to develop a low-complexity Frequency and DOA estimation method under the temporal-spatial undersampling condition.
It is well known that the most commonly-used technique is the Multiple Signal Classification (MUSIC) [6], which is based on subspace decomposition. Despite its effectiveness, a MUSIC-based estimator involves heavy computational burden when dealing with sparse arrays. Typically, it usually experiences the conversion from the observed physical covariance matrix to the virtual Nyquist-sense covariance matrix and the spatial smoothing operation [7]- [9] which inevitably yields heavy computational burden when searching for a global spatial spectrum. To lower the complexity, Zhou etc. divided the co-prime array into a pair of sparse uniform linear sub-arrays and then implemented the MUSIC algorithm on their covariance matrix respectively [10], which compromised the computation complexity and the phase ambiguity. However, this operation pays the cost of accuracy degradation. In addition, these MUSIC-based estimators only address the DOA estimation rather than the Frequency and DOA estimation.
Recently, several joint estimation algorithms along with array receiver structures have been proposed [11]- [13]. These proposal can greatly alleviate the sampling rate issue, since they adopt some sub-sampling methods. However, one common drawback of these joint estimators is that the costs of analog circuits in them are a bit high. Typically, for the two multi-coset sampling-based joint estimators (the JDFSD4MU estimator in [11] and the JDF4SA estimator in [12]), each sensor of a sparse array needs to be equipped with L(L > 2) branches of analog delayers. This arrangement not only inevitably consumes a mass of analog circuits, but also results in a difficulty in synchronization between time shift elements at quite high rates (as described in [13]). As to the CaSCADE joint estimator proposed in [13], the cost of analog circuits is also high, since it employs the modulated wideband converter (MWC), in which the impinging wideband signal at each sensor needs to be mixed by a periodic function, lowpass filtered and sampled at a low rate. Moreover, the method in [13] cannot be directly extended to the sparse arrays.
Unlike the multi-coset undersampling and the MWC, coprime undersampling is a direct means of discretization, which requires no analog circuits. Moreover, it also works well even when the channel number equals 2. Specifically, the parameterization in coprime undersampling is guided by Chinese Remainder Theorem (CRT) [14]- [16], which is a powerful approach to achieve the parameter estimation under the subsampling condition. Authors in [17]- [19] presented the CRT-based estimation methods to achieve DOA estimation with high accuracy. In [20], [21], some methods utilizing CRT were proposed to acquire the frequency estimates. The CRT-based estimators mentioned above all address a single parameter, although many signal processing applications benefit from the two combined.
To embrace the demand of practical applications, we propose a sequential estimator via combining the closedform robust CRT, spectrum correction and pattern clustering to realize the retrieval of multiple targets' frequencies and DOAs. The novelty is that we employ the combination of spectrum correction and pattern clustering as an effective solution to the difficulty in classifying the CRT remainders acquired from the DFT spectrum of the collected samples.
The three main contributions are summarized as follows: 1) Guided by the closed-form robust CRT [22], we propose a relaxed coprime array which allows a much higher array sparsity than the original one, thereby reducing mutual coupling between sensors. 2) Guided by the closed-form robust CRT, we propose a relaxed coprime temporal underampling scheme which requires only 2 low-rate analog-to-digital converters (ADCs) at each sensor. Moreover, no analog delayed units are included in this scheme, and thus this sampling method can be easily implemented from the hardware point of view.
3) Combining the closed-form robust CRT, spectrum correction and pattern clustering, we propose a sequential estimator which realizes the retrieval of multiple targets' frequencies and DOAs at a low computation complexity. Further, this combination yields a merit that the Degree-of-Freedom (DOF) can be enhanced via increasing the snapshot size rather than enlarging the array aperture, which differs greatly with the existing joint estimators [11]- [13]. The rest of this paper is organized as follows. In Section II-A, we propose the relaxed coprime array and formulate the signal model; In Section III, we elaborate on the CRT reconstruction and remainder acquisition for the frequency estimation. In Section IV, the DOA estimation principle and realization are considered based on the knowledge obtained from the process of frequency estimation. The procedures of the proposed estimator are detailed in Section V. In Section VI, we present the performance analysis on the proposed multi-target sequential estimator. In SectionVII, we show the numerical experiments before the paper is concluded in Section VIII.

II. SYSTEM MODEL FOR MULTIPLE TARGETS BASED ON THE RELAXED COPRIME ARRAY
A. RELAXED COPRIME SPARSE ARRAY Let f max be the upper frequency of the incident signal and thus the unit wavelength λ = c/f max (c refers to the light velocity). As Fig. 1 illustrates, the proposed array is a nonuniform linear array with L sensors (each sensor is equipped with 2 undersampling ADCs). Guided by the closed-form CRT (see [22] for details), an individual inter-element spacing d l (i.e. the distance between the l-th sensor and the (l + 1)-th sensor) is defined as where 1 , . . . , L−1 are pairwisely coprime integers. We infer that any selection of 1 , . . . , L−1 surely leads the proposed array to a high sparsity, i.e., d l λ/2, l = 1, . . . , L − 1. Further, due to the cascade product of k in (1), the greater the sensor number L is, the sparser this array tends to be. Thus, the proposed array is named relaxed coprime sparse array.
Accordingly, the distance between the l-th senor and the 1 st sensor is B. MULTI-TARGET DATA MODEL BASED ON THE RELAXED COPRIME ARRAY Assume that there are D uncorrelated far-field narrowband targets {s i (t)}(with the incident direction θ i and the carrier frequency f i ), i = 1, . . . , D, impinging on our proposed relaxed coprime sparse array in Fig. 1. Assume the i-th signal FIGURE 1. Relaxed coprime sparse array configuration including L sensors, with distance d l between the l -th sensor and the (l + 1)-th sensor. Each sensor includes 2 ADCs with sampling rates f s1 f s2 at each sensor.
where a i , φ i are unknown amplitude and initial phase, respectively. For the direction of arrival θ i , the observed signal at the l-th sensor can be expressed as where ζ l (t) is the zero-mean Gaussian noise. We address the recovery of frequencies and DOAs for multiple sources based on the received signals.

A. CRT RECONSTRUCTION MODEL FOR THE FREQUENCY ESTIMATION
Note that, limited by the realizable sampling rate and the affordable power consumption, at any sensor, f i can hardly be estimated by employing a common algorithm to process discretized data collected from a single high-rate ADC. Instead, as Fig. 1 illustrates, guided by the closed-form robust CRT, we substitute the single high-rate ADC with two low-rate ADCs at each sensor and design the relaxed coprime sampling method for data discretization. Specifically, the two subsampling rates of ADC1, ADC2 in Fig. are set as where η 1 , η 2 are coprime integers and M f is a positive integer. That's the point we refer to it as the relaxed coprime undersampling. Substituting t = n/f s1 , t = n/f s2 into (4), we can acquire two M f -length discretized sequences x l,1 (n), x l,2 (n), where ϕ l i refers to the i-th target's impinging phase at the l-th sensor, i.e., Due to undersampling, for any individual target, its two digital frequencies f i /f s1 , f i /f s2 in (6) surely exceed 1, i.e., where i = 1, . . . , D. Thus, x l,1 (n), x l,2 (n) can be rewritten as A i e j(2πf i,1 n+ϕ l i ) + ζ 1 (n) A i e j(2πf i,2 n+ϕ l i ) + ζ 2 (n), For convenience, denotẽ Substituting (5), (10) into (8) One can find that, Eq.(11) is fully in accordance with the model of closed-form robust CRT [22] with 2 reconstruction paths, in which 2M2 } refer to the remainders and n i,1 , n i,2 are 2 folding integers to be determined.
For an individual l-th sensor, as long as two normalized frequency estimatesf l i,1 ,f l i,2 are acquired, n i,1 , n i,2 can be determined using the closed-form robust CRT (see Appendix B), which calculates the frequency estimate aŝ Finally, averaging all L sensors' frequency estimates in (12) yields a more accurate frequency estimatef î Accordingly, the wavelength estimateλ i = c/f i can also be readily calculated. Furthermore, since that the frequency of signals can be shown as the peak of the DFT spectra, and the normalized fre- (14) refer to the peak indices of their M f -point hanning-windowed DFT spectra X l,1 (k), X l,2 (k), andδ l i,1 ,δ l i,2 are their fractional frequency offsets. At this point, it can offer an effective way to recover the parameters by combining the closed-form robust CRT and DFT analysis. As a result, the extremely urgent task is to achieve the classification of the frequency remainders. Towads this purpose, we introduce the spectrum correction and pattern clustering as follows.

B. PATTERN CLUSTERING BASED REMAINDER ACQUISITION
From the aforementioned analysis, there are several problems to be solved. 1) it is notable that theδ l i,1 ,δ l i,2 in (14) cannot be directly acquired from the DFT results influenced by spectral leakage. 2) We suppose there are D sources, consequently, there exist D spectral peaks across either X l,1 (k) or X l,2 (k), k = 0, . . . , M f − 1. In other words, D contemporary congruence equation sets (similar to (11)) can be built up. On this occasion, identifying the spectral peak index pair {k l i,1 , k l i,2 } among 2D spectral peaks in X l,1 (k), X l,2 (k) is crucial. 3) Since a received signal contains multiple tones, the interferences among them further add to the difficulties of the parameter estimation. To overcome the difficulties, we propose to combine the spectrum correction and pattern clustering.

1) SPECTRUM CORRECTION
As is known, for the signal with frequency f i , influenced by the spectral leakage, the DFT spectrum of it exhibits multiple spectral bins, which are scattered around the peak bin. The observed peak index bin differs from the ideal one, thus the parameter pair {A i , ϕ l i } in (9) cannot be accurately taken from the peak bin. To enhance the remainder acquisition accuracy, it is vital to adopt some correction measures to estimate the ideal parameter triple (i.e.{f i , A i , ϕ l i }). Firstly, we implement the windowed DFT rather than the direct DFT on each sequence pair {x l,1 (n), x l,2 (n)}. Then, we employ the ratio-interpolation based spectrum corrector [23], by which an exponential sequence's 3 parameters (the frequency offset, amplitude and phase) can be accurately corrected from its hanning windowed DFT results. Notablely, the operation only deals with two spectral bins, i.e. the peak bin and the adjacent sub-peak bin, indicating the computation complexity is significantly low (See Appendix A for detail). Implementing this procedure on x l,1 (n) and x l,2 (n), one can acquire

2) PEAK PAIR IDENTIFICATION
Since two normalized frequenciesf i,1 ,f i,2 stem from the same target with frequency f i , the exponential tone A i e j(2πf i,1 n+ϕ l i ) in x l,1 (n) and the exponential tone A i e j(2πf i,2 n+ϕ l i ) in x l,2 (n) share a common amplitude-phase parameter pair {A i , ϕ l i } (see (9)), which provides a reference for identifying the peak index pair {k l i,1 , k l i,2 }. Moreover, for an individual tone with the frequency f i in (8), it is certain that the two normalized frequenciesf i,1 =f i,2 since the two undersampling rates f s1 = f s2 . Consequently, for any i-th individual tone, it is very likely that the spectral peak indices k l i,1 , k l i,2 of X l,1 (k) and X l,2 (k) are distinct. It is crucial to identify the peak index pair Further, a parameter pair {Â,φ} recovered by the spectrum corrector can be combined into a vector quantityẑ =Âe jφ , which helps to better identify peak pairs. Particularly, one can recover a vector quantityẑ l i,1 =Â l i,1 e jφ l i,1 at an individual peak index k l i,1 of DFT spectrum X l,1 (k) of the ADC1, whereas the other vector quantityẑ l p,2 =Â l p,2 e jφ l p,2 can also be constructed at any peak index k l p,2 of the DFT spectrum X l,2 (k) of the ADC2. Hence, if the distance d l i,p = ||ẑ l i,1 −ẑ l p,2 || (|| · || refers to the Frobenius norm) is sufficiently small, it is reasonable to consider thatẑ l i,1 ,ẑ l p,2 stem from a common target, vice versa. Accordingly, in terms of (9), other recovered parameterŝ Once all the peak pairs of L sensors are collected, all corrected parameters of any individual target can be clustered together. In this way, the remainder acquisition for the CRT reconstruction can be achieved. Substituting the acquired parameters into the CRT algorithm (See Appendix B for detail), the D frequency estimationf i , i = 1, . . . , D can be acquired.

IV. DOA ESTIMATION BASED ON THE CRT AND DFT ANALYSIS A. ESTIMATION PRINCIPLE
In terms of (7), as for the i-th source, its theoretical phase difference between the (l + 1)-th sensor and the l-th sensor is denoted as Since d l λ/2, it is known that, increasing the interelement spacing in arrays could lead to array ambiguities. In other VOLUME 8, 2020 words, φ l i surely contains 2π -ambiguity, i.e., where n l is an unknown integer. Combining (15) with (16) yields Let an intermediate quantityd be formulated as where M θ is a specified positive integer. Combining (1) with (18), we have Then, substituting (19) into (17) yields For convenience, denote Hence, Eq. (20) can be rewritten as a system of simultaneous congruences One can find that, Eq. (22) is fully in accordance with the model of robust Chinese Remainder Theorem [22], in which M 1 , . . . , M L−1 refer to the moduli, r 1 i , . . . , r L−1 i refer to the remainders and n 1 i , . . . , n L−1 i refer to the folding integers to be determined.
As [22] pointed out, for the robust CRT, the upper bound of reconstruction is the least common multiple (lcm) of all the moduli. Specifically, for the model in (22), Q θ i should satisfy Since f max ≥ f i and thus λ ≤ λ i , combining (18) with (21), one can prove that (23) constantly holds.
As to the i-th target, following the procedure of the closed-form robust CRT (see Appendix B), one can acquire L − 1 folding integers n 1 i , . . . , n L−1 i and thus Q θ i can be estimated asQ Accordingly, the final DOA is estimated aŝ

B. ESTIMATION REALIZATION
What matters in the DOA estimation process is similar to that in the frequency estimation, obviously, the preceding steps in Section III can give a clear indication. According to the peak identification results, L vector quantity . . , L can be obtained. To enhance the accuracy of φ l i in (24), at the l-th sensor, it is necessary to average two vector quantitiesÂ l i,1 e jφ l i,1 ,Â l p,2 e jφ l p,2 resulting from spectrum correction, i.e., Similarly, one can also acquire the corrected phase estimatê ϕ l+1 of the arriving phase at the (l + 1)-th sensor. Note thatφ l i ,φ l+1 i are two observation phases falling in the range [−π, π). The phase-difference estimate φ l i ∈ [0, 2π ), which is required by DOA estimation (24), can be calcluated through the following modulo operation From the above analysis, one can see that the proposed relaxed coprime array's superiority over the existing sparse arrays lies in that the phase ambiguities can be removed by means of CRT and spectrum correction. On the other hand, the DOA estimation takes advantage of the accurate phase and amplitude knowledge acquired in frequency estimation, that is why we call it the sequential estimator.

V. THE PROCEDURE OF THE FREQUENCY AND DOA ESTIMATION FOR MULTIPLE TARGETS
To facilitate comprehension, the dataflow of correction operations on the sequence pair x l,1 (n), x l,2 (n) is summarized in Fig. 2(a). Further, the overall dataflow of the proposed Frequency and DOA sequential estimator is given in Fig. 2(b). To summarize, we present the procedures of the sequential Frequency and DOA estimator, which consists of 5 steps and can be divided into 2 stages: This stage mainly addresses the peak pair identification (i.e., to acquire D peak indicators c 1 , . . . , c D ), which achieves the remainder acquisition to reconstruct the CRT model.

VI. PERFORMANCE ANALYSIS OF THE MULTITARGET SEQUENTIAL ESTIMATOR A. ANALYSIS OF ANTINOISE ROBUSTNESS AND ACCURACY
Proposition 1: The successful recovery of the frequencies and DOAs for multiple sources based on the proposed estimator requires that the following two inequalities hold.
wheref i,1 ,f i,2 are the i-th target's ideal normalized frequencies and ϕ l i is the i-th target's ideal phase difference between the l-th sensor and (l + 1)-th sensor.
Proof: As [22] pointed out, the successful recovery of the robust CRT requires that each remainder error be less than a quarter of the greatest common divisor (gcd) of all moduli.
For the frequency estimation, the moduli are and thus M f = gcd{M 1 ,M 2 }, thus the remainder errors of the i-th target at the l-th sensor should satisfy Hence (30) holds. From Proposition 1, four conclusions can be drawn: 1) With the sample length M f fixed, the greater the minimum of η 1 , η 2 is (i.e., higher sampling rates give rise to shorter observation time intervals), the smaller the error tolerance 1 4η 1 (or 1 4η 2 ) of frequency detection is, thereby lowering the robustness of the frequency recovery. For the sake that these two error tolerances are independent of the snapshot length M f , we specify the undersampling rates of two ADCs as M f η 1 , M f η 2 .
2) With the sample length M f fixed, the greater the minimum of 1 , . . . , L−1 is (i.e., the array layout becomes sparser), the lower the robustness of the DOA recovery is, which can be viewed as the cost paid by our proposed array's relaxation on interelement spacings. 3) Increasing the sample length M f , which helps to improve the accuracy of spectrum correction and thus reduce all remainder errors in (29), (30), will enhance the antinoise robustness of both frequency recovery and DOA recovery. 4) Within the proposed estimation framework, the key is to detect the frequencies of the undersampled sequence at each ADC. In other words, the question whether the proposed method can succeed depends heavily on the accuracy of the frequency estimation. In fact, DFT-based spectrum correction is not the unique means of frequency detection. In recent years, several convex optimization methods [24] exploiting signal sparsity were proposed to achieve high-accuracy frequency estimation. However, these methods are not incorporated into our proposed estimator. The reasons lie in two aspects: • Other than frequency estimates, the DFT-based spectrum correction also provides phase estimates (See Appendix A), which can be further reused to estimate DOAs. In contrast, the convex optimization methods can hardly recover the phase information, which makes data-reuse based DOA estimation impossible.
• The DFT-based spectrum corrector does not involve any optimization process and thus consumes a lower computational complexity.
The 3 rd and 4 th conclusions are obvious. The former two conclusions will be verified by numerical simulations in Section VII.

B. ANALYSIS OF THE DEGREE OF FREEDOM
The degree of freedom (DOF) refers to the maximum number of targets that an estimator can resolve [25]. To address this problem, as far as the identification mechanism of multiple targets is concerned, it is necessary to make a comparison between the existing sparse arrays and the proposed relaxed array.
For the existing sparse arrays [26]- [28], their DOA estimators usually experience the following steps [25]- [27], [29], [30]: original covariance matrix calculation (i.e., in the physical undersampling sense), continual lag extraction from the difference set, covariance matrix reconstruction (i.e., in the virtual Nyquist-sampling sense), spatial smoothing and MUSIC decomposition. Hence, the DOF equals the dimension of signal subspace in the final MUSIC decomposition, which is closely related to the number of sensors (e.g, the original coprime array's DOF is MN , whereas its sensor number is N +2M −1 [27]). In other words, the enhancement of DOF can only be realized through spatial-domain ways rather than temporal-domain ways.
For our proposed sequential estimator, its multi-targetidentification mechanism is completely different with those existing sparse arrays-based estimators. From Step 3 of the procedure in Section V one can find that the DOF equals the number of obtainable peak indicators c i , i = 1, . . . , D. Since an indicator c i represents the pairing relationship between a DFT bin cluster at ADC1 and another DFT bin cluster at ADC2, the number of DFT clusters actually equals the number of indicators (i.e., the DOF). It is well known that, for a complex exponential tone, the windowed DFT operation plays the role of suppressing spectral leakage, resulting in that a DFT bin cluster covers less than 5 bins (as will be illustrated in Fig. 3). Further, recall that each ADC collects M f samples and the DFT size is also M f , from which it can be inferred that the DOF of the proposed sequential estimator is M f /5.
In other words, the DOF of the proposed estimator is enhanced through increasing the sample length rather than improving the spatial-domain structure.

C. ANALYSIS OF COMPUTATIONAL COMPLEXITY
From Fig. 2 and the retrieval procedure listed in Section V, one can notice that the phase adjustment, the averaging operation, the reverse-sine operation are no more than simple algebra calculations. Besides, the pattern clustering involves only a small number of vector quantities. Moreover, the clustering results at any individual sensor also apply to other sensors. Hence, the computational complexity consumed by pattern clustering is ignorable. Moreover, the CRT we adopt works in a closed-form way, which involves no tedious searching operations.
Thus, the proposed estimator's computational complexity mainly depends on the DFT-based spectrum correction. As the Appendix A illustrates, the calculation involved in spectrum correction (33)-(38) are no more than simple algebra operations, also. As a result, the computational complexity mainly focuses on 2L times of M f -point windowed DFT. Obviously, no optimization operations and searching operations are involved, thus the proposed estimator's computational complexity is much lower than the existing estimators.

D. MERITS AND SHORTCOMINGS OF THE PROPOSED ESTIMATOR
Based on the above performance analysis, we summarize the merits and shortcomings of the proposed estimator.
The specific merits include: 1) We present the realxed coprime array configuration, which allows a much higher sparsity than the existing sparse arrays, thereby considerably reducing the mutual coupling. The closed-form expression of the relaxed coprime array is specified by Eq. (1), which makes it easy to apply optimization in practical engineering. 2) In the coprime subsampling scheme, only two lowrate ADCs are required at each sensor, achieving the hardware-cost saving. Compared to other undersampling methods for the joint estimation (e.g. the multicoset and the interleaved approach), which are limited by the requirement of accurate time-delay synchronization, it is more applicable in practical applications. 3) Different from the conventional coprime array-based estimators, whose DOF are subject to the number of sensors, the proposed estimator can achieve a higher DOF by increasing the sampling length in temporal domain. The characteristic can be extended to reduce the sensor-number cost. Typically, the proposed estimator applies to the case that the sensor number L = 3, where the CRT remainders are provided by the phase difference and thus the number of reconstruction paths equals 2. In this case, multi-target estimation can also be realized as long as the sample length M f is long enough.

4) The proposed estimator can obtain the frequencies and
DOAs in the multi-target case, and the computation complexity mainly depends on the 2L times of M f -point windowed DFT.
The shortcomings lie in two aspects: • As aforementioned, the high array sparsity is at the cost of lowering antinoise robustness of the DOA recovery.
• Our proposed sequential estimator cannot discriminate coherent targets. The reasons are as follows. As mentioned in Section III-B.2, the peak pair identification (involved in Stage 1) is the basis for both frequency estimation and DOA estimation (involved in Stage 2). If two targets share the same frequency, their DFT bin clusters inevitably overlap, resulting in a failure of peak pair identification.

VII. SIMULATION RESULTS
In this section, the simulation is carried on the platform of MATLAB R2016b, with an Inter Core i5 2.60 GHz. Firstly, we will demonstrate the proposed sequential frequency and DOA estimator step by step. To disclose our estimator's DOA distinguishing ability, i.e. the ability to estimate the closely spaced signals, both the case that all targets have distinct DOAs and the case that several targets share the same DOA are taken into account. Secondly, we will investigate the ADC sampling rates' influence on the antinoise robustness of frequency estimation; Thirdly, we will investigate the array sparsity's influence on the antinoise robustness of DOA estimation, also. VOLUME 8, 2020

A. STEPWISE DEMO OF THE MULTITARGET SEQUENTIAL ESTIMATOR
The input parameters of the procedure in Section V are specified as follows: the upper frequency f max = 2GHz (thus the unit wavelength λ = c/f max = 0.15m), the coprime integers The signal parameters (i.e., frequency f i , amplitude A i , initial phase φ i , incident direction θ i ) in the data model (4) are specified in Table 1. The signal-to-noise ratio (SNR) at each sensor (with additive Gaussian white noise) is set as 10 dB. Substituting M f , η 1 , η 2 into (5), one can calculate two ADCs' sampling rates D), indicating that the temporal-domain sampling rates are sufficiently low. Substituting λ, 1 , 2 into (1), one can calculate 2 interelement spacings {d 1 , d 2 } = {0.9m, 0.75m} λ/2, indicating that the array sparsity is sufficiently high, also.
The following detailedly demonstrates the procedure in Section V step by step.
From Table 4, one can find that the frequency estimation errors are less than 20Hz and the DOA estimation errors are less than 0.81 • .

B. DOA ESTIMATION FOR MULTIPLE TARGETS WITH THE SAME DOA
To assess the estimation performance under the condition that several targets share the same DOA, which we can refer to as the applicability to the closely spaced signals, we still employ the model parameterization in Section VII-A except that all the DOAs θ 1 , . . . θ 4 in Table 1 are altered as 80 • . Fig. 5 illustrates the clusters of vector quantities and Table 5 lists the estimation results of frequencies and DOAs.
From Fig. 5 and Table 5, it is discovered that the proposed estimator still applies to the same DOA case, indicating that our estimator suits well for the case of the most closely spaced signals. This attributes to the fact that the spectral overlappings across multiple targets only depend on the targets' frequencies and the ADC sampling rates, irrelevant of the target DOAs. As long as spectral overapplings are bypassed, our proposed estimator works well.

C. ANTINOISE ROBUSTNESS VERIFICATION FOR FREQUENCY ESTIMATION
In order to investigate sampling rates' effect on the antinoise robustness of frequency estimation, in this simulation, we compare 2 groups of coprime integers η 1 , η 2 ([1801, 1811] and [3301, 3307]). In terms of (5), the sampling rates f s1 , f s2 equal [0.988672, 0.997888]MHz and [1.690112, 1.693184]MHz. Other parameter settings are the same as Section VII-A.
For each SNR case, 1000 Monte Carlo trials were conducted. We use the detection probability P d to evaluate the antinoise robustness. Its criterion is as follows: for any trial, if |f i − f i | < f i /1000, then it passes; otherwise it fails. Fig. 6 (a)illustrates the P d curves for these two coprime integer groups (take the case target 1 for example), showing that the P d curve for η 1 = 3301, η 2 = 3307 is 4dB right to the P d curve for η 1 = 1801, η 2 = 1811. Thus, the former antinoise robustness is inferior to the latter. This is because, with the fixed sample length M f , a higher sampling rate results in a shorter observation time interval and thus deteriorates the antinoise robustness, thereby confirming the 1 st conclusion drawn in Section VI-A.

D. ARRAY SPARSITY's EFFECT ON DOA ANTINOISE ROBUSTNESS
In order to investigate array sparsity's effect on the antinoise robustness of DOA estimation, in this simulation, we compare 2 groups of coprime integers 1 , 2 ( [5,6] and [13,17]). In terms of (1), the interelement spacings d 1    in Table 1 into λ i = c/f i , i = 1, . . . , D). Thus, the array sparsity is sufficiently high. Other parameters are the same as Section VII-A.
For each SNR case, 1000 Monte Carlo trials were conducted. We use the detection probability P d to evaluate the antinoise robustness. Its criterion is as follows: for any trial, if |θ i − θ i | < 1 • , then it passes; otherwise it fails. Fig. 6(b) illustrates the P d curves for these two coprime integer groups (take the DOA recovery of target 1 for example), which shows that the P d curve for 1 = 13, 2 = 17 is 3 dB right to the P d curve for 1 = 5, 2 = 6 (i.e., the antinoise robustness of the former case is inferior to the latter). This confirms the conclusion that the increase in the array sparsity deteriorates the antinoise robustness, as drawn in Section VI-A.

E. THE EFFECT OF THE SOURCE SIGNALS' ANGLE ON THE DOA ANTINOISE ROBUSTNESS
In the previous simulation, we compare the effect of the array sparsity and sampling rates on the estimation performance. To further explore the performance of the proposed estimator, in this subsection, we compare the DOA antinoise robustness under different angles with the identical estimation process. Specifically, we compare 3 angles θ 1 = 30 • , θ 2 = 50 • , θ 3 = 80 • based on the model parameterization in Section VII-B, i.e. all targets have the same angle θ m (m=1, 2, 3) in one test, then evaluate the antinoise robustness in the case of For each SNR case, 1000 Monte Carlo trials were conducted. We use the detection probability P d to evaluate the antinoise robustness. Its criterion is as follows: for any trial, if |θ i − θ i | < 1 • , then it passes; otherwise it fails. Fig. 6 (c)illustrates the P d curves for these three different angles (take the DOA recovery of Target 1 as an example), which shows that the estimation performance for θ = 30 • is the best, and that for θ = 50 • , θ = 80 • takes the 2nd, 3rd place, respectively. The outcome implies that the proposed estimator is more suitable for the signals with smaller angles despite the identical estimation condition. We emphasize that the proposed estimator works in a sequential way. Concretely, we first combine the DFT analysis and CRT reconstruction to acquire the frequency estimates. Then we utilize the phase values and the clustering result arising from the frequency estimation to achieve the CRT reconstruction for  DOA estimation. From (24), it is apparent that the solution to the CRT reconstruction is not the ultimate DOA estimate, as opposed to the frequency estimation. Instead, the DOA estimate is obtained via (25), which involves arcsin(·) operation. Notably, the derivative of the inverse sine function arcsin(·) is progressively increasing, indicating the same error value in (24) can result in different error values in (25), and the estimation performance difference accordingly.

F. COMPARISONS WITH ESTIMATORS BASED ON THE ORIGINAL COPRIME SPARSE ARRAY AND THE NYQUIST ULA
Compared with the original coprime array based joint estimators [30], [32], in which the parameter pairs{f i , θ i } are separable and the estimation can be achieved in a multidimensional searching way, the data model in our proposed estimator is inseparable and the estimation is acquired in a sequential way.
In particular, we aims to recovery the carrier frequency whereas the frequencies detected by some of the existing joint estimators refer to the carrier frequency offset(or Dopplereffect related frequencies). Limited by the data mode and parameter properties, the comparison for frequency and DOA can be hardly be conducted.
Instead, we compare the DOA detection performances between our proposed estimator and the original coprime array based DOA estimator [17], where the MUSIC algorithm is utilized. Besides, as a reference, the traditional Nyquist ULA based estimator was also considered.
All of these 3 estimators were specified with the same target number D = 4, the same sample length M f = 512 and the same sensor number L = 6 (i.e. the minimum sensor number of the original coprime sparse array The D = 4 targets were specified with the same signal parameters listed in Table 1. For our proposed sequential estimator, except the sensor number L and the integers 1 , . . . , 5 , other parameters were the same as those in Section VII-A. Fig. 7 illustrates the root-mean-square error (RMSE) curves of DOAs detected by these 3 estimators (For each SNR case, 1000 Monte Carlo trials were conducted), from which 2 conclusions can be drawn: 1) In the high-SNR region, for any target, the proposed sequential estimator exhibits much higher accuracy than the other two estimators. Overall speaking, the RMSE curve of the proposed estimator is about 2 ∼ 3 orders of magnitude lower than the other two RMSE curves. This improvement attributes to 3 aspects: Firstly, the hanningwindowed operation plays the role of suppressing the interferences among multiple targets, which cannot be realized by the other two subspace decomposition based estimators; Secondly, the ratio spectrum corrector can accurately provide 3-parameter estimates, ensuring that accurate remainders are fed to the subsequent CRT recovery; Thirdly, as [22] pointed out, the CRT reconstruction tool itself does not introduce any estimation errors. 2) In the low-SNR region, the proposed sequential estimator's robustness is inferior to the other two estimators. Specifically, for each RMSE curve of any estimator, there exists an SNR threshold below which the DOA estimation performance acutely deteriorates. Overall speaking, the traditional Nyquist ULA based estimator has the lowest SNR threshold, and the original coprime array based estimator, the proposed sequential estimator take the 2nd, 3rd place, respectively. For our estimator, this deficiency can be regarded as the cost paid by two aspects of performance improvements: On one hand, the degradation of antinoise robustness in the low-SNR region partly arises from the accuracy improvement in the high-SNR region; On the other hand, as theoretically analyzed in Section VI-A, our proposed array's relaxation on interelement spacings sacrifices the estimator's antinoise robustness. Also, for our proposed estimator, we can discover that Target 1 has the highest SNR threshold, indicating that it has the worst antinoise robustness. From Table 1, Target 1 has the smallest amplitude value, accordingly, the peak bins of Target 1 are the most susceptible to noise in DFT spectrum, and consequently, the robustness of DOA estimation is the worst.
Further, as analyzed in Section II-A, the greater the sensor number L is, the sparser the relaxed coprime array tends to be, indicating that the antinoise robustness becomes worse, also. Therefore, it is not suggested to deploy too many sensors in the relaxed coprime array. Practically, L = 3 is sufficiently enough, since the DOF of the proposed estimator is enhanced through increasing the sample length M f rather than increasing the sensor number L.

VIII. CONCLUSION
This paper proposes a relaxed coprime array-based multitarget sequential Frequency and DOA estimator. The proposed estimator applies to the multitarget case This array has two characteristics: 1) It allows a much higher array sparsity than the existing difference based coarrays; 2) It also works well even if only 3 sensors are included in this array, which is beyond the reach of the existing sparse arrays. Combining this array with 3 techniques (spectrum correction, closedform robust CRT and patten clustering), we devised a novel sequential estimator which possesses 4 merits: 1) Suitable for the multi-target case; 2) Suitable for highly sparse sensor arrangement; 3) Low computational complexity; 4) DOF is enhanced by increasing the sample length.
Since our proposed sequential estimator possesses the above 4 merits under the temporal-spatial undersampling condition, it suits well with the increasing tendency of spectrum utilization in higher frequency bands. Hence, our proposed estimator has a promising prospect in wide application fields such as radar, sonar, communication, etc.
where [·] represents the rounding operation such that Step 2: Calculate the remainder ofq i,1¯ i,1 modulo i : where¯ i,1 is the modular multiplicative inverse of 1 modulo i and can be previously calculated.
Step 3: Calculate an unknown folding integern 1 : where b i,1 is the modular multiplicative inverse of γ 1 / i modulo i , which can be calculated in advance also, and γ i is defined as Step 4: Calculate other J − 1 folding integersn î n i =n Step 5: Recover the final real-valued numberQ aŝ ANAN LIU (Member, IEEE) was born in 1979. He received the Ph.D. degree from Tianjin University, Tianjin, China, in 2010. He is currently a Professor with the School of Electrical and Information Engineering, Tianjin University. His research interests include artificial intelligence and biomedical signal processing.