Designing Waveform Sets with Good Correlation and Stopband Properties for MIMO Radar via the Gradient-Based Method

Waveform sets with good correlation and/or stopband properties have received extensive attention and been widely used in multiple-input multiple-output (MIMO) radar. In this paper, we aim at designing unimodular waveform sets with good correlation and stopband properties. To formulate the problem, we construct two criteria to measure the correlation and stopband properties and then establish an unconstrained problem in the frequency domain. After deducing the phase gradient and the step size, an efficient gradient-based algorithm with monotonicity is proposed to minimize the objective function directly. For the design problem without considering the correlation weights, we develop a simplified algorithm, which only requires a few fast Fourier transform (FFT) operations and is more efficient. Because both of the algorithms can be implemented via the FFT operations and the Hadamard product, they are computationally efficient and can be used to design waveform sets with a large waveform number and waveform length. Numerical experiments show that the proposed algorithms can provide better performance than the state-of-the-art algorithms in terms of the computational complexity.


Introduction
Waveform design has received considerable attention in recent years [1] and been employed in many applications, including polarimetric radar [2], multiple-input multiple-output (MIMO) radar [3,4], stealth communications [5] and the code-division multiple-access (CDMA) system. As a category of the general waveform design research [1], the design of waveform sets (i.e., multidimensional waveforms) is an important research content of MIMO radar. Generally, waveform sets are desired to have a good correlation property, which can effectively improve radar resolution, detection performance, imaging quality, the ability to obtain information and the accuracy of MIMO channel estimation [6][7][8]. In recent years, a large number of scholars has been devoted to designing waveform sets with a good correlation property. The main research covers two aspects: one is the waveform sets with good auto-and cross-correlation properties [9][10][11][12][13][14][15][16], and the other is the complementary sets of sequences (CSS) [17][18][19][20][21][22][23][24][25].
Waveform sets with good auto-and cross-correlation properties, also known as the orthogonal waveform set (OWS), have low autocorrelation sidelobes and low cross-correlation levels. In the early stage, the simulated annealing- [9] and cross entropy-based [10] methods were proposed for OWS design. However, due to the high computational complexity, these methods are not suited to design long waveforms. To improve the computational efficiency, the CAN (cyclic algorithm-new) algorithm (gradient-correlation-SFW), which is faster than the Gra-WeCorr-SFW, for the design problem without considering the correlation weights.
The rest of the paper is organized as follows. In Section 2, the design problem is formulated. In Section 3, we develop a gradient-based algorithm by deducing the phase gradient and the step size and then summarize the algorithm. In Section 4, the simplified algorithm for the design problem without considering the correlation weights is derived. Section 5 provides several numerical experiments to verify the effectiveness of the proposed algorithms. Finally, Section 6 gives the conclusions.
Notation: Boldface upper case letters denote matrices, while boldface lower case letters denote column vectors. (·) * , (·) T and (·) H denote the complex conjugate, transpose and conjugate transpose, respectively. · and · F denote the Euclidean norm and the Frobenius norm. Re(·) and Im(·) denote the real and imaginary part, respectively. Diag(x) denotes a diagonal matrix formed with the column vector x. • denotes the Hadamard product. x(m) denotes the m-th element of the vector x. x l is the l-th iteration of x. 1 N is the all-one vectors of length N. I N denotes the N × N identity matrix. F (x) and F −1 (x) denote the 2N-point FFT and IFFT (inverse FFT) operations of x, respectively. F c (X) and F −1 c (X) represent the FFT and IFFT of each column of the matrix X, respectively. In the (I)FFT operations, if the length of x is less than 2N, x is padded with trailing zeros to length 2N. e (·) is the exponent arithmetic applied to the scalar, vector or matrix.

Problem Formulation
As mentioned in the Introduction, this paper focuses on the problem of designing waveform sets with good correlation and stopband properties. Therefore, we first establish two criteria in the frequency domain to measure the correlation and stopband properties and then formulate the waveform set design problem.

The Criterion for Good Correlation Property
where N is the waveform length and M denotes the number of the waveforms. Then, the correlation function of x i and x j is defined as: Here, we consider the design problems of complementary sets of sequences and waveform sets with both good auto-and cross-correlation properties. For the complementary set of sequences, the common criterion is the complementary integrated sidelobe level (CISL) metric [14,25], which is defined as: for the waveform sets with both good auto-and cross-correlation properties, we consider the following more general weighted measure [14]: where w t (k) = w t (−k) ≥ 0, k = 0, ..., N − 1 denote the weights assigned to different time lags. In those specified intervals that are expected to have as low as possible sidelobes, w t (k) = 1, else w t (k) = 0. Let: be the correlation vector and the corresponding correlation weight vector, respectively. The Fourier transform of the correlation function r ij (i.e., power spectrum density (PSD)) can be written as: where F is the 2N × 2N discrete Fourier transform (DFT) matrix with the following expression: ∂φ mn The first term in (19) can be simplified as: where the third equality follows from the fact (Q + Q ) T = (Q + Q ) * . Let: then (20) can be expressed as: According to (A4) and (A5), we can deduce the derivative ∂ (p mm (k)) ∂φ mn as follows: By substituting (23) into (22), we have: To compute the second term of (19), we first deduce ∂ p H mi Re (Q ) p mi ∂φ mn (i = m): Similarly to the derivation of (23), it is easy to obtain that: On the basis of (25) and (26), ∂ p H mi Re (Q ) p mi ∂φ mn can be further denoted as By substituting (24) and (27) into (19), ∂J T ∂φ mn can be simplified as: then ∂J T ∂φ mn in (28) becomes: By stacking (30) in a vector, the phase gradient ∇ φ m J T is given by: where z m (1 : N) denotes the first N elements of z m . It is worth noting that y mm and y mi defined in (21) can be calculated by the Hadamard product and the FFT operations: From the derivation above, it is easy to see that the calculation of the phase gradient is a little bit cumbersome. In order to make the calculation process concise, we write the gradient in the form of matrix. Let X = [x 1 , ..., x M ] and Φ = [φ 1 , ..., φ M ] denote the matrix of waveform set and the corresponding phase matrix, respectively. Then, the spectrum matrix of X is S = F c (X) = [f 1 , ..., f M ]. Define: then the matrix form of (29) is given by: Thus, according to (31), the matrix of phase gradient G can be expressed as: where Z(1 : N, :) denotes the submatrix formed with the first N rows of Z. By defining Y and Y 1 , we can easily calculate the gradient matrix. However, it is still hard to calculate Y 1 directly by (33).
To efficiently calculate Y 1 , define: Then it is easy to verify that: where vec (Y 1 ) denotes the column vector consisting of all of the columns of Y 1 . Thus, we can obtain Y 1 via the inverse operation of vec (·), i.e., where Y can be calculated by (32). It is easy to see that y mi can be implemented by four FFT (IFFT) operations. Since y mi = y * im , the calculation of Y takes 4 M(M+1) 2 = 2M 2 + 2M FFT (IFFT) operations.

Step Size Calculation via Taylor Series Expansion
The traditional methods for obtaining the step size are the linear search methods, which require many iterations and thus are quite time consuming. To reduce the computing expense, here we propose to calculate the step size directly. Assume that Φ l is the phase matrix of the present iteration point X l = x l 1 , ..., x l M = e jΦ l , and D l = d l 1 , ..., d l M is the descent direction. Then, the new iteration point can be denoted as: i.e., where µ is the step size. Thus, the linear search problem can be formulated as the following minimization problem: By taking the derivative of (41), we have: To simplify (41) and (42), we first deduce p l+1 ij . By using the Taylor series expansion, p l+1 ij can be approximated as (see Appendix C): where: Let then (43) can be rewritten as: By replacing p l+1 ij withp l+1 ij , the approximate function of h(µ) in (41) can be denoted as: Since p l+1 ij =p l+1 ij holds when µ = 0, it is easy to verify that: which indicates that h(µ) and h 1 (µ) have the same function value and slope at µ = 0. As D l is the descent direction, the slopes of these two functions are less than zero. Thus, these two functions have at least one minimum point greater than zero. Since h(µ) is sensitive to the waveform phases, the optimal step size of h(µ) is very small and close to zero. Consequently, we can use the minimum point of h 1 (µ) to approximate the optimal step size.
To calculate the minimum point of h 1 (µ), we replace p l+1 ij in (42) withp l+1 ij , then the derivative of h 1 (µ) with respective to µ is given by: where: To simplify the calculation, let: Then (50) can be rewritten as: where the first term of each equality follows from the fact a H Qb = λa are the arbitrary vectors), and the second term of each equality follows from the fact . Let ∂h 1 ∂µ = 0, then the minimum point of h 1 (µ) can be obtained by solving the following cubic equation: It is well known that a cubic equation with real coefficients has three roots, in which there is at least a real root. Therefore, we can choose the positive root that is closest to zero as the precise estimation of the step size.
In order to facilitate the calculation, we write the above derivation in the form of matrix. Define: then according to (44) we have: Then the inverse Fourier transforms of P, C, C can be denoted as: According to (52), (56) and (57), the coefficients of the cubic equation can be reformulated as:

Algorithm Summary
After deducing the phase gradient and the step size, it is easy to solve the unconstrained problem (17) by using the conjugate gradient algorithm (CGA). Here, we apply the classical Polak-Ribiere-Polyak CGA (PRP-CGA) to deal with the waveform design problem. The searching direction of the PRP-CGA can be expressed as: where g l and d l are the gradient vector and the direction vector, respectively. For the problem here, g l and d l are defined as: where g l m = ∇ φ m J T X l , m = 1, ..., M. Since the step size is an approximate value, g l+1 T d l is not equal to zero. Thus, d l+1 may not be descendent, i.e., is not always satisfied. For guaranteeing that the searching direction is descendant, we adopt the following modified direction: Actually, (59) can also be expressed as: Since and (63) can be rewritten as: Thus, (62) can be expressed as the following matrix form: On the basis of the above derivation, the gradient-based algorithm, which we call Gra-WeCorr-SFW, is summarized in Algorithm 1.

Simplified Algorithm for the Design Problem without Considering the Correlation Weights
In Section 3, we present a gradient-based algorithm to handle Problem (17). From (9), we can see that the criterion J CF can be simplified when w t (k) = 1, k = 1 − N, ..., N − 1 (i.e., the correlation weights are not taken into account). Thus, in this section, we derive a simplified algorithm for the design problem without considering the correlation weights. Let w t (k) = 1, k = 1 − N, ..., N − 1, then the criterion J CF can be simplified as: Thus we can rewrite the problem (17) as: where λ 1 = 1−λ 2N . Like the derivation in Section 3, we deduce the gradient and the step size.
Compute coefficients a, b, c, d according to (58). (53), and then choose the positive root that is closest to zero as the step size µ l . (32).

Step Size
Similarly to (41), the linear search problem here can be denoted as: By taking the derivative of h(µ), we have: Then, the approximate derivative can be obtained by replacing p l+1 pq withp l+1 pq defined in (46): According the definition (56), we have: where P d , C d and C d can be expressed as the following equalities according to (45): Thus, (77) can be written more compactly as: where: By solving the cubic equation a 1 µ 3 + b 1 µ 2 + c 1 µ + d 1 = 0, the step size can be easily obtained.

Algorithm Summary
In the previous two subsections, the gradient and step size for Problem (68) are derived. Then, the simplified algorithm (which we call Gra-Corr-SFW) for the design problem without considering the correlation weights is summarized in Algorithm 2. It is easy to observe that both Algorithm 1 and Algorithm 2 can be easily implemented by the (I)FFT operations and the Hadamard product. Since the Hadamard product is more efficient than the (I)FFT operation, we mainly use the number of (I)FFT operations to measure the time complexity of the algorithms. From Section 3, we know that the calculation of Y needs 2M 2 + 2M (I)FFT operations. Thus, the Gra-WeCorr-SFW requires 5M 2 + 6M (I)FFT operations at each iteration. Compared to the Gra-WeCorr-SFW, the Gra-Corr-SFW is simpler and only needs 4M (I)FFT operations at each iteration. It is worth noting that the calculation of the cubic function is very simple and needs only a small amount of computation.

4: Solve the cubic equation (80), and then choose the positive root
which is closest to zero as the step size µ l .
Until convergence Table 1. The per iteration computational complexities of different algorithms.

MM-WeCorr-acc
To analyze the convergence speed, Table 1 presents the per iteration computational complexities of the proposed and existing algorithms. As shown in Table 1, the proposed Gra-WeCorr-SFW requires fewer (I)FFT operations than the MM-WeCorr-acc (MM-WeCorr-acceleration) at each iteration. Compared to these two algorithms, the per iteration computational complexities of the rest of the three algorithms (MM-Corr-acc (MM-Corr-acceleration), MDISAA-SFW and Gra-Corr-SFW), which do not consider the correlation weights, are much smaller. In addition to the per iteration computational complexity, the iteration number is also an important factor affecting the convergence speed. Thus, it is difficult to compare the convergence performance of the algorithms via the per iteration computational complexity. In the following section, several numerical experiments are provided to show the convergence performance of the proposed algorithms.

Numerical Experiments
To illustrate the effectiveness and superiority of the proposed algorithms, several numerical experiments are presented in this section. We first validate the monotonicity of the proposed algorithms and then assess the performance of the algorithms by designing three different waveform sets. The proposed algorithms are compared with the MM-Corr-acc [14], MM-WeCorr-acc [14] and MDISAA-SFW [37] algorithms, where MM-Corr-acc and MM-WeCorr-acc are the state-of-the-art algorithms for designing waveform sets with a good correlation property, and MDISAA-SFW is the state-of-the-art algorithm for designing the orthogonal waveform set (OWS) with the stopband constraint.
All of the experiments are performed on a PC with a 3.60-GHz i7-4790 CPU and 8GB RAM. The software environment is MATLAB 2012b. In the following experiments, all of the algorithms are initialized by the unimodular waveform sets with random phases.

Verification of The Monotonicity
In this subsection, we investigate the monotonicity of the proposed algorithms in three different waveform design problems, which are respectively complementary sets of sequences (CSS) design, waveform set design with zero correlation zone and orthogonal waveform set (OWS) design with the stopband constraint. In order to measure the monotonicity, define the relative error of the l-th iteration as follows: where µ l denotes the approximate step size of the l-th iteration, which is obtained via the Taylor series expansion, andμ l denotes the optimal step size obtained by the searching method. The smaller the value of ε l r is, the closer the approximate step size is to the optimal step size. Thus, ε l r can be used to express the accuracy of the approximate step size. When ε l r < 1, we have h(µ l ) < h(0), which means the objective function is decreasing at the l-th iteration. Therefore, we can measure the monotonicity of the algorithms by using the following peak relative error P re : where N I is the iteration number.  To simulate ε l r and P re , we use the Gra-WeCorr-SFW to design waveform sets with M = 3 waveforms and each waveform of length N = 256. The simulation parameters of different design problems are shown in Table 2. Additionally, we stop the algorithm after 200 iterations. Figure 1 shows the evolution curves of the relative error with respect to the iteration number. From Figure 1, we can see that the relative error is a little larger at the initial iterations. However, it decreases rapidly and is substantially below 10 −4 after several iterations. In the whole iteration process, the relative error is less than one, which means that the objective function is monotonically decreasing. Figure 2 shows the peak relative error of 100 random trials. It is easy to see that the peak relative error of all three cases is very small and less than one. This indicates that the monotonicity of the Gra-WeCorr-SFW is not affected by the initial iteration point. Since the Gra-Corr-SFW is a simplified algorithm of the Gra-WeCorr-SFW, it can also guarantee that the objective function is monotonically decreasing at each iteration.

CSS Design or OWS Design
In Appendix A, we have demonstrated that the CSS design is equivalent to the OWS design. In this subsection, we apply the Gra-Corr-SFW to design CSS (or OWS). Since the lower bound of the CISL is zero, we choose CISL ≤ ε as the stopping criterion for all of the algorithms in this experiment. The weight coefficient λ and the correlation weights are the same as Case 1 in Table 2. Figure 3 shows the normalized autocorrelation sum of the waveform sets designed by the Gra-Corr-SFW, where the normalized autocorrelation sum is defined as: The evolution curves of MM-Corr-acc and Gra-Corr-SFW are respectively shown in Figure 4a,b. From these two subfigures, we observe that the convergence speed of the proposed Gra-Corr-SFW is much faster than that of the MM-Corr-acc. Even for (N = 1280, M = 15), the proposed algorithm takes 1.35 s to coverage to CISL ≤ 10 −13 , while the MM-Corr-acc takes 48.11 s.
Further, to eliminate the randomness, we repeat the algorithms 100 times for each (M, N) pair and record the average iteration number N I , the average running time t and the average peak sidelobe level of the autocorrelation sum P sum , where P sum is defined as: Then, we choose ε = 10 −1 for the stopping criterion. The performance parameters (N I , t, P sum ) of these two algorithms are provided in Table 3. As can be seen from this Table, since the algorithms stop the iteration when CISL ≤ 10 −1 , the P sum of these two algorithms are basically the same. At the same time, the average running time in Table 3 shows that the Gra-Corr-SFW is an order of magnitude faster than the MM-Corr-acc. From Table 1, we can see that the per iteration computational complexities of both the MM-Corr-acc and the Gra-Corr-SFW are O (8MN log 2N). Therefore, the reason for the slower convergence of the MM-Corr-acc may be that the MM strategy in this algorithm makes the objective function loose, so that the iteration number of the MM-Corr-acc is larger than that of the Gra-Corr-SFW (as shown in Table 3). Table 3. The comparison of the performance parameters between MM-Corr-acc and Gra-Corr-SFW.

Waveform Set Design with Zero Correlation Zone
To verify the effectiveness of the Gra-WeCorr-SFW, we consider the problem of suppressing the correlation sidelobes at the specified intervals, i.e., designing the waveform set with the zero correlation zone (ZCZ), and compare the performance with the MM-WeCorr-acc. The experimental parameters are the same as Case 2 in Table 2. For both algorithms, we choose the value of the objective function ψ in (4) as the stopping criterion, i.e., ψ ≤ ε. The auto-and cross-correlations of the waveform sets designed by MM-WeCorr-acc and Gra-WeCorr-SFW are shown in Figure 5. From this figure, we can see that the proposed algorithm can also generate a waveform set with correlation sidelobes that are almost zero at the specified intervals. Similar to Figure 4, for each (M, N) pair, we simulate the evolution curves of the objective function ψ with respect to the running time. Additionally, the results of these two algorithms are presented in Figure 6. It is easy to see that the Gra-WeCorr-SFW is faster than the MM-WeCorr-acc. In addition, we can also find that when the objective function ψ is less than 10 −11 , the convergence speed of the MM-WeCorr-acc decreases, especially for a large (M, N) pair.  Further, Table 4 presents the comparison of performance parameters between MM-WeCorr-acc and Gra-WeCorr-SFW, where P zcz is the peak sidelobe level at the zero correlation zone defined as: Similarly, for each (M, N) pair, the algorithms are repeated 100 times. As can be seen from Table 4, the P zcz of the Gra-WeCorr-SFW is a little bit lower than that of the MM-WeCorr-acc. Moreover, due to fewer iterations, the average running time of the Gra-WeCorr-SFW is about 50% of that of the MM-WeCorr-acc, which indicates that the former is more computationally efficient. Table 4.

OWS Design with the Stopband Constraint
In this subsection, we consider the problem of designing OWS with the stopband constraint. Since the Gra-Corr-SFW is a simplified version of the Gra-WeCorr-SFW and needs much fewer (I)FFT operations, it is more efficient than the Gra-WeCorr-SFW. Thus, we only investigate the performance of the Gra-Corr-SFW and compare it with the MDISAA-SFW. Here, X l+1 − X l F ≤ 10 −4 is employed as the stopping criterion, and the experimental parameters are the same as Case 3 in Table 2. First, we apply these two algorithms to design OWS with the stopband constraint. Suppose the waveform number M is three and the waveform length N is 256. The correlation levels and spectral power of the waveform sets designed by MDISAA-SFW and Gra-Corr-SFW are shown in Figure 7 and Figure 8. From Figure 7, we can see that the waveform set designed by the MDISAA-SFW has better autocorrelation performance, while the waveform set designed by the Gra-Corr-SFW has better cross-correlation performance. This is because the MDISAA-SFW optimizes the autocorrelation explicitly and, thus, places more emphasis on suppressing the autocorrelation sidelobes. Figure 8 indicates that the stopband spectral power of the waveform set designed by the Gra-Corr-SFW is less than that of the waveform set designed by the MDISAA-SFW. In order to fully compare the algorithms, we define the peak autocorrelation sidelobe P ac , peak cross-correlation P cc , peak stopband power P sp and integrated sidelobe level (ISL) [16] as follows: where P i stop denotes the peak stopband power of the waveform x i . In the above four performance parameters, P ac and P cc are used to measure the autocorrelation and cross-correlation performances, respectively; P sp is the parameter related to the stopband performance; and ISL indicates the overall sidelobe performance. Since the decrease of the bandwidth leads to the broadening of the main lobe, r ii (1), i = 1, ..., M can be regarded as the autocorrelation main lobe. Thus, the definition of P ac does not take k = 1 into account. Here, we consider the waveform sets with M = 3 waveforms and each waveform of length N ∈ {256, 512, 1024}. For each (M, N) pair, we repeat these two algorithms 100 times and record the average values of the performance parameters. The comparison of the performance parameters between MDISAA-SFW and Gra-Corr-SFW is shown in Table 5. In addition to the parameters defined in (88), the average iteration number N I and the running time t are provided in Table 5. From this table, we obtain three findings, as follows. Above all, the P cc and P sp of the Gra-Corr-SFW are lower than that of the MDISAA-SFW, which means that the cross-correlation and stopband performances of the Gra-Corr-SFW are better than that of the MDISAA-SFW. However, in terms of autocorrelation performance, the Gra-Corr-SFW is inferior to the MDISAA-SFW. From the ISL in Table 5, it is easy to see that the overall sidelobe performance of the Gra-Corr-SFW is better than that of the MDISAA-SFW. Furthermore, although the per iteration computational complexity of the Gra-Corr-SFW is slightly higher than that of the MDISAA-SFW (see Table 1), the Gra-Corr-SFW needs fewer iterations due to the fast convergence of the gradient algorithm, which makes the proposed algorithm more efficient. Finally, with the increase of the waveform length, both the correlation and stopband performances of these two algorithms are improved, especially the cross-correlation performance. Next, we investigate the performance of the algorithms under different λ. The average values (100 trials) of P ac , P cc and P sp are shown in Figure 9. From this figure, we can see that with the increase of λ, the P ac and P cc of these two algorithms change slowly, while the P sp is sensitive to the change of λ. It can also be observed that except the autocorrelation performance, the cross-correlation and stopband performances of the Gra-Corr-SFW are better than that of the MDISAA-SFW. This means that the stopband property can be easily obtained by using the Gra-Corr-SFW. Moreover, since the correlation performance changes slowly with the increase of λ, we can choose a relatively large λ (e.g., λ ∈ (0.85, 0.95)) for the Gra-Corr-SFW to obtain both good correlation and stopband properties.

Conclusions
In this paper, we propose an efficient algorithm named Gra-WeCorr-SFW for designing the waveform set with a good correlation and stopband properties. The algorithm optimizes the objective function directly and can guarantee that the objective function is monotonically decreasing at each iteration. By changing the design parameters, the Gra-WeCorr-SFW can be used to generate different waveform sets, such as CSS, the waveform set with ZCZ and OWS with the stopband constraint.
As the main steps can be implemented by the FFT operations and the Hadamard product, the proposed algorithm is computationally efficient and can be used to design waveform sets with large M and N. Besides, the simplified version of the Gra-WeCorr-SFW (named Gra-Corr-SFW) is proposed for the design problem without considering the correlation weights. Compared to the Gra-WeCorr-SFW, the simplified algorithm requires fewer FFT operations and is faster. Numerical experiments show that the proposed algorithms are faster than the state-of-the-art algorithms (MM-WeCorr-acc and MM-Corr-acc) when designing CSS or the waveform set with ZCZ. In the case of designing OWS with the stopband constraint, the simplified algorithm has better stopband performance and computational efficiency compared with the MDISAA-SFW.