Design of multi-carrier MIMO radar array for DOA estimation

Multi-carrier (MC) multiple-input multiple-output (MIMO) radar offers an additional degree of freedom in the array optimization through the carrier frequencies. In this paper, we study the MC-MIMO array optimization with respect to the direction of arrival (DOA) estimation based on the Cramer-Rao bound (CRB). In particular, we choose the transmit and receive antenna positions as well as the carrier frequencies to minimize the single-target CRB subject to a constraint of the peak sidelobe level. A genetic algorithm is used to solve the problem and numerical examples demonstrate the superiority of our approach over both single-carrier MIMO radar and existing design rules.


I. INTRODUCTION
To improve the direction of arrival (DOA) estimation accuracy, multiple-input multiple-output (MIMO) radar proved to be a suitable technique in comparison to conventional phased array [1], [2]. The idea is a virtual array with M ×N antennas, where M is the number of transmit (Tx) and N is the number of receive (Rx) antennas. The concept of multi-carrier (MC) MIMO extends the MIMO idea by using C carrier frequencies to obtain a virtual array of C × M × N channels [3]. The steering vector of this virtual array depends on the antenna positions in wavelength. Hence the use of multiple carriers results in a scaling of the single-carrier (SC) virtual array. This enables a large aperture of the virtual non-uniform array for a high DOA estimation accuracy while satisfying the spatial Nyquist condition to avoid grating lobes at the same time. The price of such an MC-MIMO radar in comparison to SC-MIMO is however a coupling of the DOA and range estimation [4], [5], [3]. Though the coupling can be resolved in the Cramer-Rao bound (CRB) with a certain choice of the coordinate system (local decoupling) [6], the calculation of the sidelobes as well as the estimation have to be done in the joint range-DOA domain.
In the literature, the term "Multi-Carrier Radar" appears also to describe Orthogonal Frequency Multiplexing (OFDM, [7]), where the available frequency band is divided into subcarriers to simultaneously transmit symbols at different Tx of a MIMO system. In contrast, this work uses multiple carrier frequencies to improve the DOA estimation and is independent of the recent progress in OFDM. Furthermore, some of the listed literature uses the term "Frequency Diverse Radar" to describe this MC approach, however Frequency Diverse Radar is also the name of a technique to overcome frequency-dependent fluctuation of the radar cross section (RCS, [8]) or to intentionally create a coupling of range and DOA [4]. We would like to mention the difference to those techniques here.
The CRB was calculated in [9] for a SC-MIMO and in [6] for MC-MIMO radar. It expresses the achievable DOA estimation accuracy of unbiased estimators. It depends on the positions of the Tx and Rx [1] for SC-MIMO and additionally on the carrier frequencies for MC-MIMO radar [6].
To obtain a high DOA estimation accuracy, the design of the array is thus important. Existing approaches optimizes the correlation matrix of probing signals of an MIMO radar with given antenna positions [2], [10], [11]. The minimization of the deviation (in the least squares sense) from a desired beampattern (beampattern matching) is another option [2], [10], [11]. [12] gives basic design rules for the antenna positions of a linear monostatic MIMO radar. [13] suggests an optimization of the antenna positions to minimize the single-target CRB of a phased array subject to a given maximum sidelobe level using a genetic algorithm (GA). [3] gives design rules to avoid spatial aliasing in a colocated MC-MIMO system, but without taking the CRB or the sidelobes into account.
In this work, we extend the array optimization in [13] to the case of a colocated MC-MIMO radar. We want to design an array with a low single-target CRB of DOA for a high-accuracy DOA estimation while not exceeding a given maximum sidelobe level of the range-DOA ambiguity function to improve the noise and interference robustness of the DOA estimator. One major difference to [13] is the virtual array of size C × M × N of the MC-MIMO radar which is determined by the M + N antenna positions and C carrier frequencies.
Hence, the number of parameters to be optimized for a MC-MIMO system is much lower than that of a phased array with the same number of channels (e.g. only 7 parameters must be optimized for a 2 × 4 × 4 MC-MIMO linear array with 32 channels, whereas 31 parameters are required for a comparable phased array, see Sec. IV-B). For this purpose, we formulate an optimization problem which yields a high DOA estimation accuracy and a limited sidelobe-level. The optimization problem is examined (i.e. we derive the range in which the ambiguity function must be evaluated) and compared to the phased array design problem. Experiments compare the proposed array design approach to existing MC-MIMO and single-carrier design approaches.
The rest of this paper is organized as follows: We introduce the signal model in Sec. II. The ambiguity function and the CRB used as design metrics are stated in Sec. III. Sec. IV describes the array optimization. Results are discussed in Sec. V. Sec. VI concludes the paper.
We use the following notations in this work: Underlined variables denote vectors while bold variables denote matrices. j is the imaginary unit, ⊙ the elementwise Hadamard product and ⊗ the Kronecker tensor product. 1 N is a vector of N ones.

II. MULTI-CARRIER MIMO SIGNAL MODEL
We focus on the array optimization for DOA estimation in this paper and make the following assumptions: • The MC-MIMO radar has M Tx and N Rx antennas and C carrier frequencies which result in a virtual array with C × M × N channels. • We consider the single-target case. Multiple targets will be studied in another publication due to limited space. • The signals of the different channels are orthogonal using an appropriate multiplexing scheme (e.g. TDM, FDM, CDM [14] or OFDM). • All antennas are isotropic and have the same antenna gain. The same array design approach can be easily extended to known non-isotropic characteristics. • Pulse compression was already performed, i.e. symbols such as frequency-modulated continuous wave (FMCW) chirps or OFDM symbols are already removed. The signal vector contains the signal of one range-Doppler cell. We make this assumption only to provide an array design tool independent of the modulation. In our model, the coupling of range and DOA is still present and is not assumed to be resolved by the pulse compression. • The noise is an additive temporally and spatially white Gaussian random process with variance σ 2 . • For simplicity, we assume a stationary scene in this paper, i.e. we do not estimate the Doppler caused by a moving target. This will also be done in a future work. The vector x ∈ C CMN contains the baseband signals of all C × M × N channels. The vector of unknown parameters θ = [u T , r] T contains the range r and the electrical angles u ∈ R 2 of the target [3]: The rows of P Tx ∈ R M×2 and P Rx ∈ R N

A. Cramer-Rao bound
The CRB gives a lower bound for the variance of unbiased estimators. For multi-carrier arrays, the joint CRB of DOA and range was derived in [6]. Without loss of generality, we assume in the following that the origin of the coordinate system is always placed in the centroid of the virtual array. In this case, the CRBs of DOA and range are decoupled [6]. This allows us to study the CRB of DOA independent of range and simplifies the expression of the CRB.
The individual CRBs of the two electrical angles u x and u y for a general MC-MIMO planar array is given by γ f = ||f || 2 is the sample power of the carrier frequencies and is the sample correlation matrix (except for a missing normalization 1 MN ) of the (x,y)-position of the single-carrier virtual array P. It is equal to the sample covariance matrix (except for a scaling) since the origin of the coordinate system is assumed to be in the array centroid. ρ = ||α|| 2 is the signal power of the target.
In the case of a linear array, the CRB is

B. Ambiguity function
While the CRB is a local property of the likelihood function around the true parameter value, the (noise-free) ambiguity function also describes the global properties of an array such as sidelobes. The noise-free ambiguity function is the power of the correlation of the steering vector a(θ 0 ) with θ 0 = [u 0x , u 0y , r 0 ] T (or θ 0 = [u 0x , r 0 ] T for a linear array) of the true target with the steering vector of any hypothesized target a(θ) [15]. Though the choice of the coordinate system leads to a range-DOA decoupling of the CRB, the ambiguity function does not decouple in general. While for SC radar the ambiguity function of DOA (known as beampattern) and range decouple due to space-time separability [16], this is not the case for MC radar [6]. The range leads to different phase-shifts for different carrier frequencies which are superimposed by the phaseshifts caused by DOA. Hence the ambiguity function must be evaluated at all combinations of range and electrical angles to find the height of the sidelobes. The ambiguity function of θ is 3-dimensional for a planar array or 2-dimensional for a linear array.
The sidelobes reflect the correlation between the signals and hence the probability of global errors in the DOA estimation. High sidelobes are unwanted, as they generally lead to an increased SNR threshold where the estimation reaches the CRB [17] and to a performance degradation when multiple targets are present. A common measure of the sidelobes is the highest sidelobe level or peak sidelobe level (PSLL).
To determine the PSLL, the ambiguity function must be evaluated in a certain parameter range. Due to symmetry (θ 0 and θ are exchangeable in eq. 9), it is sufficient to evaluate the ambiguity function for a fixed u 0x and u 0y (e.g. u 0x = 0 and u 0y = 0) and varying u x and u y , where |u 0x − u x | ≤ 2 and |u 0y − u y | ≤ 2 to observe all possible sidelobes. Though e.g. u 0x = 0 and u x = 2 is physically not possible due to the definition of ||u|| ≤ 1, this formulation is useful in eq. 9, as β(θ 0 , θ) has the same value as for u 0x = −1 and u x = 1. If the region of interest for u can be limited, an evaluation of the ambiguity function in a smaller region is sufficient.
According to [5], the multi-carrier ambiguity function is periodic in the range direction. Extending [5] to the case of non-equidistant carrier frequencies, the uniform grid with the largest grid spacing determines the unambiguous range r amb . Numerically, we can find r amb by calculating the greatest common divisor (gcd) of the carrier frequencies differences δ f . As the gcd is defined only for integers, we have to normalize the carrier frequencies with respect to machine precision first. The unambiguous range is then given by r amb = c /gcd(δ f ). r amb can be much smaller than the largest range of the pulse compression (FMCW processing), e.g. r amb = c /|f1−f2| in the case of two carriers. Therefore we evaluate the ambiguity function between an arbitrary r 0 and r 0 + r amb . We note that the ambiguous range and hence the values of θ for which the ambiguity function is evaluated to find the PSLL depend on the optimization parameter f . Therefore the ambiguous range must be recalculated in each step of the optimization.

A. Problem formulation
For the design of a MC-MIMO array, we optimize the positions of the Tx and Rx antennas P Tx and P Rx as well as the carrier frequencies f . Our approach is the minimization of a cost function with respect to some constraints:

arg min
PTx,PRx,f η · CRB ux + (1 − η) · CRB uy (10) To make the design of the MC-MIMO array independent of the DOA estimator (e.g. Bartlett Beamformer or MUSIC [18]), we use a cost function based on the CRB. We combine both CRBs in eq. 4 and eq. 5 into a single cost function by a weighted combination: Choosing η ∈ [0, 1] allows a weighting of the DOA estimation accuracies in the x-and y-direction. η can be set to η = 0.5 for an equal weighting or to η = 1 or η = 0 if only the x-or y-direction should be considered (e.g. for linear arrays). The first line of constraints are bounds for the carrier frequencies and antenna positions. R f denotes the available region for the carrier frequencies f i . Similarly, all Tx and Rx antenna positions p Tx,i ∈ R 2 and p Rx,i ∈ R 2 are restricted to the regions R Tx and R Rx , respectively. They define in particular the maximum aperture allowed for the physical Tx and Rx array. The second line of constraints define the minimum antenna spacing due to antenna considerations like patch size, tapering and crosstalk. The last constraint defines the maximum PSLL value PSLL max which is not allowed to be exceeded in order to enhance the noise and interference robustness of the DOA estimation.
The PSLL value has to be determined from the complete sidelobe region which starts at the first minimum of the ambiguity function after the mainlobe. In this region, the ambiguity function is evaluated on a grid to get a coarse estimate of the PSLL. The location of the maximum is then used as the starting value for a local maximization of the ambiguity function. A gradient-based or downhill simplex (Nelder-Mead) method can be used for this purpose. (Such a local maximization is only applied to refine the PSLL found with the grid search; the grid search must already sufficiently resolve the ambiguity function to ensure finding the highest sidelobe level.) The result is the location and value of the PSLL.

B. Complexity consideration
In the common case of a contiguous R Tx and R Rx , e.g. R Tx and R Rx only limit the array aperture, the number of optimization parameters can be slightly reduced. The CRB and ambiguity function do not depend on absolute position of the virtual array P, hence we can assign the outermost (not to affect the constraints) position of the aperture to one Tx and Rx antenna and exclude them from the set of optimization parameters.
If R f , R Tx , R Rx and δ min should be normalized to one of the carrier frequencies and R f is contiguous, one can similarly assign the highest (not to affect the constraints) frequency of R f to one element of f . This carrier is then removed from the optimization set. A corresponding scaling of f , P and r would result in the same term Pθ in eq. 1 and hence the same PSLL and CRB.
The total number of parameters to be optimized is in general C + (M + N )d. If both of the above simplification can be applied, the complexity is reduced to (C − 1)+ (M + N − 2)d. d is the dimension of the array, i.e. d = 1 for a linear array and d = 2 for a planar array. This is much less than (CM N − 1)d parameters for a SC phased array with the same number of channels as in [13].

C. Optimization algorithm
Eq.10 poses a non-convex optimization problem. A local optimization is hence not suitable, as it can converge to a poor local optimum.
To overcome this problem, we apply a genetic algorithm which is known to be suitable for non-convex optimization tasks [19], [20] due to the random nature of the mutation step. Though genetic algorithms may not converge to the global optimum, they will most likely find a good local optimum. We apply a GA, because GA were already successfully applied for phased array design in [13], however other global optimizers such as simulated annealing might also be suitable.
The idea of genetic algorithms is a set of possible solutions (the population) which is improving in each iteration (generation). The children of a population are generated by altering the best individuals of one generation (parents), by either crossover (random combination of parameter values of two parents) or mutation (addition of a Gaussian distributed random value to each parameter of one parent). The lower the cost function of a individual, the higher the probability that it will be selected as a parent of the next generation. We use the Augmented Lagrangian Genetic Algorithm (ALGA) to solve the MC-MIMO design problem, but other genetic algorithms can be suitable as well. The idea of ALGA is to implement the nonlinear constraints as penalty terms with iteratively re-estimated weights which are added to the cost function value and then apply an unconstrained optimization [21]. In total, the virtual array has C × M × N = 32 different channels. For the optimization of a 32-channel linear phased array as in [13], 31 parameters must be optimized. By our MC-MIMO radar, only (C − 1) + (M + N − 2)d = 7 parameters need to be optimized.
The ALGA generated 20% of the child population by mutation and 80% by crossover. The population size was 10,000. We used the MATLAB implementation of ALGA from the global optimization toolbox. The result is an array with p  Table I lists the performance metrics for the design results. The achievable performance is expressed in terms of γ f · γ xx , which is inversely proportional to the CRB (but independent of the SNR). The PSLL is given as well. Fig. 1 shows the range-DOA ambiguity functions for θ 0 = [−1, 0] to illustrate the shape of the sidelobes. The ambiguity function of the SC MIMO radar is independent of range, while the MC MIMO ambiguity function must be evaluated for different ranges to find the PSLL. Fig. 2 shows a cut through the range-DOA ambiguity function at r 0 to illustrate the shape of the main beam and to better compare the ambiguity functions of the different design approaches. The sidelobes can only be partially observed in Fig. 2, as they might occur at other ranges (see Fig. 2). It    can be observed that the design rule of [3] is better than the optimized SC-MIMO array in terms of CRB, but at the cost of higher sidelobes. The optimized MC-MIMO array outperforms [3] in terms of both CRB and PSLL, and the optimized SC-MIMO array in terms of CRB at similar PSLL. Fig. 3 shows the performance of the deterministic maximum likelihood estimator (DMLE) [3] for a single target at u 0x = −0.8. The root-mean-square error (RMSE) of the DOA estimation and root-CRB are plotted over the SNR of one channel. It can be observed that all DMLE achieve their CRBs for high SNR, confirming the CRB as a good design criterion. Furthermore, the threshold where the RMSE of the estimator achieves the √ CRB of the optimized SC-MIMO and MC-MIMO array are equal, while the threshold of [3] is slightly higher. This is due to the height of the sidelobes. Below the threshold point the SC-MIMO would outperform the MC-MIMO system, but arrays are usually not used for such a low SNR. We show this region only to demonstrate the influence of the sidelobes on the threshold point. Beyond the threshold point, the SC-MIMO has the worst performance, followed by the design rule of [3] while the proposed design shows the best (lowest RMSE) results. The exploitation of MC hence enables a better array design not only in terms of technical metrics like CRB and PSLL, but also in terms of RMSE.
VI. SUMMARY In this paper we examine the optimization of an MC-MIMO array. We propose a CRB-based cost function with the PSLL as a constraint for the optimization. The optimization is solved using a genetic algorithm. Experiments demonstrate the through MC, a better array design is achievable in terms of PSLL and CRB as well as in terms of RMSE.
ACKNOWLEDGMENT This work was supported by the German Federal Ministry of Education and Research, Grant No. 13N13480.