Direction-of-Arrival Estimation for Circulating Space-Time Coding Arrays: From Beamspace MUSIC to Spatial Smoothing in the Transform Domain

As a special type of coherent collocated Multiple-Input Multiple-Output (MIMO) radar, a circulating space-time coding array (CSTCA) transmits an identical waveform with a tiny time shift. It provides a simple way to achieve a full angular coverage with a stable gain and a low sidelobe level (SLL) in the range domain. In this paper, we address the problem of direction-of-arrival (DOA) estimation in CSTCA. Firstly, we design a novel two-dimensional space-time matched filter on receiver. It jointly performs equivalent transmit beamforming in the angle domain and waveform matching in the fast time domain. Multi-beams can be formed to acquire controllable transmit freedoms. Then, we propose a beamspace multiple signal classification (MUSIC) algorithm applicable in case of small training samples. Next, since targets at the same range cell show characteristics of coherence, we devise a transformation matrix to restore the rotational invariance property (RIP) of the receive array. Afterwards, we perform spatial smoothing in element domain based on the transformation. In addition, the closed-form expression of Cramer-Rao lower bound (CRLB) for angle estimation is derived. Theoretical performance analysis and numerical simulations are presented to demonstrate the effectiveness of proposed approaches.

Space-time coding is based on the colored space-time waveform transmission principle. It allows transmitting different waveforms in different directions with a wide angular coverage [11][12][13][14][15]. The performance of four typical space-time coding waveforms for active antenna systems was assessed by ambiguity functions in [3]. The impact of mutual coupling on MIMO radar with space-time coding was investigated in [4], and a calibration procedure for the transmit beamforming was presented. Particularly, the concept of circulating space-time coding was proposed in [1]. As a special type of coherent colocated MIMO radar [16], CSTCA is simple to be implemented in engineering [17,18]. In [17], a low-cost digital beamforming technique was introduced according to circulating signal principle. In [18], the element embedded pattern error and antenna coupling were directly calibrated in operating mode by exploiting circulating space-time coding. In [19], an extended circulating space-time code with constant modulus was suggested for transmit beampattern synthesis. The joint slow-time coding with circulating codes was derived to restore the range resolution and suppress interferences in [20].
Direction-of-arrival (DOA) estimation of multiple targets is one of the most important radar applications [21][22][23]. Methods in the context of colocated MIMO radar have been developed extensively. The most popular approaches are multiple signal classification (MUSIC) and estimation of signal parameters via rotational invariance techniques (ESPRIT) with high resolution capabilities [24][25][26][27]. However, MIMO radar demands to transmit orthogonal waveforms ensuring independence between sensor channels. The completely orthogonal waveforms are hardly available for multi-channel transmissions resulting in gain fluctuations in different directions. Besides, the correlations between distinct waveforms usually lead to a high sidelobe level in range domain, which has an adverse impact on weak target detection. Additionally, the constraints of constant modulus cannot be satisfied without performance loss in such cases [28,29].
The transmit beamspace design for colocated MIMO radar has been widely explored, which is generally formulated as an optimization problem with a tedious solving process [30]. Compared with current beamspace direction finding techniques [30][31][32], the proposed MUSIC algorithm can circumvent requirements on beamspace design. Besides, DOAs can be estimated properly in case of small samples with a low sidelobe level (SLL) in range domain. Therefore, CSTCA scheme is an advisable choice in parameter estimation, which has not been suggested in the literature.
The Cramér-Rao lower bound (CRLB) as a benchmark for unbiased estimators is of pivotal importance for target parameter estimation [33][34][35]. In [36,37], the CRLB of colocated MIMO radar was derived. However, the structure of the likelihood function is changed due to the intrinsic property of coupling between pulse compression and beamforming in CSTCA. Thus, the CRLB in CSTCA is different from that in MIMO radar and deserving of further analysis.
Unlike existing DOA estimation studies that mainly concentrate on the characteristics of the covariance matrix in receive beamspace [38][39][40][41][42], we focus on the transmit beamspace to estimate DOAs in CSTCA. Our main contributions include: (1) We design a novel two dimensional space-time matched filter incorporating pulse compression with beamforming. It can form a cluster of equivalent transmit multi-beams to obtain degree-of-freedoms (DOFs) in transmit domain. (2) The beamspace MUSIC method in CSTCA is proposed. We devise the searching steering vector in beamspace and derive the spatial spectrum. Moreover, since multi-targets in the same range cell appears characteristics of coherence, we design a transformation matrix which can turn data from beamspace into element space. The desired property of rotational invariance property (RIP) at the receive array can be guaranteed through transformation. Eventually, we estimate DOAs by spatial smoothing. (3) The closed-form expressions of Cramer-Rao lower bound for angle estimation is derived for the first time.
The remainder of this paper is organized as follows: the signal model of circulating space-time coding array is formulated in Section 2. In Section 3, we present the scheme of bi-dimensional space-time matched filtering. Additionally, the beamspace MUSIC method and the spatial smoothing method in transform domain are proposed in application of different scenarios. Experimental results demonstrating the effectiveness of the proposed algorithms are given in Section 4, followed by conclusions drawn in Section 5.

Signal Model of Circulating Space-Time Coding Array
Consider a uniform linear M element CSTCA with half-wavelength spacing d. An identical linear frequency modulation (LFM) waveform circulates through every channel. The tiny time shift ∆t is employed across array elements. The transmitted signal of the mth element can be expressed as: where f 0 is the carrier frequency. The time shift is inversely proportional to signal bandwidth B. The transmitted waveform s(t) can be written as: where rect( t T p ) = 1, 0 ≤ t ≤ T p 0, else is the rectangular envelope, and µ = B/T p is the chirp rate with T p being the pulse duration. Furthermore, the time-variant waveform can be combined with the time-invariant spatial phase codes in space-time coding to improve the range resolution [2]. This kind of transmit diversity technique imposes an additional phase shift on the emitted signal. Meanwhile, the waveform itself is unchangeable. Herein, we associate LFM waveform as the time domain signal with Barker code c = [c 1 , c 2 , . . . , c M ] as the spatial domain signal. Thus, the m-th hybrid signal can be represented by: The overall transmitted signal at angle θ in range R can be modelled as: where τ = R/c is the time delay. As illustrated in Figure 1, the transmitted signal is highly overlapped along the time axis. The envelope of overall signal is similar to the trapezoid with time duration T p in −3 dB level due to the transient effect. Besides, the pulse duration of overall signal is expanded into T p + (M − 1)∆t leading to an edge effect [6]. Note that the transmitted signal of every element can't be summed in any range cell, since ∆t = 1/B exactly corresponds to the size of a range resolution cell, that is, c/2B. Thus, CSTCA is capable of full spatial coverage with a range-angle-dependent transmit beampattern as shown in Figure 2. It is shown that the mainlobe points to different directions at distinct time instants during the pulse duration.
The overall transmitted signal at angle θ in range R can be modelled as: where τ = R/c is the time delay. As illustrated in Figure 1, the transmitted signal is highly overlapped along the time axis. The envelope of overall signal is similar to the trapezoid with time duration Tp in −3 dB level due to the transient effect. Besides, the pulse duration of overall signal is expanded into Tp + (M − 1)Δt leading to an edge effect [6]. Note that the transmitted signal of every element can't be summed in any range cell, since Δt = 1/B exactly corresponds to the size of a range resolution cell, that is, c/2B. Thus, CSTCA is capable of full spatial coverage with a range-angle-dependent transmit beampattern as shown in Figure 2. It is shown that the mainlobe points to different directions at distinct time instants during the pulse duration.
The overall transmitted signal at angle θ in range R can be modelled as: where τ = R/c is the time delay. As illustrated in Figure 1, the transmitted signal is highly overlapped along the time axis. The envelope of overall signal is similar to the trapezoid with time duration Tp in −3 dB level due to the transient effect. Besides, the pulse duration of overall signal is expanded into Tp + (M − 1)Δt leading to an edge effect [6]. Note that the transmitted signal of every element can't be summed in any range cell, since Δt = 1/B exactly corresponds to the size of a range resolution cell, that is, c/2B. Thus, CSTCA is capable of full spatial coverage with a range-angle-dependent transmit beampattern as shown in Figure 2. It is shown that the mainlobe points to different directions at distinct time instants during the pulse duration.   Without loss of generality, suppose Q far-field point targets at angles θ 1 , . . . , θ Q at ranges R 1 , . . . , R Q . After backscattering, the received signal at the nth element can be concisely written as where ξ q denotes the complex scattering coefficient of the qth target, τ m,n,q = τ q − d(m − 1) sin θ q /c − d(n − 1) sin θ q /c is the round-trip time delay with τ q = 2R q /c, n(t) is the additive Gaussian white noise with zero-mean. Under the assumption of a narrowband sounding signal, we have . Therefore, (5) can be rewritten as After down-conversion, the received signal from N receive channels can be concisely formulated as where A R (θ) = a R (θ 1 ), a R (θ 2 ), . . . , a R θ Q is the N × Q receive array manifold with N × 1 receive steering vector a R θ q = 1, exp j2πd sin θ q /λ 0 , · · · , exp j2πd(N − 1) sin θ q /λ 0 , ξ Q e −j2π f 0 τ Q ) denotes a diagonal matrix with qth diagonal elements ξ q e j2π f 0 τ q , w 2 = [1, 1, . . . , 1] M×1 , n(t) is the N × 1 Gaussian white noise. stands for the Hadamard product, (·) T denotes the transposition operation.

Two-Dimensional Space-Time Matched Filtering Design
Based on the time-angle-dependent transmit beampattern in CSTCA, the conventional matched filter in fast time domain is no longer applicable. In this section, we devise a two-dimensional space-time matched filter to perform transmit beamforming and pulse compression simultaneously. The matched function can be formulated as: where θ i is the matched angle, namely, the beamforming direction. It satisfies the following condition: Under the constraint of (9), the steering vector can be turned into a discrete Fourier basis, that is, e −j2πmi/M with m = 0, . . . , M − 1, i = 0, . . . , M − 1. The M matched filters can perform ordinary beamforming to obtain M orthogonal beams for omnidirectional detection. Since we focus on the feature of transmit beamspace, we adopt non-adaptive beamforming at the receiver. The output after receive beamforming can be written as: where w R = [1, 1, . . . , 1] N×1 is the receive weight vector, τ i = 2R i /c, i = 1, . . . , Q. Then, the signal is fed back into the matched filter bank. The output of the ith matched filter channel can be written as: where ⊗ denotes the convolution operator. Note that the matched filtering can also be conducted before receive beamforming. Two processing procedures are theoretically equivalent with differences in practical application. Here, we utilise receive beamforming before matched filter to reduce the system complexity. Next, substitute (8) and (10) into (11), the output is expanded as: where the auto-correlation function ACF can be extended as The ACF is a function of both continuous time τ q and discrete time i q , m, n with i q = τ q /∆t. The envelope in (13) is similar with sinc function. Besides, it is shown in (12) that the equivalent transmit beampattern via matched filtering is determined by the exponential term, and the time shift ∆t has no influence on it. Thereby, the angular resolution in CSTCA is the same as that in phased array counterpart.
It is worth noticing that the space-time matched filter can bring an equivalent transmit beampattern gain with a value of 10log 10 M 2 . Additionally, the pulse compression gain is equal to 10log 10 BT/M, which is reduced by a factor of M in comparison with conventional pulse compression gain 10log 10 BT. The main reasons are listed as follows. At first, the power of transmitted signal is uniformly distributed in space. As a result, the power for one of the M beams is only 1/M of the total radiated power. Secondly, as shown in Figure 3, the time duration of the mainlobe (the segment between two dotted lines) is only the fraction 1/M of the whole pulse duration, since the width of mainlobe at a given probing angle is π/M in typical beampattern. Consequently, the effective coherent summation time of space-time matched filter is equal to T p /M, which means only 1/M of the bandwidth is sent in LFM signal. Thus, the pulse compression gain is degraded by a factor M. Furthermore, the range resolution is also reduced.
where ⨂ denotes the convolution operator. Note that the matched filtering can also be conducted before receive beamforming. Two processing procedures are theoretically equivalent with differences in practical application. Here, we utilise receive beamforming before matched filter to reduce the system complexity. Next, substitute (8) and (10) into (11), the output is expanded as: where the auto-correlation function ACF can be extended as The ACF is a function of both continuous time τq and discrete time iq, m, n with iq = τq/Δt. The envelope in (13) is similar with sinc function. Besides, it is shown in (12) that the equivalent transmit beampattern via matched filtering is determined by the exponential term, and the time shift Δt has no influence on it. Thereby, the angular resolution in CSTCA is the same as that in phased array counterpart.
It is worth noticing that the space-time matched filter can bring an equivalent transmit beampattern gain with a value of 10log10M 2 . Additionally, the pulse compression gain is equal to 10log10BT/M, which is reduced by a factor of M in comparison with conventional pulse compression gain 10log10BT. The main reasons are listed as follows. At first, the power of transmitted signal is uniformly distributed in space. As a result, the power for one of the M beams is only 1/M of the total radiated power. Secondly, as shown in Figure 3, the time duration of the mainlobe (the segment between two dotted lines) is only the fraction 1/M of the whole pulse duration, since the width of mainlobe at a given probing angle is π/M in typical beampattern. Consequently, the effective coherent summation time of space-time matched filter is equal to Tp/M, which means only 1/M of the bandwidth is sent in LFM signal. Thus, the pulse compression gain is degraded by a factor M. Furthermore, the range resolution is also reduced.

Beamspace MUSIC Algorithm
In this subsection, we propose a simple beamspace MUSIC algorithm. Unlike conventional beamspace DOA estimation methods, the spatial spectrum estimation in CSTCA can be directly performed in beamspace without design of beamforming matrix, since multiple beams are formed via space-time matched filtering.
Because different matched angles of filters give rise to different outputs. Peaks only appear via

Beamspace MUSIC Algorithm
In this subsection, we propose a simple beamspace MUSIC algorithm. Unlike conventional beamspace DOA estimation methods, the spatial spectrum estimation in CSTCA can be directly performed in beamspace without design of beamforming matrix, since multiple beams are formed via space-time matched filtering.
Because different matched angles of filters give rise to different outputs. Peaks only appear via the matched filter wherein the matched angle is closest to the true target angle. For other matched filters, it is mismatched. We regard the peak locations as the target ranges. Thus, the target range can be considered as a priori knowledge. The samples at peaks are picked out. They serve as training samples to construct the sample covariance matrix (SCM). Firstly, we can obtain: where To acquire a sufficient amount of training samples, we exploit the coherent pulse train in slow time domain. Received data from K pulses is arranged into the M×K SCM, written as: The data cube and matched filter configuration are portrayed in Figure 4 to visualize the processing. As can be seen, range cells of potential targets are selected and are combined with the K pulses. In practice, the SCM is basically estimated by secondary training samples referred to as sample matrix inverse (SMI) method, which is expressed as: where (·) H denotes the Hermitian operation. filters, it is mismatched. We regard the peak locations as the target ranges. Thus, the target range can be considered as a priori knowledge. The samples at peaks are picked out. They serve as training samples to construct the sample covariance matrix (SCM). Firstly, we can obtain: where [ ] To acquire a sufficient amount of training samples, we exploit the coherent pulse train in slow time domain. Received data from K pulses is arranged into the M×K SCM, written as: The data cube and matched filter configuration are portrayed in Figure 4 to visualize the processing. As can be seen, range cells of potential targets are selected and are combined with the K pulses. In practice, the SCM is basically estimated by secondary training samples referred to as sample matrix inverse (SMI) method, which is expressed as: where (·) H denotes the Hermitian operation.  The eigen-decomposition of SCM � yields � = = + , where ES and EN are the signal subspace matrix and orthogonal noise subspace matrix, ΛS = diag{λ1, …, λQ} and ΛN = diag{σ 2 , …, σ 2 } collect the eigenvalues with λ1 ≥ λ2 ≥ … ≥ λQ ≥ λQ+1 =…= σ 2 , respectively. Subsequently, the beamspace spatial spectrum is established as: where abeam(θ) is the M × 1 beamspace steering vector. It is noteworthy that abeam(θ) is distinguishable with the element space steering vector aelement(θ) = [1, exp{j2πdsinθ/λ},…, exp{j2π(M − 1)dsinθ/λ}]. We devise the beamspace steering vector by the following steps. Firstly, signal reflected from direction θ, range Rq via ith matched filter can be written as: Next, stack the outputs from M channels into a column vector: Then, the samples in target range bins are selected out to make up the beamspace searching vector, that is: where E S and E N are the signal subspace matrix and orthogonal noise subspace matrix, Subsequently, the beamspace spatial spectrum is established as: where a beam (θ) is the M × 1 beamspace steering vector. It is noteworthy that a beam (θ) is distinguishable with the element space steering vector a element (θ) = [1, exp{j2πdsinθ/λ}, . . . , exp{j2π(M − 1)dsinθ/λ}]. We devise the beamspace steering vector by the following steps. Firstly, signal reflected from direction θ, range R q via ith matched filter can be written as: Next, stack the outputs from M channels into a column vector: Then, the samples in target range bins are selected out to make up the beamspace searching vector, that is: Since the pulse compression and beamforming are implemented concurrently, data via matched filtering is decoupled in time and angle domains. It implies that the range bin corresponding to any one of targets contains the angle information to reflect characteristics of the beamspace steering vector. Therefore, we only choose an arbitrary target to construct the beamspace searching vector in (20), which can significantly reduce the computational burden compared with utilization of entire Q targets. Moreover, the amplitudes of Q spectral peaks in (20) are different. The spectral peak intensity of the chosen qth target is slightly higher than others. When we use all the target samples to obtain beamspace steering vector, the peak intensities of Q targets will tend to be the same. However, it is insignificant to performance improvements and increases the computational complexity.

Spatial Smoothing Method in Transform Domain
In this subsection, we consider multiple closely spaced targets whose spacing is less than a range resolution cell. Since the range cell consisting of different target samples can't satisfy the independent and identically distributed (IID) conditions and embodies feature of coherent signals, the proposed beamspace MUSIC method fails to estimate DOA accurately. To deal with this issue, we design a transformation matrix to apply spatial smoothing in element space.
Classical spatial smoothing technique makes use of the linear phase relationship between elements to guarantee nonsingularity of SCM. Motivated by the fact that the ordinary beamforming in beamspace via space-time matched filter is based on the discrete Fourier transform (DFT), we may conversely perform IDFT to return to element space. However, the beamforming matrix in space-time matched filter isn't strictly in accordance with the discrete Fourier basis. Correspondingly, we can't obtain the rotational invariance property in element space through IDFT directly. Remarkably, we design the transformation matrix with an approximate Vandermonde structure, which can be expressed as: Thus, the data matrix in element space can be written as: Figure 5 instantiates the regularity of amplitude-phase for data in element space. As can be seen, the amplitude is approximate to the unit value after transformation. Meanwhile, the phase is uniformly distributed at intervals of 90 • in case of a single target at θ 0 = 30 • , R 0 = 5 km. It verifies that the transformation procedure in (21) can restore the rotational invariance structure. In sequence, the range bins of targets are extracted out of data matrix X. The M × 1 data vector can be given as: where τ i = 2R i /c, i = 1, . . . , Q. Then, we can obtain the M × K data matrix by sampling from K pulses, which is expressed as: Next, we divide the M elements uniform linear array (ULA) into P = M − N 0 + 1 overlapping subarrays of size N 0 . Each subarray is separated by one element. The covariance matrix of the pth subarray can be written as:R where X p = [X (p), X (p + 1), · · · , X (p + N 0 + 1)] T is the N 0 × K data matrix of the pth subarray, X (p) denotes the pth row vector of X . The spatial smoothing matrix is the average of subarray covariance matrix, expressed as:R Eventually, we can adopt basic MUSIC method to estimate DOAs since we have transformed data into the element space. Because of limitation of space, we won't repeat it here.  (25) where = [ ( ), ( + 1), ⋯ , ( + + 1)] is the N0 × K data matrix of the pth subarray, ( ) denotes the pth row vector of . The spatial smoothing matrix is the average of subarray covariance matrix, expressed as: Eventually, we can adopt basic MUSIC method to estimate DOAs since we have transformed data into the element space. Because of limitation of space, we won't repeat it here.

Simulation Results
In this section, we consider a 13-element uniform linear CSTCA with half-wavelength spacing. The reference carrier frequency is f0 = 10 GHz and the time shift is Δt = 0.05 μs. We set SNR = −20 dB unless otherwise specified. Note that the said SNR is before matched filtering. The additive noise is zero-mean complex Gaussian distributed. We use 0.001° mesh grids to search candidate peaks, and conduct computer simulations with 200 independent Monte Carlo trials. The accuracy of DOA estimates is evaluated by the root mean square errors (RMSEs), expressed as: where Mon denotes the number of Monte Carlo trial runs, ( ) is the estimated angle of the qth target in the ith run, and θq is the true angle of the qth target. Another performance measurement metric is the resolution probability reflecting detection efficiency. In this case, we assume two closelyspaced point targets from directions θ1 = 4°and θ2 = 6°. It succeeds in resolving if and only if the following inequality holds [43]:

Simulation Results
In this section, we consider a 13-element uniform linear CSTCA with half-wavelength spacing. The reference carrier frequency is f 0 = 10 GHz and the time shift is ∆t = 0.05 µs. We set SNR = −20 dB unless otherwise specified. Note that the said SNR is before matched filtering. The additive noise is zero-mean complex Gaussian distributed. We use 0.001 • mesh grids to search candidate peaks, and conduct computer simulations with 200 independent Monte Carlo trials. The accuracy of DOA estimates is evaluated by the root mean square errors (RMSEs), expressed as: where Mon denotes the number of Monte Carlo trial runs,θ q (i) is the estimated angle of the qth target in the ith run, and θ q is the true angle of the qth target. Another performance measurement metric is the resolution probability reflecting detection efficiency. In this case, we assume two closely-spaced point targets from directions θ 1 = 4 • and θ 2 = 6 • . It succeeds in resolving if and only if the following inequality holds [43]: Figure 6 plots the one-dimensional range profile via space-time matched filtering in case of a single target at θ 1 = 0 • and R 1 = 5 km. The right side figure is the zoom-in of the range profile in range of [4.5 km, 5.5 km]. It is shown that the two peaks are situated at target locations. Note that the −3 dB width of the mainlobe for hybrid signal is evidently narrower than that for circulating LFM signal. As mentioned before in Section 3.1, the range resolution is degraded by a factor of M when transmitting circulating LFM signal in CSTCA. Particularly, the hybrid code scheme can significantly improve the range resolution. In addition, the sidelobe level for circulating LFM waveform is nearly below −40 dB in most range areas without window weighting. It is much lower than the SLL of LFM waveform itself, which is generally −13.2 dB. The reason is that transmitted signal of every element influences each other via sidelobes. As a result, the sidelobes can compensate each other for most ranges. Based on the positive impact of the sidelobes' interaction, the SLL is reduced. However, the SLL for hybrid code is increased considerably, wherein it is almost 20 dB higher than that for circulating LFM signal. Since the time and space are dependent in CSTCA, both the spatial waveform and the temporal waveform are demanded to have low range SLLs to reduce the synthesized SLL. However, the peak SLL of spatial phase code is much higher than that for pure circulating temporal waveform, since the code length is limited by the element number. Thus, the peak SLL is larger for hybrid code. We can conclude that the range resolution enhancement in hybrid code is at cost of the side lobe level. below −40 dB in most range areas without window weighting. It is much lower than the SLL of LFM waveform itself, which is generally −13.2 dB. The reason is that transmitted signal of every element influences each other via sidelobes. As a result, the sidelobes can compensate each other for most ranges. Based on the positive impact of the sidelobes' interaction, the SLL is reduced. However, the SLL for hybrid code is increased considerably, wherein it is almost 20 dB higher than that for circulating LFM signal. Since the time and space are dependent in CSTCA, both the spatial waveform and the temporal waveform are demanded to have low range SLLs to reduce the synthesized SLL. However, the peak SLL of spatial phase code is much higher than that for pure circulating temporal waveform, since the code length is limited by the element number. Thus, the peak SLL is larger for hybrid code. We can conclude that the range resolution enhancement in hybrid code is at cost of the side lobe level.

DOA Estimation Performance Analysis for Beamspace MUSIC Algorithm
In this example, we assess the performance of proposed beamspace MUSIC algorithm. The analytical expression for exact Cramer-Rao lower bounds (CRLBs) is derived to benchmark the method (see Appendix).
Assume that two targets are at θ1 = 0°, θ2 = 10°, and R1 = 3 km, R2 = 5 km. As shown in Figure 7, the proposed method can accurately estimate DOAs with two peaks at targets' locations, that is, θ1 = 0°and θ2 = 10°. The beamforming angle in space-time matched filter is uniformly distributed in form of sin(θ) in Figure 7a in contrast of a uniform distribution of θ in Figure 7b. As can be seen, the performance in Figure 7a outperforms that in Figure 7b. Since the steering vector 2 ( −1) sin / in matched filter is consistent with Discrete Fourier orthogonal basis in Figure 7a. As a result, the transformation from element space to beamspace via matched filtering is an orthogonal transformation, which can ascertain no loss of energy in the process. Therefore, the corresponding sample covariance matrix can provide better orthogonality between noise subspace and beamspace steering vector in (17). Moreover, there is little discrepancy between the two paradigms. Because the imposed Barker code alters the original phases of circulating LFM signal without any changes in the waveform itself.

DOA Estimation Performance Analysis for Beamspace MUSIC Algorithm
In this example, we assess the performance of proposed beamspace MUSIC algorithm. The analytical expression for exact Cramer-Rao lower bounds (CRLBs) is derived to benchmark the method (see Appendix A).
Assume that two targets are at θ 1 = 0 • , θ 2 = 10 • , and R 1 = 3 km, R 2 = 5 km. As shown in Figure 7, the proposed method can accurately estimate DOAs with two peaks at targets' locations, that is, θ 1 = 0 • and θ 2 = 10 • . The beamforming angle in space-time matched filter is uniformly distributed in form of sin(θ) in Figure 7a in contrast of a uniform distribution of θ in Figure 7b. As can be seen, the performance in Figure 7a outperforms that in Figure 7b. Since the steering vector e j2π(n−1)d sin θi/λ in matched filter is consistent with Discrete Fourier orthogonal basis in Figure 7a. As a result, the transformation from element space to beamspace via matched filtering is an orthogonal transformation, which can ascertain no loss of energy in the process. Therefore, the corresponding sample covariance matrix can provide better orthogonality between noise subspace and beamspace steering vector in (17). Moreover, there is little discrepancy between the two paradigms. Because the imposed Barker code alters the original phases of circulating LFM signal without any changes in the waveform itself. In this example, we consider two targets at θ1 = −3°, θ2 = 5°, respectively. The sample number K is set to be 100. The pulse compression gain is approximate to 10 dB, and the equivalent transmit beampatten gain is nearly 20 dB. Therefore, overall gain is 30 dB. The RMSEs with respect to input SNRs are plotted in Figure 8, where the SNRs vary over a range of −40 dB to 0 dB. It can be observed that the RMSEs decrease monotonically as SNRs increase. A little difference exists between the two cases (with or without Barker code). It is seen from Figure 8 that RMSEs of DOA estimation are close to the CRLBs at moderate and high SNRs regions. The accuracy of the angle estimation is enhanced with increment of SNRs. Figure 9 plots RMSEs versus the number of snapshots at SNR = −20 dB. It is shown that the angle can be estimated accurately at small number of snapshots. The direction finding performance attains the CRLBs when the snapshot number is more than 20.   In this example, we consider two targets at θ 1 = −3 • , θ 2 = 5 • , respectively. The sample number K is set to be 100. The pulse compression gain is approximate to 10 dB, and the equivalent transmit beampatten gain is nearly 20 dB. Therefore, overall gain is 30 dB. The RMSEs with respect to input SNRs are plotted in Figure 8, where the SNRs vary over a range of −40 dB to 0 dB. It can be observed that the RMSEs decrease monotonically as SNRs increase. A little difference exists between the two cases (with or without Barker code). It is seen from Figure 8 that RMSEs of DOA estimation are close to the CRLBs at moderate and high SNRs regions. The accuracy of the angle estimation is enhanced with increment of SNRs. Figure 9 plots RMSEs versus the number of snapshots at SNR = −20 dB. It is shown that the angle can be estimated accurately at small number of snapshots. The direction finding performance attains the CRLBs when the snapshot number is more than 20. In this example, we consider two targets at θ1 = −3°, θ2 = 5°, respectively. The sample number K is set to be 100. The pulse compression gain is approximate to 10 dB, and the equivalent transmit beampatten gain is nearly 20 dB. Therefore, overall gain is 30 dB. The RMSEs with respect to input SNRs are plotted in Figure 8, where the SNRs vary over a range of −40 dB to 0 dB. It can be observed that the RMSEs decrease monotonically as SNRs increase. A little difference exists between the two cases (with or without Barker code). It is seen from Figure 8 that RMSEs of DOA estimation are close to the CRLBs at moderate and high SNRs regions. The accuracy of the angle estimation is enhanced with increment of SNRs. Figure 9 plots RMSEs versus the number of snapshots at SNR = −20 dB. It is shown that the angle can be estimated accurately at small number of snapshots. The direction finding performance attains the CRLBs when the snapshot number is more than 20.   In this example, we consider two targets at θ1 = −3°, θ2 = 5°, respectively. The sample number K is set to be 100. The pulse compression gain is approximate to 10 dB, and the equivalent transmit beampatten gain is nearly 20 dB. Therefore, overall gain is 30 dB. The RMSEs with respect to input SNRs are plotted in Figure 8, where the SNRs vary over a range of −40 dB to 0 dB. It can be observed that the RMSEs decrease monotonically as SNRs increase. A little difference exists between the two cases (with or without Barker code). It is seen from Figure 8 that RMSEs of DOA estimation are close to the CRLBs at moderate and high SNRs regions. The accuracy of the angle estimation is enhanced with increment of SNRs. Figure 9 plots RMSEs versus the number of snapshots at SNR = −20 dB. It is shown that the angle can be estimated accurately at small number of snapshots. The direction finding performance attains the CRLBs when the snapshot number is more than 20.    Figure 10 displays the resolution probability against SNRs at K = 50. It is shown that the resolution probability becomes larger when SNRs increase. It achieves 100% at SNR = −18 dB. In Figure 11, the variation of resolution probability with number of snapshots is depicted at SNR = −20 dB. It is observed that the resolution probability increases with a larger sample number. In Figure 12, the dependence between resolution probability and the angle separation is plotted at SNR = −20 dB and K = 50. The angles of two targets are set as θ 1 = θ cen − ∆θ/2 and θ 2 = θ cen + ∆θ/2 with central angle θ cen = 5 • . The angle separation is from 1 • to 6 • . The minimum angular separation for a specific resolution performance can be determined from the curves. It is shown that the resolution probability increases monotonically with increased angle separations in the range of [2.5 • , 4 • ]. In addition, measurements of two transmit waveforms show slightly differences in DOA estimation performance in Figures 10-12. Sensors 2018, 18, x FOR PEER REVIEW 11 of 18 Figure 10. displays the resolution probability against SNRs at K = 50. It is shown that the resolution probability becomes larger when SNRs increase. It achieves 100% at SNR = −18 dB. In Figure 11, the variation of resolution probability with number of snapshots is depicted at SNR = −20 dB. It is observed that the resolution probability increases with a larger sample number. In Figure 12, the dependence between resolution probability and the angle separation is plotted at SNR = −20 dB and K = 50. The angles of two targets are set as θ1 = θcen − Δθ/2 and θ2 = θcen + Δθ/2 with central angle θcen = 5°. The angle separation is from 1° to 6°. The minimum angular separation for a specific resolution performance can be determined from the curves. It is shown that the resolution probability increases monotonically with increased angle separations in the range of [2.5°, 4°]. In addition, measurements of two transmit waveforms show slightly differences in DOA estimation performance in Figures 10-12.     Figure 10. displays the resolution probability against SNRs at K = 50. It is shown that the resolution probability becomes larger when SNRs increase. It achieves 100% at SNR = −18 dB. In Figure 11, the variation of resolution probability with number of snapshots is depicted at SNR = −20 dB. It is observed that the resolution probability increases with a larger sample number. In Figure 12, the dependence between resolution probability and the angle separation is plotted at SNR = −20 dB and K = 50. The angles of two targets are set as θ1 = θcen − Δθ/2 and θ2 = θcen + Δθ/2 with central angle θcen = 5°. The angle separation is from 1° to 6°. The minimum angular separation for a specific resolution performance can be determined from the curves. It is shown that the resolution probability increases monotonically with increased angle separations in the range of [2.5°, 4°]. In addition, measurements of two transmit waveforms show slightly differences in DOA estimation performance in Figures 10-12.      Figure 10. displays the resolution probability against SNRs at K = 50. It is shown that the resolution probability becomes larger when SNRs increase. It achieves 100% at SNR = −18 dB. In Figure 11, the variation of resolution probability with number of snapshots is depicted at SNR = −20 dB. It is observed that the resolution probability increases with a larger sample number. In Figure 12, the dependence between resolution probability and the angle separation is plotted at SNR = −20 dB and K = 50. The angles of two targets are set as θ1 = θcen − Δθ/2 and θ2 = θcen + Δθ/2 with central angle θcen = 5°. The angle separation is from 1° to 6°. The minimum angular separation for a specific resolution performance can be determined from the curves. It is shown that the resolution probability increases monotonically with increased angle separations in the range of [2.5°, 4°]. In addition, measurements of two transmit waveforms show slightly differences in DOA estimation performance in Figures 10-12.

DOA Estimation Performance Analysis for Spatial Smoothing in Transform Domain
In this example, we demonstrate the effectiveness of the spatial smoothing technique in transform domain in case of two targets in the same range cell under test (CUT). The CSTCA is divided into 3 overlapped subarrays. The other parameters are identical with those in beamspace MUSIC simulations. Figure 13 plots the spatial spectrum of two transmit waveforms. It is observed that the spectral peaks point to true target angles at θ 1 = −5 • , θ 2 = 15 • and R 1 = R 2 = 5 km, respectively. Moreover, the spatial spectrum of pure circulating LFM waveform has lower background noise level and more evident peak than the hybrid code. The reason is that Barker code introduces random phases across array elements, which has adverse influence on spatial phase relationship in element domain. In this regard, pure circulating LFM waveform is more beneficial to angle estimation.

DOA Estimation Performance Analysis for Spatial Smoothing in Transform Domain
In this example, we demonstrate the effectiveness of the spatial smoothing technique in transform domain in case of two targets in the same range cell under test (CUT). The CSTCA is divided into 3 overlapped subarrays. The other parameters are identical with those in beamspace MUSIC simulations. Figure 13 plots the spatial spectrum of two transmit waveforms. It is observed that the spectral peaks point to true target angles at θ1 = −5°, θ2 = 15° and R1 = R2 = 5 km, respectively. Moreover, the spatial spectrum of pure circulating LFM waveform has lower background noise level and more evident peak than the hybrid code. The reason is that Barker code introduces random phases across array elements, which has adverse influence on spatial phase relationship in element domain. In this regard, pure circulating LFM waveform is more beneficial to angle estimation. The RMSEs curves versus SNRs at K = 50 are displayed in Figure 14. It is shown that the angular RMSEs descend linearly when SNR > −30 dB demonstrating the excellent performance in angular estimation. Figure 15 plots RMSEs curves versus snapshots, wherein two curves exhibit similar variation trends in comparison with Figure 9. Furthermore, the CRLB curve in element space is analogous to that in beamspace in Figures 8 and 9 for the following two reasons. Firstly, the transformation from beamspace to element space through (21) doesn't introduce additional information such as new independent training samples. Thereby, the structure of SCM isn't changed. Secondly, the total energy is lossless in the orthogonal transformation step by (21). The reduced array aperture in subarray partition would cause a little gain loss. For the sake of simplicity, the CRLB curve in element space isn't derived here.  The RMSEs curves versus SNRs at K = 50 are displayed in Figure 14. It is shown that the angular RMSEs descend linearly when SNR > −30 dB demonstrating the excellent performance in angular estimation. Figure 15 plots RMSEs curves versus snapshots, wherein two curves exhibit similar variation trends in comparison with Figure 9. Furthermore, the CRLB curve in element space is analogous to that in beamspace in Figures 8 and 9 for the following two reasons. Firstly, the transformation from beamspace to element space through (21) doesn't introduce additional information such as new independent training samples. Thereby, the structure of SCM isn't changed. Secondly, the total energy is lossless in the orthogonal transformation step by (21). The reduced array aperture in subarray partition would cause a little gain loss. For the sake of simplicity, the CRLB curve in element space isn't derived here.

DOA Estimation Performance Analysis for Spatial Smoothing in Transform Domain
In this example, we demonstrate the effectiveness of the spatial smoothing technique in transform domain in case of two targets in the same range cell under test (CUT). The CSTCA is divided into 3 overlapped subarrays. The other parameters are identical with those in beamspace MUSIC simulations. Figure 13 plots the spatial spectrum of two transmit waveforms. It is observed that the spectral peaks point to true target angles at θ1 = −5°, θ2 = 15° and R1 = R2 = 5 km, respectively. Moreover, the spatial spectrum of pure circulating LFM waveform has lower background noise level and more evident peak than the hybrid code. The reason is that Barker code introduces random phases across array elements, which has adverse influence on spatial phase relationship in element domain. In this regard, pure circulating LFM waveform is more beneficial to angle estimation. The RMSEs curves versus SNRs at K = 50 are displayed in Figure 14. It is shown that the angular RMSEs descend linearly when SNR > −30 dB demonstrating the excellent performance in angular estimation. Figure 15 plots RMSEs curves versus snapshots, wherein two curves exhibit similar variation trends in comparison with Figure 9. Furthermore, the CRLB curve in element space is analogous to that in beamspace in Figures 8 and 9 for the following two reasons. Firstly, the transformation from beamspace to element space through (21) doesn't introduce additional information such as new independent training samples. Thereby, the structure of SCM isn't changed. Secondly, the total energy is lossless in the orthogonal transformation step by (21). The reduced array aperture in subarray partition would cause a little gain loss. For the sake of simplicity, the CRLB curve in element space isn't derived here.  The resolution probability curves versus input SNRs, the number of snapshots and the angular separation are depicted in Figures 16-18, respectively. As can be seen in Figure 16, the resolution probability for pure circulating space-time LFM signal is equal to 100% when SNRs ≥ −12 dB with K = 200. While it is −9 dB for the Barker code. It indicates that the circulating space-time LFM signal is more sensitive to the variation of SNRs. As illustrated in Figure 17, the probability of resolution is improved as the increase of snapshot number with SNRs = −25 dB. It is worth noting that the coherent pulse trains as sample data can also be applied to Doppler processing subsequently in future work. As shown in Figure 18, a good performance can be achieved when the angular separation is in excess of 4° in both curves.   The resolution probability curves versus input SNRs, the number of snapshots and the angular separation are depicted in Figures 16-18, respectively. As can be seen in Figure 16, the resolution probability for pure circulating space-time LFM signal is equal to 100% when SNRs ≥ −12 dB with K = 200. While it is −9 dB for the Barker code. It indicates that the circulating space-time LFM signal is more sensitive to the variation of SNRs. As illustrated in Figure 17, the probability of resolution is improved as the increase of snapshot number with SNRs = −25 dB. It is worth noting that the coherent pulse trains as sample data can also be applied to Doppler processing subsequently in future work. As shown in Figure 18, a good performance can be achieved when the angular separation is in excess of 4 • in both curves. The resolution probability curves versus input SNRs, the number of snapshots and the angular separation are depicted in Figures 16-18, respectively. As can be seen in Figure 16, the resolution probability for pure circulating space-time LFM signal is equal to 100% when SNRs ≥ −12 dB with K = 200. While it is −9 dB for the Barker code. It indicates that the circulating space-time LFM signal is more sensitive to the variation of SNRs. As illustrated in Figure 17, the probability of resolution is improved as the increase of snapshot number with SNRs = −25 dB. It is worth noting that the coherent pulse trains as sample data can also be applied to Doppler processing subsequently in future work. As shown in Figure 18, a good performance can be achieved when the angular separation is in excess of 4° in both curves.   The resolution probability curves versus input SNRs, the number of snapshots and the angular separation are depicted in Figures 16-18, respectively. As can be seen in Figure 16, the resolution probability for pure circulating space-time LFM signal is equal to 100% when SNRs ≥ −12 dB with K = 200. While it is −9 dB for the Barker code. It indicates that the circulating space-time LFM signal is more sensitive to the variation of SNRs. As illustrated in Figure 17, the probability of resolution is improved as the increase of snapshot number with SNRs = −25 dB. It is worth noting that the coherent pulse trains as sample data can also be applied to Doppler processing subsequently in future work. As shown in Figure 18, a good performance can be achieved when the angular separation is in excess of 4° in both curves.

Conclusions
The circulating space-time coding array (CSTCA) radar employs a tiny time shift across array elements. It can provide a simple way to acquire a wide angular coverage with a stable gain by transmitting an identical waveform, and obtains a low sidelobe level in range domain. As a special type of coherent colocated MIMO radar, CSTCA is much simpler to be implemented in practice. In this paper, we solve the issue of DOA estimation in CSTCA. At first, we designed a two-dimensional space-time matched filter to form multi-beams affording controllable transmit freedom, and the pulse compression is performed simultaneously. Then, we proposed the beamspace MUSIC method to estimate DOAs. We devised the beamspace searching vector to derive the spatial spectrum. Afterwards, we designed a transformation matrix to map the received data cube from beamspace into element space in case of closely-spaced targets. Based on the transformation, the RIP can be restored and spatial smoothing is performed. Theoretical performance analysis and numerical simulations on the RMSEs, the probability of resolution as well as CRLB curve can demonstrate the effectiveness of proposed approaches. Our future work will focus on the practical obstacles including DOA estimation under additive colored Gaussian noise scenarios. Multiple parameters optimization problem for performance enhancement is also one of ongoing investigations.
where z(k) denotes the kth column of Z in (15), y(k) denotes the signal components of Z . Then, the log-likelihood function is: First, we calculate the derivatives of (30) with respect to θ, ξ and ξ. We have: ∂ ln L ∂θ