Efficient Low-PAR Waveform Design Method for Extended Target Estimation Based on Information Theory in Cognitive Radar

Abstract

This paper addresses the waveform design problem of cognitive radar for extended target estimation in the presence of signal-dependent clutter, subject to a peak-to-average power ratio (PAR) constraint. Owing to this kind of constraint and the convolution operation of the waveform in the time domain, the formulated optimization problem for maximizing the mutual information (MI) between the target and the received signal is a complex non-convex problem. To this end, an efficient waveform design method based on minimization–maximization (MM) technique is proposed. First, by using the MM approach, the original non-convex problem is converted to a convex problem concerning the matrix variable. Then a trick is used for replacing the matrix variable with the vector variable by utilizing the properties of the Toeplitz matrix. Based on this, the optimization problem can be solved efficiently combined with the nearest neighbor method. Finally, an acceleration scheme is used to improve the convergence speed of the proposed method. The simulation results illustrate that the proposed method is superior to the existing methods in terms of estimation performance when designing the constrained waveform.

1. Introduction

Cognitive radar (CR) is a new intelligent closed-loop radar system that can perceive the surrounding complicated electromagnetism environment in real time and make reasoning decisions on this basis [1,2]. In CR, adaptive transmitted waveform design based on the perceived prior knowledge of environment and target is one of the key technologies which can significantly improve the performance of target detection, parameter estimation, recognition, and tracking in complicated environments [3].

For each of these missions, there are corresponding valid waveform design methods [4]. When designing waveform for target estimation, minimum mean squared error (MMSE) criterion [5], the minimum Cramer-Rao lower bound (CRLB) criterion [6,7] and maximum mutual information (MI) [5,8,9] criterion are usually selected. However, the CRLB criterion is only suitable for situations in which the target information is unknown. Meanwhile, maximization of MI and the minimization of the MMSE lead to the same solution when the target information is known.

The authors in [10] proposed an estimation waveform design method based on MI in noise, in which the original non-convex problem was converted to a convex problem by using the convex optimization method. Considering the signal-dependent clutter, the waveform design methods for target estimation based on MI and MMSE were proposed in [11,12], respectively. It is worth noting that the envelope constraint on the transmitted waveform was not considered in any of those studies, which made it difficult to meet the hardware constraints and maximize the power efficiency [13]. For this reason, unimodular or low peak-to-average power ratio (PAR) waveform is always applied in radar systems [13,14]. Nevertheless, unimodular waveform may lead to the degradation of waveform performance [15]. To tackle this problem, some researchers had used a more general low-PAR constraint to replace the unimodular constraint to further improve the waveform performance [6,7,14,15,16,17,18].

In [6], subject to the constraint, a frequency domain based PAR waveform design method with the MI criterion was proposed. However, the algorithm in [6] can only find the envelope of optimal waveform spectrum. Since the spectral phase cannot be determined, the number of time domain signals that satisfy the unique waveform spectrum magnitude is infinite. Therefore, it may result in a decline of the waveform performance when the frequency domain waveform spectrum is transformed to the time domain waveform [19]. To solve this problem, an algorithm based on the sequence linear programming (SLP) in time domain was proposed in [20], and the original non-convex problem was converted a convex problem which could be solved efficiently. It is worth noting that the optimization problem in [20] is the relaxation of the original problem so the synthesized waveform may be still the suboptimal solution.

As is well known, the minorization–maximization (MM) method is a powerful optimization technique to solve the hard problem that is difficult to tackle directly [21]. The core principle of MM is to transform the original problem into a series of simple problems which can be tackled efficiently and converge to the stationary optimal solution of the original problem [22]. Motivated by the ascent property and superior convergence of MM, it has been applied in many fields [23,24,25,26,27].

In this paper, we propose an efficient low-PAR cognitive waveform design method based on the MM approach for target estimation, which is directly studied in the time domain. Based on the MM approach, the original non-convex problem is converted to a convex problem with respect to (w.r.t.) a matrix variable. To reduce computation cost, the convex problem is further converted to quadratic programming (QP) problem w.r.t. a vector variable by utilizing the properties of the Toeplitz matrix. Based on this, the QP problem is converted to a simple convex problem which can be tackled efficiently by using the nearest neighbor method. Finally, the squared iterative methods (SQUAREM) is used to improve the convergence speed of the proposed method. The simulation results demonstrate that the synthesized waveform can be obtained efficiently within the given low-PAR range and the proposed method has better estimation performance than the existing methods.

The remainder of the paper is organized as follows: Section 2 gives the baseband radar signal model. In Section 3, the optimal criterion based MI is formulated, and an efficient low-PAR waveform design method based on MM is proposed. A detailed performance analysis of the proposed method is provided in Section 4. Section 5 presents our simulation results. Finally, the conclusion is summarized in Section 6.

Notation: Scalars are represented by italic letters, vectors and matrices are denoted by boldface lowercase and uppercase letters, respectively. The superscripts in ${(·)}^{\mathrm{T}}$ and ${(·)}^{\mathrm{H}}$ represent the transpose and Hermitian transpose operations, respectively. $A(m,n)$ denotes the element located in the $m$th row and $n$th column of $A$.$\mathrm{F}(·)$ denotes the Toeplitz matrix mapping function of a vector, $\Re (·)$, $\Im (·)$, $\left|\hspace{0.17em}·\hspace{0.17em}\right|$ and $\Vert \hspace{0.17em}·\hspace{0.17em}\Vert$ represent the real part, imaginary part, modulus and 2-norm of a complex scalar/vector/matrix, respectively. $ℂ$ is the set of complex-valued number. The symbol ‘$\otimes$’ and ‘*’ denote the Kronecker product and the convolution operation, respectively. Finally, $\mathcal{C}\mathcal{N}(0,A)$ denotes a circular symmetric complex Gaussian distribution with zero mean and the covariance matrix $A$.

2. Signal Model

In this paper, we consider the waveform design of cognitive radar for target estimation in the presence of signal-dependent clutter. The scattering characteristic of the target is represented by the target impulse response (TIR) [28], and the signal-dependent clutter is represented by the clutter impulse response (CIR) [29]. Generally, the prior knowledge of the target and environment (noise and clutter) can be obtained by some cognitive methods [30,31] and is assumed to be known when designing waveform for simplicity. It is assumed that the influence of sidelobes has been mitigated by sidelobe blanking technology in front of the receiver. Meanwhile, we focus on the analysis of single-input single-output radar in this paper which can be straightforwardly extended to multiple-input multiple-output radar case. Then, the discrete baseband signal model is shown in Figure 1.

As illustrated in Figure 1, the target and clutter are modelled by the finite impulse response filter, and the waveform is assumed to be energy-limited. $\mathit{s}\in {ℂ}^{N_{\mathit{s}}\times 1}$ denotes a transmitted waveform with length $N_{\mathit{s}}$ $\mathit{t}\in {ℂ}^{N_{\mathit{t}}\times 1}$ and $c\in {ℂ}^{N_c\times 1}$ denote the TIR and CIR, respectively. According to the radar signal model in [8], it is assumed that $N_{\mathit{t}}=N_c$ to simplify the derivations. If $N_c>N_{\mathit{t}}$, it is necessary to apply a zero-filling operation to the TIR to make the TIR and CIR sampling points equal. $n\in {ℂ}^{N_n\times 1}$ denotes the sum of the noise and the interference, $N_n=N_{\mathit{s}}+N_{\mathit{t}}-1$. Let $N=N_n$, and x is the echo with length $N_x=N$. Then, the model can be described as:

$$x\hspace{0.17em}=\mathit{t}\ast \mathit{s}+c\ast \mathit{s}+n=T\mathit{s}+C\mathit{s}+n=S\mathit{t}+Sc+n$$

(1)

where $S\mathit{t}=T\mathit{s}$ and $Sc=C\mathit{s}$ can be obtained due to the reciprocity of the convolution operation. The convolution matrices $\mathit{S}$ and $\mathit{T}$ are Toeplitz matrices corresponding to $\mathit{s}$ and $\mathit{t}$, respectively. We use the function ‘$\mathrm{F}(·)$’ represents their mapping relationship in this paper, i.e., $\mathit{T}=\mathrm{F}\left(\mathit{t}\right)$, $\mathit{S}=\mathrm{F}\left(\mathit{s}\right)$. Taking transmitted waveform as an example, the convolution matrix $\mathit{S}$ can be written as:

$$S=\left[\begin{array}{cccc}\mathit{s}\left(1\right)& 0& \cdots & 0\\ ⋮& \mathit{s}\left(1\right)& \ddots & ⋮\\ \mathit{s}(N_{\mathit{s}})& ⋮& \ddots & 0\\ 0& \mathit{s}(N_{\mathit{s}})& \ddots & \mathit{s}\left(1\right)\\ ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & 0& \mathit{s}(N_{\mathit{s}})\end{array}\right]\in {ℂ}^{\left(N_{\mathit{s}}+N_{\mathit{t}}-1\right)\times N_{\mathit{t}}}.$$

(2)

3. Waveform Design Method

In this section, we utilize the MM technique to solve the estimation waveform design problem based on information theory.

3.1. MM Method

The MM method refers to the minorization–maximization method, which can transform the original complex problem into a series of simple problems that can be tackled efficiently and converge to the stationary optimal solution of the original problem [22]. Now we first give a brief description of MM. Consider a general maximization problem

$$\begin{array}{l}\underset{x}{\max }\hspace{0.17em}\hspace{0.17em}f\left(\mathit{x}\right)\\ \mathrm{s}.\mathrm{t}.\hspace{1em}x\in \mathsf{\Theta }\end{array}$$

(3)

where $f\left(\mathit{x}\right)$ is a function which is difficult to solve directly. Then, the approximate function $Y\left(\mathit{x};{\mathit{x}}^k\right)$ is commonly used to replace the original function $f\left(\mathit{x}\right)$. More precisely, MM can get the optimal solution ${\mathit{x}}^{k+1}$ of the $\left(k+1\right)\mathrm{th}$ iteration based on the known ${\mathit{x}}^k$ according to the following criterion.

$${\mathit{x}}^{k+1}\in \underset{\mathit{x}\in \mathsf{\Theta }}{\arg \hspace{0.17em}\max }\hspace{0.17em}\hspace{0.17em}Y\left(\mathit{x};{\mathit{x}}^k\right)$$

(4)

where $Y\left(\mathit{x};{\mathit{x}}^k\right)$ is said to minorize the function $f\left(\mathit{x}\right)$ at the point ${\mathit{x}}^k$, which satisfies

$$f\left(\mathit{x}\right)\ge Y\left(\mathit{x};{\mathit{x}}^k\right)\hspace{0.17em}\mathrm{for}\hspace{0.17em}\forall \mathit{x}\in \mathsf{\Theta },$$

(5)

$$f\left({\mathit{x}}^k\right)=Y\left({\mathit{x}}^k;{\mathit{x}}^k\right)$$

(6)

Then, it can be seen that the objective value is increased monotonically at every iteration, i.e.,

$$f\left({\mathit{x}}^{k+1}\right)\ge Y\left({\mathit{x}}^{k+1};{\mathit{x}}^k\right)\ge Y\left({\mathit{x}}^k;{\mathit{x}}^k\right)=f\left({\mathit{x}}^k\right)$$

(7)

The first inequality and the third equality hold due to the properties of (5) and (6), respectively. The second inequality holds according to (4). Next, MM is utilized to solve the estimation waveform based on information theory.

3.2. Problem Formulation

In this paper, the maximization of MI between the received signal and target is used as the optimization criterion for waveform design. According to [28], supposing t, $c$ and $n$ are mutually independent and $\mathit{t}\sim \mathcal{C}\mathcal{N}(0,R_{\mathit{t}})$, $\mathit{c}\sim \mathcal{C}\mathcal{N}(0,R_{\mathit{c}})$, and $\mathit{n}\sim \mathcal{C}\mathcal{N}(0,R_{\mathit{n}})$. Then, the MI of received signal $\mathit{x}$ and target $\mathit{t}$ can be formulated as [32]:

$$I\left(\mathit{x};\mathit{t}|S\right)=h\left(\mathit{x}|S\right)-h\left(\mathit{x}|\mathit{t},S\right)$$

(8)

where $h\left(\mathit{x}|S\right)$ denotes the entropy of received signal $\mathit{x}$ when the Toeplitz matrix of transmit waveform $S$ is known, and $h\left(\mathit{x}|\mathit{t},S\right)$ denotes the entropy of $\mathit{x}$ when $S$ and TIR are known.

For the given $S$, $\mathit{t}$ and $\mathit{x}$ obey the joint Gaussian distribution which can be expressed as:

$$\left(\begin{array}{c}\mathit{t}\\ \mathit{x}\end{array}\right)\sim \mathcal{C}\mathcal{N}\left[\left(\begin{array}{c}{0}_{{\mathit{N}}_{\mathit{t}}}\\ 0_{{\mathit{N}}_{\mathit{x}}}\end{array}\right),\left(\begin{array}{cc}{R}_{\mathit{t}}& R_{\mathit{t}}{S}^{\mathrm{H}}\\ S{R}_{\mathit{t}}& S{R}_{\mathit{t}}{S}^{\mathrm{H}}+S{R}_{\mathit{c}}{S}^{\mathrm{H}}+R_{\mathit{n}}\end{array}\right)\right]$$

(9)

Let $R_{\mathit{x}}=S{R}_{\mathit{t}}{S}^{\mathrm{H}}+S{R}_{\mathit{c}}{S}^{\mathrm{H}}+R_{\mathit{n}}$, we can get

$$p\left(\mathit{x}|S\right)=\frac{1}{{\pi }^{N_{\mathit{x}}}\det \left(R_{\mathit{x}}\right)}\exp \left[-{\mathit{x}}^{\mathrm{H}}{R}_{\mathit{x}}^{-1}\mathit{x}\right]$$

(10)

Then the entropy $h\left(\mathit{x}|S\right)$ and $h\left(\mathit{x}|\mathit{t},S\right)$ can be expressed as:

$$\begin{array}{ll}h\left(\mathit{x}|S\right)& =-{\displaystyle \int \mathit{p}\left(\mathit{x}|S\right)\ln \mathit{p}\left(\mathit{x}|S\right)d\mathit{x}}\\ & =\ln \det \left(S{R}_{\mathit{t}}{S}^{\mathrm{H}}+S{R}_{\mathit{c}}{S}^{\mathrm{H}}+R_{\mathit{n}}\right)+N_{\mathit{x}}\ln \pi +N_{\mathit{x}}\end{array}$$

(11)

$$\begin{array}{ll}h\left(\mathit{x}|\mathit{t},S\right)& =-{\displaystyle \int p\left(\mathit{x}|\mathit{t},S\right)}\ln \mathit{p}\left(\mathit{x}|\mathit{t},S\right)d\mathit{x}\\ & =\ln \det \left(S{R}_{\mathit{c}}{S}^{\mathrm{H}}+R_{\mathit{n}}\right)+N_{\mathit{x}}\ln \pi +N_{\mathit{x}}\end{array}$$

(12)

Bring (11) and (12) into (8), the objective function can be written as:

$$\mathit{I}\left(\mathit{x},\mathit{t}|\mathit{s}\right)=\ln \det \left[I_{{\mathit{N}}_{\mathit{x}}}+S{R}_{\mathit{t}}{S}^{\mathrm{H}}{\left(S{R}_{\mathit{c}}{S}^{\mathrm{H}}+R_{\mathit{n}}\right)}^{-1}\right]$$

(13)

To meet the hardware constraints and maximize the power efficiency, the PAR constraint must be considered. Let the total energy of the transmitted waveform be $E_{\mathit{s}}$. Without any loss of generality, it can be assumed that $E_{\mathit{s}}=N_{\mathit{s}}$. Then, PAR can be defined as:

$$\mathrm{PAR}(\mathit{s})\triangleq \frac{{\max }_j{\left|\mathit{s}(j)\right|}^2}{\frac{1}{N_{\mathit{s}}}{{\displaystyle {\sum }_{n=1}^{N_{\mathit{s}}}\left|\mathit{s}(j)\right|}}^2}=\underset{j}{\max }{\left|\mathit{s}\left(j\right)\right|}^2\le \eta ,\hspace{1em}\eta \in \left[1,N_{\mathit{s}}\right]$$

(14)

where $\mathit{s}\left(j\right)$ is the $j\mathrm{th}$ element of $\mathit{s}$, and $\eta$ is a predefined parameter that denotes the maximum allowed PAR. Note that the PAR constraint is equivalent to a unimodular constraint when $\eta =1$, while it becomes a redundant constraint when $\eta =N_{\mathrm{s}}$.

Then, the optimization problem can be formulated as:

$$\mathit{P}\{\begin{array}{l}\underset{S}{\max \hspace{0.17em}}\ln \det \left[I_{N_{\mathit{x}}}+S{R}_{\mathit{t}}{S}^{\mathrm{H}}{\left(S{R}_{\mathit{c}}{S}^{\mathrm{H}}+R_{\mathit{n}}\right)}^{-1}\right]\\ \mathrm{s}.\mathrm{t}.\hspace{0.17em}\hspace{0.17em}{\mathit{s}}^{\mathrm{H}}\mathit{s}\le E_{\mathit{s}}\\ \hspace{1em}\hspace{1em}{\left|\mathit{s}\left(j\right)\right|}^2\hspace{0.17em}\le \eta ,\hspace{1em}j=1,2,\cdots ,N_{\mathit{s}}\\ \hspace{1em}\hspace{1em}\mathit{S}=\mathrm{F}\left(\mathit{s}\right).\end{array}$$

(15)

It can be seen that the objective function in problem $\mathit{P}$ is non-convex, the two quadratic inequality constraints are nonhomogeneous [33]. So $\mathit{P}$ is a non-convex problem which is difficult to solve. Therefore, we need to transform $\mathit{P}$ into a convex problem.

3.3. Waveform Design

The key to solving $\mathit{P}$ is to convert the non-convex objective function into a convex function. First, the objective function in (13) can be reformulated as:

$$\ln \det \left[I_{{\mathit{N}}_{\mathrm{x}}}+S{R}_{\mathrm{t}}{S}^{\mathrm{H}}{\left(S{R}_{\mathrm{c}}{S}^{\mathrm{H}}+R_{\mathrm{n}}\right)}^{-1}\right]=\ln \det \left[I_{{\mathit{N}}_{\mathrm{t}}}+R_{\mathrm{t}}^{1/2}{S}^{\mathrm{H}}{\left(S{R}_{\mathrm{c}}{S}^{\mathrm{H}}+R_{\mathrm{n}}\right)}^{-1}S{R}_{\mathrm{t}}^{1/2}\right]$$

(16)

According to the Woodbury identity [34], we can have

$$\begin{array}{l}\ln \det \left[I_{{\mathit{N}}_{\mathrm{t}}}+R_{\mathrm{t}}^{1/2}{S}^{\mathrm{H}}{\left(S{R}_{\mathrm{c}}{S}^{\mathrm{H}}+R_{\mathrm{n}}\right)}^{-1}S{R}_{\mathrm{t}}^{1/2}\right]\\ =-\ln \det \left[I_{{\mathit{N}}_{\mathrm{t}}}-R_{\mathrm{t}}^{1/2}{S}^{\mathrm{H}}{\left[\left(S{R}_{\mathrm{c}}{S}^{\mathrm{H}}+R_{\mathrm{n}}\right)+S{R}_{\mathrm{t}}{S}^{\mathrm{H}}\right]}^{-1}S{R}_{\mathrm{t}}^{1/2}\right]\\ =-\ln \det \left[I_{{\mathit{N}}_{\mathrm{t}}}-R_{\mathrm{t}}^{1/2}{S}^{\mathrm{H}}{\left(S{R}_{\mathrm{tc}}{S}^{\mathrm{H}}+R_{\mathrm{n}}\right)}^{-1}S{R}_{\mathrm{t}}^{1/2}\right]\end{array}$$

(17)

where $R_{\mathrm{tc}}=R_{\mathit{t}}+R_{\mathit{c}}$. Let $\mathit{J}=\left[\begin{array}{l}{\mathit{I}}_{N_{\mathit{t}}}\\ 0_{N\times N_{\mathit{t}}}\end{array}\right]\in {ℂ}^{\left(N_{\mathit{t}}+N\right)\times N_{\mathit{t}}}$ and

$$\mathit{V}=\left[\begin{array}{cc}{\mathit{I}}_{N_{\mathit{t}}}& {\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{S}}^{\mathrm{H}}\\ \mathit{S}{\mathit{R}}_{\mathit{t}}^{1/2}& \mathit{S}{\mathit{R}}_{\mathit{t}\mathit{c}}^{1/2}{\mathit{S}}^{\mathrm{H}}+{\mathit{R}}_{\mathit{n}}\end{array}\right]\in {ℂ}^{\left(N_{\mathit{t}}+N\right)\times \left(N_{\mathit{t}}+N\right)}$$

(18)

According to the inversion identity of block matrix [34], the expression in (17) can be reformulated as:

$$I_{{\mathit{N}}_{\mathit{t}}}-R_{\mathit{t}}^{1/2}{S}^{\mathrm{H}}{\left(S{R}_{\mathit{t}\mathit{c}}{S}^{\mathrm{H}}+R_{\mathit{n}}\right)}^{-1}S{R}_{\mathit{t}}^{1/2}={\left({\mathit{J}}^{\mathrm{H}}{\mathit{V}}^{-1}\mathit{J}\right)}^{-1}$$

(19)

So, we can recast the objective function of $\mathit{P}$ as follows:

$$\ln \hspace{0.17em}\det \left[I_{{\mathit{N}}_{\mathit{t}}}+R_{\mathit{t}}^{1/2}{S}^{\mathrm{H}}{\left(S{R}_{\mathit{c}}{S}^{\mathrm{H}}+R_{\mathit{n}}\right)}^{-1}S{R}_{\mathit{t}}^{1/2}\right]=\ln \hspace{0.17em}\det \left[{\mathit{J}}^{\mathrm{H}}{\mathit{V}}^{-1}\mathit{J}\right]$$

(20)

The following lemma provides a way to solve the non-convex design problem by utilizing the MM approach.

Lemma 1:

For any full-column rank matrix$\mathit{J}\in {ℂ}^{n\times m}\left(m\le n\right)$, if$\mathit{V}\in {ℂ}^{n\times n}$is a positive definite matrix, then$\ln \hspace{0.17em}\det \left[{\mathit{J}}^{\mathrm{H}}{\mathit{V}}^{-1}\mathit{J}\right]$is convex w.r.t.$\mathit{V}$.

Then the proof of Lemma 1 can be found in [35]. Based on Lemma 1, we can find that $\mathit{J}$ is a full-column rank matrix and $\mathit{V}$ is a positive definite matrix, so we can know that $\ln \hspace{0.17em}\det \left[{\mathit{J}}^{\mathrm{H}}{\mathit{V}}^{-1}\mathit{J}\right]$ is convex w.r.t. $\mathit{V}$. Therefore, by using its tangent plane [33] with a given $\mathit{V}$, this term can be minorized as:

$$\ln \det \left[{\mathit{J}}^{\mathrm{H}}{\mathit{V}}^{-1}\mathit{J}\right]\ge \ln \det \left[{\mathit{J}}^{\mathrm{H}}{\left({\mathit{V}}^k\right)}^{-1}\mathit{J}\right]+\mathrm{tr}\left[{\mathit{Q}}^{\left(k\right)}\left(\mathit{V}-{\mathit{V}}^k\right)\right]$$

(21)

where ${\mathit{Q}}^k=-{\left({\mathit{V}}^k\right)}^{-1}\mathit{J}\left[{\mathit{J}}^{\mathrm{H}}{\left({\mathit{V}}^k\right)}^{-1}\mathit{J}\right]{\mathit{J}}^{\mathrm{H}}{\left({\mathit{V}}^k\right)}^{-1}$, the right-hand side of (21) is the first-order approximation of $\ln \det \left[{\mathit{J}}^{\mathrm{H}}{\mathit{V}}^{-1}\mathit{J}\right]$ for a given ${\mathit{V}}^k$ at $k\mathrm{th}$ iteration. Let ${\mathit{Q}}^k=\left[\begin{array}{cc}{\mathit{Q}}_{11}^k& {\mathit{Q}}_{12}^k\\ {\left({\mathit{Q}}_{12}^k\right)}^{\mathrm{H}}& {\mathit{Q}}_{22}^k\end{array}\right]$ with the same partitioning as that of $\mathit{V}$ in (18), where ${\mathit{Q}}_{11}^k\in {ℂ}^{N_{\mathit{t}}\times N_{\mathit{t}}}$, ${\mathit{Q}}_{12}^k\in {ℂ}^{N_{\mathit{t}}\times N}$, and ${\mathit{Q}}_{22}^k\in {ℂ}^{N\times N}$. Now the MM is applied to minorize the function of the left side in (21). Then, the right-hand side of (21) can be rewritten as:

$$Q_0^k+\mathrm{tr}\left[{\mathit{Q}}_{12}^k\mathit{S}{\mathit{R}}_{\mathit{t}}^{1/2}\right]+\mathrm{tr}\left[{\left({\mathit{Q}}_{12}^k\right)}^{\mathrm{H}}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{S}}^{\mathrm{H}}\right]+\mathrm{tr}\left[{\mathit{Q}}_{22}^k\mathit{S}{\mathit{R}}_{\mathit{t}\mathit{c}}{\mathit{S}}^{\mathrm{H}}\right]$$

(22)

where $Q_0^k=\ln \det \left[{\mathit{J}}^{\mathrm{H}}{\left({\mathit{V}}^k\right)}^{-1}\mathit{J}\right]+\mathrm{tr}\left[{\mathit{Q}}_{11}^k+{\mathit{Q}}_{22}^k{\mathit{R}}_{\mathit{n}}\right]$ is the constant term.

Ignoring the constant term, and using the identity that $\mathrm{tr}\left[\mathit{A}\mathit{B}\right]=\mathrm{tr}\left[\mathit{B}\mathit{A}\right]$ [34], the problem $\mathit{P}$ can be recast as:

$${\mathit{P}}_1\{\begin{array}{l}\underset{S}{\max \hspace{0.17em}}\mathrm{tr}\left[\mathit{S}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{Q}}_{12}^k+{\left({\mathit{Q}}_{12}^k\right)}^{\mathrm{H}}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{S}}^{\mathrm{H}}+{\mathit{Q}}_{22}^k\mathit{S}{\mathit{R}}_{\mathit{t}\mathit{c}}{\mathit{S}}^{\mathrm{H}}\right]\\ \mathrm{s}.\mathrm{t}.\hspace{0.17em}{\mathit{s}}^{\mathrm{H}}\mathit{s}\le E_{\mathit{s}}\\ \hspace{1em}\hspace{1em}{\left|\mathit{s}\left(j\right)\right|}^2\hspace{0.17em}\le \eta ,\hspace{1em}j=1,2,\cdots ,N_{\mathit{s}}\\ \hspace{1em}\hspace{1em}\mathit{S}=\mathrm{F}\left(\mathit{s}\right).\end{array}$$

(23)

It can be seen that the objective function of (23) is still non-convex [33]. Then, let ${\lambda }_{\min }$ denotes the smallest eigenvalue of $\mathit{S}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{Q}}_{12}^k+{\mathit{Q}}_{21}^k{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{S}}^{\mathrm{H}}+{\mathit{Q}}_{22}^k\mathit{S}{\mathit{R}}_{\mathit{t}\mathit{c}}{\mathit{S}}^{\mathrm{H}}$, so the problem ${\mathit{P}}_1$ can be rewritten as:

$${\mathit{P}}_2\{\begin{array}{l}\underset{\mathit{S},{\lambda }_{\min }}{\max \hspace{0.17em}}\hspace{1em}{\lambda }_{\min }\\ \mathrm{s}.\mathrm{t}.\hspace{1em}\mathit{S}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{Q}}_{12}^k+{\left({\mathit{Q}}_{12}^k\right)}^{\mathrm{H}}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{S}}^{\mathrm{H}}+{\mathit{Q}}_{22}^k\mathit{S}{\mathit{R}}_{\mathit{t}\mathit{c}}{\mathit{S}}^{\mathrm{H}}\succcurlyeq {\lambda }_{\min }{\mathit{I}}_N\\ \hspace{1em}\hspace{1em}{\mathit{s}}^{\mathrm{H}}\mathit{s}\le E_{\mathit{s}}\\ \hspace{1em}\hspace{1em}{\left|\mathit{s}\left(j\right)\right|}^2\hspace{0.17em}\le \eta ,\hspace{1em}j=1,2,\cdots ,N_{\mathit{s}}\\ \hspace{1em}\hspace{1em}\mathit{S}=\mathrm{F}\left(\mathit{s}\right).\end{array}$$

(24)

Then, the first constraint can be converted to a convex set by utilizing the Schur complement theorem [36] which is defined as follows:

Lemma 2.

(Schur complement theorem):Let$\mathit{A}=\left[\begin{array}{cc}{\mathit{A}}_{11}& {\mathit{A}}_{12}\\ {\mathit{A}}_{12}^{\mathrm{H}}& {\mathit{A}}_{22}\end{array}\right]$, then we can get$\mathit{A}\succ 0$if and only if${\mathit{A}}_{22}\succ 0$and${\mathit{A}}_{11}-{\mathit{A}}_{12}{\mathit{A}}_{22}^{-1}{\mathit{A}}_{12}^{\mathrm{H}}\succ 0$.

According to the Schur complement theorem, the first constraint of this problem is equivalent to

$$\left[\begin{array}{cc}-{\lambda }_{\min }{\mathit{I}}_N+\mathit{S}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{Q}}_{12}^k+{\left({\mathit{Q}}_{12}^k\right)}^{\mathrm{H}}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{S}}^{\mathrm{H}}& {\left({\mathit{Q}}_{22}^k\right)}^{1/2}\mathit{S}\\ {\mathit{S}}^{\mathrm{H}}{\left({\mathit{Q}}_{22}^k\right)}^{1/2}& -{\mathit{R}}_{\mathit{t}\mathit{c}}^{-1}\end{array}\right]\succcurlyeq 0$$

(25)

To make the optimization problem more intuitional, ${\mathit{P}}_2$ can be recast as:

$${\mathit{P}}_3\{\begin{array}{l}\underset{\mathit{S},\hspace{0.17em}{\lambda }_{\min }}{\max \hspace{0.17em}}\hspace{0.17em}{\lambda }_{\min }\\ \mathrm{s}.\mathrm{t}.\hspace{1em}\left[\begin{array}{cc}-{\lambda }_{\min }{\mathit{I}}_N+\mathit{S}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{Q}}_{12}^k+{\left({\mathit{Q}}_{12}^k\right)}^{\mathrm{H}}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{S}}^{\mathrm{H}}& {\left({\mathit{Q}}_{22}^k\right)}^{1/2}\mathit{S}\\ {\mathit{S}}^{\mathrm{H}}{\left({\mathit{Q}}_{22}^k\right)}^{1/2}& -{\mathit{R}}_{\mathit{t}\mathit{c}}^{-1}\end{array}\right]\succcurlyeq \mathbf{0}\\ \hspace{1em}\hspace{1em}{\left[\mathit{S}(m:m+N_{\mathit{s}}-1,m)\right]}^{\mathrm{H}}\mathit{S}(m:m+N_{\mathit{s}}-1,m)\le E_{\mathit{s}},\hspace{1em}m=1,2,\cdots ,N_{\mathit{t}}\\ \hspace{1em}\hspace{1em}{\left|\mathit{S}\left(i,j\right)\right|}^2\hspace{0.17em}\le \eta ,\hspace{1em}i=1,2,\cdots ,N,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,2,\cdots ,N_{\mathit{t}}.\end{array}$$

(26)

where ‘$\mathit{S}(m:m+N_{\mathit{s}}-1,m)$’ represents the elements in the $m\mathrm{th}$ column, and the $m\mathrm{th}$ to $\left(m+N_{\mathit{s}}-1\right)\mathrm{th}$ rows of $\mathit{S}$. We can see that the first inequality constraint is the linear inequality with regard to matrix variable $\mathit{S}$. According to [33], any line is affine, so the first inequality constraint is a convex set. Then the second constraint belongs to a Euclidean ball, and the third constraint belongs to a norm ball, and both of them are convex set [33]. Therefore, the constraints of ${\mathit{P}}_3$ are convex sets with regard to matrix $\mathit{S}$.

In addition, it is obvious that the objective function of ${\mathit{P}}_3$ is an auxiliary variable which is also convex. Hence, ${\mathit{P}}_3$ is a convex problem with regard to matrix variable $\mathit{S}$ and it can be solved by applying the interior point method [33] with CVX toolbox [37]. However, it has an approximate computational complexity of $O\left(N_{\mathit{s}}^{3.5}\right)$ [38] at each iteration, which may bring a high computation cost especially when $N_{\mathit{s}}$ is large. Therefore, a fast optimization method is needed.

3.4. A Fast Optimization Method

To reduce the cost of computation, we should convert the original problem ${\mathit{P}}_1$ to a form that is easier to solve. First, we can convert the matrix variable $\mathit{S}$ to the form of vector. Then, the identities that $\mathrm{tr}\left(\mathit{A}\mathit{B}\right)=\mathrm{ve}{\mathrm{c}}^{\mathrm{T}}\left({\mathit{B}}^{\mathrm{T}}\right)\mathrm{vec}\left(\mathit{A}\right)$ and $\mathrm{tr}\left(\mathit{A}\mathit{B}\mathit{C}\mathit{D}\right)=\mathrm{ve}{\mathrm{c}}^{\mathrm{T}}\left({\mathit{D}}^{\mathrm{T}}\right)\left({\mathit{C}}^{\mathrm{T}}\otimes \mathit{A}\right)\mathrm{vec}\left(\mathit{B}\right)$ [39] can be used to recast the objective function in (23), which can be rewritten as:

$$\begin{array}{l}\mathrm{tr}\left[\mathit{S}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{Q}}_{12}^k+{\left({\mathit{Q}}_{12}^k\right)}^{\mathrm{H}}{\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{S}}^{\mathrm{H}}+{\mathit{Q}}_{22}^k\mathit{S}{\mathit{R}}_{\mathit{t}\mathit{c}}{\mathit{S}}^{\mathrm{H}}\right]\\ =\mathrm{ve}{\mathrm{c}}^{\mathrm{T}}\left({\left({\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{Q}}_{12}^k\right)}^{\mathrm{T}}\right)\mathrm{vec}\left(\mathit{S}\right)+\mathrm{ve}{\mathrm{c}}^{\mathrm{T}}\left({\left({\mathit{S}}^{\mathrm{H}}\right)}^{\mathrm{T}}\right)\mathrm{vec}\left({\left({\mathit{Q}}_{12}^k\right)}^{\mathrm{H}}{\mathit{R}}_{\mathit{t}}^{1/2}\right)+\mathrm{ve}{\mathrm{c}}^{\mathrm{T}}\left({\left({\mathit{S}}^{\mathrm{H}}\right)}^{\mathrm{T}}\right)\left({\mathit{R}}_{\mathit{t}\mathit{c}}^{\mathrm{T}}\otimes {\mathit{Q}}_{22}^k\right)\mathrm{vec}\left(\mathit{S}\right)\\ ={\left({\mathit{u}}^k\right)}^{\mathrm{T}}\mathrm{vec}\left(\mathit{S}\right)+\mathrm{ve}{\mathrm{c}}^{\mathrm{T}}\left({\left({\mathit{S}}^{\mathrm{H}}\right)}^{\mathrm{T}}\right)v^k+\mathrm{ve}{\mathrm{c}}^{\mathrm{T}}\left({\left({\mathit{S}}^{\mathrm{H}}\right)}^{\mathrm{T}}\right){\mathit{G}}^k\mathrm{vec}\left(\mathit{S}\right)\end{array}$$

(27)

where

$$\begin{array}{l}{\mathit{u}}^k=\mathrm{vec}\left({\left({\mathit{R}}_{\mathit{t}}^{1/2}{\mathit{Q}}_{12}^k\right)}^{\mathrm{T}}\right)\in {ℂ}^{N_{\mathit{t}}N\times 1},\\ v^k=\mathrm{vec}\left({\left({\mathit{Q}}_{12}^k\right)}^{\mathrm{H}}{\mathit{R}}_{\mathit{t}}^{1/2}\right)\in {ℂ}^{N_{\mathit{t}}N\times 1},\\ {\mathit{G}}^k=\left({\mathit{R}}_{\mathit{t}\mathit{c}}^{\mathrm{T}}\otimes {\mathit{Q}}_{22}^k\right)\in {ℂ}^{N_{\mathit{t}}N\times N_{\mathit{t}}N}\end{array}$$

(28)

For further simplification, considering that $\mathit{S}$ is a convolution matrix with Toeplitz structure (shown in (2)) which consists of $\mathit{s}$, we can recast the right hand side of (27) as:

$${\left({\tilde{\mathit{u}}}^k\right)}^{\mathrm{T}}\mathit{s}+{\mathit{s}}^{\mathrm{H}}{\tilde{\mathit{v}}}^k+{\mathit{s}}^{\mathrm{H}}{\tilde{\mathit{G}}}^k\mathit{s}$$

(29)

where

$$\begin{array}{l}{\tilde{\mathit{u}}}^k={\displaystyle {\sum }_{i=1}^{N_{\mathit{t}}}{\tilde{\mathit{u}}}_i^k}\in {ℂ}^{N_{\mathit{s}}\times 1},\\ {\tilde{\mathit{u}}}_i^k=\left[{\mathit{u}}^k\left(\left(i-1\right)\ast N+i\right),{\mathit{u}}^k\left(\left(i-1\right)\ast N+i+1\right),\cdots ,{\mathit{u}}^k\left(\left(i-1\right)\ast N+i+N_{\mathit{s}}-1\right)\right]\end{array}$$

(30)

$$\begin{array}{l}{\tilde{\mathit{v}}}^k={\displaystyle {\sum }_{i=1}^{N_{\mathit{t}}}{\tilde{\mathit{v}}}_i^k}\in {ℂ}^{N_{\mathit{s}}\times 1},\\ {\tilde{\mathit{v}}}_i^k=\left[{\mathit{v}}^k\left(\left(i-1\right)\ast N+i\right),{\mathit{v}}^k\left(\left(i-1\right)\ast N+i+1\right),\cdots ,{\mathit{v}}^k\left(\left(i-1\right)\ast N+i+N_{\mathit{s}}-1\right)\right]\end{array}$$

(31)

and ${\tilde{\mathit{G}}}^{\left(k\right)}={\displaystyle {\sum }_{i=1}^{N_{\mathit{t}}}{\displaystyle {\sum }_{j=1}^{N_{\mathit{t}}}{\mathit{B}}_{i,j}^k}}$, where ${\mathit{B}}_{i,j}^k\in {ℂ}^{N_{\mathit{s}}\times N_{\mathit{s}}}$ can be expressed as:

$${\mathit{B}}_{i,\hspace{0.17em}j}^k=\left[\begin{array}{ccc}{\mathit{G}}^k\left[\left(i-1\right)N+j\hspace{0.17em},\left(\hspace{0.17em}i-1\right)N+j\right]& \cdots & {\mathit{G}}^k\left[\left(i-1\right)N+j\hspace{0.17em},\left(\hspace{0.17em}i-1\right)N+j+N_{\mathit{s}}-1\right]\\ ⋮& \ddots & ⋮\\ {\mathit{G}}^k\left[\left(i-1\right)N+j+N_{\mathit{s}}-1\hspace{0.17em},\left(\hspace{0.17em}i-1\right)N+j\right]& \cdots & {\mathit{G}}^k\left[\left(i-1\right)N+j+N_{\mathit{s}}-1\hspace{0.17em},\left(\hspace{0.17em}i-1\right)N+j+N_{\mathit{s}}-1\right]\end{array}\right]$$

(32)

Then the optimization problem ${\mathit{P}}_1$ can be further rewritten as:

$${\mathit{P}}_4\{\begin{array}{l}\underset{\mathit{s}}{\max \hspace{0.17em}}\mathrm{Re}\left({\mathit{s}}^{\mathrm{H}}{\mathit{r}}^k\right)+{\mathit{s}}^{\mathrm{H}}{\tilde{\mathit{G}}}^k\mathit{s}\\ \mathrm{s}.\mathrm{t}.\hspace{1em}{\mathit{s}}^{\mathrm{H}}\mathit{s}\le E_{\mathit{s}}\\ \hspace{1em}\hspace{1em}{\left|\mathit{s}\left(j\right)\right|}^2\hspace{0.17em}\le \eta ,\hspace{1em}j=1,2,\cdots ,N_{\mathit{s}}.\end{array}$$

(33)

where ${\mathit{r}}^k={\tilde{\mathit{u}}}^k+{\tilde{\mathit{v}}}^k$. It can be seen that ${\mathit{P}}_4$ is a quadratically constrained quadratic program (QCQP) problem, which can be solved by the power method-like in [26]. More precisely, ${\mathit{P}}_4$ can be reformulated as:

$${\mathit{P}}_5\{\begin{array}{l}\underset{\mathit{s}}{\min \hspace{0.17em}}{\Vert \mathit{s}-a^k\Vert }_2^2\\ \mathrm{s}.\mathrm{t}.\hspace{1em}{\mathit{s}}^{\mathrm{H}}\mathit{s}\le E_{\mathit{s}}\\ \hspace{1em}\hspace{1em}{\left|\mathit{s}\left(j\right)\right|}^2\hspace{0.17em}\le \eta ,\hspace{1em}j=1,2,\cdots ,N_{\mathit{s}}\end{array}$$

(34)

where ${\mathit{a}}^k$ consists of the first $N_{\mathit{s}}$ entries of ${\tilde{\mathit{W}}}^k{\tilde{\mathit{s}}}^k$, ${\tilde{\mathit{s}}}^k={\left[{\left({\mathit{s}}^k\right)}^{\mathrm{T}},\hspace{0.17em}1\right]}^{\mathrm{T}}$, let

$${\mathit{W}}^k=\left[\begin{array}{cc}{\tilde{\mathit{G}}}^k& {\mathit{r}}^k\\ {\left({\mathit{r}}^k\right)}^{\mathrm{H}}& 0\end{array}\right]$$

(35)

Then ${\tilde{\mathit{W}}}^k={\mu }^k{\mathit{I}}_{N_{\mathit{s}}+1}-{\mathit{W}}^k$, where ${\mu }^k$ is a constant that is larger than the maximum eigenvalue of ${\mathit{W}}^k$ to make sure that ${\tilde{\mathit{W}}}^k$ is positive definite. It can be seen that ${\mathit{P}}_5$ is a convex problem which can be solved by using the interior point method. However, we can find that the form of objective function and constraints of ${\mathit{P}}_5$ are the same as the nearest neighbor method with a lower complexity of $O\left(N_{\mathit{s}}^2\right)$ [40]. Hence, ${\mathit{P}}_5$ can be solved efficiently.

Then, let $I^k$ denote the value of MI at $k\mathrm{th}$ iteration, $\tau$ the termination tolerance and $\gamma$ the maximum iterative number. According to the above steps, the proposed MM-based method is summarized in Box 1.

Box 1. The proposed minimization–maximization (MM)-based method for low-PAR estimation waveform design.

Step 0:

Set $k=0$, generate a random waveform ${\mathit{s}}^k$, initialize the $\tau$ and $\gamma$.

Step 1:

Use (28) to update ${\mathit{u}}^k$, ${\mathit{v}}^k$ and ${\mathit{G}}^k$;

Step 2:

Use (30) and (31) to update ${\tilde{\mathit{u}}}^k$ and ${\tilde{\mathit{v}}}^k$;

Step 3:

Use (32) to update ${\mathit{B}}_{i,j}^k$, ${\tilde{\mathit{G}}}^k={\displaystyle {\sum }_{i=1}^{N_{\mathit{t}}}{\displaystyle {\sum }_{i=1}^{N_{\mathit{t}}}{\mathit{B}}_{i,j}^k}}$;

Step 4:

${\mathit{r}}^k={\tilde{\mathit{u}}}^k+{\tilde{\mathit{v}}}^k$, use (35) to update ${\mathit{W}}^k$, ${\tilde{\mathit{W}}}^k={\mu }^k{\mathit{I}}_{N_{\mathit{s}}+1}-{\mathit{W}}^k$;

Step 5:

Get ${\mathit{a}}^k$ from the first $N_{\mathit{s}}$ entries of ${\tilde{\mathit{W}}}^k{\tilde{\mathit{s}}}^k$;

Step 6:

Solve ${\mathit{P}}_5$ to update ${\mathit{s}}^{k+1}$, set $k=k+1$;

Step 7:

Go back to step 1 until $\left|I^k-I^{k-1}\right|/I^k\le \tau$ or the iteration number is larger than $\gamma$.

3.5. Acceleration Scheme

For the MM method, the convergence speed depends mainly on the minorized function. It is worth noting that the minorized function of the proposed method (as shown in the right-hand side of (21)) may be relatively loose as a lower bound of the original function. Therefore, although the computational cost of the proposed method is low at each iteration, the convergence speed may still be slow. For the purpose of accelerating the convergence speed, the acceleration scheme (SQUAREM) in [41] is adopted in this paper. Then, we give the modified version of SQUAREM according to the optimization problem we meet.

Let ${\mathrm{L}}_{\mathrm{MM}}(·)$ denote the fixed-point map of the proposed MM-based method which can be described as:

$${\mathit{s}}^{k+1}={\mathrm{L}}_{\mathrm{MM}}\left({\mathit{s}}^k\right)$$

(36)

However, SQUAREM is not applicable to the case of the limited-energy and PAR constraints, and the monotonicity of the proposed method cannot be guaranteed by using SQUAREM. To this end, the first problem can be solved by using the nearest neighbor method to deal with ${\mathit{P}}_5$, and the second problem can be tackled by utilizing a backtracking strategy. Then, the acceleration scheme based on SQUAREM is summarized in Box 2.

Box 2. The Acceleration scheme based on squared iterative methods (SQUAREM).

Step 0:

Set k = 0, generate a random waveform s^k, initialize the τ and γ;

Step 1:

${\mathit{s}}_1={\mathrm{L}}_{\mathrm{MM}}\left({\mathit{s}}^k\right)$;

Step 2:

${\mathit{s}}_2={\mathrm{L}}_{\mathrm{MM}}\left({\mathit{s}}_1\right)$;

Step 3:

$\mathit{e}={\mathit{s}}_1-{\mathit{s}}^k$;

Step 4:

$\mathit{q}={\mathit{s}}_2-{\mathit{s}}_1-\mathit{e}$;

Step 5:

$\alpha =-\Vert \mathit{e}\Vert /\Vert \mathit{q}\Vert$;

Step 6:

${\mathit{a}}^k={\mathit{s}}^k-2\alpha \mathit{e}+{\alpha }^2\mathit{q}$;

Step 7:

Solve to P₅ update s^k+1;

Step 8:

while$I^{k+1}<I^k$ do

Step 9:

  $\alpha =\left(\alpha -1\right)/2$;

Step 10:

 ${\mathit{a}}^k={\mathit{s}}^k-2\alpha \mathit{e}+{\alpha }^2\mathit{q}$;

Step 11:

 Solve P₅ to update s^k+1;

Step 12:

end while

Step 13:

Set k = k + 1;

Step 14:

Go back to step 1 until $\left|I^k-I^{k-1}\right|/I^k\le \tau$ or the iteration number is larger than γ.

4. Performance Analysis

4.1. Convergence

Let $U\left(\mathit{s};{\mathit{s}}^k\right)$ represent the minorized function (i.e., the objective function of ${\mathit{P}}_4$) which can be written as:

$$U\left(\mathit{s};{\mathit{s}}^k\right)=\mathrm{Re}\left({\mathit{s}}^{\mathrm{H}}{\mathit{r}}^k\right)+{\mathit{s}}^{\mathrm{H}}{\tilde{\mathit{G}}}^k\mathit{s}$$

(37)

As is well known, ${\mathit{P}}_4$ is a QCQP problem. According to [42], it had been proved that $U\left({\mathit{s}}^{k+1};{\mathit{s}}^k\right)>U\left({\mathit{s}}^k;{\mathit{s}}^k\right)$ under the PAR and transmitted energy constraints. Then, we can have

$$I\left({\mathit{s}}^k\right)=U\left({\mathit{s}}^k;{\mathit{s}}^k\right)\le U\left({\mathit{s}}^{k+1};{\mathit{s}}^k\right)\le I\left({\mathit{s}}^{k+1}\right)$$

(38)

where $I\left({\mathit{s}}^k\right)$ denotes MI value of the original function at $k\mathrm{th}$ iteration. The first equality and the third inequality hold due to the properties of (6) and (5), respectively. Hence, we can know that the proposed method is monotonically increasing.

In addition, the waveform $\left\{{\mathit{s}}^k\right\}$ is energy-limited and its every point bounded with $\left|{\mathit{s}}^k\left(j\right)\right|\hspace{0.17em}\le \sqrt{\eta }\hspace{0.17em}\hspace{0.17em}\left(j=1,2,\cdots ,N_{\mathit{s}}\right)\hspace{0.17em}$. Therefore, according to the Theorem 2.17 in [42] and the fact that $I\left(\mathit{s}\right)$ and $U\left(\mathit{s};{\mathit{s}}^k\right)$ have the same gradient value when $\mathit{s}={\mathit{s}}^k$, we can know that at least one limit point exists and the MI of the synthesized waveform has its upper bound.

4.2. Computational Complexity

The proposed MM-based method converts the original problem into a simple problem which can be solved efficiently. In every iteration, updating ${\mathit{Q}}^k$, performing eigenvalue decomposition to obtain ${\mu }^k$, and solving ${\mathit{P}}_5$ with the nearest neighbor method, which have complexities of $O\left({\mathit{N}}_{\mathrm{s}}^2\right)$, $O\left({\mathit{N}}_{\mathrm{s}}^3\right)$ and $O\left({\mathit{N}}_{\mathrm{s}}^2\right)$, respectively. Therefore, the total computation complexity of the proposed method is $O\left({\mathit{N}}_{\mathrm{s}}^3+2{\mathit{N}}_{\mathrm{s}}^2\right)$ in every iteration.

5. Simulation Results

In this section, several numerical simulations are performed to demonstrate the performance of the proposed method. Assuming that the length of the transmitted waveform is $N_{\mathit{s}}=20$, the initial waveform ${\mathit{s}}_0$ is generated by a random phase-coded signal. The length of TIR and CIR are $N_{\mathit{t}}=N_{\mathit{c}}=30$. Meanwhile, both the target and clutter are mutually independent circular symmetric complex Gaussian random vector, i.e., $\mathit{t}\sim \mathcal{C}\mathcal{N}(0,{\mathit{R}}_{\mathit{t}})$, $\mathit{c}\sim \mathcal{C}\mathcal{N}(0,{\mathit{R}}_{\mathit{c}})$ [28], where ${\mathit{R}}_{\mathit{t}}={\sigma }_{\mathit{t}}^2{\mathit{R}}_{\mathit{t}0}$ and ${\mathit{R}}_{\mathit{c}}={\sigma }_{\mathit{c}}^2{\mathit{R}}_{\mathit{c}0}$. According to [5], ${\mathit{R}}_{\mathit{t}0}=U_{\mathit{t}}{\Lambda }_{\mathit{t}}{U}_{\mathit{t}}^{\mathrm{H}}$ and ${\mathit{R}}_{\mathit{c}0}=U_{\mathit{c}}{\Lambda }_{\mathit{c}}{U}_{\mathit{c}}^{\mathrm{H}}$ are normalized covariance matrix, where ${\Lambda }_{\mathit{t}}$ and ${\Lambda }_{\mathit{c}}$ have the same structure. Hence, taking the ${\Lambda }_{\mathit{t}}$ as an example, ${\Lambda }_{\mathit{t}}\in {ℂ}^{N_{\mathit{t}}\times N_{\mathit{t}}}$ is a diagonal matrix and $U_{\mathit{t}}$ is the $N_{\mathit{t}}\times N_{\mathit{t}}$ unitary discrete Fourier transform (DFT) matrix with its $\left(i,j\right)\mathrm{th}$ entry given by

$$\frac{1}{\sqrt{N_{\mathit{t}}}}\exp \left[\frac{-j2\pi \left(i-1\right)\left(j-1\right)}{N_{\mathit{t}}}\right],\forall i,j\in \left[1,N_{\mathit{t}}\right]$$

(39)

The noise is white Gaussian with zero mean and covariance matrix ${\mathit{R}}_{\mathit{n}}={\sigma }_{\mathit{n}}^2{\mathit{I}}_N$, where ${\sigma }_{\mathit{n}}^2=0.5$ denotes the variance of noise. Then, we perform 300 Monte Carlo trials for each combination of parameters and the termination tolerance $\tau ={10}^{-6}$. The MATLAB 2013b version is used to perform the simulations with a standard PC (CPU Core i5-3230M 2.6GHz and 4GB RAM).

5.1. Effectiveness Verification

In this subsection, we demonstrate the effectiveness of the proposed method. First, we give the typical set of eigenvalues for the normalized matrices ${\mathit{R}}_{\mathit{t}0}$ and ${\mathit{R}}_{\mathit{c}0}$ as shown in Figure 2.

Let the total energy of waveform $E_{\mathit{s}}=N_{\mathit{s}}$, ${\sigma }_{\mathrm{t}}^2=1$ and ${\sigma }_{\mathit{c}}^2=1$. Figure 3 shows the convergence of the proposed MM-based method. In addition, the upper bound is obtained by using the Lagrange multipliers method to solve ${\mathit{P}}_4$ under only the energy constraint. From Figure 3a ($\eta =N_{\mathit{s}}$ which is equivalent to the energy constraint) and Figure 3b ($\eta =1$ which is equivalent to the constant-modulus constraint), it can be seen that the accelerated case is much faster than the case without acceleration. Meanwhile, the MI of the synthesized waveform with $\eta =N_{\mathit{s}}$ can get the upper bound and the case with $\eta =1$ is about 0.07 away from the upper bound, this is because the large PAR value ($\eta =N_{\mathit{s}}$) has larger feasible set region than small PAR value ($\eta =1$).

Let the transmitted energy $E_{\mathit{s}}$ range from 1 to 30, ${\sigma }_{\mathrm{t}}^2=1$ and ${\sigma }_{\mathit{c}}^2=1$. Then, Figure 4a,b shows the estimation performance comparison of the proposed method, the Sequence Linear Programming-based Waveform Design algorithm (SLPWD) in [20] and the frequency domain-based Cognitive REceiver and Waveform design algorithm (CREW(fre)) in [6] versus the transmitted energy with $\eta =E_{\mathit{s}}$ and $\eta ={E_{\mathit{s}}/N}_{\mathit{s}}$ (i.e., constant-modulus constraint), respectively.

From Figure 4 we can see that the MI of proposed method is larger than that of other methods. Meanwhile, the MI of SLPWD algorithm in [20] is larger than CREW(fre) algorithm in [6], this is because CREW(fre) addresses the waveform design in frequency domain and the waveform spectrum does not contain the phase information, which may result in a decline of the waveform performance when the waveform spectrum is converted to time domain.

Let $E_{\mathit{s}}=N_{\mathit{s}}$, ${\sigma }_{\mathrm{t}}^2=1$, the clutter-noise-ration (CNR) ranges from −10 dB to 10 dB. Figure 5a,b shows the estimation performance of the proposed method, SLPWD in [20] and CREW(fre) in [6] versus CNR with $\eta =E_{\mathit{s}}$ and $\eta =1$, respectively. It can be seen that the MI of the proposed method is better than the other methods. Hence, it demonstrates the effectiveness of the proposed method.

Figure 6 shows the performance assessment for the estimation of target $\mathit{t}$ with the proposed method. For the optimal waveform based on maximizing MI, mean squared error (MSE) will be used for performance assessment. In addition, we also compared with the CRLB. The MSE and CRLB for transmitted waveform are derived in Appendix A. From Figure 6 we can see that the estimated performance gets closer to the CRLB as the transmitted energy increases. So, it verifies that the waveform generated by the proposed method has good estimation performance.

5.2. Influence of PAR

In this subsection, we discuss the influence of the PAR on the synthesized waveform. Figure 7 shows the MI of the synthesized waveforms under different PAR values. We can see that the curves can converge to their respective stationary values which become larger as $\eta$ increases, and this is because the feasible set region in ${\mathit{P}}_5$ becomes larger as $\eta$ increases. However, since the energy of the transmitted waveform is limited, the waveform performance has its upper bound. Then, we can also see that the curves can be monotonically convergent to the same stationary value and the curves almost overlap when $\eta \ge 1.5$. Figure 8 shows the real and imaginary parts of the waveforms under different PAR constraints. When $\eta =N_{\mathit{s}}$, the distribution radii of the corresponding points are large, which is not favorable for practical applications. In contrast, the results obtained with $\eta =1$ are unimodular and lie on the unit circle. Meanwhile, the distribution radii of the waveform with $\eta =2$ are close to those of the waveform with $\eta =1$, and the performance is very close to that of the waveform with $\eta =N_{\mathit{s}}$, as shown in Figure 7. This result indicates that the low-PAR waveform (for example $\eta =2$) not only meet the hardware constraints but also have better estimation performance than a unimodular waveform. Hence, the low-PAR waveform is more suitable for practical applications.

6. Conclusions

In this paper, we proposed an efficient low-PAR waveform design method of cognitive radar for the extended target estimation in the presence of signal-dependent clutter. To tackle the original non-convex problem in the time domain, an efficient method is proposed by using the MM technique. Meanwhile, to improve the convergence speed, an acceleration approach is given based on the SQUAREM. Numerical experiments demonstrate the effectiveness of the proposed method for the given PAR. Compared with the existing method, the proposed method demonstrates the advantage w.r.t. estimation performance. Moreover, the proposed method can be used in the waveform design of cognitive radar systems since the high computational efficiency will enable real-time waveform changes. Possible future research tracks include the extension to cases with low autocorrelation sidelobes and spectral constraints with imprecise prior information.