LPI Sequences Optimization Method against Summation Detector Based on FFT Filter Bank

Waveform design is a crucial factor in electronic surveillance (ES) systems. In this paper, we introduce an algorithm that designs a low probability of intercept (LPI) radar waveform. Our approach directly minimizes the detection probability of summation detectors based on FFT filter banks. The algorithm is derived from the general quadratic optimization framework, which inherits the monotonic properties of such methods. To expedite overall convergence, we have integrated acceleration schemes based on the squared iterative method (SQUAREM). Additionally, the proposed algorithm can be executed through fast Fourier transform (FFT) operations, enhancing computational efficiency. With some modifications, the algorithm can be adjusted to incorporate spectral constraints, increasing its flexibility. Numerical experiments indicate that our proposed algorithm outperforms existing ones in terms of both intercept properties and computational complexity.

Keywords: LPI; waveform design; unimodular sequences; propagation distance

During radar detection and tracking, the radar waveform is susceptible to detection, posing a significant threat to the radar's operation and survival [[1]]. As a result, radar systems continually demand advanced low-intercept technology. The objective is to develop a low probability of intercept (LPI) radar systems that can effectively detect targets while minimizing the likelihood of detection by ES systems [[3]]. This ensures enhanced security for both the radar system and its associated platforms.

Indeed, the concept of LPI radar was first developed in the 1980s [[4]]. Over the years, through iterative advancements, a plethora of low-interception technologies have progressively been implemented. Currently, the development of LPI radar primarily targets the two critical interception stages of the ES system: the power interception phase and the radar signature interception phase [[5]]. The power interception phase involves the ES's broad reception of all electromagnetic signals, detecting signals when their amplitude surpasses a certain threshold [[6]]. On the other hand, the radar signature interception phase refers to the process where the ES system extracts and identifies features from the detected signals [[7]].

The use of low-power continuous wave technology to diminish the chance of power interception is a superior low-interception method. It minimizes the peak of instantaneous power in the waveform by spreading the energy over an extended duration, thereby reducing the reconnaissance system's probability of interception [[8]]. Low-side-lobe antenna technology deliberately prevents the ES system from acquiring excess radiated energy, effectively preventing energy leakage and decreasing the interception likelihood [[10]]. Similarly, the dual/multi-station radar system reduces the interception probability by positioning the radar transmitter as far from the reconnaissance receiver as possible, thereby limiting the receiver's energy acquisition [[11]].

The frequency shortcut [[12]], pseudo-random [[14]], and complex waveform modulation techniques [[15]] are employed with the aim of reducing the interception probability of radar signatures. However, these techniques also inadvertently increase the processing complexity of the reconnaissance system. Furthermore, the issue of designing a low probability of identification (LPID) waveforms is being examined, with a specific focus on the aspect of time complexity [[16]].

The aforementioned tools for low-interception do not delve into the potential influence that waveform design could have on low-interception performance. In reality, waveform design is a crucial method for impacting ES systems [[18]]. The design of LPI radar waveforms has been employed as a countermeasure against ES systems as far back as the 1990s [[19]]. Among the various techniques, shortening the parameters of the waveform is an effective low-intercept method that targets the reduction in feature interception probability [[21]].

The swift fluctuation in radar waveform parameters makes it difficult for the receiver to accurately identify the relevant parameters of the waveform, which in turn obstructs effective classification and identification [[23]]. Furthermore, the advent of innovative radar systems has led to a heightened degree of versatility in waveform design [[24]]. The development of these cutting-edge radar systems primarily focuses on MIMO radar waveform design, utilizing waveform diversity to boost the adaptability of LPI radar waveform design. The incorporation of transmit beamforming further enhances the radar's stealth detection capabilities [[25]].

Nonetheless, as processing power improves, an increasing number of reconnaissance systems can intercept complex waveforms. They achieve this by analyzing and extracting suspicious features from the received signals [[26]]. Even under low signal-to-noise ratio conditions, contemporary low-interception systems are capable of identifying the potential characteristics of radar waveforms [[28]].

In recent years, a number of scholars have utilized Kullback–Leibler (KL) divergence for the measurement to represent the feature interception and identification process of radar missions and ES systems, which takes into account the detection characteristics under the joint entropy constraint for LPI waveform design, aiming to make the reconnaissance system lose more information [[29]]. However, simply utilizing the clustering and dispersion characteristics of the frequency domain distribution is not well adapted to the characteristics of modern wide-band ES systems [[31]].

To better obscure waveform features and enhance detection capabilities, a waveform design method was developed. This method conceals the spectral features of the waveform within the spectral features of the background noise. Under a specific received signal-to-noise ratio, the discrepancy between the periodic spectrum of the designed waveform and that of the Gaussian white noise is minimized. This minimization increases the complexity of feature extraction for reconnaissance receivers, thereby boosting the waveform's low-interception performance.

Waveform design methods aimed at reducing the probability of power interception in LPI radar initially distribute signal power over a wider frequency band using pulse compression theory [[33]]. The advent of new regime radars then amplified the signal's bandwidth, resulting in a wider signal bandwidth and a larger frequency range for energy distribution. Furthermore, MIMO radar waveform design minimizes the peak power captured by the interception system through the transmission beamforming technique [[35]]. The emergence of advanced frequency diversity array (FDA) radars [[37]], which utilize electronically scanned methods, further broadened the technical scope of LPI radar waveform design research.

Undeniably, diminishing the peak radiated power of passive reconnaissance systems is an efficacious method for low interception. However, this leads to a decline in detection performance as the energy radiated at the target diminishes. Therefore, target detection with low radiated power is typically accomplished through energy accumulation over the duration of the pulse train, which unfortunately compromises radar sensitivity.

Detecting signals that are embedded in the radio spectrum is a common challenge in various applications, including spectrum detectors in cognitive radio systems [[38]], cognitive radar systems [[39]], and ES systems [[40]]. More specifically, the challenge involves detecting signals within a spectrum using a receiver that has an instantaneous bandwidth (IBW).

One possible solution for this challenge is to use the total power radiometer [[41]], also known as the energy detector. This detector does not require a calculated frequency spectrum and compresses the entire band (i.e., the samples) into a single detection variable. This variable is then compared to a detection threshold. However, the frequency resolution of these detectors is limited to the entire IBW.

Another type of detector involves calculating the detection variable for each band in the digital spectrum, which is then compared to a threshold [[42]]. The frequency resolution of these detectors is limited to the bandwidth of the band.

We chose to use the FFT filter bank-based summation detector, which is a modified version of the basic FFT filter bank detector [[43]]. The FFT summation detector groups FFT bins to correspond with the desired channelization and estimates the power contained in each channel by summing the power computed for the individual bins.

In this paper, we present the MPI (monotonic minimizer for the probability of intercept) algorithm, which is designed to efficiently reduce the detection probability of summation detectors based on FFT filter banks in a monotonic manner. The MPI algorithm is derived by applying a general quadratic optimization framework method to the probability of intercept minimization problem, resulting in a simple closed-form solution at each iteration. The power method-like iteration (PMLI) formulation, as mentioned in [[44]], provides valuable guarantees regarding the convergence of the waveform itself. Furthermore, the algorithm's implementation is computationally efficient due to the utilization of fast Fourier transform (FFT). However, the nature of the unimodular constraint may lead to the slow convergence of the MPI algorithm, particularly for large lengths. To address this issue, we propose the application of acceleration schemes based on the squared iterative method (SQUAREM). Additionally, for scenarios involving spectral and similarity constraints, we developed a modified algorithm of a similar form.

Figure 1 shows the LPI radar and ES system framework. Conventional LPI radars employ a generalized low-intercept design scheme, including power control, complex modulation, and ultra-low sidelobe, to make it more difficult to process the radar signal at the ES side. This achieves LPI potency to some extent.

The problem of interest is to design a complex unimodular sequence ${\left\{s_n\right\}}_{n=1}^N$ , where N is the length of the sequence, minimizing the probability of detection $P_d$ , i.e.,

(1) $$\begin{array}{c}min\mathit{s}{P}_d\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

Next, we will restate Equation (1) from the perspective of a summation detector based on an FFT filter bank detection algorithm.

Assuming that there are M channels in the sampled bandwidth, which are uniformly distributed in frequency, an integer number of K FFT bins are assigned per channel for power estimation. Therefore, an FFT of length $M\times K$ is required to compute the power for the M channels. Let $\mathbf{W}$ be the diagonal matrix consisting of a linear-phase FIR filter of length $MK$ , which can be expressed as $\mathbf{W}=\mathrm{Diag}{[w_1,...,w_{MK}]}^{\mathrm{T}}$ .

Under the assumption of a complex Gaussian white noise background and a distribution where $\mathit{n}$ follows a mean of 0 and a variance of ${\sigma }_n^2$ , the received signal consists of a combination of $\mathit{n}$ and $\mathit{s}$ . After down-conversion, the received vector $\mathit{r}=[r_1,r_2,...,r_{(P-1)(1-\gamma )\mathrm{MK}+\mathrm{MK}}]$ is divided into $\mathrm{P}$ overlapping sample vectors ${\mathit{r}}_1,...,{\mathit{r}}_P$ as follows:

(2) $${\mathit{r}}_m={r_{m(1-\gamma )MK+1}}^{,}$$

where each vector ${\mathit{r}}_m$ has $\gamma MK$ samples in common with the preceding vector ${\mathit{r}}_{m-1}$ and $0\le \gamma <1$ is overlapping ratio. The step length can be expressed as $L=(1-\gamma )MK$ . The vector ${\mathit{r}}_m$ is windowed by the FIR filter matrix $\mathbf{W}$ , resulting in the windowed sample vector ${\mathit{x}}_m$ :

(3) $${\mathit{x}}_m={w_1}^{r_{mL+1}}$$

For the simplicity of expression, define the selection matrix ${\mathbf{J}}_p$ , which can be expressed as:

(4) $${\mathbf{J}}_p=\begin{array}{ccccc}{\mathbf{0}}^{MK\times (p-1)L}& |& {\mathbf{I}}_{MK}& |& {\mathbf{0}}^{MK\times (P-p)L}\end{array}$$

where ${\mathbf{I}}_{MK}$ is a unit diagonal matrix of dimension $MK$ . The detection variable ${\kappa }_m$ for the mth channel of the FFT filter bank-based summation detector is expressed as [[45]]:

(5) $$\begin{array}{cc}\hfill {\kappa }_m=& \frac{1}{P}\sum _{p=1}^P{\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{r}\right)\hfill \\ \hfill =& \frac{p_s}{P}\sum _{p=1}^P{\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{s}\right)+c_m\hfill \end{array}$$

where ${\mathbf{C}}_m={[{\mathit{\alpha }}_{(m-1)K},...,{\mathit{\alpha }}_{mK-1}]}^{\mathrm{T}}$ and ${\mathit{\alpha }}_k=1/\left(MK\right)·{[1,{\mathrm{e}}^{j\frac{2\pi k}{MK}},...,{\mathrm{e}}^{j\frac{2\pi k(MK-1)}{MK}}]}^{\mathrm{T}}$ . In addition, $p_s$ is the received power of $\mathit{s}$ and the scalar $c_m$ represents the detection scalar associated with noise. When the influence of filters is disregarded, $\mathbb{E}\left\{c_m\right\}=\frac{1}{M}{\sigma }_n^2$ , where $\mathbb{E}\{·\}$ is the operator for mean. A signal is declared to exist in the mth channel only if ${\kappa }_m\ge {\mu }_{th}$ , where ${\mu }_{th}$ is the detection threshold. It should be noted that adjusting the value of ${\mu }_{th}$ can achieve the desired false alarm probability $P_{fa}$ .

In the case where the detection threshold ${\mu }_{th}$ and ${\sigma }_n^2$ is fixed, the waveform design problem for LPI in the mth channel can be expressed as $min\mathit{s}{\kappa }_m$ . The detection vector $\kappa$ is defined as a column vector comprising the elements ${\kappa }_1,{\kappa }_2,...,{\kappa }_M$ . Therefore, the output signal-to-noise ratio (SNR) of the mth channel can be expressed as:

(6) $$\begin{array}{c}\hfill SN{R}_m=\frac{p_s{\sum }_{p=1}^P{\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{s}\right)}{P{c}_m}\end{array}$$

The maximum output SNR of all channels is equal to the total output SNR of the FFT filter bank-based summation detector. Assuming uniform noise levels across all channels, the waveform design problem for LPI scenarios can be succinctly formulated as $min\mathit{s}maxm{\kappa }_m,m=1,...,M$ . Consequently, the LPI problem (1) becomes

(7) $$\begin{array}{c}min\mathit{s}{\parallel \kappa \parallel }_{\infty }\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

where ${\parallel ·\parallel }_{\infty }$ is the infinity norm operator.

Assuming that the optimization variable $\mathit{s}$ is an N-dimensional vector, where N equals $MK+(P-1)L$ , without loss of generality. In this case, we can disregard the energy distortion caused by FIR filters and noise, and set $\mathbf{W}$ equal to the identity matrix $\mathbf{I}$ . Then, the sufficiency problem of Equation (7) can be formulated as follows (as proven in Appendix A):

(8) $$\begin{array}{c}min\mathit{s}maxm,p{\left({\mathbf{C}}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

According to Parseval's theorem, ${\sum }_{m=1}^M{\left({\mathbf{C}}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)=MK$ and it can be further obtained that $0\le maxm,p{\left({\mathbf{C}}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)\le MK$ . The optimal value of Equation (8) must be greater than K; hence, the equivalence problem for Equation (8) can be formulated as follows:

(9) $$\begin{array}{c}min\mathit{s}{\sum }_{m=1}^M{\sum }_{p=1}^P{\parallel {\left({\mathbf{C}}_m^{\mathrm{H}}{\mathit{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)-K\parallel }^2\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

It is important to note that the number of channels is not a priori information, and therefore, Equation (9) can be formulated more appropriately as follows:

(10) $$\begin{array}{c}min\mathit{s}{\sum }_{m=1}^{MK}{\sum }_{p=1}^P{\parallel {\left({\alpha }_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\alpha }_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)-1\parallel }^2\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

Expanding the square in the objective function yields:

(11) $$\begin{array}{c}min\mathit{s}{\sum }_{m=1}^{MK}{\sum }_{p=1}^P({({({\mathit{\alpha }}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s})}^{\mathrm{H}}({\mathit{\alpha }}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}))}^2\hfill \\ -2{({\mathit{\alpha }}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s})}^{\mathrm{H}}({\mathit{\alpha }}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s})+1)\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

Using Parseval's theorem, it can be shown that the second term in the objective function is a constant, i.e., ${\sum }_{m=1}^{MK}{\sum }_{p=1}^P{\left({\mathit{\alpha }}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathit{\alpha }}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)\equiv PMK$ . Therefore, by ignoring the constant terms, the problem can be simplified as follows:

(12) $$\begin{array}{c}min\mathit{s}{\sum }_{m=1}^{MK}{\sum }_{p=1}^P{\left({\left({\mathit{\alpha }}_m^{\mathrm{H}}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathit{\alpha }}_m^{\mathrm{H}}{\mathit{J}}_p\mathit{s}\right)\right)}^2\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

Let us denote the sliding window matrix as $\mathbf{S}$ .

(13) $$\mathbf{S}=\begin{array}{ccccc}{\mathbf{J}}_1\mathit{s}& |& \cdots & |& {\mathbf{J}}_P\mathit{s}\end{array}$$

Equation (12) can be rewritten as

(14) $$\begin{array}{c}min\mathit{s}{\sum }_{m=1}^{MK}{\left({\left({\mathit{\alpha }}_m^{\mathrm{H}}\mathbf{S}\right)}^{\mathrm{H}}\left({\mathit{\alpha }}_m^{\mathrm{H}}\mathbf{S}\right)\right)}^2\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

For the sake of better expression, let us define $\mathbf{X}={\mathbf{SS}}^{\mathrm{H}}$ and ${\mathbf{A}}_m={\alpha }_m{\alpha }_m^{\mathrm{H}}$ . With these definitions, Equation (12) can be rewritten as:

(15) $$\begin{array}{c}min\mathit{s}{\sum }_{m=1}^{MK}\mathrm{Tr}{\left(\mathbf{X}{\mathbf{A}}_m\right)}^2\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

It is worth noting that Equation (15) changes the quartic function on $\mathbf{S}$ into a quadratic function on $\mathbf{X}$ .

The optimization problem for LPI radar waveform design against a summation detector based on FFT filter bank, as denoted by (15), lacks the capability to construct a fuzzy function for the waveform. This limitation prevents the guarantee of certain desirable properties, such as a lower integrated sidelobe level (ISL).

Similarity constraints refer to the requirement that the designed waveform has a certain degree of similarity or correlation among its different components. Several studies have shown that enforcing a similarity constraint on a waveform can control the shape of its ambiguity function [[46]], partially circumventing issues such as significant modulus variation, poor range resolution, and/or high peak sidelobe levels. In Equation (15), we consider the following similarity constraint:

(16) $$\parallel \widehat{\mathit{s}}-{\mathit{s}}_0{\parallel }_{\infty }\le \xi$$

where ${\mathit{s}}_0$ denotes the reference waveforms, $\xi$ is a use-specific parameter which rules the size of similarity region. For a given $\xi >0$ , (16) can be transformed into the following optimization problem:

(17) $$\begin{array}{c}min\mathit{s}{\sum }_{m=1}^{MK}\mathrm{Tr}{\left(\mathbf{X}{\mathbf{A}}_m\right)}^2\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \\ {\parallel \widehat{\mathit{s}}-{\mathit{s}}_0\parallel }_{\infty }\le \xi \hfill \end{array}$$

Furthermore, in practical applications, radar systems must satisfy spectral constraints in addition to having good autocorrelation properties. For example, certain spectral powers must be below specified levels.

Spectral constraints, such as limiting the power in a certain frequency band, denoted by the set of indices $\Omega \subset [1,MK]$ , can be expressed as follows: the power in $\Omega$ should be lower than a certain threshold $ϵ$ .

(18) $$\sum k\in \Omega {\parallel {\mathit{\alpha }}_k^{\mathrm{H}}\mathbf{S}\parallel }_2^2\le ϵ$$

From (16) and (18), the optimization problem (15) can be rewritten as:

(19) $$\begin{array}{c}min\mathit{s}{\sum }_{m=1}^{MK}\mathrm{Tr}{\left(\mathbf{X}{\mathbf{A}}_m\right)}^2\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \\ {\parallel \widehat{\mathit{s}}-{\mathit{s}}_0\parallel }_{\infty }\le \xi \hfill \\ {\sum }_{k\in \Omega }{\parallel {\mathit{\alpha }}_k^{\mathrm{H}}\mathbf{S}\parallel }_2^2\le ϵ\hfill \end{array}$$

In this section, we develop an algorithm for LPI waveform design based on the general quadratic optimization framework.

As $\mathrm{Tr}\left(\mathbf{X}{\mathbf{A}}_m\right)=\mathrm{vec}{\left(\mathbf{X}\right)}^{\mathrm{H}}·\mathrm{vec}\left({\mathbf{A}}_m\right)$ , Equation (15) can be considered as a quadratic optimization problem of generic form as follows:

(20) $$\begin{array}{c}min\mathit{s}\mathrm{vec}{\left(\mathbf{X}\right)}^{\mathrm{H}}\mathbf{\Phi }\mathrm{vec}\left(\mathbf{X}\right)\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

where $\mathbf{\Phi }$ is a positive definite and $\mathbf{\Phi }={\sum }_{m=1}^{MK}\mathrm{vec}\left({\mathbf{A}}_m\right)·\mathrm{vec}{\left({\mathbf{A}}_m\right)}^{\mathrm{H}}$ . Because of $\mathbf{\Phi }⪯\mathbf{I}$ , Equation (20) is NP-hard for the constraint set [[48]]. A monotonically decreasing objective of (20) can be achieved by iteratively updating $\mathrm{vec}\left(\mathbf{X}\right)$ by solving the following nearest problem at each iteration, where $\mathbf{G}=\mathbf{\Phi }-\mathbf{I}$ :

(21) $$\begin{array}{c}max\mathit{s}{\parallel \mathrm{vec}\left(\mathbf{X}\right)-\mathbf{G}\mathrm{vec}{\left(\mathbf{X}\right)}^{\left(i\right)}\parallel }_2\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

where ${\mathbf{X}}^{\left(i\right)}={\mathbf{S}}^{\left(i\right)}{\left({\mathbf{S}}^{\left(i\right)}\right)}^{\mathrm{H}}$ at iteration k. However, according to (21), $\mathrm{vec}\left(\mathbf{X}\right)$ is a quadratic form of $\mathrm{vec}\left(\mathbf{S}\right)$ , unimodular constraints are equivalent to making the $(p-1)MK+p$ th element of $\mathrm{vec}\left(\mathbf{X}\right)$ equal to MK. To solve the problem given in Equation (21), an efficient way is to expand it as follows:

(22) $$\begin{array}{c}\begin{array}{cc}\hfill max\mathit{s}& \mathrm{vec}{\left(\mathbf{X}\right)}^{\mathrm{H}}\mathrm{vec}\left(\mathbf{X}\right)-2\mathrm{Re}\left(\mathrm{vec}{\left(\mathbf{X}\right)}^{\mathrm{H}}\mathbf{G}\mathrm{vec}\left({\mathbf{X}}^{\left(i\right)}\right)\right)\hfill \\ \hfill & +\mathrm{vec}{\left({\mathbf{X}}^{\left(i\right)}\right)}^{\mathrm{H}}\left({\mathbf{G}}^{\mathrm{H}}\mathbf{G}\right)\mathrm{vec}\left({\mathbf{X}}^{\left(i\right)}\right)\hfill \\ \hfill s.t.& \left|s_n\right|=1,n=0,...,N-1\hfill \end{array}\hfill \end{array}$$

Additionally, we have $MK{P}^2\le \mathrm{vec}{\left(\mathbf{X}\right)}^{\mathrm{H}}·\mathrm{vec}\left(\mathbf{X}\right)\le {\left(MKP\right)}^2$ (proof in Appendix B), which is just a constant. After ignoring the constant term in Equation (23), the majorized problem of Equation (22) is given by:

(23) $$\begin{array}{c}min\mathit{s}\mathrm{Re}\left(\mathrm{vec}{\left(\mathbf{X}\right)}^{\mathrm{H}}\mathbf{G}\mathrm{vec}\left({\mathbf{X}}^{\left(i\right)}\right)\right)\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

which can be rewritten as

(24) $$\begin{array}{c}min\mathit{s}{\sum }_{m=1}^{MK}\mathrm{Re}(\mathrm{Tr}\left({\mathbf{X}}^{\left(i\right)}{\mathbf{A}}_m\right)·\mathrm{Tr}\left(\mathbf{X}{\mathbf{A}}_m\right))-\mathrm{Tr}\left({\mathbf{X}}^{\left(i\right)}\mathbf{X}\right)\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

We define ${\mathbf{T}}^{\left(i\right)}={\mathbf{A}}^{\mathrm{H}}{\mathbf{S}}^{\left(i\right)}$ , where $\mathbf{A}$ is the discrete Fourier transform (DFT) matrix given by $\mathbf{A}=\frac{1}{\sqrt{MK}}{[{\mathit{\alpha }}_1,...,{\mathit{\alpha }}_{MK}]}^{\mathrm{T}}$ , and ${\mathit{q}}^{\left(i\right)}\in {\mathbb{C}}^{MK}$ . The m-th element of ${\mathit{q}}^{\left(i\right)}$ can be expressed as $q_m^{\left(i\right)}={\mathit{t}}_m^{\left(i\right)}·{\left({\mathit{t}}_m^{\left(i\right)}\right)}^{\mathrm{H}}$ , where $\left({\mathit{t}}_m^{\left(i\right)}\right)$ is the m-th row vector of ${\mathbf{T}}^{\left(i\right)}$ . Let $\mathbf{Q}=\mathrm{Diag}\left\{{\mathit{q}}^{\left(i\right)}\right\}$ , then Equation (24) can be written as:

(25) $$\begin{array}{c}min\mathit{s}\mathrm{Re}\left(\mathrm{Tr}\left({\mathbf{S}}^{\mathrm{H}}(\mathbf{A}{\mathbf{Q}}^{\left(i\right)}{\mathbf{A}}^{\mathrm{H}}-{\mathbf{S}}^{\left(i\right)}{\left({\mathbf{S}}^{\left(i\right)}\right)}^{\mathrm{H}}){\mathbf{S}}^{\left(i\right)}\right)\right)\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

It is worth noting that $\mathbf{Q}$ is essentially a spectral power sum matrix that encompasses all bands and Equation (25) is a quadratic function on $\mathbf{S}$ . Due to ${\mathbf{A}}^{\mathrm{H}}\mathbf{A}=\mathbf{I}$ , Equation (25) can be formulated as:

(26) $$\begin{array}{c}min\mathit{s}\mathrm{Re}\left(\mathrm{Tr}\left({\mathbf{S}}^{\mathrm{H}}\left({\mathbf{A}}^{\mathrm{H}}{\mathbf{\Psi }}^{\left(i\right)}\mathbf{A}\right){\mathbf{S}}^{\left(i\right)}\right)\right)\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

where ${\mathbf{\Psi }}^{\left(i\right)}={\mathbf{Q}}^{\left(i\right)}-{\mathbf{T}}^{\left(i\right)}{\left({\mathbf{T}}^{\left(i\right)}\right)}^{\mathrm{H}}$ , let ${\mathbf{M}}^{\left(i\right)}={\mathbf{\Psi }}^{\left(i\right)}-q_{\max }^{\left(i\right)}\mathbf{I}$ , where $q_{\max }^{\left(i\right)}=maxp\{q_p^{\left(i\right)}:p=1,...,MK\}$ . Equation (27) can be written as:

(27) $$\begin{array}{c}min\mathit{s}\mathrm{Re}\left(\mathrm{Tr}\left({\mathbf{S}}^{\mathrm{H}}\left(\mathbf{A}{\mathbf{M}}^{\left(i\right)}{\mathbf{A}}^{\mathrm{H}}\right){\mathbf{S}}^{\left(i\right)}\right)\right)\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

By defining ${\mathbf{Z}}^{\left(i\right)}=-\mathbf{A}{\mathbf{M}}^{\left(i\right)}{\mathbf{T}}^{\left(i\right)}$ , when we only consider the unimodular constraint, the majorized problem of Equation (27) can be simplified as follows:

(28) $$\begin{array}{c}min\mathit{s}{\parallel \mathbf{S}-{\mathbf{Z}}^{\left(i\right)}\parallel }_{\mathrm{Fro}}\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \end{array}$$

${\parallel ·\parallel }_{\mathrm{Fro}}$ in (28) represents the Frobenius norm operator. It is easy to see that the problem given in Equation (27) has a closed-form solution, which is given by:

(29) $$\mathbf{S}=e^{j\arg \left({\mathbf{Z}}^{\left(i\right)}\right)}$$

It is important to note that the solution of (29) considers only unimodular constraints, specifically $|s_n|=1,n=0,...,N-1$ , and the solution of the similarity constraint and the spectral constraints in (19) will be presented in Section 4.2. The overall algorithm is summarized in Algorithm 1. While local optimization algorithms typically ensure a monotonic behavior of the optimization objective and eventual convergence through the optimization process, the monotonic minimizer for probability of intercept (MPI) formulation itself provides useful insights, as demonstrated in Appendix C.

The sliding window matrix $\mathbf{S}$ , defined by Equation (13), can be uniquely obtained from $\mathit{s}$ . However, the converse is not true. According to the mapping relation, $\mathit{s}$ can be mapped by several columns of $\mathbf{S}$ , with the number of mappings depending on the overlap time of the receiver window. To better extract $\mathit{s}$ from $\mathbf{S}$ , waveform extraction approaches are introduced based on the number-first and length-first principles, respectively.

In this subsection, we present a sequence extraction approach based on the time-domain circular shift theorem, which yields a sequence characterized by a constant envelope in the time–frequency domain.

Let ${\mathbf{S}}^{*}$ be an optimized sliding window matrix, where ${\mathit{s}}_m$ represents the mth column of ${\mathbf{S}}^{*}$ . Assuming that $\parallel \mathbf{A}{\mathit{s}}_m{\parallel }_2^2\le \varrho$ , where ϱ is a positive number, there exists a unique transformation ${\widehat{s}}_m\left(n\right)=s_m\left(mod(n,MK)\right)$ . This transformation guarantees that the resulting sliding window matrix ${\widehat{\mathbf{S}}}_m$ satisfies the inequality $\parallel \mathbf{A}{\widehat{\mathbf{S}}}_m{\parallel }_{\mathrm{Fro}}^2\le P\varrho$ .

From Equation (13), the sliding window matrix of ${\widehat{\mathit{s}}}_m$ is

(30) $${\widehat{\mathbf{S}}}_m=\begin{array}{ccccc}{\mathbf{J}}_1{\widehat{\mathit{s}}}_m& |& \cdots & |& {\mathbf{J}}_P{\widehat{\mathit{s}}}_m\end{array}$$

Let ${\widehat{\mathbf{U}}}_m=\mathbf{A}{\widehat{\mathbf{S}}}_m$ and ${\mathit{u}}_p^m,p=1,...,P$ is the pth colomn of ${\widehat{\mathbf{U}}}_m$ . According to the time-domain cyclic shift theorem [[50]]

(31) $${\parallel \mathbf{A}{\mathit{s}}_m\parallel }_2^2={\parallel {\mathit{u}}_p^m\parallel }_2^2\le \varrho ,p=1,...,P$$

$\parallel {\widehat{\mathbf{U}}}_m{\parallel }_{\mathrm{Fro}}^2$ can be expanded as

(32) $$\begin{array}{cc}\hfill {\parallel {\widehat{\mathbf{U}}}_m\parallel }_{\mathrm{Fro}}^2=& \sum _{p=1}^P\sum _{i=1}^{MK}{\left|u_{i,p}^m\right|}^2\hfill \\ \hfill =& \sum _{p=1}^P{\parallel {\mathit{u}}_p^m\parallel }_2^2\hfill \\ \hfill =& P·{\parallel {\mathit{u}}_p^m\parallel }_2^2\le P·\varrho ,p=1,...,P\hfill \end{array}$$

where $u_{i,p}^m$ is the $(i,p)$ th element of ${\widehat{\mathbf{U}}}_m$ . Proof is complete. □

From Theorem 1, we can obtain the expression for $\mathit{s}$ from the sliding window matrix $\mathbf{S}$ as follows:

(33) $${\widehat{\mathit{s}}}_i\left(n\right)={\mathit{s}}_i\left(mod(n,MK)\right),n=1,...,N,i=1,...,P$$

where ${\mathit{s}}_i$ denotes the i-th column vector of $\mathbf{S}$ . It should be noted that any column in $\mathbf{S}$ is convergent. Therefore, we can obtain at most P distinct $\widehat{\mathit{s}}$ from $\mathbf{S}$ . However, the cyclic nature of $\widehat{\mathit{s}}$ leads to a deterioration in the correlation property. When N is relatively large, this approach is not suitable for constructing the sequence.

In this subsection, we introduce a length-first extraction method. As mentioned in the previous section, the number-first extraction approach, while better conforming to (14), is not suitable for constructing long sequences when N is relatively large due to the cyclic nature of the time-domain. In such cases, we can extract $\widehat{\mathit{s}}$ as follows:

(34) $$\widehat{\mathit{s}}=\mathrm{vec}\left(\mathbf{S}\right)$$

Clearly, the $\widehat{\mathit{s}}$ extracted in this way has a maximum length of $PMK$ . It is worth mentioning that the sequence obtained by this extraction will be somewhat distorted for $\gamma >0$ . However, due to Theorem 1, the sequence is somewhat robust. For ease of discussion, in this paper, we only consider the length-first extracted method.

In this subsection, we provide optimized solutions for similarity constraints and spectral constraints for MPI.

The similarity constrain in Equation (16) is equivalent to $\mathrm{Re}\left({\widehat{s}}_n^{*}{s}_0\left(n\right)\right)\ge 1-{\xi }^2/2$ for $n=1,2,...,N$ , implying that $\arg \left({\widehat{s}}_n\right)\in [{\gamma }_n,{\gamma }_n+\delta ]$ , where ${\gamma }_n=\arg \left(s_0\left(n\right)\right)-\arccos (1-{\xi }^2/2)$ , $\delta =2\arccos (1-{\xi }^2/2)$ . Letting ${\phi }_n$ be the phase of ${\widehat{s}}_n$ , the optimal solution ${\phi }_n^{*}$ can be given as [[51]]:

(35) $${\phi }_n^{*}=\begin{array}{cc}{\gamma }_n+\delta ,\hfill & ifcos({\phi }_n+{\gamma }_n+\delta )\ge cos({\phi }_n+{\gamma }_n)\hfill \\ {\gamma }_n,\hfill & ifcos({\phi }_n+{\gamma }_n+\delta )<cos({\phi }_n+{\gamma }_n)\hfill \end{array}$$

As shown in Equation (18), the optimization problem is utilized to enforce similarity constraints and spectral constraints, ensuring that the designed waveforms satisfy these requirements.

For a given $ϵ>0$ in (19), we can always find a $\rho$ such that problem (14) can be transformed into the following equivalent problem:

(36) $$\begin{array}{c}min\mathbf{S}{\sum }_{m=1}^{MK}{\left({\left({\mathit{\alpha }}_m^{\mathrm{H}}\mathbf{S}\right)}^{\mathrm{H}}\left({\mathit{\alpha }}_m^{\mathrm{H}}\mathbf{S}\right)\right)}^2+\rho \sum k\in \Omega {\left({\mathit{\alpha }}_k^{\mathrm{H}}\mathbf{S}\right)}^{\mathrm{H}}\left({\mathit{\alpha }}_k^{\mathrm{H}}\mathbf{S}\right)\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \\ {\parallel \widehat{\mathit{s}}-{\mathit{s}}_0\parallel }_{\infty }\le \xi \hfill \end{array}$$

where $\Omega$ is the low-power frequency band. From here, we can follow the derivation of MPI to obtain the spectral-MPI algorithm for problem (26), and (26) can be rewritten as:

(37) $$\begin{array}{c}min\mathbf{S}\mathrm{Re}(\mathrm{Tr}({\mathbf{T}}^{\mathrm{H}}({\widehat{\mathbf{Q}}}^{(i)}-{\mathbf{T}}^{(i)}{({\mathbf{T}}^{(i)})}^{\mathrm{H}}){({\mathbf{T}}^{(i)})}^{\mathrm{H}}))\hfill \\ s.t.\left|s_n\right|=1,n=0,...,N-1\hfill \\ {\parallel \widehat{\mathit{s}}-{\mathit{s}}_0\parallel }_{\infty }\le \xi \hfill \end{array}$$

where ${\widehat{\mathbf{Q}}}^{\left(i\right)}=\mathrm{Diag}\left\{{\widehat{\mathit{q}}}^{\left(i\right)}\right\}$ , and the m-th element of ${\widehat{\mathit{q}}}^{\left(i\right)}$ can be rewritten as

(38) $${\widehat{q}}_m^{\left(i\right)}=\begin{array}{cc}{\mathit{t}}_p^{\left(i\right)}{\left({\mathit{t}}_p^{\left(i\right)}\right)}^{\mathrm{H}}+\rho /2,\hfill & m\in \Omega \hfill \\ {\mathit{t}}_p^{\left(i\right)}{\left({\mathit{t}}_p^{\left(i\right)}\right)}^{\mathrm{H}},\hfill & otherwise.\hfill \end{array}$$

By defining ${\widehat{q}}_{\max }^{\left(i\right)}=maxm\{{\widehat{q}}_m^{\left(i\right)}:m=1,...,n_0\}$ , ${\widehat{\mathbf{\Psi }}}^{\left(i\right)}={\widehat{\mathbf{Q}}}^{\left(i\right)}-{\mathbf{T}}^{\left(i\right)}{\left({\mathbf{T}}^{\left(i\right)}\right)}^{\mathrm{H}}$ and ${\widehat{\mathbf{M}}}^{\left(i\right)}={\widehat{\mathbf{\Psi }}}^{\left(i\right)}-{\widehat{q}}_{\max }^{\left(i\right)}\mathbf{I}$ , the problem can be recast as shown in (28):

(39) $$\begin{array}{c}min\mathbf{S}\mathrm{Re}(\mathrm{Tr}({\mathbf{S}}^{\mathrm{H}}(\mathbf{A}{\widehat{\mathbf{M}}}^{(i)}{\mathbf{A}}^{\mathrm{H}}){\mathbf{S}}^{(i)}))\hfill \\ s.t.|s_n|=1,n=0,...,N-1\hfill \\ {\parallel \widehat{\mathit{s}}-{\mathit{s}}_0\parallel }_{\infty }\le \xi \hfill \end{array}$$

This problem has the same form as problem (28), for which the similarity constraint solution is shown in Equation (35). Therefore, we can follow the same steps as before to derive the spectral-MPI algorithm, the main steps of which are listed in Algorithm 2.

The computational complexity of each iteration in MPI is primarily dominated by two matrix multiplications involving $\mathbf{A}$ , as shown in Algorithm 1. These computations can be easily performed using FFT and IFFT operations. Specifically, in the i-th iteration, the i-th column of ${\mathbf{T}}^{\left(i\right)}$ corresponds to the $MK$ -FFT operation of the i-th column of ${\mathbf{S}}^{\left(i\right)}$ , with a computational complexity of $\mathcal{O}(MKPlog(MK\left)\right)$ . The matrix product ${\mathbf{T}}^{\left(i\right)}{\left({\mathbf{T}}^{\left(i\right)}\right)}^{\mathrm{H}}$ needs to be computed only once, and $\mathbf{Q}$ can be obtained from a diagonal vector, with a computational complexity of $\mathcal{O}(MK+P{\left(MK\right)}^2)$ . Similarly, the matrix ${\mathbf{Z}}^{\left(i\right)}$ can be obtained by P $MK$ -IFFT operation, with a computational complexity of $\mathcal{O}(MKPlog(MK\left)\right)$ . Additionally, the computational complexity of performing complex angle operations on ${\widehat{\mathbf{Z}}}^{\left(i\right)}$ is denoted by $\mathcal{O}\left(MKP\right)$ . In summary, the computational complexity of each iteration can be expressed as $\mathcal{O}(2MKPlog(MK)+MK(MKP+P+1\left)\right)$ . The MPI algorithm exhibits high computational efficiency by using P FFT and IFFT operations, making it suitable for designing long sequences.

We also propose an improved version of the MPI algorithm to speed up the convergence process.

The MPI algorithm is derived based on the majorization–minimization principle, and the nature of the majorization functions determines the convergence rate of the algorithm. However, numerical simulations have shown that MPI converges very slowly when N is large. In this subsection, we address this issue and propose an acceleration scheme to improve the convergence speed of the algorithm. The acceleration scheme is called the squared iterative method (SQUAREM), which was originally proposed in [[52]] to accelerate any EM algorithm. Here, we outline the main steps of SQUAREM and modify it accordingly to address the waveform design problem we encounter.

Let $\mathbb{F}$ denote the fixed-point mapping during the iterations of the proposed MPI algorithm (noting that the mapping varies with the enforced constrains), i.e.,

(40) $${\mathbf{S}}^{(i+1)}=\mathbb{F}\left({\mathbf{S}}^{\left(i\right)}\right)$$

Motivated by the Cauchy–Barzilai–Borwein (CBB) method for accelerating the convergence of the classical Cauchy method for solving linear equations, the SQUAREM method iterates as follows:

(41) $${\mathbf{S}}^{(i+1)}=e^{jarg({\mathbf{S}}^{\left(i\right)}-2{\alpha }^{\left(i\right)}{\mathbf{R}}^{\left(i\right)}+{\left({\alpha }^{\left(i\right)}\right)}^2{\mathbf{V}}^{\left(i\right)})}$$

where ${\mathbf{R}}^{\left(i\right)}=\mathbb{F}\left({\mathbf{S}}^{\left(i\right)}\right)-{\mathbf{S}}^{\left(i\right)}$ , ${\mathbf{V}}^{\left(i\right)}=\mathbb{F}\left(\mathbb{F}\left({\mathbf{S}}^{\left(i\right)}\right)\right)-\mathbb{F}\left({\mathbf{S}}^{\left(i\right)}\right)-{\mathbf{R}}^{\left(i\right)}$ , and ${\alpha }^{\left(i\right)}$ is the step-length in i-th iteration. Follow the step update approach as follows:

(42) $${\alpha }^{\left(i\right)}=-\frac{{\parallel \mathbf{R}\parallel }_{\mathrm{Fro}}}{{\parallel \mathbf{V}\parallel }_{\mathrm{Fro}}}$$

where ${\parallel ·\parallel }_{\mathrm{Fro}}$ is the Frobenius norm. To ensure the descent property, a backtracking strategy is adopted. This strategy repeatedly halves the distance between ${\alpha }^{\left(i\right)}$ and $-1$ (i.e., ${\alpha }^{\left(i\right)}⟵({\alpha }^{\left(i\right)}-1)/2$ ) until the descent property is maintained. Monotonicity can be determined by checking whether ${\parallel {\mathbf{A}}^{\mathrm{H}}{\mathbf{S}}^{(i+1)}\parallel }_{\mathrm{Fro}}\le {\parallel {\mathbf{A}}^{\mathrm{H}}{\mathbf{S}}^{\left(i\right)}\parallel }_{\mathrm{Fro}}$ holds or not. The accelerated MPI based on SQUAREM is summarized in Algorithm 3.

It is worth noting that in each iteration, the acceleration scheme requires the computation of ${\mathbf{S}}^{(i+1)}$ and ${\mathbf{S}}^{(i+2)}$ , as well as the Frobenius norm and angle operations associated with $MKP$ , which increases the computational complexity of the acceleration scheme. However, the acceleration scheme can significantly reduce the number of iterations, thereby achieving a faster convergence time.

In this section, we provide several numerical examples to demonstrate the performance of the proposed algorithm. All experiments were performed on a PC with a 2.9 GHz AMD-Ryzen7-4800HS CPU and 24 GB RAM.

We first compared the convergence speed of the MPI algorithm with the accelerated MPI algorithm, measured by the square root of (14), denoted by ${\Lambda }^{\left(i\right)}$ . For all algorithms in this experiment, the stopping criterion was set to $|{\Lambda }^{(i+1)}-{\Lambda }^{\left(i\right)}\le {10}^{-8}|$ , and the initial sequence ${\left\{x_n^{\left(0\right)}\right\}}_{n=1}^N$ was chosen to be ${\left\{e^{j2\pi {\theta }_n}\right\}}_{n=1}^N$ , where ${\left\{{\theta }_n\right\}}_{n=1}^N$ are independent random variables uniformly distributed in $[0,1]$ . In the spectrum-MPI of this subsection, the low power frequency band $\Omega$ is set to $[0.2,0.8]$ and $\rho =1$ . To reduce the amount of computation, the sequence extraction approach in this paper is based on the principle of length-first.

Figure 2 shows the convergence properties of the MPI algorithm and the spectral MPI algorithm with their respective acceleration schemes for different parameters. We observe that the proposed MPI and spectral-MPI algorithms always converge monotonically, and the accelerated algorithm is effective in increasing the convergence rate. Additionally, Figure 3 compares the convergence time of the MPI, accelerated MPI, and CAN [[53]] algorithms with respect to the sequence length N. The results show that the accelerated MPI algorithm has a faster convergence time than the MPI algorithm. Moreover, since the size of the FFT operation is related to the size of $MK$ , the convergence time of the proposed algorithm does not significantly increase with the sequence length, unlike the case of CAN.

Subsequently, an experiment is conducted to illustrate the superior interference performance of the proposed sequence in comparison to commonly utilized sequences like SIN, LFM, CAN, and Golomb [[54]]. Among these, SIN and LFM can be respectively defined as a sine signal represented by $exp\left(j0.2\pi n\right)$ and a linear frequency modulated signal represented by $exp(j0.2\pi n^2/N)$ . Assuming a free-space propagation model and neglecting gain and loss, the relationship between the maximum propagation distance $r_{max}$ and the minimum received power $p_{smin}$ can be expressed as [[55]]:

(43) $$r_{max}={\frac{P_t{G}_t{G}_r{\lambda }^2}{4^{\pi }{p}_{smin}}}^{1/2}$$

where $P_t$ , $\lambda$ , $G_t$ , and $G_r$ represent the transmit power, wavelength, transmit gain, and receive gain, respectively. In this experiment, we only consider the cases where $G_t=1$ and $G_r=1$ . On one hand, $P_t={10}^2$ watts and $\lambda =0.1$ meters. The sequence is parameterized by $N=512,M=16,K=16$ . The minimum output power $p_{omin}$ after processing with the FFT filter bank-based summation detector can be represented as $p_{omin}=p_{smin}·\frac{maxm{\sum }_{p=1}^P{\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{s}\right)}{PK}$ .

Figure 4 compares the intercept distances between common sequences and the proposed sequence for $p_{smin}\in [-90,-70]$ dBm and ES system setting: $M=8,K=16$ . In this configuration, if $\frac{maxm{\sum }_{p=1}^P{\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{s}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{W}{\mathbf{J}}_p\mathit{s}\right)}{PK}=1$ , the theoretical lower bound is computed. Obviously, the proposed sequences have a smaller maximum propagation distance than the normal sequence under the same $p_{smin}$ , especially for SIN with LFM sequences. Moreover, the proposed sequence almost reaches the theoretical lower bound.

Figure 5 compares the maximum intercept distance of the proposed waveform for a minimum received power $p_{smin}$ of −80 dBm and window lengths $MK$ ranging from 80 to 256. The results demonstrate that the proposed algorithm does not require prior knowledge of the receiver window length $MK$ and maintains a low interception property across different window lengths.

To investigate the degree to which the similarity constraint improves the autocorrelation properties of the sequence, we conducted an experiment to evaluate the merit factor (MF) under a different $\xi$ . The MF is defined as the ratio of the central lobe energy to the total energy of all other lobes, as proposed by Golay in 1972 [[56]].

(44) $$\mathrm{MF}=\frac{|r_0{|}^2}{2{\sum }_{k=1}^{k=N-1}{\left|r_k\right|}^2}$$

where ${\mathit{r}}_k$ denotes aperiodic autocorrelations which can be defined as:

(45) $$r_k=\sum _{n=1}^{N-k}{x}_n{x}_{n+k}^{*}=r_{-k}^{*},k=0,...,N-1$$

The accelerated-MPI algorithm initialized by Golomb or random (RN) sequence are applied to the following lengths: $N=2^7,2^8,2^9,2^{10},2^{11},2^{12}$ , with $\xi =0.01$ . We chose the CAN sequence with superior MF properties as the reference sequence. The MF of the resulting sequences is shown in Figure 6. We observe that the MF of the proposed algorithm is almost independent of the sequence length. Figure 7 compares the MF of the accelerated-MPI initialized by Golomb and RN, respectively, versus $\xi$ from 0 to 1, and the proposed sequence is parameterized by $N=512,M=8,K=16$ . The results demonstrate that the similarity constraint can improve certain properties as $\xi$ increases, which increases the flexibility of the design.

In addition, we considered the design of spectral constraints with low spectral power in frequency bands $[0,0.2)\bigcup (0.3,0.5)\bigcup (0.8,1]$ . Figure 8 shows the time–frequency energy flow of the output sequence generated by spectral-MPI when initialized with a Golomb sequence and $\rho =1$ . To observe the change in time–frequency energy more intuitively due to spectral constraints, Figure 9 shows the time–frequency energy flow without spectral constraints. We observe that the power in the pre-specified frequency bands has been suppressed, and the time–frequency flow is similar to the constant-envelope smooth.

We have introduced the MPI algorithm, an efficient method for minimizing the probability of intercept. This algorithm is derived from the general quadratic optimization framework. In scenarios involving spectral and similarity constraints, we also developed the spectral-MPI algorithm to achieve enhanced properties within the same framework as MPI. Additionally, we have explored acceleration schemes to expedite the MPI algorithm. Numerical results demonstrate that, when applied to FFT filter bank-based summation detectors, the proposed MPI algorithm can generate sequences with a lower probability of intercept compared to common sequences. Moreover, it can be utilized to design sequences with suppressed spectral power in arbitrary frequency bands. Importantly, all of the proposed algorithms can be efficiently implemented using FFT and have demonstrated computational efficiency in practical applications.

Graph: Figure 1 The sketch of the LPI radar and ES system framework.

Graph: Figure 2 Comparison of the convergence between MPI and spectral-MPI with different parameters and their respective accelerated versions.

Graph: Figure 3 Comparison of the convergence time of the spectral-MPI and accelerate-spectral-MPI.

Graph: Figure 4 Comparison of the maximum propagation distance of common sequence and accelerated-MPI initialized by different sequences.

Graph: Figure 5 Comparison of the maximum propagation distance between normal waveforms and accelerated-MPI with different MK.

Graph: Figure 6 Comparison of the merit factor of common sequences and accelerated-MPI initialized by different sequences.

Graph: Figure 7 Comparison of the merit factor of accelerated-MPI initialized by different sequences versus ξ.

Graph: Figure 8 Time–frequency energy flow of spectral-MPI with low spectral power in frequency bands [0,0.2)⋃(0.3,0.5)⋃(0.8,1].

Graph: Figure 9 Time–frequency energy flow of spectral-MPI with low spectral power in frequency bands [0,1].

Conceptualization, methodology, software, and writing—original draft preparation, Q.L.; supervision, F.G.; software and validation, K.X.; writing—review and editing, Z.L.; funding acquisition, W.H. All authors have read and agreed to the published version of the manuscript.

The original contributions presented in this study are included in the article, whilst further inquiries can be directed to the corresponding author.

The authors declare no conflicts of interest.

To prove that (8) is a sufficient condition for (7), we need to prove

$$\underset{m,p}{max}{\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right)\ge {\parallel \kappa \parallel }_{\infty }$$ (A1)

Assuming that there exist $m_0$ and $p_0$ such that

$$\begin{array}{cc}\hfill & \underset{m,p}{max}{\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_m^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right)\hfill \\ \hfill & ={\left({\mathbf{C}}_{m_0}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_0}\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_{m_0}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_0}\mathit{r}\right)\hfill \end{array}$$ (A2)

When $P=1$ , it is obvious that

$${\left({\mathbf{C}}_{m_0}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_0}\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_{m_0}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_0}\mathit{r}\right)\equiv {\parallel \kappa \parallel }_{\infty }$$ (A3)

When $P>1$ , assuming there exists $m_1$ such that

$${\parallel \kappa \parallel }_{\infty }=\frac{1}{P}\sum _{p=1}^P{\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right)$$ (A4)

Let $p=p_1$ be the maximal term in the summation on the right-hand side of (A4), i.e.,

$$\begin{array}{cc}\hfill & {\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_1}\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_1}\mathit{r}\right)\hfill \\ \hfill & \ge {\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right),\forall p\in [1,P]\hfill \end{array}$$ (A5)

Therefore, we have:

$$\begin{array}{cc}\hfill & {\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_1}\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_1}\mathit{r}\right)\hfill \\ \hfill & \ge \frac{1}{P}\sum _{p=1}^P{\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_{m_1}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_p\mathit{r}\right)\hfill \\ \hfill & ={\parallel \kappa \parallel }_{\infty }\hfill \end{array}$$

According to the definition,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\left({\mathbf{C}}_{m_0}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_0}\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_{m_0}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_0}\mathit{r}\right)\hfill \\ \hfill & \ge {\left({\mathbf{C}}_{m_0}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_0}\mathit{r}\right)}^{\mathrm{H}}\left({\mathbf{C}}_{m_0}^{\mathrm{H}}\mathbf{H}{\mathit{J}}_{p_0}\mathit{r}\right)\hfill \\ \hfill & \ge {\parallel \kappa \parallel }_{\infty }\hfill \end{array}\end{array}$$

This proves that (8) is a sufficient condition for (7). The proof is complete.     □

of upper and lower bounds for $\mathrm{vec}{\left(\mathbf{X}\right)}^{\mathrm{H}}·\mathrm{vec}\left(\mathbf{X}\right)$ . The matrix $\mathbf{S}$ can be represented as:

$$\mathbf{S}=\left[\begin{array}{cccc}{s}_1& s_{L+1}& \cdots & s_{(P-1)L+1}\\ s_2& s_{L+2}& \cdots & s_{(P-1)L+2}\\ ⋮& ⋮& \ddots & ⋮\\ s_{MK}& s_{L+MK}& \cdots & s_{(P-1)L+MK}\end{array}\right]$$ (A6)

Since $\mathbf{X}=\mathbf{S}{\mathbf{S}}^{\mathrm{H}}$ , the elements of $\mathbf{X}$ , denoted by $x_{ij}$ , can be expressed as:

$$x_{ij}=\sum _{p=1}^P{s}_{i+(p-1)L}^{*}{s}_{j+(p-1)L}$$ (A7)

Given that ${\mathit{s}}^{\mathrm{H}}{\mathit{e}}_n^{\mathrm{H}}{\mathit{e}}_n\mathit{s}=1$ , for $n=0,...,N-1$ , we can express $\mathrm{vec}{\left(\mathbf{X}\right)}^{\mathrm{H}}·\mathrm{vec}\left(\mathbf{X}\right)$ as:

$$\begin{array}{cc}\hfill \mathrm{vec}{\left(\mathbf{X}\right)}^{\mathrm{H}}·\mathrm{vec}\left(\mathbf{X}\right)& =\sum _{i=1}^{MK}\sum _{j=1}^{MK}{\left|\sum _{p=1}^P{s}_{p+(i-1)L}^{*}{s}_{p+(j-1)L}\right|}^2\hfill \\ \hfill & =\sum _{i=1}^{MK}\sum _{j=1j\ne i}^{MK}{\left|{\tilde{r}}_{i,j}\right|}^2+MK{P}^2\hfill \\ \hfill & =\sum _{i=1}^{MK}\sum _{j=1j\ne i}^{MK}{\left|{\mathit{s}}_i^{\mathrm{H}}·{\mathit{s}}_j\right|}^2+MK{P}^2\hfill \end{array}$$ (A8)