Perturbation Transmit Beamformer Based Fast Constant Modulus MIMO Radar Waveform Design

Abstract

In this paper, a fast method to generate a constant-modulus (CM) waveform for a multiple-input, multiple-output, (MIMO) radar is proposed. To simplify the optimization process, the design of the transmit waveform is decoupled from the design of transmit beamformers (TBs) and subpulses. To further improve the computational efficiency, the TBs’ optimization is conducted in parallel, and a linear programming model is proposed to match the desired beampattern. Additionally, we incorporate the perturbation vectors into the TBs’ optimization so that the TBs can be adjusted to satisfy the CM constraint. To quickly generate the CM subpulses with the desired range-compression (RC) performance, the classical linear frequency modulation (LFM) signal and non-LFM (NLFM) are adopted as subpulses. Meanwhile, to guarantee the RC performance of the final angular waveform, the selection of LFM signal parameters is analyzed to achieve a low cross-correlation between subpulses. Numerical simulations verify the transmit beampattern performance, RC performance, and computational efficiency of the proposed method.

1. Introduction

In recent years, due to the increasing complexity of electromagnetic environments and the growing demand for detection, how to utilize the waveform-diversity capability of a radar’s transmitting end has attracted much attention [1]. Each transmitting antenna of a multiple-input, multiple-output, (MIMO) radar can transmit different waveforms. Therefore, it is a typical waveform-diversity radar. Compared to a traditional phased-array radar, a MIMO radar has advantages, such as the flexible design of the transmit beampattern, a high spatial resolution, and a strong multiple-target tracking ability [2,3,4,5,6,7,8,9].

A MIMO radar is usually divided into two types, based on the distance between radar antennas, namely a distributed MIMO radar and a colocated MIMO radar. In this paper, we focus on the study of a colocated MIMO radar, which can improve spatial resolution, flexibility in the transmit beampattern design, and multiple-target tracking capability. The transmit waveform is critical for realizing its advantages over a conventional phased-array radar [10]. For a colocated MIMO radar, the waveform design is usually divided into a quasi-orthogonal waveform design [11,12,13,14,15] and a partially correlated waveform design [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. The quasi-orthogonal waveform is defined as waveforms approaching orthogonality at all mutual delays so that the target echoes caused by different waveforms can be separated at the receiving end. Many works have been dedicated to lowering the auto-correlation and cross-correlation levels of transmit waveforms to obtain ideal quasi-orthogonal waveforms [11,12,13,14,15]. The MIMO radar-transmitting quasi-orthogonal waveforms can achieve the highest spatial resolution [37]. However, the cost of this high spatial resolution is an omnidirectionally transmit beampattern and a low transmission gain [35].

In order to meet the detection demand for different scenarios, such as multiple-target tracking [23,24,25,26], interference mitigation [27,28,29], and region-of-interest searches [20,30,36], a compromise between the transmit beampattern and transmission gain should be achieved. Therefore, some work has begun to focus on a partially correlated waveform design. The partially correlated waveform design can be classified into two sub-classes. One approach aims at obtaining the maximization output of the signal-to-interference-plus-noise ratio [17,18,19]. Another approach is to focus on the design of the transmit beampattern and optimize the temporal characteristic of the transmit waveform simultaneously [16,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,36]. The first sub-class not only utilizes the optimization resources of the transmitting end but also uses the optimization resources of the receiving end. However, in these works, prior knowledge of interference is assumed. Additionally, when considering the constant-modulus (CM) constraint that guarantees the radio frequency amplifier’s operating at the maximum efficiency, joint optimization is an intractable problem, and it cannot be solved within the acceptable computation time for a radar. For the second sub-class, some work has only involved the transmit beampattern design and ignored the temporal characteristics of the transmit waveform [28,31,36]. In [31,38], the classical two-step scheme was adopted, which involves optimizing the covariance matrix first and then synthesizing the CM waveform to approximate the optimal covariance matrix. It should be noted that covariance matrix optimization is a semi-definite programming (SDP) problem with a heavy computational burden [39]. The methods proposed in [36] achieve high computational efficiency, but the generated waveform does not meet the CM constraint. In [28], a parallel design strategy was adopted to improve computational efficiency. However, the sidelobe control was ignored.

It is well known that radar detection performance is closely related to the temporal characteristics of waveforms. Therefore, some work has investigated the joint design of waveforms’ spatial and temporal characteristics for a MIMO radar. Because of their good Doppler tolerance and ease of generation, in [20,30], linear frequency modulation (LFM) waveforms were used directly. However, to match a desired beampattern, both the methods in [20,30] involved nonlinear programming to optimize the parameter of LFM, which was solved using a sequential quadratic programming algorithm with high computational complexity [30]. In addition to adopting the LFM signal directly, the similarity constraint was also adopted to optimize the ambiguity function (AF) of the transmit waveform [21,32]. All the above work has focused on the temporal characteristics of each waveform out of each transmitting antenna. Due to the different signal-processing procedures, some work has been dedicated to designing the temporal characteristics of angular waveforms (AWs), which refer to the synthesized waveforms of the transmit waveform in different directions [22,26,33,34]. The work [22,26,34] applied advanced, non-convex optimization algorithms to the joint design of the transmit beampattern and AWs. However, under the CM constraint, the computational burden of these algorithms is heavy, and the convergence speed is slow. Time-diversity technology was utilized in the work [16] so that the CM transmit waveform could be generated directly. However, in [16], the flexibility of a MIMO radar transmit beampattern was not fully explored.

Motivated by the above issues, in this paper, we focus on the fast method to generate CM waveforms, which makes sense in rapidly changing detection environments. Additionally, the joint design of a transmit beampattern and AWs is considered under the proposed signal model. Specifically, the main contributions of this paper can be summarized as follows:

(1) 

Direct CM waveform synthesis based on proposed transmission strategy:

In the proposed transmission strategy, the transmit waveform can be obtained via a product of a transmit beamformer (TB) matrix and a subpulse vector. Therefore, under the proposed transmission strategy, once the TBs and subpulses that satisfy the CM constraints are obtained, the final CM transmit waveform is achieved.

(2) 

Fast and flexible transmit beampattern design based on perturbation TB:

In order to quickly realize the desired transmit beampattern, we formulated the beampattern design as a linear programming (LP) problem and a set of parallel TBs’ optimization. Both the above problems can be solved efficiently. Additionally, the perturbation vectors are introduced in the TBs’ design so that the TBs can be adjusted to meet CM constraints. Combining the LP and TB design, the final beampattern is flexible, and most scenarios are included, such as multiple beams, a wide beam, sidelobe control, and notches’ generation.

(3) 

Fast and low cross-correlation CM subpulses’ generation based on analytical expression:

In order to quickly generate the CM subpulses, we directly adopted the classical radar waveforms as the subpulses, such as the LFM signal and the non-LFM (NLFM) signal. Both of them have analytical expressions, and they are easy to generate. Moreover, both of them have an excellent AF feature. In the proposed signal model, we show that the temporal characteristic of AWs is consistent with that of subpulses. In order to distinguish between the different subpulses, the relation between cross-correlation and LFM parameters was analyzed. Then, a proper parameter setting was proposed to produce low-cross-correlation subpulses.

The rest of the paper is organized as follows. Section 2 introduces the signal model and clarifies the transmission strategy and signal processing scheme. In Section 3, the problem model is formulated, and the proposed method is shown in Section 4. The simulation results are presented in Section 5. Conclusions are drawn in Section 6.

Notation: Standard case letters stand for scalars; bold, uppercase letters denote matrices; and bold, lowercase letters denote vectors. The notations ${\left(•\right)}^{*}$, ${\left(•\right)}^{\mathrm{T}}$, and ${\left(•\right)}^{\mathrm{H}}$ denote the conjugate, the transpose, and the conjugate transpose, respectively. ${\mathbb{C}}^{L\times L}$ denotes the $L\times L$ complex space. ⊙ denotes the Hadamard product. The subscript ${∥•∥}_2$ denotes the Euclidean vector norm. $\left|•\right|$ denotes the modulus operation. For a vector, the operation $\left|•\right|$ is imposed in an element-wise way. $angle(•)$ denotes the phase of each element in a vector. $diag\left(\mathbf{x}\right)$ represents a diagonal matrix with its main diagonal filled with $\mathbf{x}$, and j denotes an imaginary unit.

2. Signal Model

2.1. Transmission Strategy

A colocated MIMO radar with N transmitting antennas in a uniform linear array (ULA) was considered. Let $s_n\left(t\right)$ denote the waveform emitted via the n-th transmitting antenna, and the transmit waveform can be written as

$$\mathbf{S}\left(\begin{array}{c}t\end{array}\right)={\left[\begin{array}{c}{s}_1\left(\begin{array}{c}t\end{array}\right),s_2\left(\begin{array}{c}t\end{array}\right),\cdots ,s_N\left(\begin{array}{c}t\end{array}\right)\end{array}\right]}^{\mathrm{T}}$$

(1)

where $0\le t\le T_p,T_p$ is the pulse duration.

The AW at direction $\theta$ can be expressed as

$$y\left(\begin{array}{c}t,\theta \end{array}\right)={\mathbf{a}}^{\mathrm{T}}\left(\theta \right)\mathbf{S}\left(t\right)$$

(2)

where

$$\mathbf{a}\left(\theta \right)={\left[e^{-j2\pi \left({\displaystyle \frac{N-1}{2}}\right)dsin\left(\theta \right)/\lambda },e^{-j2\pi \left({\displaystyle \frac{N-1}{2}}-1\right)dsin\left(\theta \right)/\lambda },\cdots ,e^{j2\pi \left({\displaystyle \frac{N-1}{2}}\right)dsin\left(\theta \right)/\lambda }\right]}^{\mathrm{T}}/\sqrt{N}$$

denotes the N-dimensional transmission-steering vector. d is the inter-transmitter spacing, and $\lambda$ is the signal wavelength.

The average power received at direction $\theta$ is expressed as

$$P\left(\begin{array}{c}\theta \end{array}\right)={\displaystyle \frac{1}{T_p}}{\int }_0^{T_p}{\mathbf{a}}^{\mathrm{T}}\left(\begin{array}{c}\theta \end{array}\right)\mathbf{S}\left(\begin{array}{c}t\end{array}\right){\mathbf{S}}^{\mathrm{H}}\left(\begin{array}{c}t\end{array}\right){\mathbf{a}}^{*}\left(\begin{array}{c}\theta \end{array}\right)dt$$

(3)

where $P\left(\theta \right)$ is the transmit beampattern.

To simplify the optimization complexity, the transmission strategy in [28] is introduced here. The transmit waveform is resolved into a product of a TB matrix, $\mathbf{X}\in {\mathbb{C}}^{N\times K}$, and a subpulse vector, $\mathbf{f}\left(t\right)\in {\mathbb{C}}^{K\times 1}:$

$$\begin{array}{c}\mathbf{S}\left(t\right)=\mathbf{Xf}\left(t\right)\hfill \\ \hspace{1em}\hspace{1em}=\left[{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_K\right]\left[\begin{array}{c}{f}_1\left(t\right)\\ f_2\left(t-T_1\right)\\ ⋮\\ f_K\left(t-{\displaystyle \sum _{k=1}^{K-1}}{T}_k\right)\end{array}\right]\hfill \end{array}$$

(4)

where ${\mathbf{x}}_k={\left[\begin{array}{cccc}{x}_k\left(\begin{array}{c}1\end{array}\right)& x_k\left(\begin{array}{c}2\end{array}\right)& \cdots & x_k\left(\begin{array}{c}\mathit{N}\end{array}\right)\end{array}\right]}^{\mathrm{T}}$ is an N-dimensional TB, and the subpulse vector is composed of the time shift of subpulses $f_k\left(t\right)$ with support $\left[\begin{array}{c}0,T_k\end{array}\right]$. Obviously, there is ${\displaystyle \sum _{k=1}^K}{T}_k=T_p$.

In this way, the transmit waveform can be constructed as follows:

$$\mathbf{S}\left(\begin{array}{c}t\end{array}\right)={\mathbf{x}}_1{f}_1\left(\begin{array}{c}t\end{array}\right)+\sum _{k=2}^K{\mathbf{x}}_k{f}_k\left(t-\sum _{\tilde{k}=1}^{k-1}{T}_{\tilde{k}}\right)$$

(5)

Similarly, the transmit beampattern can be rewritten as

$$P\left(\theta \right)={\displaystyle \frac{1}{T_p}}\sum _{k=1}^K{P}_k\left(\theta \right){\int }_0^{T_k}|f_k\left(t\right){|}^{2}dt$$

(6)

where

$$P_k\left(\theta \right)={\mathbf{a}}^{\mathrm{T}}\left(\theta \right){\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{H}}{\mathbf{a}}^{*}\left(\theta \right)$$

(7)

represents the beampattern of the k-th TB.

2.2. Signal Processing Scheme

As we know, the signal processing scheme of the receiving end determines which properties of the transmit waveform we need to optimize. Here, the signal processing scheme in [33] is introduced. In this scheme, the spatial receiving channel (SRC) consists of a receive beamformer and a range compressor (RC). The details of each SRC and the introduced signal processing scheme are given in Figure 1.

In Figure 1, $\mathit{M}$ is the number of receiving antennas, and $\mathbf{w}\left(\begin{array}{c}\theta \end{array}\right)={\left[\begin{array}{c}{w}_1\left(\begin{array}{c}\theta \end{array}\right),w_2\left(\begin{array}{c}\theta \end{array}\right),\cdots ,w_M\left(\begin{array}{c}\theta \end{array}\right)\end{array}\right]}^{\mathrm{T}}$ denotes the receiving beamformer for the direction $\theta$. Review (2) and (5); under the proposed transmission strategy, the AW at direction $\theta$ can be rewritten as

$$\begin{array}{cc}\hfill y(t,\theta )& \hfill ={\mathbf{a}}^{\mathrm{T}}\left(\theta \right)\left({\mathbf{x}}_1{f}_1\left(t\right)+\sum _{k=2}^K{\mathbf{x}}_k{f}_k\left(t-\sum _{\tilde{k}=1}^{k-1}{T}_{\tilde{k}}\right)\right)\hfill \\ \hfill & \hfill ={\alpha }_1\left(\begin{array}{c}\theta \end{array}\right)f_1\left(\begin{array}{c}t\end{array}\right)+\sum _{k=2}^K{\alpha }_k\left(\begin{array}{c}\theta \end{array}\right)f_k\left(t-\sum _{\tilde{k}=1}^{k-1}{T}_{\tilde{k}}\right)\hfill \end{array}$$

(8)

where ${\alpha }_k\left(\begin{array}{c}\theta \end{array}\right)={\mathbf{a}}^{\mathrm{T}}\left(\begin{array}{c}\theta \end{array}\right){\mathbf{x}}_k$. Then, at the receiving array, we have

$$\tilde{\mathbf{y}}(t,\theta )=\tilde{\mathbf{a}}\left(\theta \right)y(t,\theta )$$

(9)

where

$$\tilde{\mathbf{a}}\left(\theta \right)={\left[e^{-j2\pi \left({\displaystyle \frac{M-1}{2}}\right)dsin\left(\theta \right)/\lambda },e^{-j2\pi \left({\displaystyle \frac{M-1}{2}}-1\right)dsin\left(\theta \right)/\lambda },\cdots ,e^{j2\pi \left({\displaystyle \frac{M-1}{2}}\right)dsin\left(\theta \right)/\lambda }\right]}^{\mathrm{T}}/\sqrt{M}$$

denotes the M-dimensional receiving steering vector.

The output of the SRC for the direction ${\theta }_0$ can be given by

$$\begin{array}{cc}\hfill {\tilde{y}}_{rb}\left(t,{\theta }_0\right)& \hfill ={\mathbf{w}}^{\mathrm{H}}\left({\theta }_0\right)\tilde{\mathbf{y}}\left(t,\theta \right)\hfill \\ & \hfill =\tilde{\alpha }\left({\theta }_0,\theta \right)y\left(t,\theta \right)\hfill \end{array}$$

(10)

where $\tilde{\alpha }\left({\theta }_0,\theta \right)={\mathbf{w}}^{\mathrm{H}}\left({\theta }_0\right)\tilde{\mathbf{a}}\left(\theta \right)$. Let ${\theta }_N$ denote the natural first null beamwidth after the received beamforming. It can be noticed that, when $\theta \notin \left(\begin{array}{c}{\theta }_0-{\theta }_N/2,{\theta }_0+{\theta }_N/2\end{array}\right)$, we have $\left|\tilde{\alpha }({\theta }_0,\theta )\right|\ll 1$. In the mainlobe region, the signals have similar temporal characteristics [40]. Therefore, after the received beamforming, only the signals from direction ${\theta }_0$ are considered [26], and the output of the SRC for the direction ${\theta }_0$ is approximated as

$${\tilde{y}}_{rb}\left(\begin{array}{c}t,{\theta }_0\end{array}\right)=\tilde{\alpha }\left(\begin{array}{c}{\theta }_0,{\theta }_0\end{array}\right)y\left(\begin{array}{c}t,{\theta }_0\end{array}\right)$$

(11)

It can be seen that after, the received beamforming, the signal’s temporal characteristic entirely depends on the AW at the beam-pointing direction. If ${\alpha }_{\overline{k}}\left({\theta }_0\right)/{\alpha }_k\left({\theta }_0\right)\gg 1,$ $k=1$, $\cdots ,K,k\ne \overline{k}$, the AW in (8) can be approximated as

$$y\left(t,{\theta }_0\right)={\alpha }_{\overline{k}}\left({\theta }_0\right)f_{\overline{k}}\left(t-\sum _{\tilde{k}=1}^{\overline{k}-1}{T}_{\tilde{k}}\right)$$

(12)

In this situation, the subpulses $f_{\overline{k}}\left(\begin{array}{c}t\end{array}\right)$, which were carried via beamformer ${\mathbf{x}}_k$, decide the characteristics of AW at direction ${\theta }_0$, and the corresponding RC can adopt $f_{\overline{k}}\left(t\right)$.

If ${\alpha }_{\overline{k}}\left({\theta }_0\right)/{\alpha }_{\widehat{k}}\left({\theta }_0\right)\approx 1,\overline{k}\ne \widehat{k}$, the AW in (8) can be approximated as

$$y\left(t,{\theta }_0\right)={\alpha }_{\overline{k}}\left({\theta }_0\right)f_{\overline{k}}\left(t-\sum _{\tilde{k}=1}^{\overline{k}-1}{T}_{\tilde{k}}\right)+{\alpha }_{\widehat{k}}\left({\theta }_0\right)f_{\widehat{k}}\left(t-\sum _{\tilde{k}=1}^{\widehat{k}-1}{T}_{\tilde{k}}\right)$$

(13)

In this situation, the cross-correlation between $f_{\overline{k}}\left(\begin{array}{c}t\end{array}\right)$  and $f_{\widehat{k}}\left(t\right)$ should be designed so that the subpulses carried via different TBs can be separated. In this way, the RC can use $f_{\overline{k}}\left(\begin{array}{c}t\end{array}\right)$ or $f_{\widehat{k}}\left(t\right)$, and the characteristics of AW in direction ${\theta }_0$ are decided via  $f_{\overline{k}}\left(\begin{array}{c}t\end{array}\right)$   or  $f_{\widehat{k}}\left(t\right)$.

3. Problem Formulation

According to (6), it can be noticed that the spatial properties of a transmit waveform are decided according to the beampattern of TBs and the corresponding subpulses’ energy. Moreover, from (12) and (13), we can find that the temporal properties of a transmit waveform are determined according to subpulses. Through the above analysis, the transmit waveform design is divided into two parts: the transmit beampattern design and the subpulses’ correlation design.

3.1. The Transmit Beampattern Design

Let $P_d\left(\theta \right)$ denote the desired transmit beampattern, and let ${\beta }_k={\int }_0^{T_k}{\left|f_k\left(t\right)\right|}^{2}dt/T_p$ denote the proportion of the power of the k-th subpulse. When considering the CM constraint, without losing generality, it is supposed that $\left|\begin{array}{c}{f}_k\left(\begin{array}{c}t\end{array}\right)\end{array}\right|=1,t\in \left[\begin{array}{c}0,T_k\end{array}\right]$, $k=1$, $\cdots ,K$. Then, it is natural that ${\displaystyle \sum _{k=1}^K}{\beta }_k=1$, and the transmit beampattern (6) can be rewritten as

$$P\left(\begin{array}{c}\theta \end{array}\right)=\sum _{k=1}^K{P}_k\left(\begin{array}{c}\theta \end{array}\right){\beta }_k$$

(14)

In summary, the transmit beampattern design can be formulated as the following optimization problem:

$$\begin{array}{cc}\hfill & \hfill \underset{{\beta }_k,P_k\left(\theta \right)}{min}\sum _{l=1}^L{u}_l\left|P_d\left({\theta }_l\right)-\sum _{k=1}^K{\beta }_k{P}_k\left({\theta }_l\right)\right|\hfill \\ \hfill & \hfill s.t.{\beta }_k\ge 0,k=1,...,K\hfill \\ \hfill & \hfill \sum _{k=1}^K{\beta }_k=1\hfill \\ \hfill & \hfill P_k\left(\theta \right)={\mathbf{a}}^{\mathrm{T}}\left(\theta \right){\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{H}}{\mathbf{a}}^{*}\left(\theta \right)\hfill \end{array}$$

(15)

where L is the spatial grid number, $u_l$ is the weight for the l-th grid point, and $P_d\left(\theta \right)$ is the desired transmit beampattern. The objective function of problem (15) denotes the absolute error between the optimized transmit beampattern and the desired transmit beampattern. The constraints in problem (15) can be concluded according to the above derivation about TB beampatterns. In problem (15), the CM constraint is ignored, and it will be considered in the following solution process. Nonetheless, the problem (15) is still a quartic programming problem, and it is difficult to solve efficiently. In the following subsection, a simplified method is proposed to solve the problem (15).

3.2. The Subpulses Correlation Design

There are two parts that are considered in the subpulses’ correlation design, which are auto-correlation and cross-correlation. Auto-correlation is related to RC performance, and cross-correlation determines the isolation degree between different subpulses. Therefore, subpulse correlation design can be specifically considered according to the radar’s working mode. When the radar is in tracking mode, the different TBs are designed for targets at different directions, and we have ${\alpha }_{\overline{k}}\left(\theta \right)/{\alpha }_k\left(\theta \right)\gg 1,k=1,\cdots ,K,k\ne \overline{k}$. In this situation, according to (11), only the auto-correlation of subpulses can be considered. The efficiency of waveform generation was considered, and the classical radar signals with analytical expressions, such as LFM and NLFM, could be adopted here.

When a radar is in search mode, the scanning beams are close to each other, and the correlation between TBs cannot be ignored. In this situation, in addition to auto-correlation, the cross-correlation properties between subpulses should be considered so that the subpulses carried via different TBs can be separated. Accordingly, the following problem is proposed to optimize the cross-correlation between subpulses:

$$\begin{array}{l}\underset{\begin{array}{c}{f}_{\overline{k}}\left(t\right),f_{\widehat{k}}\left(t\right)\\ \overline{k}=1,\cdots K\\ \widehat{k}=1,\cdots K\\ \overline{k}\ne \widehat{k}\end{array}}{min}\left|{\int }_{-\infty }^{\infty }{f}_{\overline{k}}\left(t\right)f_{\widehat{k}}^{*}\left(t+\tau \right)dt\right|\\ s.t.\left|f_{\overline{k}}\left(t\right)\right|=1,t\in \left[0,T_{\overline{k}}\right],\overline{k}=1,\cdots K\\ \left|f_{\widehat{k}}\left(t\right)\right|=1,t\in \left[0,T_{\widehat{k}}\right],\widehat{k}=1,\cdots K,\overline{k}\ne \widehat{k}\end{array}$$

(16)

Some existing work has proposed methods to solve the above problem [11,12,13,14,15]. However, because the above problem is NP-hard, most of these methods cannot obtain a solution within the accepted computation time for a radar. For example, the proposed method in [15] involved singular value decomposition, and the computation complexity was around $O\left(I{N}^3\right)$. In [14], though the majorization–minimization algorithm was used, the computation complexity was around $O\left(I{N}^2\right)$. In the above, I represents the number of iterations. In the following section, we analyze the cross-correlation between different-parameter LFM signals, and we propose a method to obtain LFM signals with a low cross-correlation level and without complicated computations.

4. Problem Analysis and Proposed Method

In Section 3, the transmit waveform design is divided into the transmit beampattern design and the subpulses’ correlation design. In this section, problem (15) and problem (16) are analyzed successively, and the corresponding waveform design methods are proposed.

4.1. Analysis of Transmit Beampattern Design Problem (15)

Problem (15) is a quartic programming problem about the power proportion of subpulses and the TB beampattern. It is hard to get a satisfactory result directly with an accepted computation cost. However, it can be noticed that, when $P_k\left(\theta \right),k=1,\cdots ,K$ are fixed, problem (15) can be simplified as an LP problem concerning ${\beta }_k,k=1,\cdots ,K$:

$$\begin{array}{l}\underset{{\beta }_k,{\eta }_l}{min}\sum _{l=1}^L{u}_l{\eta }_l\\ s.t.{\beta }_k\ge 0,k=1,\cdots ,K\\ \sum _{k=1}^K{\beta }_k=1\\ P_d\left({\theta }_l\right)-\sum _{k=1}^K{\beta }_k{P}_k\left({\theta }_l\right)\le {\eta }_l\\ \sum _{k=1}^K{\beta }_k{P}_k\left({\theta }_l\right)-P_d\left({\theta }_l\right)\le {\eta }_l\end{array}$$

(17)

where ${\eta }_l,l=1,\cdots ,L$ are auxiliary variables. Problem (17) can be solved efficiently via some classical algorithms, e.g., the simplex method and dual simplex method. Based on the above analysis, an alternating optimization algorithm was employed for problem (15); in particular, we optimized ${\beta }_k,k=1,\cdots ,K$ for fixed TB beampatterns first, and then the TBs were optimized for the fixed ${\beta }_k,k=1,\cdots ,K$. For clarity, the proposed strategy is summarized in Algorithm 1.

Algorithm 1 Proposed strategy for transmit beampattern design.

Input: Initial TBs ${\mathbf{x}}_k^{\left(0\right)},k=1,\cdots ,K,\ell =0$ Output: The power proportion of subpulses ${\beta }_k^{(\ell )},k=1,\cdots K$ and TBs ${\mathbf{x}}_k^{(\ell )},k=1,\cdots ,K$ 1: For fixed TBs ${\mathbf{x}}_k^{(\ell )},k=1,\cdots ,K$, solve problem (17) using dual simplex method, and obtain the optimized ${\beta }_k^{(\ell )},k=1,\cdots K$; 2: According to (14), calculate the transmit beampattern $P^{(\ell )}\left(\theta \right)$, and obtain ${\theta }_{max}=\underset{{\theta }_l}{max}\left\{\left|P^{(\ell )}\left({\theta }_l\right)-P_d\left({\theta }_l\right)\right|,l=1,\cdots ,L\right\}$; 3: Let $\ell =\ell +1$; 4: Optimize ${\mathbf{x}}_k^{(\ell )},k=1,\cdots ,K$ to satisfy $P_k^{(\ell )}\left({\theta }_{max}\right)=P_d\left({\theta }_{max}\right),k=1,\cdots ,K$; 5: Repeat steps 1–4 until some termination conditions are satisfied.

4.2. Perturbation TB-Based Transmit Beampattern Design

It can be seen that the remaining problem of Algorithm 1 is how to optimize ${\mathbf{x}}_k,k=1,\cdots ,K$ to satisfy $P_k\left({\theta }_{max}\right)=P_d\left({\theta }_{max}\right),k=1,\cdots ,K$. Additionally, another key issue is that the TBs should satisfy the following CM constraint:

$$\left|\begin{array}{c}{\mathbf{x}}_k\end{array}\right|=\mathsf{\xi },k=1,\cdots ,K$$

(18)

where $\mathsf{\xi }$ denotes a constant vector, and all adjusted TB modulus values will be equal to it. Combining it with $\left|\begin{array}{c}{f}_k\left(\begin{array}{c}t\end{array}\right)\end{array}\right|=1,t\in \left[\begin{array}{c}0,T_k\end{array}\right],k=1,\cdots ,K$, we actually require the transmit waveform of each transmitting antenna to meet the CM constraint.

To satisfy Equation (18), the idea of perturbation TBs is proposed. In array signal processing, the array error cannot be ignored. Therefore, many works have proposed a robust beamforming algorithm under an array error case [41,42,43,44,45,46,47]. Almost all of the above works described the array error as a norm-bounded uncertainty set, and they proposed methods such as the worst-case optimization technique, convex optimization algorithm, and closed-form solution iteration to solve this problem. In this paper, we do not consider the array error, but we regard this error as an artificial perturbation to force the TB satisfying the CM constraint. Suppose that there is a set of perturbation vectors, ${\mathsf{\lambda }}_k,k=1,\cdots ,K$ and ${∥{\mathsf{\lambda }}_k∥}_2=\epsilon {∥{\mathbf{x}}_k∥}_2,k=1,\cdots ,K,\epsilon \ll 1$. The key of the proposed idea is adding different perturbation vectors for each TB to satisfy the following equation:

$$\left|{\mathbf{x}}_k+{\mathsf{\lambda }}_k\right|=\mathsf{\xi },k=1,\cdots ,K$$

(19)

Let $V_{kr}\left(\begin{array}{c}\theta \end{array}\right)$ and $V_{ku}\left(\begin{array}{c}\theta \end{array}\right)$ stand for the normalized transmit magnitude response and its upper boundary, respectively, and they can be expressed as

$$V_{kr}\left(\theta \right)=\left|{\left({\mathbf{x}}_k+{\mathsf{\lambda }}_k\right)}^{\mathrm{H}}{\mathbf{a}}^{*}\left(\theta \right)\right|/\left|{\left({\mathbf{x}}_k+{\mathsf{\lambda }}_k\right)}^{\mathrm{H}}{\mathbf{a}}^{*}\left({\theta }_k\right)\right|$$

(20)

$$V_{ku}\left(\theta \right)={\displaystyle \frac{\left|{\mathbf{x}}_k^{\mathrm{H}}{\mathbf{a}}^{*}\left(\theta \right)\right|/\left|{\mathbf{x}}_k^{\mathrm{H}}{\mathbf{a}}^{*}\left({\theta }_k\right)\right|+\epsilon {∥{\mathbf{x}}_k∥}_2/\left|{\mathbf{x}}_k^{\mathrm{H}}{\mathbf{a}}^{*}\left({\theta }_k\right)\right|}{1-\epsilon {∥{\mathbf{x}}_k∥}_2/\left|{\mathbf{x}}_k^{\mathrm{H}}{\mathbf{a}}^{*}\left({\theta }_k\right)\right|}}$$

(21)

where ${\theta }_k$ denotes the mainlobe direction of ${\mathbf{x}}_k$ and the derivation of the upper boundary; see Appendix A.

Review Algorithm 1; after perturbation, when ${\theta }_{max}$ is in the sidelobe region of $P_d\left(\begin{array}{c}\theta \end{array}\right)$, $P_k\left({\theta }_{max}\right)=P_d\left({\theta }_{max}\right),k=1,\cdots ,K$ is equivalent to

$$V_{kr}^2\left({\theta }_{max}\right)=P_d\left({\theta }_{max}\right),k=1,\cdots ,K$$

(22)

where it is assumed that $P_d\left(\begin{array}{c}\theta \end{array}\right)$ and $P_k\left(\theta \right)$ have been normalized. For the sidelobe region, (22) can be relaxed as

$$V_{ku}^2\left({\theta }_{max}\right)=P_d\left({\theta }_{max}\right),k=1,\cdots ,K$$

(23)

It is assumed that the TBs are designed to generate a “fan” of beams jointly covering the target region [28], as shown in Figure 2. Based on this assumption, for the mainlobe region of the desired beampattern, the beampattern-matching problem is converted to maximize the mainlobe gain of each TB and allocate the power proportion for them to approach the desired mainlobe shape. This simplification cannot guarantee the optimal solution for problem (17), but it can decouple the optimization between TBs and improve the computational efficiency. After perturbation, the gain of the TBs can be written as

$${\gamma }_k={\displaystyle \frac{{\left|{\left({\mathbf{x}}_k+{\mathsf{\lambda }}_k\right)}^{\mathrm{H}}{\mathbf{a}}^{*}\left({\theta }_k\right)\right|}^2}{{∥{\mathbf{x}}_k+{\mathsf{\lambda }}_k∥}_2^2}},k=1,\cdots ,K$$

(24)

Using the triangle inequality, we have the lower boundary of the perturbation TB gain:

$${\gamma }_k\ge {\displaystyle \frac{1}{{\left(1+\epsilon \right)}^2}}{\left({\displaystyle \frac{\left|{\mathbf{x}}_k^{\mathrm{H}}{\mathbf{a}}^{*}\left({\theta }_k\right)\right|}{{∥{\mathbf{x}}_k∥}_2}}-\epsilon \right)}^2,k=1,\cdots ,K$$

(25)

When ${\theta }_{max}$ is in the mainlobe region of $P_d\left(\begin{array}{c}\theta \end{array}\right)$, we only focus on the transmission gain of each perturbation TB. According to inequality (25), the problem for the k-th TB can be relaxed as

$$\underset{{\mathbf{x}}_k}{max}{\displaystyle \frac{\left|{\mathbf{x}}_k^{\mathrm{H}}{\mathbf{a}}^{*}\left({\theta }_k\right)\right|}{{∥{\mathbf{x}}_k∥}_2}}$$

(26)

In summary, for the k-th TB, the following problem is proposed:

$$\begin{array}{l}\underset{{\mathbf{x}}_k}{max}{\displaystyle \frac{\left|{\mathbf{x}}_k^{\mathrm{H}}{\mathbf{a}}^{*}\left({\theta }_k\right)\right|}{{∥{\mathbf{x}}_k∥}_2}}\\ s.t.V_{ku}\left({\theta }_{max}\right)=\sqrt{P_d\left({\theta }_{max}\right)}\end{array}$$

(27)

Problem (27) can be solved by using $C^2-WORD$ efficiently [47]. Compared with the traditional adaptive beamforming algorithms, which need covariance matrix inversion, the computational complexity of $C^2-WORD$ to obtain a closed-form solution is only $O\left(N^2\right)$.

After all TBs are obtained, the CM constraint is considered. Suppose that ${\tilde{\mathbf{x}}}_k$, $k=1$, $\cdots ,K$ is the solution to (27), and we have

$${\overline{\mathbf{x}}}_c={\displaystyle \frac{{\displaystyle \sum _{k=1}^K}\left|{\tilde{\mathbf{x}}}_k\right|}{K}}$$

(28)

Then, the final TBs are set as ${\mathbf{x}}_k={\overline{\mathbf{x}}}_c\odot angle\left({\tilde{\mathbf{x}}}_k\right),k=1,\cdots ,K$. It can be proven that the modulus difference in solutions obtained from $C^2-WORD$ is small. The details are seen in Appendix B. Finally, the algorithm concerning perturbation TBs’ design is summarized in Algorithm 2.

After the above analysis, the proposed strategy in Algorithm 1 can be refined into the above Algorithm 3, which is called the perturbation TB-based transmit beampattern design (PTB-TBD) here.

Algorithm 2 Perturbation vector-based CM TBs’ design.

Input: Initial TBs ${\mathbf{x}}_k^{(\ell )},k=1,\cdots ,K,$ the mainlobe directions of TBs ${\theta }_k,k=1,\cdots ,K$, the maximum error direction ${\theta }_{max}$, the desired beampattern $P_d\left(\begin{array}{c}\theta \end{array}\right)$, and the perturbation parameter $\epsilon$; Output: Optimized TBs ${\mathbf{x}}_k^{(\ell +1)},k=1,\cdots ,K$; 1: In parallel, solve problem (27), and obtain the solution, ${\tilde{\mathbf{x}}}_k,k=1,\cdots ,K$; 2: Obtain ${\overline{\mathbf{x}}}_c$ via (28); 3: Calculate ${\mathbf{x}}_k^{(\ell +1)}=angle\left({\tilde{\mathbf{x}}}_k\right)\odot {\overline{\mathbf{x}}}_c$; 4: Calculate the distance $d_k=\left|\right|{\overline{\mathbf{x}}}_c-|{\tilde{\mathbf{x}}}_k{\left|\right||}_2,k=1,\cdots ,K$; 5: If $d_k\le \epsilon {\parallel {\mathbf{x}}_k\parallel }_2,k=1,\cdots ,K$, output $d_k={∥{\overline{\mathbf{x}}}_c-\left|{\tilde{\mathbf{x}}}_k\right|∥}_2,k=1,\cdots ,K$; otherwise, increase $\epsilon$, and repeat above steps.

Algorithm 3 PTB-TBD.

Input: Initial TBs ${\mathbf{x}}_k^{\left(0\right)},k=1,\cdots ,K,\ell =0$, and ${\theta }_{\kappa },k=1,\cdots ,K$; Output: The power proportion of subpulses ${\beta }_k^{(\ell )},k=1,\cdots K$ and TBs ${\mathbf{x}}_k^{(\ell )},k=1,\cdots ,K$; 1: For fixed TBs ${\mathbf{x}}_k^{(\ell )},k=1,\cdots ,K$, solve problem (17) using the dual simplex method, and obtain the optimized ${\beta }_k^{(\ell )},k=1,\cdots K$; 2: According to (14), calculate the transmit beampattern $P^{(\ell )}\left(\theta \right)$, and obtain ${\theta }_{max}=\underset{{\theta }_l}{max}\left\{P^{(\ell )}\left({\theta }_l\right)-P_d\left({\theta }_l\right),{\theta }_l\in {\Omega }_s\right\}$, where ${\Omega }_s$ denotes the sidelobe region; 3: Let $\ell =\ell +1$; 4: Optimize ${\mathbf{x}}_k^{(\ell )},k=1,\cdots ,K$ using Algorithm 2; 5: Repeat steps 1–4 until some termination conditions are satisfied.

4.3. Subpulses’ Cross-Correlation Design and CM Waveform Synthesis

As described in Section 3.2, the cross-correlation between subpulses should be optimized when the radar is in search mode. However, problem (16) is difficult to obtain a solution for within the accepted computation time. Given that the strong Doppler tolerance signal is more proper in search mode, two LFM signals with different parameters are analyzed in this section.

The LFM signals with different carrier frequencies are written as

$$s_i\left(t\right)=exp\left(j2\pi f_{i}t+j\pi \mu t^2\right),\left|t\right|\le {\displaystyle \frac{T_k}{2}},i=1,2$$

(29)

where $\mu$ denotes the chirp rate. Without a loss of generality, it is assumed that the LFM signals have an equal pulse width.

When $0\le \tau \le T_k$, the cross-correlation of the above signals is

$$\begin{array}{l}{R}_{12}\left(\tau \right)={\int }_{-\infty }^{\infty }{s}_1\left(t\right)s_2^{*}\left(t+\tau \right)dt\\ \hspace{1em}\hspace{1em}\hspace{1em}={\int }_{-{\displaystyle \frac{T_k}{2}}}^{{\displaystyle \frac{T_k}{2}}-\tau }exp\left(j2\pi f_{1}t+j\pi \mu t^2\right)exp\left(-j2\pi f_2(t-\tau )-j\pi \mu {(t-\tau )}^2\right)dt\\ \hspace{1em}\hspace{1em}\hspace{1em}=exp\left(3j\pi f_2\tau -2j\pi \mu {\tau }^2-j\pi f_1\tau \right)\times sinc\left((f_1-f_2+\mu \tau )(T_k-\tau )\right)\times (T_k-\tau )\end{array}$$

(30)

When $-T_k\le \tau \le 0$, the cross-correlation of the above signals is

$$\begin{array}{l}{R}_{12}\left(\tau \right)={\int }_{-\infty }^{\infty }{s}_1\left(t\right)s_2^{*}\left(t+\tau \right)dt\\ \hspace{1em}\hspace{1em}\hspace{1em}={\int }_{-{\displaystyle \frac{T_k}{2}}-\tau }^{{\displaystyle \frac{T_k}{2}}}exp\left(j2\pi f_{1}t+j\pi \mu t^2\right)exp\left(-j2\pi f_2(t-\tau )-j\pi \mu {(t-\tau )}^2\right)dt\\ \hspace{1em}\hspace{1em}\hspace{1em}=exp\left(3j\pi f_2\tau -2j\pi \mu {\tau }^2-j\pi f_1\tau \right)\times sinc\left((f_1-f_2+\mu \tau )(T_k+\tau )\right)\times (T_k+\tau )\end{array}$$

(31)

In summary, the cross-correlation between the above LFM signals can be expressed as

$$\left|R_{12}\left(\tau \right)\right|=\left|sinc\left(\left(f_1-f_2+\mu \tau \right)\left(T_k-\left|\tau \right|\right)\right)\times \left(T_k-\left|\tau \right|\right)\right|,-T_k\le \tau \le T_k$$

(32)

From (32), it can be seen that $\left|R_{12}\left(\tau \right)\right|$ decreases as $\left|\begin{array}{c}{f}_1-f_2\end{array}\right|$ increases. Based on this conclusion, we can generate low cross-correlation LFM signals by allocating different carrier frequencies. In this way, though the optimal solution of (16) is not obtained, the expected cross-correlation level can still be obtained by expanding the difference between carrier frequencies.

Finally, the overall fast CM waveform synthesis process is summarized in the flowchart Figure 3. After a review of the entire CM waveform synthesis process, it can be seen that the key to achieving fast performance through the proposed method lies in the following four points:

• Matching the transmit beampattern design via LP solving;
• Obtaining TBs via a closed-form solution without matrix inversion;
• The parallel optimization of TBs and the CM operation based on perturbation vectors;
• Generating low-cross-correlation CM subpulses based on analytical expressions.

5. Simulation Results

This section consists of three parts: transmit beampattern performance, RC performance, and a comparison with existing methods. In the transmit beampattern part, different desired transmit beampatterns are adopted, and the effectiveness of the proposed LP is shown. In the RC performance part, the auto- and cross-correlation of the designed subpulses are shown, and the final RC performance is presented. Finally, in the comparison part, we compare the waveform performance and computation time of the proposed method with those of other existing methods, and the numerical results verify the fast generation performance of the proposed method.

In all the simulations, a colocated MIMO ULA with $N=16$ transmitting antennas, separated by a half-wavelength, was considered. For Algorithm 3, the termination condition was set such that the iteration number reached 500 or $\forall \theta \in {\Omega }_s,P\left(\theta \right)-P_d\left(\theta \right)\le 0$, and the weights for the spatial grid point were set to 1. Unless otherwise specified, the simulation parameters remained unchanged. All the simulations were performed in the MATLAB 2017b/Windows 10 environment on a computer with a 2.1 GHz CPU and 32 GB of RAM.

5.1. Transmit Beampattern Design

In this subsection, we present the results of numerical simulations that were conducted to evaluate the proposed algorithm for the transmit beampattern design. The search mode was considered first, and the desired beampattern is expressed as

$$P_d\left(\theta \right)=\left\{\begin{array}{cc}\hfill 1,\hfill & \hfill \theta \in \left[-{30}^{\circ },{30}^{\circ }\right]\hfill \\ \hfill 0.01,\hfill & \hfill \theta \notin \left[-{30}^{\circ },{30}^{\circ }\right]\hfill \end{array}\right.$$

(33)

where a wide beam with a −20 dB sidelobe level was adopted. In order to avoid missing the effective TBs, the spacing between ${\theta }_k$ was set to half of the 3-dB beamwidth. The specific value for each ${\theta }_k$ is listed in

Table 1. The mainlobe directions of TBs and the power proportion of the corresponding subpulse.

	${\mathit{x}}_1$	${\mathit{x}}_2$	${\mathit{x}}_3$	${\mathit{x}}_4$	${\mathit{x}}_5$	${\mathit{x}}_6$	${\mathit{x}}_7$	${\mathit{x}}_8$	${\mathit{x}}_9$
${\theta }_k$	$-30°$	$-26.3°$	$-22.5°$	$-18.8°$	$-15°$	$-11.3°$	$-7.5°$	$-3.8°$	0°
${\beta }_k$	0	0.1488	0	0	0.1194	0.0350	0.1038	0.0455	0.0981
	${\mathit{x}}_{\mathbf{10}}$	${\mathit{x}}_{\mathbf{11}}$	${\mathit{x}}_{\mathbf{12}}$	${\mathit{x}}_{\mathbf{13}}$	${\mathit{x}}_{\mathbf{14}}$	${\mathit{x}}_{\mathbf{15}}$	${\mathit{x}}_{\mathbf{16}}$	${\mathit{x}}_{\mathbf{17}}$
${\theta }_k$	3.8°	7.5°	11.3°	15°	18.8°	22.5°	26.3°	30°
${\beta }_k$	0.0443	0.1038	0.0342	0.1206	0	0	0.1464	0

. In Figure 4, the initial TBs’ beampattern and the desired beampattern are shown. It can be seen that the mainlobe region is covered with the initial transmit beams. Let the perturbation parameter $\epsilon =0.03$, and let us conduct the proposed PTB-TBD. In Table 1, we represent the final power proportion of each subpulse. It can be found that some TBs were not used. Figure 5 shows the optimized TBs beampattern and the final transmit beampattern, respectively. It can be seen that the beams of TBs were concentrated in the mainlobe region, and the sidelobe level was lower than −20 dB. After power proportion allocation, we can see that the final transmit beampattern could approach the desired mainlobe shape, and the sidelobe satisfied the expected level. To verify the influence of the perturbation parameter $\epsilon$, the final transmit beampatterns under different $\epsilon$ are shown in Figure 6. It can be seen that the impact of the perturbation parameter on the transmit beampattern performance was minimal. This result may imply that the optimal CM solution is included in the uncertain set $\epsilon \le 0.03$. Finally, Figure 7 depicts the amplitude of each TB obtained using the proposed method. It can be found that, for each transmitting antenna, the TBs satisfied the CM constraint.

Next, the multiple-target tracking scene was considered. The desired beampattern was obtained by adding three traditional single beams whose beam directions were $-{35}^{\circ }$, $0^{\circ }$, and ${35}^{\circ }$, respectively, and the desired sidelobe level was set to −16 dB. Because the target directions are known in the tracking mode, the ${\theta }_k,k=1,2,3$ were set to $-{35}^{\circ }$, $0^{\circ }$ and ${35}^{\circ }$, respectively. The perturbation parameter remained unchanged. Figure 8 shows the transmit beampattern performance. We can see that the performance was similar to that in Figure 5. In Figure 9, we show the influence of the perturbation parameter on the transmit beampattern with multiple beams. It can be seen that the result is similar to that in search mode. In Figure 10, the iterative process is shown. It can be found that the direction with the maximum sidelobe level was optimized in each iteration. Additionally, the sidelobe level gradually decreased as the iteration progressed.

Finally, the transmit beampattern design in the interference scenario was considered. It makes sense to form a transmit beampattern with notches at the interference directions [48]. Therefore, an interference region was added to the desired beampattern mentioned above, respectively. For the search mode, the interference region was $\left[{35}^{\circ },{40}^{\circ }\right]$, and the notches’ level was −30 dB. For the tracking mode, the interference region was $\left[{15}^{\circ },{20}^{\circ }\right]$, and the notches’ level was −26 dB. Figure 11 shows the transmit beampattern performance. It can be observed that the transmit beampattern produced using the proposed method could simultaneously form sidelobes and notches lower than the expected levels.

5.2. Subpulses’ Design and RC Performance

In this part, some simulations are presented that were conducted to verify the effectiveness of the proposed subpulse design method and the proposed signal processing scheme. Search mode was considered first, and the transmit beampattern result in Figure 5 was adopted here. From Table 1, we can see the useful TBs and the corresponding power proportion. It was assumed that the pulse width was 100 $\mathsf{\mu }$s, and the width of each subpulse could be obtained via the power proportion. For search mode, the LFM signal was adopted, and the signal expression for each subpulse was

$$s_{\tilde{k}}\left(t\right)=exp\left(j2\pi f_{\tilde{k}}t+j\pi {\mu }_{\tilde{k}}{t}^2\right),\left|t\right|\le {\displaystyle \frac{T_{\tilde{k}}}{2}},\tilde{k}=1,\cdots ,\tilde{K}$$

(34)

where $f_{\tilde{k}}$, ${\mu }_{\tilde{k}}$, and $T_{\tilde{k}}$ denote the start frequency, chirp rate, and subpulse width of the $\tilde{k}$th subpulse, respectively. $\tilde{K}$ is the number of useful TBs whose power proportion is non-zero. From Table 1, we can see that $\tilde{K}=11$. To ensure the same range resolution for each subpulse, we set ${\mu }_{\tilde{k}}=B/T_{\tilde{k}}$ and $B=5\mathrm{MHz}$. Finally, in order to obtain a low cross-correlation, according to (32), the start frequency was set to

$$f_{\tilde{k}+1}-f_{\tilde{k}}=B,\tilde{k}=1,\cdots ,\tilde{K}-1$$

(35)

Figure 12 shows the frequency spectrum of each subpulse. In Figure 12, subpulses 1 to 11 are carried via the useful TBs ${\mathbf{x}}_{\tilde{k}},\tilde{k}=1,\cdots ,11$, respectively. It can be observed that the interval between each start frequency is B, and the frequency spectrum is non-overlapping between subpulses. Figure 13 shows the correlation performance of subpulses. It can be seen that the cross-correlation level between subpulses is lower than −25 dB, which can satisfy the detection requirement of a radar.

In Figure 14, we verified the RC performance of the proposed signal processing scheme. The RC process was implemented using

$$\mathrm{RC}\left(\tau \right)=\left|{\int }_{-\infty }^{\infty }{\tilde{y}}_{rb}\left(t,{\theta }_k\right)f_k^{*}\left(t+\tau \right)dt\right|$$

where ${\theta }_k$ denotes the mainlobe direction of the k-th TB, and $f_k\left(t\right)$ is the subpulse carried via the k-th TB. According to Table 1, it can be understood that the corresponding RC for the AW at $-26.3^{\circ }$ is subpulse 1. Figure 14a shows the RC result of the AW at $-26.3^{\circ }$ and the auto-correlation of subpulse 1. It can be seen that the RC performance was almost the same as that of the expected LFM signal. It can be noticed that some TBs were not used, for example, the TB pointing to $-{30}^{\circ }$. Given that the signals within the mainlobe region were highly correlated, subpulse 1 was also adopted as the RC for the AW at $-{30}^{\circ }$. Figure 14b shows the RC result for the AW at $-{30}^{\circ }$ and the auto-correlation of subpulse 1. It can be found that results similar to those of $-26.3^{\circ }$ could be obtained.

In the tracking mode, only the auto-correlation of subpulses was considered, and the transmit beampattern result in Figure 8 was adopted. For the result in Figure 8, the power proportion for each TB was 0.3330, 0.3334, and 0.3336, respectively. It was assumed that the pulse width was 10 $\mathsf{\mu }s$, and the width of each subpulse could be obtained using the power proportion. For the tracking mode, the NLFM signal was adopted, and the phase expression of each subpulse was

$${\phi }_{\tilde{k}}\left(t\right)={\displaystyle \frac{\pi {\xi }_{\tilde{k}}{t}^2}{T_{\tilde{k}}}}-{\displaystyle \frac{\pi {\zeta }_{\tilde{k}}{T}_{\tilde{k}}\sqrt{1-4t^2/T_{\tilde{k}}^2}}{2}},\tilde{k}=1,2,3$$

(36)

where

$$\begin{array}{cc}\hfill & \hfill {\xi }_1=-30/T_1,{\xi }_2=30/T_2,{\xi }_3=60/T_3\hfill \\ & \hfill {\zeta }_1=-8/T_1,{\zeta }_2=13/T_2,{\zeta }_3=18/T_3\hfill \end{array}$$

and, for the TB pointing to $-{35}^{\circ }$, $0^{\circ }$ and ${35}^{\circ }$, the width of the subpulse was $T_1=3.33$ $\mathsf{\mu }$s, $T_2=3.334$ $\mathsf{\mu }$s, and $T_3=3.336$ $\mathsf{\mu }$s, respectively. Figure 15 shows the auto-correlation performance of subpulses. It can be seen that, because of the NLFM form, the subpulses had a sidelobe level lower than −20 dB. With the AW at $-{30}^{\circ }$ taken as an example, Figure 16 shows the RC result and the auto-correlation of subpulse 1. It can be seen that the RC performance was almost the same as that of the expected NLFM signal.

5.3. Comparison with Existing Methods

For comparison, the alternating-direction method of multipliers (ADMM) in [26] is considered here. For the transmit beampattern, the work [26] aimed at concentrating the transmission power in the mainlobe regions and maximizing the transmission gain, but it ignored the sidelobe level control. For temporal characteristics, the work [26] focused on the correlation of the beam directions’ AW, and the AW was adopted as the RC. In the proposed method, the subpulse was adopted as the RC, and the correlation of subpulses was considered. It is difficult to solve problem (15) directly. In order to verify the improvement in the solving efficiency using the proposed method, we converted problem (15) into the following SDP problem (37), and we solved it using CVX. According to the optimal solution to the SDP problem, the CM waveform was synthesized via the cyclic algorithm (CA) [38].

$$\begin{array}{l}\underset{\mathbf{R}}{min}{\displaystyle \sum _{l=1}^L}{u}_l\left|P_d\left({\theta }_l\right)-{\mathbf{a}}^{\mathrm{T}}\left({\theta }_l\right)\mathbf{R}{\mathbf{a}}^{*}\left({\theta }_l\right)\right|\\ s.t.\mathbf{R}\ge 0\end{array}$$

(37)

where $\mathbf{R}={\displaystyle \sum _{k=1}^K}{\beta }_k{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{H}}$.

In the following simulations, the multiple-target tracking scene was considered. For the proposed method and SDP, the desired beampattern was obtained by adding two traditional single beams whose beam directions were $-{20}^{\circ }$ and ${30}^{\circ }$, respectively, and the desired sidelobe level was set to −16 dB. For the ADMM, the approaching points were calculated according to the description in [26]. For the temporal design, the proposed method adopted the NLFM signal as subpulses, and the parameter in (36) was changed as follows:

$$\begin{array}{cc}\hfill & \hfill {\xi }_1=45/T_1,{\xi }_2=45/T_2\hfill \\ \hfill & \hfill {\zeta }_1=13/T_1,{\zeta }_2=13/T_2\hfill \end{array}$$

and the pulse width $T_p=T_1+T_2=3.2$ $\mathsf{\mu }$s. For the ADMM, the same NLFM signal was adopted as the desired AW. It should be noted that the pulse width in ADMM was 1.6 $\mathsf{\mu }$s in order to maintain the same range resolution as that in the proposed method. After temporal sampling, the pulse code lengths of all methods were 320.

In Figure 17, we set $N=10$ and show the transmit beampattern result of the proposed method, SDP, CA, and ADMM. Additionally, we also show the result of the proposed method without considering the CM constraint. It can be seen that the proposed method and ADMM achieved similar mainlobe performance, but the mainlobe generated via SDP and CA was slightly wider. This is because SDP and CA use the criterion of minimizing the matching error, while the proposed method and ADMM maximize the mainlobe gain during waveform generation. It can be found that the proposed method achieved a lower sidelobe level than the other methods because only the proposed method focuses on controlling sidelobe levels. Comparing the proposed method with and without a CM constraint, we can find that the presence or absence of a CM constraint has almost no effect on the transmit beampattern performance of the proposed method. Figure 18 shows the RC performance of different methods. In Figure 18, the RC results of AW in beam directions are shown. It can be seen that the proposed method has a lower temporal sidelobe level than the other methods. This is because the proposed method directly uses NLFM signals, while ADMM optimizes phase-coded signals to approximate the desired NLFM signals, and CA does not consider AW properties.

In

Table 2. The ETPI (s) versus the transmitting antenna number.

$\mathit{N}$	10	30	50	70	90
Proposed	0.0533	0.0587	0.0604	0.0621	0.0743
ADMM	1.3164	2.8288	7.3455	12.4647	18.6786
SDP	0.1281	4.2662	10.5617	18.4455	33.3160
CA	0.0012	0.0143	0.0717	0.4978	1.3980

, we show the execution time of per iteration (ETPI) under the scenario in Figure 17 to compare the computational efficiency of all methods. In Table 2, except for the transmitting antenna number N, all the other settings are the same as those in Figure 17. The results of each method are the average of 100 trials with different random initial points. From Table 2, we can find that the ETPI of the proposed method is much less than that of ADMM and SDP, and the gap in the ETPI became larger as the transmitting antenna number increased. Moreover, the simulation experience revealed that the typical iteration number of the proposed method is about 10∼20, while the typical iteration number of ADMM is about 1000∼2000. Although CA is more efficient with a small number of elements, it is a subsequent step of SDP, and the combination of the two is necessary to obtain a CM waveform. In

Table 3. The ETPI (s) versus the code length.

$\mathit{Code}\mathit{Length}$	32	64	128	256	512
Proposed	0.0516	0.0524	0.0527	0.0530	0.0539
ADMM	0.0360	0.0865	0.2126	0.8668	4.4349
SDP	0.1235	0.1247	0.1293	0.1287	0.1286
CA	0.0002	0.0003	0.0006	0.0011	0.0015

, the temporal sampling frequency was changed to evaluate the influence of the code length on the computational efficiency. In Table 3, except for the code length, all the other settings are the same as those in Figure 17. Because the subpulses are produced via analytical expressions, it can be seen that the ETPI of the proposed method remained almost unchanged with a code-length increase. Additionally, due to the fact that SDP optimization only involves the number of array elements, its ETPI is almost unaffected. However, the ETPI of ADMM increases significantly with the code length. According to the analysis in [26], we can understand that the computational complexity of ADMM is about $\mathrm{O}\left(N{L}^2\right)$, where L denotes the code length.

6. Conclusions

In this paper, we have proposed a fast method to produce a CM MIMO radar waveform. In order to improve the computational efficiency, a series of reasonable, simplified operations were proposed. The final simulations showed that the simplified operations can guarantee the beampattern and RC performance of the transmit waveform while improving the computational efficiency. In the proposed method, the key idea is the perturbation vector, which guarantees the CM property without a heavy computational burden. However, the boundary of perturbation parameter $\epsilon$ was not clarified, which may need further research. Moreover, for the transmit beampattern design, the improper desired beampattern may lead to no solution. Therefore, the transmit beampattern boundary of the proposed method should be investigated in the future.