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

Abstract

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.

1. 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 [11] based on fast Fourier transform (FFT) is proposed to minimize the autocorrelation sidelobes and the cross-correlation. This algorithm is computationally efficient and can be used for the design of long waveforms. As it is impossible to design completely orthogonal (i.e., both the autocorrelation sidelobes and the cross-correlation are zeroes) waveform sets [11,16], [11] proposes to design the waveform sets that are orthogonal only at the specified intervals and extends the classical WeCAN (weighted cyclic algorithm-new) [26] algorithm to MIMO radar. To solve the same problem, [12] develops the LBFGS (limited-memory Broyden, Fletcher, Goldfarb and Shanno) iterative algorithm, which is more efficient than the WeCAN algorithm. However, because of the complicated linear search rule for determining the step size, the LBFGS iterative algorithm is still time consuming. Recently, the majorization-minimization (MM)-based algorithms (i.e., MM-Corr (MM-correlation) and MM-WeCorr (MM-weighted correlation)) are proposed in [14]. These two algorithms are also based on FFT operations and much faster than the CAN and WeCAN algorithms [11].

Another waveform set with a good correlation property is the complementary sets of sequences (CSS). A waveform set is called CSS if and only if the autocorrelation sum of the waveforms is a delta function [17]. CSS design is proposed to overcome the difficulty of generating a single unimodular waveform with ideal (impulse-like) autocorrelation. A common application of the CSS is the pulse compression [18,19,20]. In pulse compression radar, the complementary sequences are used to modulate consecutive pulses in a coherent pulse train. Then, the autocorrelation sidelobes can be reduced via the coherent [19] or noncoherent [18] accumulation, which can be regarded as the process of obtaining the autocorrelation sum of the complementary sequences. Moreover, due to the good correlation property, CSS has been widely applied to the CDMA system [21], ISI (intersymbol interference) channel estimation [22], orthogonal frequency division multiplexing (OFDM) [23], and many other areas. The main methods of designing CSS are the analytical construction methods, which have great limitation in generating long waveforms. To overcome this problem, [24] introduces a computational framework based on an iterative twisted approximation (ITROX) for periodically complementary sets of sequences design. Subsequently, [25] extends the CAN algorithm [26] and develops a fast algorithm named CANARY (CAN complementary). Additionally, [14] applies the MM method to the design of CSS.

In addition to the good correlation property, waveform sets are expected to have a good stopband property when the radar systems work in a crowded electromagnetic environment. Waveforms with the stopband property, also known as the sparse frequency waveforms (SFW) in many literature works, are a kind of waveforms with several frequency stopbands. The applications of SFW include ultra-wide bandwidth (UWB) systems [27], high frequency surface wave radar (HFSWR) [28,29] and cognitive radar [30]. By designing waveforms with the stopband property, it can effectively overcome the narrowband interference in the congested frequency bands. At present, there are many research works on the design of a single waveform with the stopband property [31,32,33,34,35,36], but they cannot be used in MIMO systems. Therefore, [36] proposes an iterative algorithm combined with the steepest descent (SD) method for MIMO waveform design. By searching along the gradient direction, the convergence speed of this algorithm is improved. However, the computation of the step size along the gradient direction is complicated, which makes the algorithm still costly. In order to improve the computational efficiency, [37] proposes an algorithm named MDISAA-SFW (multi-dimensional iterative spectral approximation algorithm-SFW) based on alternating projection and phase retrieval.

In this paper, we consider the problem of designing unimodular waveform sets with good correlation and stopband properties and propose a gradient-based algorithm, i.e., Gra-WeCorr-SFW (gradient-weighted correlation-SFW). By using the relationship between the correlation function and the power spectrum density (PSD), the design problem is formulated as an unconstrained minimization problem in the frequency domain. Then, the phase gradient is deduced, and its matrix form is given. In order to avoid searching the step size, we use the Taylor series expansion to derive the step size, which is more efficient than the traditional searching methods. Since both the gradient and the step size can be implemented via the FFT operations and the Hadamard product, the proposed algorithm has high computational efficiency. We also deduce the simplified algorithm named Gra-Corr-SFW (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. $\mathrm{Re}(·)$ and $\mathrm{Im}(·)$ denote the real and imaginary part, respectively. $\mathrm{Diag}(\mathbf{x})$ denotes a diagonal matrix formed with the column vector $\mathbf{x}$. ∘ denotes the Hadamard product. $x(m)$ denotes the m-th element of the vector $\mathbf{x}$. ${\mathbf{x}}^l$ is the l-th iteration of $\mathbf{x}$. ${\mathbf{1}}_N$ is the all-one vectors of length N. ${\mathbf{I}}_N$ denotes the $N\times N$ identity matrix. $\mathcal{F}(\mathbf{x})$ and ${\mathcal{F}}^{-1}(\mathbf{x})$ denote the $2N$-point FFT and IFFT (inverse FFT) operations of $\mathbf{x}$, respectively. ${\mathcal{F}}_c(\mathbf{X})$ and ${\mathcal{F}}_c^{-1}(\mathbf{X})$ represent the FFT and IFFT of each column of the matrix $\mathbf{X}$, respectively. In the (I)FFT operations, if the length of $\mathbf{x}$ is less than $2N$, $\mathbf{x}$ is padded with trailing zeros to length $2N$. $e^{(·)}$ is the exponent arithmetic applied to the scalar, vector or matrix.

2. 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.

2.1. The Criterion for Good Correlation Property

Let ${\{{\mathbf{x}}_m\}}_{m=1}^M$ be the complex waveform set to be designed. The vector form of each waveform can be expressed as:

$${\mathbf{x}}_m={\left[x_m(1),x_m(2),...,x_m(N)\right]}^T,m=1,...,M,$$

(1)

where N is the waveform length and M denotes the number of the waveforms. Then, the correlation function of ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$ is defined as:

$$\begin{array}{c}{r}_{ij}(k)={\displaystyle \sum _{n=k+1}^N}{x}_i(n)x_j^{*}(n-k)=r_{ji}^{*}(-k),i,j=1,...,M,k=1-N,...,N-1.\hfill \end{array}$$

(2)

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:

$$\mathrm{CISL}={\displaystyle \sum _{\genfrac{}{}{0pt}{}{k=1-N}{k\ne 0}}^{N-1}}{\left|{\displaystyle \sum _{m=1}^M}{r}_{mm}(k)\right|}^2.$$

(3)

for the waveform sets with both good auto- and cross-correlation properties, we consider the following more general weighted measure [14]:

$$\psi ={\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\displaystyle \sum _{k=1-N}^{N-1}}{w}_t(k){\left|r_{ij}(k)\right|}^2-w_t(0)M{N}^2,$$

(4)

where $w_t(k)=w_t(-k)\ge 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:

$$\begin{array}{cc}\hfill {\mathbf{r}}_{ij}& ={\left[r_{ij}(0),r_{ij}(1),...,r_{ij}(N-1),0,r_{ij}(1-N),...,r_{ij}(-1)\right]}^T,\hfill \\ \hfill {\mathbf{w}}_t& ={\left[w_t(0),w_t(1),...,w_t(N-1),1,w_t(1-N),...,w_t\left(-1\right)\right]}^T\hfill \end{array}$$

(5)

be the correlation vector and the corresponding correlation weight vector, respectively. The Fourier transform of the correlation function ${\mathbf{r}}_{ij}$ (i.e., power spectrum density (PSD)) can be written as:

$$\begin{array}{cc}\hfill {\mathbf{p}}_{ij}& =\mathcal{F}\left({\mathbf{r}}_{ij}\right)=\mathbf{F}{\mathbf{r}}_{ij},\hfill \\ \hfill {\mathbf{r}}_{ij}& ={\mathcal{F}}^{-1}\left({\mathbf{p}}_{ij}\right)=\frac{1}{2N}{\mathbf{F}}^H{\mathbf{p}}_{ij},\hfill \end{array}$$

(6)

where $\mathbf{F}$ is the $2N\times 2N$ discrete Fourier transform (DFT) matrix with the following expression:

$$\mathbf{F}(m,n)=e^{-j\frac{2mn\pi }{2N}},1\le m,n\le 2N.$$

(7)

According to (5) and (6), (4) can be expressed in the frequency domain as:

$$\psi ={\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\mathbf{r}}_{ij}^H{\mathbf{D}}_t{\mathbf{r}}_{ij}-w_t(0)M{N}^2=\frac{1}{4N^2}{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{ij}^H\mathbf{F}{\mathbf{D}}_t{\mathbf{F}}^H{\mathbf{p}}_{ij}-w_t(0)M{N}^2,$$

(8)

where ${\mathbf{D}}_t=\mathrm{Diag}\left({\mathbf{w}}_t\right)$. It is worth noting that the value of the $N+1$-th element of ${\mathbf{w}}_t$ has no effect on the objective function. In order to facilitate the derivation below, the $N+1$-th element of ${\mathbf{w}}_t$ is set to be one.

Actually, when $w_t(k)=1,k=1-N,...,N-1$ (i.e., the correlation weights are not taken into account), Criterion (8) is equivalent to Criterion (3) (see Appendix A), which means that Criterion (3) is a special case of Criterion (8). Thus, here, we just consider Criterion (8). By ignoring the constant term in (8), the criterion related to the correlation property is given by:

$$J_{CF}=\frac{1}{4N^2}{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{ij}^H\mathbf{F}{\mathbf{D}}_t{\mathbf{F}}^H{\mathbf{p}}_{ij}.$$

(9)

2.2. The Criterion for the Good Stopband Property

The waveform set with a good stopband property means that the PSD of each waveform has several frequency stopbands. Without loss of generality, we consider that the frequency is normalized. Define the set of frequency stopbands as:

$${\mathsf{\Omega }}_f={\displaystyle \underset{k=1}{\overset{n_s}{\cup }}}(f_{k1},f_{k2})\subset [0,1],$$

(10)

where $(f_{k1},f_{k2})$ denotes one stopband and $n_s$ denotes the number of the stopbands. Considering the $2N$-point FFT operations, the corresponding point set of the stopband set ${\mathsf{\Omega }}_f$ can be expressed as:

$${\overline{\mathsf{\Omega }}}_f={\displaystyle \underset{k=1}{\overset{n_s}{\cup }}}(2N{f}_{k1},2N{f}_{k2})\subset [0,2N].$$

(11)

Define the frequency weight vector as:

$$\begin{array}{c}{\mathbf{w}}_f={\left[w_f(1),w_f(2),...,w_f(2N)\right]}^T,\hfill \\ w_f(p)=\left\{\begin{array}{cc}1,\hfill & p\in {\overline{\mathsf{\Omega }}}_f\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right..\hfill \end{array}$$

(12)

As the PSD of each waveform is nonnegative, the criterion related to the stopband property (or PSD) can be formulated as:

$$J_{PSD}={\displaystyle \sum _{i=1}^M}{∥{\mathbf{w}}_f\circ {\mathbf{p}}_{ii}∥}^2={\displaystyle \sum _{i=1}^M}{\mathbf{p}}_{ii}^H{\mathbf{D}}_f{\mathbf{p}}_{ii},$$

(13)

where ${\mathbf{D}}_f=\mathrm{Diag}\left({\mathbf{w}}_f\right)$.

2.3. The Minimization Problem

In order to obtain good correlation and stopband properties, both $J_{CF}$ and $J_{PSD}$ should be minimized. Thus, the waveform design problem can be regarded as a multi-objective optimization (also known as Pareto optimization) problem, i.e., finding the Pareto optimal solutions that satisfy the constraints. Here, we apply the traditional weighting method to solve the problem and, thus, formulate a single-objective function as follows:

$$J_T=\lambda J_{PSD}+\left(1-\lambda \right)J_{CF}={\displaystyle \sum _{i=1}^M}{\mathbf{p}}_{ii}^H\mathbf{Q}{\mathbf{p}}_{ii}+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{ij}^H{\mathbf{Q}}^{\prime }{\mathbf{p}}_{ij},$$

(14)

where $\lambda \in [0,1]$ is a weight coefficient by which we can balance the relative weight between $J_{CF}$ and $J_{PSD}$, and:

$$\mathbf{Q}=\lambda {\mathbf{D}}_f,{\mathbf{Q}}^{\prime }=\frac{1}{4N^2}\left(1-\lambda \right)\mathbf{F}{\mathbf{D}}_t{\mathbf{F}}^H.$$

(15)

When $\lambda =0$, (14) is the criterion for designing waveform sets with a good correlation property. Additionally, when $\lambda \in (0,1]$, (14) becomes the criterion for SFW design. Generally, for maximizing the transmitter efficiency and reducing the requirement to the hardware, the unimodular constraint [38] is required in the waveform design. Therefore, the design problem can be formulated as the following minimization problem:

$$\begin{array}{c}\underset{{\{{\mathbf{x}}_m\}}_{m=1}^M}{min}{J}_T\hfill \\ \mathrm{s}.\mathrm{t}.\left|x_m(n)\right|=1,n=1,...,N,m=1,...,M.\hfill \end{array}$$

(16)

Let ${\mathsf{\varphi }}_m$ denote the phase vector of ${\mathbf{x}}_m$, i.e., ${\mathbf{x}}_m={\left[e^{j{\varphi }_m(1)},e^{j{\varphi }_m(2)},...,e^{j{\varphi }_m(N)}\right]}^T$. Then, the problem (16) can be reformulated as the following unconstrained problem with regard to ${\{{\mathsf{\varphi }}_m\}}_{m=1}^M$:

$$\underset{{\{{\mathsf{\varphi }}_m\}}_{m=1}^M}{min}{J}_T={\displaystyle \sum _{i=1}^M}{\mathbf{p}}_{ii}^H\mathbf{Q}{\mathbf{p}}_{ii}+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{ij}^H{\mathbf{Q}}^{\prime }{\mathbf{p}}_{ij}.$$

(17)

3. Problem Optimization via the Gradient Method

In this section, we optimize Problem (17) by using the gradient-based algorithm, which is able to guarantee that the objective function is monotonically decreasing at each iteration. In the following, we first derive the phase gradient and the step size and then briefly summarize the algorithm.

3.1. Phase Gradient

To facilitate the derivation, let ${\varphi }_m(n)={\varphi }_{mn}$. Before deriving the phase gradient ${\nabla }_{{\mathsf{\varphi }}_m}{J}_T$, we first deduce the derivative of $J_T$ with respect to the phase ${\varphi }_{mn}$:

$$\begin{array}{cc}\hfill \frac{\partial J_T}{\partial {\varphi }_{mn}}& =\partial \left({\displaystyle \sum _{i=1}^M}{\mathbf{p}}_{ii}^H\mathbf{Q}{\mathbf{p}}_{ii}+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{ij}^H{\mathbf{Q}}^{\prime }{\mathbf{p}}_{ij}\right)/\partial {\varphi }_{mn}\hfill \\ & =\partial \left({\mathbf{p}}_{mm}^H\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mm}\right)/\partial {\varphi }_{mn}+\partial \left({\displaystyle \sum _{\genfrac{}{}{0pt}{}{j=1}{j\ne m}}^M}{\mathbf{p}}_{mj}^H{\mathbf{Q}}^{\prime }{\mathbf{p}}_{mj}{\displaystyle \sum _{\genfrac{}{}{0pt}{}{i=1}{i\ne m}}^M}{\mathbf{p}}_{im}^H{\mathbf{Q}}^{\prime }{\mathbf{p}}_{im}\right)/\partial {\varphi }_{mn}.\hfill \end{array}$$

(18)

Since ${\displaystyle \sum _{\genfrac{}{}{0pt}{}{i=1}{i\ne m}}^M}{\mathbf{p}}_{im}^H{\mathbf{Q}}^{\prime }{\mathbf{p}}_{im}={\displaystyle \sum _{\genfrac{}{}{0pt}{}{i=1}{i\ne m}}^M}{\mathbf{p}}_{mi}^H{{\mathbf{Q}}^{\prime }}^{*}{\mathbf{p}}_{mi}$ (see Appendix B), (18) can be rewritten as:

$$\frac{\partial J_T}{\partial {\varphi }_{mn}}=\partial \left({\mathbf{p}}_{mm}^H\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mm}\right)/\partial {\varphi }_{mn}+2{\displaystyle \sum _{\genfrac{}{}{0pt}{}{i=1}{i\ne m}}^M}\partial \left({\mathbf{p}}_{mi}^H\mathrm{Re}\left({\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mi}\right)/\partial {\varphi }_{mn}.$$

(19)

The first term in (19) can be simplified as:

$$\begin{array}{cc}\hfill \frac{\partial \left({\mathbf{p}}_{mm}^H\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mm}\right)}{\partial {\varphi }_{mn}}& =\frac{\partial \left({\mathbf{p}}_{mm}^T\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mm}\right)}{\partial {\varphi }_{mn}}\hfill \\ & =\frac{\partial {\mathbf{p}}_{mm}^T}{\partial {\varphi }_{mn}}\left(\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right)+{\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right)}^T\right){\mathbf{p}}_{mm}\hfill \\ & =2\frac{\partial {\mathbf{p}}_{mm}^T}{\partial {\varphi }_{mn}}\mathrm{Re}\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mm},\hfill \end{array}$$

(20)

where the third equality follows from the fact ${\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right)}^T={\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right)}^{*}$. Let:

$$\begin{array}{cc}\hfill {\mathbf{y}}_{mm}& =\mathrm{Re}\left(\mathbf{Q}\right){\mathbf{p}}_{mm}\hfill \\ \hfill {{\mathbf{y}}^{\prime }}_{mi}& =\mathrm{Re}\left({\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mi},i=1,...,M,\hfill \end{array}$$

(21)

then (20) can be expressed as:

$$\begin{array}{cc}\hfill \frac{\partial \left({\mathbf{p}}_{mm}^H\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mm}\right)}{\partial {\varphi }_{mn}}& =2\frac{\partial {\mathbf{p}}_{mm}^T}{\partial {\varphi }_{mn}}\left({\mathbf{y}}_{mm}+{{\mathbf{y}}^{\prime }}_{mm}\right)\hfill \\ & =2{\displaystyle \sum _{k=1}^{2N}}\frac{\partial \left(p_{mm}(k)\right)}{\partial {\varphi }_{mn}}\left(y_{mm}(k)+{y^{\prime }}_{mm}(k)\right).\hfill \end{array}$$

(22)

According to (A4) and (A5), we can deduce the derivative $\partial \left(p_{mm}(k)\right)/\partial {\varphi }_{mn}$ as follows:

$$\frac{\partial \left(p_{mm}(k)\right)}{\partial {\varphi }_{mn}}=\frac{\partial \left(f_m(k)f_m^{*}(k)\right)}{\partial {\varphi }_{mn}}=2\mathrm{Re}\left(\frac{\partial f_m^{*}(k)}{\partial {\varphi }_{mn}}{f}_m(k)\right)=2\mathrm{Re}\left(-j{x}_m^{*}(n)e^{j\frac{2\pi }{2N}nk}{f}_m(k)\right).$$

(23)

By substituting (23) into (22), we have:

$$\frac{\partial \left({\mathbf{p}}_{mm}^H\left(\mathbf{Q}+{\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mm}\right)}{\partial {\varphi }_{mn}}=4\mathrm{Re}\left(-j{x}_m^{*}(n){\displaystyle \sum _{k=1}^{2N}}{f}_m(k)\left(y_{mm}(k)+{y^{\prime }}_{mm}(k)\right)e^{j\frac{2\pi }{2N}nk}\right).$$

(24)

To compute the second term of (19), we first deduce $\partial \left({\mathbf{p}}_{mi}^H\mathrm{Re}\left({\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mi}\right)/\partial {\varphi }_{mn}(i\ne m)$:

$$\frac{\partial \left({\mathbf{p}}_{mi}^H\mathrm{Re}\left({\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mi}\right)}{\partial {\varphi }_{mn}}=2\mathrm{Re}\left(\frac{\partial \left({\mathbf{p}}_{mi}^H\right)}{\partial {\varphi }_{mn}}\mathrm{Re}\left({\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mi}\right)=2\mathrm{Re}\left({\displaystyle \sum _{k=1}^{2N}}\frac{\partial \left(p_{mi}^{*}(k)\right)}{\partial {\varphi }_{mn}}{y^{\prime }}_{mi}(k)\right).$$

(25)

Similarly to the derivation of (23), it is easy to obtain that:

$$\frac{\partial \left(p_{mi}(k)\right)}{\partial {\varphi }_{mn}}=j{x}_m(n)e^{-j\frac{2\pi }{2N}nk}{f}_i^{*}(k),i\ne m.$$

(26)

On the basis of (25) and (26), $\partial \left({\mathbf{p}}_{mi}^H\mathrm{Re}\left({\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mi}\right)/\partial {\varphi }_{mn}$ can be further denoted as

$$\frac{\partial \left({\mathbf{p}}_{mi}^H\mathrm{Re}\left({\mathbf{Q}}^{\prime }\right){\mathbf{p}}_{mi}\right)}{\partial {\varphi }_{mn}}=2\mathrm{Re}\left(-j{x}_m^{*}(n){\displaystyle \sum _{k=1}^{2N}}{f}_i(k){y^{\prime }}_{mi}(k)e^{j\frac{2\pi }{2N}nk}\right).$$

(27)

By substituting (24) and (27) into (19), $\partial J_T/\partial {\varphi }_{mn}$ can be simplified as:

$$\frac{\partial J_T}{\partial {\varphi }_{mn}}=4\mathrm{Re}\left(-j{x}_m^{*}(n){\displaystyle \sum _{k=1}^{2N}}\left(f_m(k)y_{mm}(k)+{\displaystyle \sum _{i=1}^M}{f}_i(k){y^{\prime }}_{mi}(k)\right)e^{j\frac{2\pi }{2N}nk}\right).$$

(28)

Let ${\mathbf{z}}_m$ be the inverse FFT (IFFT) of ${\mathbf{f}}_m\circ {\mathbf{y}}_{mm}+{\displaystyle \sum _{i=1}^M}{\mathbf{f}}_i\circ {{\mathbf{y}}^{\prime }}_{mi}$, i.e.,

$$\begin{array}{cc}\hfill {\mathbf{z}}_m& ={\mathcal{F}}^{-1}\left({\mathbf{f}}_m\circ {\mathbf{y}}_{mm}+{\displaystyle \sum _{i=1}^M}{\mathbf{f}}_i\circ {{\mathbf{y}}^{\prime }}_{mi}\right)\hfill \\ \hfill z_m(n)& =\frac{1}{2N}{\displaystyle \sum _{k=1}^{2N}}\left(f_m(k)y_{mm}(k)+{\displaystyle \sum _{i=1}^M}{f}_i(k){y^{\prime }}_{mi}(k)\right)e^{j\frac{2\pi }{2N}nk},\hfill \end{array}$$

(29)

then $\partial J_T/\partial {\varphi }_{mn}$ in (28) becomes:

$$\frac{\partial J_T}{\partial {\varphi }_{mn}}=8N\mathrm{Re}\left(-j{x}_m^{*}(n)z_m(n)\right).$$

(30)

By stacking (30) in a vector, the phase gradient ${\nabla }_{{\mathsf{\varphi }}_m}{J}_T$ is given by:

$${\nabla }_{{\mathsf{\varphi }}_m}{J}_T={\left[\partial J_T/\partial {\varphi }_{m1},...,\partial J_T/\partial {\varphi }_{mN}\right]}^T=8N\mathrm{Re}\left(-j{\mathbf{x}}_m^{*}\circ {\mathbf{z}}_m(1:N)\right),$$

(31)

where ${\mathbf{z}}_m(1:N)$ denotes the first N elements of ${\mathbf{z}}_m$. It is worth noting that ${\mathbf{y}}_{mm}$ and ${{\mathbf{y}}^{\prime }}_{mi}$ defined in (21) can be calculated by the Hadamard product and the FFT operations:

$$\begin{array}{cc}\hfill {\mathbf{y}}_{mm}& =\lambda {\mathbf{D}}_f{\mathbf{p}}_{mm}=\lambda {\mathbf{w}}_t\circ {\mathbf{p}}_{mm},\hfill \\ \hfill {{\mathbf{y}}^{\prime }}_{mi}& =\mathrm{Re}\left(\frac{1}{4N^2}\left(1-\lambda \right)\mathbf{F}{\mathbf{D}}_t{\mathbf{F}}^H\right){\mathbf{p}}_{mi}\hfill \\ & =\frac{1}{8N^2}\left(1-\lambda \right)\left(\mathbf{F}{\mathbf{D}}_t{\mathbf{F}}^H+{\mathbf{F}}^H{\mathbf{D}}_t\mathbf{F}\right){\mathbf{p}}_{mi}\hfill \\ & =\frac{1}{4N}\left(1-\lambda \right)\left(\mathcal{F}\left({\mathbf{w}}_t\circ {\mathcal{F}}^{-1}\left({\mathbf{p}}_{mi}\right)\right)+{\mathcal{F}}^{-1}\left({\mathbf{w}}_t\circ \mathcal{F}\left({\mathbf{p}}_{mi}\right)\right)\right).\hfill \end{array}$$

(32)

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 $\mathbf{X}=[{\mathbf{x}}_1,...,{\mathbf{x}}_M]$ and $\mathsf{\Phi }=[{\mathsf{\varphi }}_1,...,{\mathsf{\varphi }}_M]$ denote the matrix of waveform set and the corresponding phase matrix, respectively. Then, the spectrum matrix of $\mathbf{X}$ is $\mathbf{S}={\mathcal{F}}_c(\mathbf{X})=[{\mathbf{f}}_1,...,{\mathbf{f}}_M]$. Define:

$$\begin{array}{cc}\hfill \mathbf{Y}& =\left[{\mathbf{y}}_{11},...,{\mathbf{y}}_{MM}\right]=\lambda \left({\mathbf{w}}_f{\mathbf{1}}_M^T\right)\circ \mathbf{S}\circ {\mathbf{S}}^{*},\hfill \\ \hfill {\mathbf{Y}}_1& =\left[{\displaystyle \sum _{i=1}^M}{\mathbf{f}}_i\circ {{\mathbf{y}}^{\prime }}_{1i},...,{\displaystyle \sum _{i=1}^M}{\mathbf{f}}_i\circ {{\mathbf{y}}^{\prime }}_{Mi}\right],\hfill \end{array}$$

(33)

then the matrix form of (29) is given by:

$$\mathbf{Z}=[{\mathbf{z}}_1,...,{\mathbf{z}}_M]={\mathcal{F}}_c^{-1}\left(\mathbf{S}\circ \mathbf{Y}+{\mathbf{Y}}_1\right).$$

(34)

Thus, according to (31), the matrix of phase gradient $\mathbf{G}$ can be expressed as:

$$\mathbf{G}=\left[{\nabla }_{{\mathsf{\varphi }}_1}{J}_T,...,{\nabla }_{{\mathsf{\varphi }}_M}{J}_T\right]=8N\mathrm{Re}\left(-j{\mathbf{X}}^{*}\circ \mathbf{Z}(1:N,:)\right),$$

(35)

where $\mathbf{Z}(1:N,:)$ denotes the submatrix formed with the first N rows of $\mathbf{Z}$. By defining $\mathbf{Y}$ and ${\mathbf{Y}}_1$, we can easily calculate the gradient matrix. However, it is still hard to calculate ${\mathbf{Y}}_1$ directly by (33). To efficiently calculate ${\mathbf{Y}}_1$, define:

$${\mathbf{Y}}^{\prime }=\left[\begin{array}{ccc}{{\mathbf{y}}^{\prime }}_{11}& \cdots & {{\mathbf{y}}^{\prime }}_{1M}\\ ⋮& \ddots & ⋮\\ {{\mathbf{y}}^{\prime }}_{M1}& \cdots & {{\mathbf{y}}^{\prime }}_{MM}\end{array}\right]\in {\mathbb{C}}^{2NM\times M},{\mathbf{S}}_{rep}={\left[\begin{array}{c}\mathbf{S}\\ ⋮\\ \mathbf{S}\end{array}\right]}_{2NM\times M}.$$

(36)

Then it is easy to verify that:

$$\left({\mathbf{S}}_{rep}\circ {\mathbf{Y}}^{\prime }\right){\mathbf{1}}_M=\mathrm{vec}\left({\mathbf{Y}}_1\right),$$

(37)

where $\mathrm{vec}\left({\mathbf{Y}}_1\right)$ denotes the column vector consisting of all of the columns of ${\mathbf{Y}}_1$. Thus, we can obtain ${\mathbf{Y}}_1$ via the inverse operation of $\mathrm{vec}\left(·\right)$, i.e.,

$${\mathbf{Y}}_1=\mathrm{ve}{\mathrm{c}}^{-1}\left(\left({\mathbf{S}}_{rep}\circ {\mathbf{Y}}^{\prime }\right){\mathbf{1}}_M\right),$$

(38)

where ${\mathbf{Y}}^{\prime }$ can be calculated by (32). It is easy to see that ${{\mathbf{y}}^{\prime }}_{mi}$ can be implemented by four FFT (IFFT) operations. Since ${{\mathbf{y}}^{\prime }}_{mi}={{\mathbf{y}}^{\prime }}_{im}^{*}$, the calculation of ${\mathbf{Y}}^{\prime }$ takes $4\frac{M(M+1)}{2}=2M^2+2M$ FFT (IFFT) operations.

3.2. 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 ${\mathsf{\Phi }}^l$ is the phase matrix of the present iteration point ${\mathbf{X}}^l=\left[{\mathbf{x}}_1^l,...,{\mathbf{x}}_M^l\right]=e^{j{\mathsf{\Phi }}^l}$, and ${\mathbf{D}}^l=\left[{\mathbf{d}}_1^l,...,{\mathbf{d}}_M^l\right]$ is the descent direction. Then, the new iteration point can be denoted as:

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}^{l+1}& ={\mathsf{\Phi }}^l+\mu {\mathbf{D}}^l,\hfill \\ \hfill {\mathbf{X}}^{l+1}& =e^{j{\mathsf{\Phi }}^{l+1}}={\mathbf{X}}^l\circ e^{j\mu {\mathbf{D}}^l},\hfill \end{array}$$

(39)

i.e.,

$${\mathbf{x}}_i^{l+1}={\mathbf{x}}_i^l\circ e^{j\mu {\mathbf{d}}_i^l},i=1,...,M,$$

(40)

where $\mu$ is the step size. Thus, the linear search problem can be formulated as the following minimization problem:

$$\underset{\mu >0}{min}h(\mu )=J_T\left({\mathbf{X}}^{l+1}\right)={\displaystyle \sum _{i=1}^M}{\left({\mathbf{p}}_{ii}^{l+1}\right)}^H{\mathbf{Qp}}_{ii}^{l+1}+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\left({\mathbf{p}}_{ij}^{l+1}\right)}^H{{\mathbf{Q}}^{\prime }\mathbf{p}}_{ij}^{l+1}.$$

(41)

By taking the derivative of (41), we have:

$$\frac{\partial h}{\partial \mu }=2\mathrm{Re}\left({\displaystyle \sum _{i=1}^M}\frac{\partial {\left({\mathbf{p}}_{ii}^{l+1}\right)}^H}{\partial \mu }{\mathbf{Qp}}_{ii}^{l+1}+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}\frac{\partial {\left({\mathbf{p}}_{ij}^{l+1}\right)}^H}{\partial \mu }{{\mathbf{Q}}^{\prime }\mathbf{p}}_{ij}^{l+1}\right).$$

(42)

To simplify (41) and (42), we first deduce ${\mathbf{p}}_{ij}^{l+1}$. By using the Taylor series expansion, ${\mathbf{p}}_{ij}^{l+1}$ can be approximated as (see Appendix C):

$${\mathbf{p}}_{ij}^{l+1}\approx {\mathbf{f}}_i\circ {\mathbf{f}}_j^{*}+j\left({{\mathbf{f}}^{\prime }}_i\circ {\mathbf{f}}_j^{*}-{\mathbf{f}}_i\circ {{\mathbf{f}}^{\prime }}_j^{*}\right)\mu -\frac{1}{2}\left({{\mathbf{f}}^{\prime \prime }}_i\circ {\mathbf{f}}_j^{*}+{\mathbf{f}}_i\circ {{\mathbf{f}}^{\prime \prime }}_j^{*}-2{{\mathbf{f}}^{\prime }}_i\circ {{\mathbf{f}}^{\prime }}_j^{*}\right){\mu }^2,$$

(43)

where:

$${\mathbf{f}}_i=\mathcal{F}\left({\mathbf{x}}_i^l\right),{{\mathbf{f}}^{\prime }}_i=\mathcal{F}\left({\mathbf{x}}_i^l\circ {\mathbf{d}}_i^l\right),{{\mathbf{f}}^{\prime \prime }}_i=\mathcal{F}\left({\mathbf{x}}_i^l\circ {\mathbf{d}}_i^l\circ {\mathbf{d}}_i^l\right).$$

(44)

Let

$$\begin{array}{cc}\hfill {\mathbf{p}}_{ij}^l& ={\mathbf{f}}_i\circ {\mathbf{f}}_j^{*},\hfill \\ \hfill {\mathbf{c}}_{ij}& =j\left({{\mathbf{f}}^{\prime }}_i\circ {\mathbf{f}}_j^{*}-{\mathbf{f}}_i\circ {{\mathbf{f}}^{\prime }}_j^{*}\right),\hfill \\ \hfill {{\mathbf{c}}^{\prime }}_{ij}& =-\frac{1}{2}\left({{\mathbf{f}}^{\prime \prime }}_i\circ {\mathbf{f}}_j^{*}+{\mathbf{f}}_i\circ {{\mathbf{f}}^{\prime \prime }}_j^{*}-2{{\mathbf{f}}^{\prime }}_i\circ {{\mathbf{f}}^{\prime }}_j^{*}\right),\hfill \end{array}$$

(45)

then (43) can be rewritten as:

$${\mathbf{p}}_{ij}^{l+1}\approx {\mathbf{p}}_{ij}^l+{\mathbf{c}}_{ij}\mu +{{\mathbf{c}}^{\prime }}_{ij}{\mu }^2={\overline{\mathbf{p}}}_{ij}^{l+1}.$$

(46)

By replacing ${\mathbf{p}}_{ij}^{l+1}$ with ${\overline{\mathbf{p}}}_{ij}^{l+1}$, the approximate function of $h(\mu )$ in (41) can be denoted as:

$$h_1(\mu )={\displaystyle \sum _{i=1}^M}{\left({\overline{\mathbf{p}}}_{ii}^{l+1}\right)}^H{\mathbf{Q}\overline{\mathbf{p}}}_{ii}^{l+1}+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\left({\overline{\mathbf{p}}}_{ij}^{l+1}\right)}^H{{\mathbf{Q}}^{\prime }\overline{\mathbf{p}}}_{ij}^{l+1}.$$

(47)

Since ${\mathbf{p}}_{ij}^{l+1}={\overline{\mathbf{p}}}_{ij}^{l+1}$ holds when $\mu =0$, it is easy to verify that:

$$h(0)=h_1(0),{\left.\frac{\partial h}{\partial \mu }\right|}_{\mu =0}={\left.\frac{\partial h_1}{\partial \mu }\right|}_{\mu =0},$$

(48)

which indicates that $h(\mu )$ and $h_1(\mu )$ have the same function value and slope at $\mu =0$. As ${\mathbf{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(\mu )$ is sensitive to the waveform phases, the optimal step size of $h(\mu )$ is very small and close to zero. Consequently, we can use the minimum point of $h_1(\mu )$ to approximate the optimal step size.

To calculate the minimum point of $h_1(\mu )$, we replace ${\mathbf{p}}_{ij}^{l+1}$ in (42) with ${\overline{\mathbf{p}}}_{ij}^{l+1}$, then the derivative of $h_1(\mu )$ with respective to $\mu$ is given by:

$$\begin{array}{cc}\hfill \frac{\partial h_1}{\partial \mu }& =2\mathrm{Re}\left({\displaystyle \sum _{i=1}^M}\left({\mathbf{c}}_{ii}^H+2{{\mathbf{c}}^{\prime }}_{ii}^H\mu \right)\mathbf{Q}\left({\mathbf{p}}_{ii}^l+{\mathbf{c}}_{ii}\mu +{{\mathbf{c}}^{\prime }}_{ii}{\mu }^2\right)\right.\hfill \\ & +\left.{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}\left({\mathbf{c}}_{ij}^H+2{{\mathbf{c}}^{\prime }}_{ij}^H\mu \right){\mathbf{Q}}^{\prime }\left({\mathbf{p}}_{ij}^l+{\mathbf{c}}_{ij}\mu +{{\mathbf{c}}^{\prime }}_{ij}{\mu }^2\right)\right)\hfill \\ & =a{\mu }^3+b{\mu }^2+c\mu +d,\hfill \end{array}$$

(49)

where:

$$\begin{array}{cc}\hfill a& =2\mathrm{Re}\left(2{\displaystyle \sum _{i=1}^M}{{\mathbf{c}}^{\prime }}_{ii}^H\mathbf{Q}{{\mathbf{c}}^{\prime }}_{ii}+2{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{{\mathbf{c}}^{\prime }}_{ij}^H{\mathbf{Q}}^{\prime }{{\mathbf{c}}^{\prime }}_{ij}\right),\hfill \\ \hfill b& =2\mathrm{Re}\left({\displaystyle \sum _{i=1}^M}\left(2{{\mathbf{c}}^{\prime }}_{ii}^H\mathbf{Q}{\mathbf{c}}_{ii}+{\mathbf{c}}_{ii}^H\mathbf{Q}{{\mathbf{c}}^{\prime }}_{ii}\right)+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}\left(2{{\mathbf{c}}^{\prime }}_{ij}^H{\mathbf{Q}}^{\prime }{\mathbf{c}}_{ij}+{\mathbf{c}}_{ij}^H{\mathbf{Q}}^{\prime }{{\mathbf{c}}^{\prime }}_{ij}\right)\right),\hfill \\ \hfill c& =2\mathrm{Re}\left({\displaystyle \sum _{i=1}^M}\left(2{{\mathbf{c}}^{\prime }}_{ii}^H{\mathbf{Qp}}_{ii}^l+{\mathbf{c}}_{ii}^H\mathbf{Q}{\mathbf{c}}_{ii}\right)+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}\left(2{{\mathbf{c}}^{\prime }}_{ij}^H{{\mathbf{Q}}^{\prime }\mathbf{p}}_{ij}^l+{\mathbf{c}}_{ij}^H{\mathbf{Q}}^{\prime }{\mathbf{c}}_{ij}\right)\right),\hfill \\ \hfill d& =2\mathrm{Re}\left({\displaystyle \sum _{i=1}^M}{\mathbf{c}}_{ii}^H{\mathbf{Qp}}_{ii}^l+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\mathbf{c}}_{ij}^H{{\mathbf{Q}}^{\prime }\mathbf{p}}_{ij}^l\right).\hfill \end{array}$$

(50)

To simplify the calculation, let:

$$\begin{array}{cc}\hfill {\mathbf{h}}_{ij}& ={\mathcal{F}}^{-1}\left({\mathbf{c}}_{ij}\right)=\frac{1}{2N}{\mathbf{F}}^H{\mathbf{c}}_{ij},\hfill \\ \hfill {{\mathbf{h}}^{\prime }}_{ij}& ={\mathcal{F}}^{-1}\left({{\mathbf{c}}^{\prime }}_{ij}\right)=\frac{1}{2N}{\mathbf{F}}^H{{\mathbf{c}}^{\prime }}_{ij}.\hfill \end{array}$$

(51)

Then (50) can be rewritten as:

$$\begin{array}{cc}\hfill a& =4\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T{\displaystyle \sum _{i=1}^M}\left({{\mathbf{c}}^{\prime }}_{ii}^{*}\circ {{\mathbf{c}}^{\prime }}_{ii}\right)+\left(1-\lambda \right){\mathbf{w}}_t^T{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}\left({{\mathbf{h}}^{\prime }}_{ij}^{*}\circ {{\mathbf{h}}^{\prime }}_{ij}\right)\right),\hfill \\ \hfill b& =2\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T{\displaystyle \sum _{i=1}^M}\left(2{{\mathbf{c}}^{\prime }}_{ii}^{*}\circ {\mathbf{c}}_{ii}+{\mathbf{c}}_{ii}^{*}\circ {{\mathbf{c}}^{\prime }}_{ii}\right)+\left(1-\lambda \right){\mathbf{w}}_t^T{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}\left(2{{\mathbf{h}}^{\prime }}_{ij}^{*}\circ {\mathbf{h}}_{ij}+{\mathbf{h}}_{ij}^{*}\circ {{\mathbf{h}}^{\prime }}_{ij}\right)\right),\hfill \\ \hfill c& =2\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T{\displaystyle \sum _{i=1}^M}\left(2{{\mathbf{c}}^{\prime }}_{ii}^{*}\circ {\mathbf{p}}_{ii}^l+{\mathbf{c}}_{ii}^{*}\circ {\mathbf{c}}_{ii}\right)+\left(1-\lambda \right){\mathbf{w}}_t^T{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}\left(2{{\mathbf{h}}^{\prime }}_{ij}^{*}\circ {\mathbf{r}}_{ij}^l+{\mathbf{h}}_{ij}^{*}\circ {\mathbf{h}}_{ij}\right)\right),\hfill \\ \hfill d& =2\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T{\displaystyle \sum _{i=1}^M}\left({\mathbf{c}}_{ii}^{*}\circ {\mathbf{p}}_{ii}^l\right)+\left(1-\lambda \right){\mathbf{w}}_t^T{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}\left({\mathbf{h}}_{ij}^{*}\circ {\mathbf{r}}_{ij}^l\right)\right),\hfill \end{array}$$

(52)

where the first term of each equality follows from the fact ${\mathbf{a}}^H\mathbf{Qb}=\lambda {\mathbf{a}}^H{\mathbf{D}}_f\mathbf{b}={\mathbf{w}}_f^T({\mathbf{a}}^{*}\circ \mathbf{b})$ ($\mathbf{a},\mathbf{b}$ are the arbitrary vectors), and the second term of each equality follows from the fact ${\mathbf{a}}^H{\mathbf{Q}}^{\prime }\mathbf{b}=\frac{1}{4N^2}\left(1-\lambda \right){\mathbf{a}}^H\mathbf{F}{\mathbf{D}}_t{\mathbf{F}}^H\mathbf{b}=\left(1-\lambda \right){\mathbf{w}}_t^T\left({\left({\mathcal{F}}^{-1}(\mathbf{a})\right)}^{*}\circ {\mathcal{F}}^{-1}(\mathbf{b})\right)$. Let $\frac{\partial h_1}{\partial \mu }=0$, then the minimum point of $h_1(\mu )$ can be obtained by solving the following cubic equation:

$$a{\mu }^3+b{\mu }^2+c\mu +d=0$$

(53)

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:

$$\mathbf{S}=\left[{\mathbf{f}}_1,...,{\mathbf{f}}_M\right],{\mathbf{S}}^{\prime }=\left[{{\mathbf{f}}^{\prime }}_1,...,{{\mathbf{f}}^{\prime }}_M\right],{\mathbf{S}}^{\prime \prime }=\left[{{\mathbf{f}}^{\prime \prime }}_1,...,{{\mathbf{f}}^{\prime \prime }}_M\right],$$

(54)

then according to (44) we have:

$$\mathbf{S}={\mathcal{F}}_c\left({\mathbf{X}}^l\right),{\mathbf{S}}^{\prime }={\mathcal{F}}_c\left({\mathbf{X}}^l\circ {\mathbf{D}}^l\right),{\mathbf{S}}^{\prime \prime }={\mathcal{F}}_c\left({\mathbf{X}}^l\circ {\mathbf{D}}^l\circ {\mathbf{D}}^l\right).$$

(55)

Let:

$$\begin{array}{cc}\hfill \mathbf{P}& =\left[{\mathbf{p}}_{11}^l,...,{\mathbf{p}}_{M1}^l,{\mathbf{p}}_{12}^l,...,{\mathbf{p}}_{M2}^l,......,{\mathbf{p}}_{1M}^l,...,{\mathbf{p}}_{MM}^l\right],\hfill \\ \hfill \mathbf{C}& =\left[{\mathbf{c}}_{11},...,{\mathbf{c}}_{M1},{\mathbf{c}}_{12},...,{\mathbf{c}}_{M2},......,{\mathbf{c}}_{1M},...,{\mathbf{c}}_{MM}\right],\hfill \\ \hfill {\mathbf{C}}^{\prime }& =\left[{{\mathbf{c}}^{\prime }}_{11},...,{{\mathbf{c}}^{\prime }}_{M1},{{\mathbf{c}}^{\prime }}_{12},...,{{\mathbf{c}}^{\prime }}_{M2},......,{{\mathbf{c}}^{\prime }}_{1M},...,{{\mathbf{c}}^{\prime }}_{MM}\right],\hfill \\ \hfill {\mathbf{P}}_d& =\left[{\mathbf{p}}_{11}^l,{\mathbf{p}}_{22}^l,...,{\mathbf{p}}_{MM}^l\right],\hfill \\ \hfill {\mathbf{C}}_d& =\left[{\mathbf{c}}_{11},{\mathbf{c}}_{22},...,{\mathbf{c}}_{MM}\right],\hfill \\ \hfill {{\mathbf{C}}^{\prime }}_d& =\left[{{\mathbf{c}}^{\prime }}_{11},{{\mathbf{c}}^{\prime }}_{22},...,{{\mathbf{c}}^{\prime }}_{MM}\right],\hfill \end{array}$$

(56)

Then the inverse Fourier transforms of $\mathbf{P},\mathbf{C},{\mathbf{C}}^{\prime }$ can be denoted as:

$$\mathbf{R}={\mathcal{F}}_c^{-1}\left(\mathbf{P}\right),\mathbf{H}={\mathcal{F}}_c^{-1}\left(\mathbf{C}\right),{\mathbf{H}}^{\prime }={\mathcal{F}}_c^{-1}\left({\mathbf{C}}^{\prime }\right).$$

(57)

According to (52), (56) and (57), the coefficients of the cubic equation can be reformulated as:

$$\begin{array}{cc}\hfill a& =4\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T\left({{\mathbf{C}}^{\prime }}_d^{*}\circ {{\mathbf{C}}^{\prime }}_d\right){\mathbf{1}}_M+\left(1-\lambda \right){\mathbf{w}}_t^T\left({{\mathbf{H}}^{\prime }}^{*}\circ {\mathbf{H}}^{\prime }\right){\mathbf{1}}_{M^2}\right),\hfill \\ \hfill b& =2\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T\left(2{{\mathbf{C}}^{\prime }}_d^{*}\circ {\mathbf{C}}_d+{\mathbf{C}}_d^{*}\circ {{\mathbf{C}}^{\prime }}_d\right){\mathbf{1}}_M+\left(1-\lambda \right){\mathbf{w}}_t^T\left(2{{\mathbf{H}}^{\prime }}^{*}\circ \mathbf{H}+{\mathbf{H}}^{*}\circ {\mathbf{H}}^{\prime }\right){\mathbf{1}}_{M^2}\right),\hfill \\ \hfill c& =2\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T\left(2{{\mathbf{C}}^{\prime }}_d^{*}\circ {\mathbf{P}}_d+{\mathbf{C}}_d^{*}\circ {\mathbf{C}}_d\right){\mathbf{1}}_M+\left(1-\lambda \right){\mathbf{w}}_t^T\left(2{{\mathbf{H}}^{\prime }}^{*}\circ \mathbf{R}+{\mathbf{H}}^{*}\circ \mathbf{H}\right){\mathbf{1}}_{M^2}\right),\hfill \\ \hfill d& =2\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T\left({\mathbf{C}}_d^{*}\circ {\mathbf{P}}_d\right){\mathbf{1}}_M+\left(1-\lambda \right){\mathbf{w}}_t^T\left({\mathbf{H}}^{*}\circ \mathbf{R}\right){\mathbf{1}}_{M^2}\right).\hfill \end{array}$$

(58)

3.3. 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:

$${\mathbf{d}}^{l+1}=-{\mathbf{g}}^{l+1}+\frac{{\left({\mathbf{g}}^{l+1}-{\mathbf{g}}^l\right)}^T{\mathbf{g}}^{l+1}}{{∥{\mathbf{g}}^l∥}^2}{\mathbf{d}}^l,$$

(59)

where ${\mathbf{g}}^l$ and ${\mathbf{d}}^l$ are the gradient vector and the direction vector, respectively. For the problem here, ${\mathbf{g}}^l$ and ${\mathbf{d}}^l$ are defined as:

$${\mathbf{g}}^l={\left[{\left({\mathbf{g}}_1^l\right)}^T,...,{\left({\mathbf{g}}_M^l\right)}^T\right]}^T,{\mathbf{d}}^l={\left[{\left({\mathbf{d}}_1^l\right)}^T,...,{\left({\mathbf{d}}_M^l\right)}^T\right]}^T,$$

(60)

where ${\mathbf{g}}_m^l={\nabla }_{{\varphi }_m}{J}_T\left({\mathbf{X}}^l\right),m=1,...,M$. Since the step size is an approximate value, ${\left({\mathbf{g}}^{l+1}\right)}^T{\mathbf{d}}^l$ is not equal to zero. Thus, ${\mathbf{d}}^{l+1}$ may not be descendent, i.e.,

$${\left({\mathbf{g}}^{l+1}\right)}^T{\mathbf{d}}^{l+1}=-{∥{\mathbf{g}}^{l+1}∥}^2+\frac{{\left({\mathbf{g}}^{l+1}-{\mathbf{g}}^l\right)}^T{\mathbf{g}}^{l+1}}{{∥{\mathbf{g}}^l∥}^2}{\left({\mathbf{g}}^{l+1}\right)}^T{\mathbf{d}}^l<0$$

(61)

is not always satisfied. For guaranteeing that the searching direction is descendant, we adopt the following modified direction:

$$\begin{array}{cc}\hfill {\tilde{\mathbf{d}}}^{l+1}& =-{\mathbf{g}}^{l+1}+\frac{{\left({\mathbf{g}}^{l+1}-{\mathbf{g}}^l\right)}^T{\mathbf{g}}^{l+1}}{{∥{\mathbf{g}}^l∥}^2}{\mathbf{d}}^l,\hfill \\ \hfill {\mathbf{d}}^{l+1}& =\left\{\begin{array}{cc}{\tilde{\mathbf{d}}}^{l+1},\hfill & {\left({\mathbf{g}}^{l+1}\right)}^T{\tilde{\mathbf{d}}}^{l+1}<0\\ -{\mathbf{g}}^{l+1},\hfill & {\left({\mathbf{g}}^{l+1}\right)}^T{\tilde{\mathbf{d}}}^{l+1}\ge 0\end{array}\right..\hfill \end{array}$$

(62)

Actually, (59) can also be expressed as:

$${\mathbf{D}}^{l+1}=-{\mathbf{G}}^{l+1}+\frac{{\left({\mathbf{g}}^{l+1}-{\mathbf{g}}^l\right)}^T{\mathbf{g}}^{l+1}}{{∥{\mathbf{g}}^l∥}^2}{\mathbf{D}}^l.$$

(63)

Since

$$\begin{array}{cc}\hfill {∥{\mathbf{g}}^l∥}^2& ={\displaystyle \sum _{m=1}^M}{∥{\mathbf{g}}_m^l∥}^2={∥{\mathbf{G}}^l∥}_F^2,\hfill \\ \hfill {\left({\mathbf{g}}^{l+1}-{\mathbf{g}}^l\right)}^T{\mathbf{g}}^{l+1}& ={∥{\mathbf{G}}^{l+1}∥}_F^2-{\mathbf{1}}_N^T\left({\mathbf{G}}^{l+1}\circ {\mathbf{G}}^l\right){\mathbf{1}}_M,\hfill \end{array}$$

(64)

and (63) can be rewritten as:

$${\mathbf{D}}^{l+1}=-{\mathbf{G}}^{l+1}+\frac{{∥{\mathbf{G}}^{l+1}∥}_F^2-{\mathbf{1}}_N^T\left({\mathbf{G}}^{l+1}\circ {\mathbf{G}}^l\right){\mathbf{1}}_M}{{∥{\mathbf{G}}^l∥}_F^2}{\mathbf{D}}^l.$$

(65)

Thus, (62) can be expressed as the following matrix form:

$$\begin{array}{cc}\hfill {\tilde{\mathbf{D}}}^{l+1}& =-{\mathbf{G}}^{l+1}+\frac{{∥{\mathbf{G}}^{l+1}∥}_F^2-{\mathbf{1}}_N^T\left({\mathbf{G}}^{l+1}\circ {\mathbf{G}}^l\right){\mathbf{1}}_M}{{∥{\mathbf{G}}^l∥}_F^2}{\mathbf{D}}^l,\hfill \\ \hfill {\mathbf{D}}^{l+1}& =\left\{\begin{array}{c}{\tilde{\mathbf{D}}}^{l+1},{\mathbf{1}}_N^T\left({\mathbf{G}}^{l+1}\circ {\tilde{\mathbf{D}}}^{l+1}\right){\mathbf{1}}_M<0\\ -{\mathbf{G}}^{l+1},{\mathbf{1}}_N^T\left({\mathbf{G}}^{l+1}\circ {\tilde{\mathbf{D}}}^{l+1}\right){\mathbf{1}}_M\ge 0\end{array}\right..\hfill \end{array}$$

(66)

On the basis of the above derivation, the gradient-based algorithm, which we call Gra-WeCorr-SFW, is summarized in Algorithm 1.

4. 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:

$$J_{CF}=\frac{1}{2N}{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{ij}^H{\mathbf{p}}_{ij}=\frac{1}{2N}{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{ii}^H{\mathbf{p}}_{jj}=\frac{1}{2N}{\displaystyle \sum _{i=1}^M}{\mathbf{p}}_{ii}^H{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{jj}.$$

(67)

Thus we can rewrite the problem (17) as:

$$J_T=\lambda {\displaystyle \sum _{i=1}^M}{\mathbf{p}}_{ii}^H{\mathbf{D}}_f{\mathbf{p}}_{ii}+{\lambda }_1{\displaystyle \sum _{i=1}^M}{\mathbf{p}}_{ii}^H{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{jj},$$

(68)

where ${\lambda }_1=\frac{1-\lambda }{2N}$. Like the derivation in Section 3, we deduce the gradient and the step size.

Algorithm 1: Gra-WeCorr-SFW.

Initialization: $l=0,\lambda ,M,N,{\mathbf{w}}_t,{\mathbf{w}}_f,{\mathbf{X}}^0=\left[{\mathbf{x}}_1^0,...,{\mathbf{x}}_M^0\right],$

$\mathbf{S}={\mathcal{F}}_c\left({\mathbf{X}}^0\right),{\mathbf{G}}^0,{\mathbf{D}}^0=-{\mathbf{G}}^0.$

Repeat

1: ${\mathbf{S}}^{\prime }={\mathcal{F}}_c\left({\mathbf{X}}^l\circ {\mathbf{D}}^l\right),{\mathbf{S}}^{\prime \prime }={\mathcal{F}}_c\left({\mathbf{X}}^l\circ {\mathbf{D}}^l\circ {\mathbf{D}}^l\right).$

2: Compute $\mathbf{P},\mathbf{C},{\mathbf{C}}^{\prime },{\mathbf{P}}_d,{\mathbf{C}}_d,{{\mathbf{C}}^{\prime }}_d$ according to (45).

3: $\mathbf{R}={\mathcal{F}}_c^{-1}\left(\mathbf{P}\right),\mathbf{H}={\mathcal{F}}_c^{-1}\left(\mathbf{C}\right),{\mathbf{H}}^{\prime }={\mathcal{F}}_c^{-1}\left({\mathbf{C}}^{\prime }\right).$

4: Compute coefficients $a,b,c,d$ according to (58).

5: Solve the cubic Equation (53), and then choose the

positive root that is closest to zero as the step size ${\mu }^l$.

6: ${\mathbf{X}}^{l+1}={\mathbf{X}}^l\circ e^{j{\mu }^l{\mathbf{D}}^l},\mathbf{S}={\mathcal{F}}_c\left({\mathbf{X}}^{l+1}\right).$

7: $\mathbf{Y}=\lambda \left({\mathbf{w}}_f{\mathbf{1}}_M^T\right)\circ \mathbf{S}\circ {\mathbf{S}}^{*}.$

8: Compute ${\mathbf{Y}}^{\prime }$ according to (32).

9: ${\mathbf{S}}_{rep}={\left[{\mathbf{S}}^T,...,{\mathbf{S}}^T\right]}_{2NM\times M}^T,{\mathbf{Y}}_1=\mathrm{ve}{\mathrm{c}}^{-1}\left(\left({\mathbf{S}}_{rep}\circ {\mathbf{Y}}^{\prime }\right){\mathbf{1}}_M\right).$

10: $\mathbf{Z}={\mathcal{F}}_c^{-1}\left(\mathbf{S}\circ \mathbf{Y}+{\mathbf{Y}}_1\right).$

11: ${\mathbf{G}}^{l+1}=8N\mathrm{Re}\left(-j{\left({\mathbf{X}}^{l+1}\right)}^{*}\circ \mathbf{Z}(1:N,:)\right).$

12: Compute the searching direction ${\mathbf{D}}^{l+1}$ according to (66).

13: $l=l+1.$

Until convergence

4.1. Phase Gradient

According to (68), we deduce the derivative $\partial J_T/\partial {\varphi }_{mn}$ as:

$$\frac{\partial J_T}{\partial {\varphi }_{mn}}=2\mathrm{Re}\left(\frac{\partial {\mathbf{p}}_{mm}^H}{\partial {\varphi }_{mn}}\left(\lambda {\mathbf{w}}_f\circ {\mathbf{p}}_{mm}+{\lambda }_1{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{jj}\right)\right).$$

(69)

Let:

$${{\mathbf{y}}^{\prime \prime }}_m=\lambda {\mathbf{w}}_f\circ {\mathbf{p}}_{mm}+{\lambda }_1{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{jj},$$

(70)

then according to (23), (69) can be rewritten as:

$$\frac{\partial J_T}{\partial {\varphi }_{mn}}=2\mathrm{Re}\left({\displaystyle \sum _{k=1}^{2N}}\frac{\partial \left(p_{mm}^{*}(k)\right)}{\partial {\varphi }_{mn}}{y^{\prime \prime }}_m(k)\right)=4\mathrm{Re}\left(-j{x}_{mm}^{*}{\displaystyle \sum _{k=1}^{2N}}{f}_m(k){y^{\prime \prime }}_m(k)e^{j\frac{2\pi }{2N}nk}\right).$$

(71)

By stacking (71) in a vector, the gradient is given by:

$${\nabla }_{{\mathsf{\varphi }}_m}{J}_T=8N\mathrm{Re}\left(-j{\mathbf{x}}_m^{*}\circ {{\mathbf{z}}^{\prime }}_m\left(1:N\right)\right),$$

(72)

where ${{\mathbf{z}}^{\prime }}_m={\mathcal{F}}^{-1}({\mathbf{f}}_m\circ {{\mathbf{y}}^{\prime \prime }}_m)$. Define ${\mathbf{Y}}_2=\left[{{\mathbf{y}}^{\prime \prime }}_1,...,{{\mathbf{y}}^{\prime \prime }}_M\right],{\mathbf{Z}}_1=\left[{{\mathbf{z}}^{\prime }}_1,...,{{\mathbf{z}}^{\prime }}_M\right]$, then it is easy to verify that:

$$\begin{array}{cc}\hfill {\mathbf{Y}}_2& =\lambda \left({\mathbf{w}}_f{\mathbf{1}}_M^T\right)\circ \mathbf{S}\circ {\mathbf{S}}^{*}+{\lambda }_1\left(\mathbf{S}\circ {\mathbf{S}}^{*}\right){\mathbf{1}}_M{\mathbf{1}}_M^T,\hfill \\ \hfill {\mathbf{Z}}_1& ={\mathcal{F}}_c^{-1}\left(\mathbf{S}\circ {\mathbf{Y}}_2\right).\hfill \end{array}$$

(73)

Thus, the gradient matrix can be expressed as:

$$\mathbf{G}=8N\mathrm{Re}\left(-j{\mathbf{X}}^{*}\circ {\mathbf{Z}}_1\left(1:N,:\right)\right).$$

(74)

4.2. Step Size

Similarly to (41), the linear search problem here can be denoted as:

$$\underset{\mu >0}{min}h(\mu )=\lambda {\displaystyle \sum _{i=1}^M}{\left({\mathbf{p}}_{ii}^{l+1}\right)}^H{\mathbf{D}}_f{\mathbf{p}}_{ii}^{l+1}+{\lambda }_1{\displaystyle \sum _{i=1}^M}{\left({\mathbf{p}}_{ii}^{l+1}\right)}^H{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{jj}^{l+1}.$$

(75)

By taking the derivative of $h(\mu )$, we have:

$$\frac{\partial h}{\partial \mu }=2\mathrm{Re}\left(\lambda {\displaystyle \sum _{i=1}^M}{\left(\frac{\partial {\mathbf{p}}_{ii}^{l+1}}{\partial \mu }\right)}^H{\mathbf{D}}_f{\mathbf{p}}_{ii}^{l+1}+{\lambda }_1{\left({\displaystyle \sum _{i=1}^M}\frac{\partial {\mathbf{p}}_{ii}^{l+1}}{\partial \mu }\right)}^H{\displaystyle \sum _{j=1}^M}{\mathbf{p}}_{jj}^{l+1}\right).$$

(76)

Then, the approximate derivative can be obtained by replacing ${\mathbf{p}}_{pq}^{l+1}$ with ${\overline{\mathbf{p}}}_{pq}^{l+1}$ defined in (46):

$$\begin{array}{cc}\hfill \frac{\partial h}{\partial \mu }\approx & 2\mathrm{Re}\left(\lambda {\displaystyle \sum _{i=1}^M}{\left({\mathbf{c}}_{ii}+2{{\mathbf{c}}^{\prime }}_{ii}\mu \right)}^H{\mathbf{D}}_f\left({\mathbf{p}}_{ii}^l+{\mathbf{c}}_{ii}\mu +{{\mathbf{c}}^{\prime }}_{ii}{\mu }^2\right)\right.\hfill \\ & +\left.{\lambda }_1{\displaystyle \sum _{i=1}^M}{\left({\mathbf{c}}_{ii}+2{{\mathbf{c}}^{\prime }}_{ii}\mu \right)}^H{\displaystyle \sum _{j=1}^M}\left({\mathbf{p}}_{jj}^l+{\mathbf{c}}_{jj}\mu +{{\mathbf{c}}^{\prime }}_{jj}{\mu }^2\right)\right).\hfill \end{array}$$

(77)

According the definition (56), we have:

$${\displaystyle \sum _{i=1}^M}{\mathbf{p}}_{ii}^l={\mathbf{P}}_d{\mathbf{1}}_M,{\displaystyle \sum _{i=1}^M}{\mathbf{c}}_{ii}={\mathbf{C}}_d{\mathbf{1}}_M,{\displaystyle \sum _{i=1}^M}{{\mathbf{c}}^{\prime }}_{ii}={{\mathbf{C}}^{\prime }}_d{\mathbf{1}}_M,$$

(78)

where ${\mathbf{P}}_d,{\mathbf{C}}_d$ and ${{\mathbf{C}}^{\prime }}_d$ can be expressed as the following equalities according to (45):

$${\mathbf{P}}_d=\mathbf{S}\circ {\mathbf{S}}^{*},{\mathbf{C}}_d=-2\mathrm{Im}\left({\mathbf{S}}^{\prime }\circ {\mathbf{S}}^{*}\right),{{\mathbf{C}}^{\prime }}_d=-\left(\mathrm{Re}\left({\mathbf{S}}^{\prime \prime }\circ {\mathbf{S}}^{*}\right)-{\mathbf{S}}^{\prime }\circ {{\mathbf{S}}^{\prime }}^{*}\right).$$

(79)

Thus, (77) can be written more compactly as:

$$\frac{\partial h}{\partial \mu }\approx a_1{\mu }^3+b_1{\mu }^2+c_1\mu +d_1,$$

(80)

where:

$$\begin{array}{cc}\hfill a_1& =4\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T\left({{\mathbf{C}}^{\prime }}_d^{*}\circ {{\mathbf{C}}^{\prime }}_d\right){\mathbf{1}}_M+{\lambda }_1{\left({{\mathbf{C}}^{\prime }}_d{\mathbf{1}}_M\right)}^H\left({{\mathbf{C}}^{\prime }}_d{\mathbf{1}}_M\right)\right),\hfill \\ \hfill b_1& =2\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T\left(2{{\mathbf{C}}^{\prime }}_d^{*}\circ {\mathbf{C}}_d+{\mathbf{C}}_d^{*}\circ {{\mathbf{C}}^{\prime }}_d\right){\mathbf{1}}_M+{\lambda }_1\left(2{\left({{\mathbf{C}}^{\prime }}_d{\mathbf{1}}_M\right)}^H\left({\mathbf{C}}_d{\mathbf{1}}_M\right)+{\left({\mathbf{C}}_d{\mathbf{1}}_M\right)}^H\left({{\mathbf{C}}^{\prime }}_d{\mathbf{1}}_M\right)\right)\right),\hfill \\ \hfill c_1& =2\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T\left(2{{\mathbf{C}}^{\prime }}_d^{*}\circ {\mathbf{P}}_d+{\mathbf{C}}_d^{*}\circ {\mathbf{C}}_d\right){\mathbf{1}}_M+{\lambda }_1\left(2{\left({{\mathbf{C}}^{\prime }}_d{\mathbf{1}}_M\right)}^H\left({\mathbf{P}}_d{\mathbf{1}}_M\right)+{\left({\mathbf{C}}_d{\mathbf{1}}_M\right)}^H\left({\mathbf{C}}_d{\mathbf{1}}_M\right)\right)\right),\hfill \\ \hfill d_1& =2\mathrm{Re}\left(\lambda {\mathbf{w}}_f^T\left({\mathbf{C}}_d^{*}\circ {\mathbf{P}}_d\right){\mathbf{1}}_M+{\lambda }_1{\left({\mathbf{C}}_d{\mathbf{1}}_M\right)}^H\left({\mathbf{P}}_d{\mathbf{1}}_M\right)\right).\hfill \end{array}$$

(81)

By solving the cubic equation $a_1{\mu }^3+b_1{\mu }^2+c_1\mu +d_1=0$, the step size can be easily obtained.

4.3. 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 ${\mathbf{Y}}^{\prime }$ 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.

Algorithm 2: Gra-Corr-SFW.

Initialization: $l=0,\lambda ,{\lambda }_1,M,N,{\mathbf{w}}_t,{\mathbf{w}}_f,{\mathbf{X}}^0=\left[{\mathbf{x}}_1^0,...,{\mathbf{x}}_M^0\right],$

$\mathbf{S}={\mathcal{F}}_c\left({\mathbf{X}}^0\right),{\mathbf{G}}^0,{\mathbf{D}}^0=-{\mathbf{G}}^0.$

Repeat

1: ${\mathbf{S}}^{\prime }={\mathcal{F}}_c\left({\mathbf{X}}^l\circ {\mathbf{D}}^l\right),{\mathbf{S}}^{\prime \prime }={\mathcal{F}}_c\left({\mathbf{X}}^l\circ {\mathbf{D}}^l\circ {\mathbf{D}}^l\right).$

2: ${\mathbf{P}}_d=\mathbf{S}\circ {\mathbf{S}}^{*},{\mathbf{C}}_d=-2\mathrm{Im}\left({\mathbf{S}}^{\prime }\circ {\mathbf{S}}^{*}\right),{{\mathbf{C}}^{\prime }}_d=-\left(\mathrm{Re}\left({\mathbf{S}}^{\prime \prime }\circ {\mathbf{S}}^{*}\right)-{\mathbf{S}}^{\prime }\circ {{\mathbf{S}}^{\prime }}^{*}\right).$

3: Compute the coefficients $a_1,b_1,c_1,d_1$ according to (81).

4: Solve the cubic Equation (80), and then choose the positive root

which is closest to zero as the step size ${\mu }^l$.

5: ${\mathbf{X}}^{l+1}={\mathbf{X}}^l\circ e^{j{\mu }^l{\mathbf{D}}^l},\mathbf{S}={\mathcal{F}}_c\left({\mathbf{X}}^{l+1}\right)$

6: ${\mathbf{Y}}_2=\lambda \left({\mathbf{w}}_f{\mathbf{1}}_M^T\right)\circ \mathbf{S}\circ {\mathbf{S}}^{*}+{\lambda }_1\left(\mathbf{S}\circ {\mathbf{S}}^{*}\right){\mathbf{1}}_M{\mathbf{1}}_M^T.$

7: ${\mathbf{Z}}_1={\mathcal{F}}_c^{-1}\left(\mathbf{S}\circ {\mathbf{Y}}_2\right).$

8: ${\mathbf{G}}^{l+1}=8N\mathrm{Re}\left(-j{\left({\mathbf{X}}^{l+1}\right)}^{*}\circ {\mathbf{Z}}_1\left(1:N,:\right)\right).$

9: Compute the searching direction ${\mathbf{D}}^{l+1}$ according to (66).

10: $l=l+1.$

Until convergence

To analyze the convergence speed,

Table 1. The per iteration computational complexities of different algorithms.

Algorithm	Number of (I)FFT	Complexity
MM-WeCorr-acc	$6M^2+8M$	$\mathcal{O}((12M+16)MNlog2N)$
MM-Corr-acc	$4M$	$\mathcal{O}(8MNlog2N)$
MDISAA-SFW	$3M$	$\mathcal{O}(6MNlog2N)$
Gra-WeCorr-SFW	$5M^2+6M$	$\mathcal{O}((10M+12)MNlog2N)$
Gra-Corr-SFW	$4M$	$\mathcal{O}(8MNlog2N)$

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.

5. 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.

5.1. 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:

$${\epsilon }_r^l=\frac{h({\mu }^l)-h({\dot{\mu }}^l)}{h(0)-h({\dot{\mu }}^l)},$$

(82)

where ${\mu }^l$ denotes the approximate step size of the l-th iteration, which is obtained via the Taylor series expansion, and ${\dot{\mu }}^l$ denotes the optimal step size obtained by the searching method. The smaller the value of ${\epsilon }_r^l$ is, the closer the approximate step size is to the optimal step size. Thus, ${\epsilon }_r^l$ can be used to express the accuracy of the approximate step size. When ${\epsilon }_r^l<1$, we have $h({\mu }^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}$:

$$P_{re}=max\left({\epsilon }_r^l,l=1,...,N_I\right),$$

(83)

where $N_I$ is the iteration number.

To simulate ${\epsilon }_r^l$ 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. Simulation parameters of different design problems.

	$\mathit{\lambda }$	Correlation weights ${\{{\mathit{w}}_{\mathit{k}}\}}_{\mathit{k}=1-\mathit{N}}^{\mathit{N}-1}$	Stopband ${\mathbf{\Omega }}_{\mathit{f}}$
Case 1${}^a$	0	$w_k=1,k=1-N,...,N-1$
Case 2${}^b$	0	$w_k=\left\{\begin{array}{cc}1,\hfill & k=0\hfill \\ 1,\hfill & k\in [-80,-51]\cup [51,80]\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$
Case 3${}^c$	0.9	$w_k=1,k=1-N,...,N-1$	$\begin{array}{c}(0.04,0.21),(0.23,0.25),(0.28,0.37),\hfill \\ (0.39,0.49),(0.52,0.56)\hfill \end{array}$

a: CSS design; b: waveform design with zero correlation zone; c: OWS design with stopband constraint.

. 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.

5.2. 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 $\mathrm{CISL}\le \epsilon$ as the stopping criterion for all of the algorithms in this experiment. The weight coefficient $\lambda$ 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:

$$r_{sum}(k)=20{log}_{10}\frac{\left|{\displaystyle \sum _{m=1}^M}{r}_{mm}(k)\right|}{{\displaystyle \sum _{m=1}^M}{r}_{mm}(0)},k=1-N,...,N-1.$$

(84)

From Figure 3, we can see that the autocorrelation sum of the CSS is a Delta function. Additionally, with the decrease of $\epsilon$, the sidelobes of the CSS are getting lower and lower.

Next, we compare the performance between the proposed Gra-Corr-SFW and the state-of-the-art MM-Corr-acc by designing waveform sets with different waveform number M and waveform length N. For both algorithms, we choose M and N as follows:

$$M=3k,N=256k,k=1,...,5.$$

(85)

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 $\mathrm{CISL}\le 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:

$$P_{sum}=max\left(r_{sum}(k),k=1,...,N-1\right).$$

(86)

Then, we choose $\epsilon ={10}^{-1}$ for the stopping criterion. The performance parameters $\left(N_I,t,P_{sum}\right)$ of these two algorithms are provided in

Table 3. The comparison of the performance parameters between MM-Corr-acc and Gra-Corr-SFW.

	MM-Corr-acc			Gra-Corr-SFW
$(M,N)$	$t(\mathrm{s})$	$N_I$	$P_{sum}(\mathrm{dB})$	$t(\mathrm{s})$	$N_I$	$P_{sum}(\mathrm{dB})$
$M=3,N=256$	1.249	1687	−63.8	0.071	239	−65.9
$M=6,N=512$	0.914	496	−68.7	0.110	132	−71.4
$M=9,N=768$	1.558	460	−72.0	0.178	131	−74.7
$M=12,N=1024$	2.531	484	−74.4	0.265	132	−77.5
$M=15,N=1280$	4.483	587	−76.1	0.343	138	−79.1

. As can be seen from this Table, since the algorithms stop the iteration when $\mathrm{CISL}\le 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 $\mathcal{O}(8MNlog2N)$. 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).

5.3. 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 $\psi$ in (4) as the stopping criterion, i.e., $\psi \le \epsilon$. 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 $\psi$ 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 $\psi$ is less than ${10}^{-11}$, the convergence speed of the MM-WeCorr-acc decreases, especially for a large $(M,N)$ pair.

Further,

Table 4. The comparison of the performance parameters between MM-WeCorr-acc and Gra-WeCorr-SFW $(\epsilon ={10}^{-10})$.

	MM-WeCorr-acc			Gra-WeCorr-SFW
$(M,N)$	$t(\mathrm{s})$	$N_I$	$P_{zcz}(\mathrm{dB})$	$t(\mathrm{s})$	$N_I$	$P_{zcz}(\mathrm{dB})$
$M=3,N=256$	0.584	335	−167.5	0.248	164	−168.2
$M=6,N=512$	2.691	407	−178.1	1.348	176	−179.5
$M=9,N=768$	10.167	445	−184.5	5.149	178	−186.3
$M=12,N=1024$	37.760	486	−189.7	15.106	179	−191.5
$M=15,N=1280$	84.759	549	−193.5	30.683	182	−195.1

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:

$$P_{zcz}=max\left(20{log}_{10}\left|\frac{r_{ij}(k)}{N}\right|,i,j=1,...,M,k\in [51,80]\right).$$

(87)

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.

5.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, ${∥{\mathbf{X}}^{l+1}-{\mathbf{X}}^l∥}_F\le {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:

$$\begin{array}{cc}\hfill P_{ac}& =max\left(20{log}_{10}\left|\frac{r_{ii}(k)}{N}\right|,i=1,...,M,k=2,...,N-1\right),\hfill \\ \hfill P_{cc}& =max\left(20{log}_{10}\left|\frac{r_{ij}(k)}{N}\right|,i,j=1,...,M,i\ne j,k=0,...,N-1\right),\hfill \\ \hfill P_{sp}& =max\left(P_{stop}^i,i=1,...,M\right),\hfill \\ \hfill ISL& ={\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{\genfrac{}{}{0pt}{}{k=-N+1}{k\ne 0}}^{N-1}}{\left|r_{ii}(k)\right|}^2+{\displaystyle \sum _{i=1}^M}{\displaystyle \sum _{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^M}{\displaystyle \sum _{k=-N+1}^{N-1}}{\left|r_{ij}(k)\right|}^2,\hfill \end{array}$$

(88)

where $P_{stop}^i$ denotes the peak stopband power of the waveform ${\mathbf{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\in \{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. The comparison of the performance parameters between MDISAA-SFW and Gra-Corr-SFW.

Algorithm	N	${\mathbf{P}}_{\mathbf{ac}}(\mathbf{dB})$	${\mathbf{P}}_{\mathbf{cc}}(\mathbf{dB})$	${\mathbf{P}}_{\mathbf{sp}}(\mathbf{dB})$	$\mathbf{ISL}$	${\mathbf{N}}_{\mathbf{I}}$	$\mathbf{t}(\mathbf{s})$
MDISAA-SFW	$N=256$	−13.6	−13.9	−17.5	893,573	11932	2.512
	$N=512$	−13.9	−16.6	−17.9	3,532,447	18,602	6.045
	$N=1024$	−14.1	−19.3	−18.2	14,013,417	27,268	14.080
Gra−Corr−SFW	$N=256$	−12.5	−15.9	−20.3	832,482	1216	0.274
	$N=512$	−12.9	−18.6	−22.2	3,321,729	2100	1.013
	$N=1024$	−13.3	−21.1	−23.8	13,237,903	3577	3.033

. 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 $\lambda$. 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 $\lambda$, the $P_{ac}$ and $P_{cc}$ of these two algorithms change slowly, while the $P_{sp}$ is sensitive to the change of $\lambda$. 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 $\lambda$, we can choose a relatively large $\lambda$ (e.g., $\lambda \in (0.85,0.95)$) for the Gra-Corr-SFW to obtain both good correlation and stopband properties.

6. 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.