Unified Integrated Sensing and Communication Signal Design: A Sphere Packing Perspective

Shuaishuai Guo, Senior Member, IEEE, and Kaiqian Qu

S. Guo and K. Qu are with the School of Control Science and Engineering, and also with Shandong Key Laboratory of Wireless Communication Technologies, Shandong University, Jinan 250062, China (e-mail<EMAIL_ADDRESS>qukaiqian@mail.sdu.edu.cn).

The design of communication signal sets is fundamentally a sphere packing problem. It aims to identify a set of $M$ points in an $N$-dimensional space, with the objective of maximizing the separability of points that represent different bits. In contrast, signals used for sensing targets should ideally be as deterministic as possible. This paper explores the inherent conflict and trade-off between communication and sensing when these functions are combined within the same signal set. We present a unified approach to signal design in the time, frequency, and space domains for integrated sensing and communication (ISAC), framing it as a modified sphere packing problem. Through adept formula manipulation, this problem is transformed into a large-scale quadratic constrained quadratic programming (QCQP) challenge. We propose an augmented Lagrangian and dual ascent (ALDA) algorithm for iterative problem-solving. The computational complexity of this approach is analyzed and found to be daunting for large, high-dimensional signal set designs. To address this, we introduce a bit-dimension-power splitting (BDPS) method. This method decomposes the large-scale QCQP into a series of smaller-scale problems that can be solved more efficiently and in parallel, significantly reducing the overall computational load. Extensive simulations have been conducted to validate the effectiveness of our proposed signal design methods in the context of ISAC.

Integrated sensing and communications, signal set design, sphere packing, large-scale quadratic constrained quadratic programming

§ INTRODUCTION

Integrated sensing and communication (ISAC) is poised to be a key technology in sixth-generation (6G) communication networks, offering ubiquitous sensing capabilities  <cit.>. ISAC enhances spectrum and energy efficiency by sharing the same signal waveform for both communication and sensing. However, achieving the optimal balance between these functions presents a significant challenge. For communication, signals representing different bits must be highly separable, whereas for sensing, signals should be as deterministic as possible. This paper delves into the conflict and trade-off inherent in integrating communication and sensing within the same signal set.

In communication, signal set design is essentially a sphere packing problem, focusing on densely packing spheres without overlap in an $N$-dimensional space. This is mathematically formulated as finding $M$ points (representing $m=\log_2 M$ bits) in an $N$-dimensional space, maximizing the Euclidean distance between the closest pair under a given average power constraint $P$, i.e., <cit.>

\begin{align}
(\textbf{P1}):~\mathrm{Given}: &~M,N,P\notag\\
\mathrm{Find}:&~\mathcal{X}=\{\mathbf{x}_1,\mathbf{x}_2,\cdots,\mathbf{x}_M\}\notag\\
\mathrm{Maximize}:&~\min_{\forall k,l,k\neq l} ||\mathbf{x}_k-\mathbf{x}_l||_2\\
\mathrm{Subject~to}:&~\frac{1}{M}\sum_{k=1}^{M}||\mathbf{x}_k||_2^2\leq P\notag.
\end{align}

For sensing applications, it is crucial for signals to be deterministic, with desirable properties like minimum integrated sidelobe level (ISL), maximum merit factor (MF), and optimal ambiguity function (AF) <cit.>. While there are various signal waveform designs for active sensing systems <cit.>, this paper assumes the availability of an optimal sensing signal[Generally, the optimal waveform for sensing is designed without considering communication performance, such as linear frequency modulation (LFM) signals and optimized spatial beampatterns towards desired directions. ], represented by $\mathbf{x}_0$. The task is to find a signal set that fulfills both communication and sensing functions, formulated as a sphere-packing problem with similarity constraint:

\begin{align}
(\textbf{P2}):~\mathrm{Given}: &~M,N,P,\mathbf{x}_0,\epsilon\notag\\
\mathrm{Find}:&~\mathcal{X}=\{\mathbf{x}_1,\mathbf{x}_2,\cdots,\mathbf{x}_M\}\notag\\
\mathrm{Maximize}:&~\min_{\forall k,l,k\neq l} ||\mathbf{x}_k-\mathbf{x}_l||_2\\
\mathrm{Subject~to}:&~\frac{1}{M}\sum_{k=1}^{M}||\mathbf{x}_k||_2^2\leq P,\notag\\
&~||\mathbf{x}_k-\mathbf{x}_0||_2\leq \epsilon, \forall k\notag,
\end{align}

where $\epsilon$ is the sensing waveform similarity tolerance. This paper aims to develop universal solutions for this fundamental problem, facilitating its application in ISAC systems

§.§ Related Works on ISAC Signal Design

The design of ISAC signals, crucial to ISAC technology, has seen extensive study. Current research divides ISAC waveform/signal design into two categories: designs using non-overlapping resources and those employing fully unified signal designs that occupy overlapping resources <cit.>.

Earlier studies have predominantly focused on ISAC waveform designs that utilize non-overlapping resources, primarily due to their simplicity and ease of implementation <cit.>. This approach typically involves allocating orthogonal resources across multiple dimensions such as time <cit.>, frequency <cit.>, and space <cit.>. A notable example in the time domain is presented in <cit.>, where a duplex radar-communication waveform design is proposed. This design divides the transmission duration into distinct radar and communication cycles. Radar cycles utilize a trapezoidal frequency modulated continuous wave (FMCW), comprising an up-chirp, a constant-frequency period, and a down-chirp. Communication cycles, on the other hand, can employ any modulation method. These cycles alternate to prevent time-domain interference, a concept that has been standardized in IEEE 802.11ad <cit.>. In the frequency domain, studies like <cit.> allocate sensing and communication functionalities to different subcarriers, considering channel conditions, functional requirements, and transmitter power budgets. Meanwhile, spatial-domain approaches <cit.> transmit sensing and communication waveforms over different spatial beams or into each other's null space to avoid interference. While these methods are straightforward and readily implementable due to their decoupling of communication and sensing functions, they exhibit a key drawback: inefficient utilization of time, spectrum, and space, leading to reduced energy efficiency.

The challenge of inefficient resource utilization in ISAC has led to a growing interest in fully unified ISAC signal designs that utilize overlapping resources. These can be broadly categorized into three approaches:

* Communication-centric design (CCD): This approach primarily enhances conventional communication waveforms with sensing capabilities. For instance, <cit.> explore integrated signal designs in multicarrier communication systems by modifying orthogonal frequency division multiplexing (OFDM) signals. Sturm and Wiesbeck in <cit.> proposed a radar processing scheme for OFDM communication waveforms, while Johnston et al. <cit.> focused on jointly designing transmission waveforms and reception filters in multiple-input multiple-output (MIMO-OFDM) systems to meet specific requirements like bit error rate (BER) and directional radiated power.

* Radar-centric design (RCD): RCD prioritizes sensing performance, incorporating communication data into radar waveforms. A common method is using radar waveforms as carriers for modulating communication symbols. For example, <cit.> employed frequency modulated continuous wave (FMCW) and chirp waveforms as carriers for communication information. Another approach, seen in <cit.>, involves selecting radar waveforms from a set to form a communication symbol set. Li <cit.> suggested using a set of radar waveform filters as a communication codebook.

* Joint design (JD): JD aims to find a flexible balance between communications and radar through optimization-based signal designs <cit.>. The work in <cit.> introduced a compromise factor for minimizing radar waveform mismatch and reducing multi-user interference in communications. This concept has influenced subsequent research, such as <cit.>. Liu et al. in <cit.> proposed a constructive interference (CI) symbol-level precoding to satisfy both radar beampattern mismatch and communication quality of service (QoS) requirements. To simplify the implementation, <cit.> introduced block-level ISAC precoding for finite alphabet inputs. Meanwhile, <cit.> considered ISAC signal design with hybrid beamforming and sub-connected structures to reduce hardware complexity and cost. Guo et al. <cit.> proposed constructing an ISAC beamformer candidate set for scenarios with limited RF chains to enhance spectral efficiency. Lastly, Liu et al. in <cit.> explored minimizing the estimation Cramér-Rao bound (CRB) for sensing while meeting communication signal-to-interference-and-noise ratio (SINR) requirements.

It's important to note that all the aforementioned studies introduce certain inherent constraints in comparison to the pure ISAC signal design problem (P2). For instance, research aiming to maximize spectral efficiency under Shannon's theorem inherently assumes a Gaussian input constraint <cit.>. Similarly, studies focused on minimizing BER or managing symbol-level interferences typically operate under the assumption that transmitted or shaped data symbols are of $M$-QAM or $M$-PSK formats <cit.>. These additional constraints can potentially compromise the effectiveness of ISAC signals. As a result, the signals or waveforms derived under these conditions may be sub-optimal compared to those that would be obtained by directly solving the pure ISAC signal design problem (P2)

§.§ Contributions

As introduced earlier, the design of communication signals inherently resembles a sphere packing problem. The integration of sensing functionality adds an additional constraint to this problem, forming (P2), a concept not fully explored in previous studies. This paper seeks to tackle this fundamental problem without imposing any extraneous constraints. The primary contributions of this work are outlined as follows:

* We present a unified approach to formulating the ISAC signal design problem across time, frequency, and space domains. This culminates in a modified sphere packing problem aimed at maximizing the minimum Euclidean distance (MMED) between high-dimensional signals, considering both sensing waveform constraints and power limitations.

* Through meticulous formula manipulation, we transform the generalized ISAC signal design problem into a large-scale quadratic constrained quadratic programming (QCQP) problem. We introduce an augmented Lagrangian and dual ascent (ALDA) algorithm to address this challenge, with a detailed analysis of the computational complexity.

* To mitigate the computational complexity, we propose a bit-dimension-power splitting (BDPS) method. This approach divides the objective function and constraints, converting the overarching QCQP problem into a series of smaller-scale QCQP sub-problems. These can be solved in parallel, significantly reducing computational demands.

* Extensive simulations validate the efficacy of our algorithms. Through thorough comparisons with existing works, we demonstrate the superiority of our ISAC signal designs. Our results indicate that our proposed designs achieve a lower symbol error rate (SER) while maintaining robust sensing performance. Moreover, the BDPS algorithm effectively lowers complexity with minimal performance degradation.

§.§ Organization

The remainder of this paper is organized as follows.
In Section II, we delve into the formulation of ISAC signal design problems across time, frequency, and spatial domains. We then unify these problems to establish a comprehensive framework. Section III is dedicated to the transformation of the unified ISAC signal design problem. In Section IV, we introduce the ALDA solution, a method we've developed to solve the transformed ISAC signal design problem. In Section V, we discuss the BDPS solution. This method is proposed as a strategy to reduce the computational complexity of the ISAC signal design problem. In Section VI, we present the results of our simulations. This section includes a thorough discussion of the outcomes, providing insights into the effectiveness and efficiency of the proposed methods. Finally, Section VII concludes the paper.

§.§ Notations

The following notations are used in this paper. Boldface lower-case letter $\mathbf{s}$, upper-case letter $\mathbf{S}$, and non-bolded letter $s$ or $S$ indicate column
vectors, matrices, and scalar, respectively. $\mathcal{S}$ stands for a set and $|\mathcal{S}|$ is the size of the set. $(\cdot)^T$ and $(\cdot)^H$ denote the
transpose and the transpose-conjugate operations, respectively. $\Im\{\cdot\}$ and $\Re\{\cdot\}$ denote the imaginary and real part of a complex number. $\mathbb{R}$ and $\mathbb{C}$ stand for the real-number space and complex-number space, respectively. $\|\cdot\|_2$ denotes the $l_2$ norm. $\otimes$ and $\odot$ represent the Kronecker product and Hadamard product of two matrices or vectors, respectively. $\mathbf{I}_N$ represents a $N\times N$ identity matrix and $\mathbf{1}_N$ represents all ones vector of length $N$. $[x]^+$ stand for the function $\max\{x,0\}$. Finally, $\diag(\mathbf{s})$ stands for a diagonal matrix with
diagonal elements being $\mathbf{s}$.

§ PROBLEM DESCRIPTION AND GENERALIZATION

[Time signal]
[Frequency signal]
[Spatial signal]
Schematic of the proposed unified ISAC signaling design.

In this paper, we proceed under the assumption that the optimal sensing signal, denoted as $\mathbf{x}_0$, has already been identified. The goal for optimal ISAC signals is to ensure they are closely aligned with $\mathbf{x}_0$ to maintain satisfactory sensing performance. The core challenge, as defined in (P2), involves identifying a set of signals, $\mathcal{X}$, comprising $M$ elements, with each signal in the set having a dimensionality of $N$. These signals, corresponding to $N/2$ complex symbols, are intended for transmission over $N/2$ time, frequency, or spatial resources, as illustrated in Fig. <ref>.

To comprehensively address the ISAC signal design, we consider the following three scenarios without loss of generality:

* The first scenario involves designing $N/2$-dimensional complex ISAC signals without any channel state information (CSI). This scenario directly tackles the original problem (P2).

* The second scenario focuses on designing $N/2$-dimensional complex ISAC signals with the knowledge of CSI at the transmitter (CSIT). In this case, the $N/2$-dimensional signal undergoes independent fading, free from mutual interference. This could occur, for instance, when $N/2$-dimensional complex symbols within a signal are transmitted over orthogonal resources, such as different time slots or orthogonal sub-carriers. The optimization problem in this scenario can be described as a modified version of (P2), taking into account the CSIT. The modified problem formulation could be represented as follows:

\begin{align}
(\textbf{P3}):~\mathrm{Given}: &~M,N,P,\mathbf{A},\mathbf{x}_0,\epsilon\notag\\
\mathrm{Find}:&~\mathcal{X}=\{\mathbf{x}_1,\mathbf{x}_2,\cdots,\mathbf{x}_M\}\notag\\
\mathrm{Maximize}:&~\min_{\forall k,l,k\neq l} ||\mathbf{A}(\mathbf{x}_k-\mathbf{x}_l)||_2\\
\mathrm{Subject~to}:&~\frac{1}{M}\sum_{k=1}^{M}||\mathbf{x}_k||_2^2\leq P,\notag\\
&~||\mathbf{x}_k-\mathbf{x}_0||_2\leq \epsilon, \forall k\notag,
\end{align}

where $\mathbf{A}=\diag([\Re(\mathbf{h})^T,\Im(\mathbf{h})^T])\in\mathbb{R}^{N\times N}$ and $\mathbf{h}\in \mathbb{C}^{N/2}$ denotes the channel vector of $N/2$ independent fading channels without mutual interference.

* The third scenario corresponds to designing $N/2$-dimensional complex ISAC signals with CSIT, where the signal experiences channels with crosstalk. An example is the transmission of $N/2$-dimensional complex symbols through MIMO channels with inter-channel interferences (ICI). Let $N_r$ represent the degrees of freedom of the received signal, and $\mathbf{H}$ be a matrix in $\mathbb{C}^{N_r\times N/2}$, characterizing the channel. The optimization problem in this scenario involves a further modification of (P2), taking into account both CSIT and the unique channel characteristics due to crosstalk. It can be expressed as

\begin{align}
(\textbf{P4}):~\mathrm{Given}: &~M,N,P,\mathbf{B},\mathbf{x}_0,\epsilon\notag\\
\mathrm{Find}:&~\mathcal{X}=\{\mathbf{x}_1,\mathbf{x}_2,\cdots,\mathbf{x}_M\}\notag\\
\mathrm{Maximize}:&~\min_{\forall k,l,k\neq l} ||\mathbf{B}(\mathbf{x}_k-\mathbf{x}_l)||_2\\
\mathrm{Subject~to}:&~\frac{1}{M}\sum_{k=1}^{M}||\mathbf{x}_k||_2^2\leq P,\notag\\
&~||\mathbf{x}_k-\mathbf{x}_0||_2\leq \epsilon, \forall k\notag,
\end{align}

\begin{equation}
\mathbf{B}= \left[ \begin{array}{ccc}
\Re{(\mathbf{H})} & -\Im{(\mathbf{H})} \\
\Im{(\mathbf{H})} & \Re{(\mathbf{H})}\end{array} \right]\in\mathbb{R}^{2N_r\times N}.
\end{equation}

We make the singular value composition (SVD) of $\mathbf{B}$ as $\mathbf{B}=\mathbf{U}\pmb{\Lambda}\mathbf{V}^T$, where $\mathbf{U}$ is a $2N_r\times 2N_r$ complex unitary matrix,  $\pmb{\Lambda}$ is a $2N_r\times N$ rectangular diagonal matrix with non-negative real numbers on the diagonal, $\mathbf{V}$ is a $N\times N$ real unitary matrix. By introducing $\mathbf{s}_i\in\mathbb{R}^{N}$ and letting $\mathbf{x}_i=\mathbf{V}\mathbf{s}_i$, $i=1,2,\cdots,M$, we can express the term $||\mathbf{B}(\mathbf{x}_k-\mathbf{x}_l)||_2$ in (P4) as

\begin{equation}
\end{equation}

and the constraint of (P4)  as

\begin{equation}
\frac{1}{M}\sum_{k=1}^{M}||\mathbf{s}_k||_2^2\leq P,
\end{equation}

\begin{equation}
||\mathbf{s}_k-\mathbf{s}_0||_2\leq \epsilon.
\end{equation}

Then, we can reformulate $(\textbf{P4})$ as

\begin{align}
(\textbf{P5}):~\mathrm{Given}: &~M,N,P,\pmb{\Lambda},\mathbf{s}_0,\epsilon\notag\\
\mathrm{Find}:&~\mathcal{S}=\{\mathbf{s}_1,\mathbf{s}_2,\cdots,\mathbf{s}_M\}\notag\\
\mathrm{Maximize}:&~\min_{k,l\in [1,M],k\neq l} ||\pmb{\Lambda}(\mathbf{s}_k-\mathbf{s}_l)||_2\\
\mathrm{Subject~to}:&~\frac{1}{M}\sum_{k=1}^{M}||\mathbf{s}_k||_2^2\leq P,\notag\\
&~||\mathbf{s}_k-\mathbf{s}_0||_2\leq \epsilon, \forall k\notag.
\end{align}

Solving (P5) is equivalently solving (P4) as $\mathbf{x}_i=\mathbf{V}\mathbf{s}_i$, $i=1,2,\cdots,M$.

Upon careful examination of the formulations of (P2), (P3), and (P5), it becomes apparent that these ISAC signal shaping problems share notable similarities. Notably, (P2) can be viewed as a specific instance of (P3) when we set $\mathbf{A}$ to be the identity matrix $\mathbf{I}_{N}$. Similarly, (P3) can be considered a particular case of (P5) under the condition that $N_r$ is equal to $N/2$. This hierarchical relationship among the problems suggests that a solution to the more general (P5) would inherently encompass solutions to the more specific scenarios presented in (P2) and (P3). Therefore, our subsequent analysis will primarily focus on addressing and solving the general problem (P5). Through this approach, we aim to develop a comprehensive solution strategy that is applicable across the varying contexts of (P2) and (P3).

§ PROBLEM TRANSFORMATION

Problem (P5) is a modified sphere packing problem with weighted distance maximization and additional sensing constraints. We aim to propose a generalized solution to this problem. To fulfill this goal, we transform the original problem into a simplified version. To do so, we introduce several auxiliary matrices and vectors. In detail,  we stack $M$ $\pmb{\Lambda}$ matrices into

\begin{equation}
\mathbf{G}=[\overbrace{\pmb{\Lambda},\pmb{\Lambda},\cdots,\pmb{\Lambda}}^{M}]\in\mathbb{R}^{2N_r\times MN},
\end{equation}

and all variable vectors into

\begin{equation}
\mathbf{z}=[\mathbf{s}_1^T,\mathbf{s}_2^T,\cdots,\mathbf{s}_M^T]^T\in\mathbb{R}^{MN}.
\end{equation}

By introducing

\begin{equation}
\mathbf{D}_{\mathbf{z}}=\diag(\mathbf{z})\in\mathbb{R}^{MN\times MN},
\end{equation}

\begin{equation}
\mathbf{e}_k=\mathbf{g}_k\otimes \mathbf{1}_{N}\in\mathbb{R}^{MN},
\end{equation}

where $\mathbf{g}_k\in \mathbb{R}^{M}$ is a standard basis vector with all zeros except a one at the $k$th position, we can express the term $||\pmb{\Lambda}(\mathbf{s}_k-\mathbf{s}_l)||_2$ as

\begin{equation}
\end{equation}

Taking the square of the term $||\pmb{\Lambda}(\mathbf{s}_k-\mathbf{s}_l)||_2$, we have

\begin{equation}\label{eq15}
\begin{split}
||\pmb{\Lambda}(\mathbf{s}_k-\mathbf{s}_l)||_2^2&=(\mathbf{e}_k-\mathbf{e}_l)^T \mathbf{D}_{\mathbf{z}}^T\mathbf{G}^T\mathbf{G}\mathbf{D}_{\mathbf{z}}(\mathbf{e}_k-\mathbf{e}_l)\\
&=\tr (\mathbf{D}_{\mathbf{z}}^T \mathbf{R}_{\mathbf{G}}\mathbf{D}_{\mathbf{z}} \mathbf{W}_{kl})
\end{split}
\end{equation}

where $\mathbf{R}_G=\mathbf{G}^T\mathbf{G}$ and $\mathbf{W}_{kl}=(\mathbf{e}_k-\mathbf{e}_l)(\mathbf{e}_k-\mathbf{e}_l)^T$. According to $\tr(\mathbf{D}_\mathbf{u}\tilde{\mathbf{A}}{\mathbf{D}}_{\mathbf{v}}\tilde{\mathbf{B}}^T)=\mathbf{u}^T(\tilde{\mathbf{A}}\odot\tilde{\mathbf{B}})\mathbf{v}$ <cit.>,  with $\mathbf{D}_{\mathbf{u}}=\diag (\mathbf{u})$ and $\mathbf{D}_{\mathbf{v}}=\diag (\mathbf{v})$, we can rewrite (<ref>) as

\begin{equation}
\end{equation}

where $\mathbf{Q}_{kl}=\mathbf{R}_{G}\odot \mathbf{W}_{kl}^T\in\mathbb{R}^{MN\times MN}$.

The first constraint of (P5) can be written as

\begin{equation}
\mathbf{z}^T\mathbf{z}\leq MP.
\end{equation}

By introducing $\mathbf{E}_k=\diag(\mathbf{e}_k)$, we can rewrite the sensing constraints in (P5) as

\begin{equation}
\begin{split}
&\leq \epsilon^2, \forall k,
\end{split}
\end{equation}

\begin{equation}
\mathbf{z}_0=[\overbrace{\mathbf{s}_0^T,\mathbf{s}_0^T,\cdots,\mathbf{s}_0^T}^{M}]^T\in\mathbb{R}^{MN},
\end{equation}

and the equation holds due to the fact $\mathbf{E}_k=\mathbf{E}_k^T\mathbf{E}_k$.

By above skillful formula manipulations, we can derive a simplified but equivalent version of (P5) as

\begin{align}
(\textbf{P6}):~\mathrm{Given}: &~M,N,P,\pmb{\Lambda},\mathbf{z}_0,\epsilon\notag\\
\mathrm{Find}:&~\mathbf{z}\notag\\
\mathrm{Maximize}:&~\min_{\forall k,l,k\neq l} \mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}\\
\mathrm{Subject~to}:&~\mathbf{z}^T\mathbf{z}\leq MP,\notag\\
&~(\mathbf{z}-\mathbf{z}_0)^T\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\leq \epsilon^2, \forall k\notag.
\end{align}

By introducing an auxiliary variable $t$, we have the equivalent epigraph form <cit.> of (P6) as

\begin{align}
(\textbf{P7}):~\mathrm{Given}: &~M,N,P,\pmb{\Lambda},\mathbf{z}_0,\epsilon\notag\\
\mathrm{Find}:&~\mathbf{z}\notag\\
\mathrm{Maximize}:&~t\\
\mathrm{Subject~to}:&~\mathbf{z}^T\mathbf{z}\leq MP,\notag\\
&~\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}\geq t, \forall k,l, k\neq l,\notag\\
&~(\mathbf{z}-\mathbf{z}_0)^T\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\leq \epsilon^2, \forall k\notag.
\end{align}

By introducing a target minimum distance $d$, we can introduce a new problem closely related to (P7) as

\begin{align}
(\textbf{P8}):~\mathrm{Given}: &~M,N,P,\pmb{\Lambda},\mathbf{z}_0,\epsilon,d\notag\\
\mathrm{Find}:&~\mathbf{z}\notag\\
\mathrm{Minimize}:&~\mathbf{z}^T\mathbf{z}\\
\mathrm{Subject~to}:&~\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}\geq d, \forall k,l, k\neq l,\notag\\
&~(\mathbf{z}-\mathbf{z}_0)^T\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\leq \epsilon^2, \forall k\notag.
\end{align}

Given $\mathbf{z}^*$ as the solution to (P8), and letting ${P}^*$ represents  the average power, i.e., ${P}^*=\frac{1}{M}(\mathbf{z}^*)^T\mathbf{z}^*$, there are three cases.

* If ${P}^*$ is  equal to ${P}$, the solution to (P8) is exactly the solution to (P7) and we have $t=d$.

* If ${P}^*$ is smaller than ${P}$, the original (P7) is feasible with an objective value $t>d$. However, $z^*$ is not the optimal solution to (P7), because we can still increase the communication and sensing performance with an additional power budget (i.e., $P-P^*$).

* If ${P}^*$ is greater than ${P}$, the original (P7) is infeasible.

From our analysis, we identify an optimal $d^*$ using one-dimensional optimization methods, such as the bisection method, to ensure that $P^*=P$. Consequently, solving problem (P8) with this $d^*$ becomes equivalent to solving (P7). The crux of the issue lies in efficiently solving (P8), a QCQP problem with $MN$ variables and $\frac{M(M+1)}{2}$ constraints.

The QCQP challenge in (P8) can be approached using various optimization strategies, such as semidefinite relaxation (SDR) techniques <cit.>. SDR approximates a solution by transforming the original problem into a semi-definite programming (SDP) problem. Specifically for (P8), this involves setting $\mathbf{C} = \mathbf{z}\mathbf{z}^H$ and relaxing the rank-1 constraint of matrix $\mathbf{C}$. A feasible solution to the original problem is then derived by reinstating the rank-1 constraint, using methods like Gaussian randomization. However, the direct application of SDR is hindered by the sheer volume of quadratic constraints and optimization variables, such as $1024$ variables and $136$ constraints for $N=32, M=16$. The extensive number of quadratic constraints significantly reduces the likelihood of obtaining rank-1 solutions.

To overcome these challenges, this paper introduces an augmented Lagrangian and dual ascent (ALDA)-based algorithm. This algorithm is specifically designed to address the complexities associated with the large number of constraints and variables in (P8).

§ SIGNAL SET DESIGN USING ALDA ALGORITHM

In this section, we present the development of the ALDA algorithm, which employs augmented Lagrangian and dual ascent techniques to address (P8). Initially, we provide a concise overview of the fundamental principles underlying the dual ascent method. Following this, we delve into a detailed exposition of the ALDA algorithm.

§.§ Dual Ascent

Dual ascent is a strategy for solving convex optimization problems<cit.>. Let's consider a general constrained convex optimization problem with equality constraints and non-equality constraints as follows:

\begin{equation}
\begin{split}
\label{4-1}
\min_{\mathbf{x}\in \mathbb{R}} &~f(\mathbf{x})\\
\mathrm{subject~to:}&~c_i(\mathbf{x})=0,~ i \in \mathcal{I},\\
&~c_i(\mathbf{x})\leq0,~i \in \mathcal{M}.
\end{split}
\end{equation}

We can construct the Lagrangian function

\begin{equation}
\begin{split}
\label{4-2}
L(\mathbf{x},\pmb{\lambda},\mathbf{v})=&f(\mathbf{x})+\sum_{i\in \mathcal{I}}\lambda_ic_i(\mathbf{x})+\sum_{i\in \mathcal{M}}v_ic_i(\mathbf{x}),
\end{split}
\end{equation}

where $\pmb{\lambda}\triangleq\{\lambda_i\},\mathbf{v}\triangleq\{v_i\}$ are the dual variables. In the feasible domain of (<ref>) $\mathcal{D}=\{\mathbf{x}|c_i(\mathbf{x})=0, i \in \mathcal{I},c_i(\mathbf{x})\leq0,i \in \mathcal{M}\}$, the original function $f(\mathbf{x})$ is the upper bound of Lagrangian function (<ref>), i.e.,

\begin{equation}
\label{4-3}
\max_{\pmb{\lambda},\mathbf{v}}L(\mathbf{x},\pmb{\lambda},\mathbf{v})=f(\mathbf{x}), ~\mathbf{x}\in \mathcal{D}.
\end{equation}

Therefore, the optimization of the original problem can be equivalently written as

\begin{equation}
\label{4-4}
\min_\mathbf{x} f(\mathbf{x})= \min_\mathbf{x} \max_{\pmb{\lambda},\mathbf{v}}L(\mathbf{x},\pmb{\lambda},\mathbf{v}),~\mathbf{x}\in \mathcal{D}.
\end{equation}

The dual problem
with respect to the original problem in (<ref>) can be expressed as

\begin{equation}
\label{4-5}
\max_{\pmb{\lambda},\mathbf{v}}\min_\mathbf{x}L(\mathbf{x},\pmb{\lambda},\mathbf{v}).
\end{equation}

when the convex optimization problem satisfies strong duality, the dual problem has the same optimal solution as the original problem. Based on this, one can update $\mathbf{x}$ and $\pmb{\lambda},\mathbf{v}$ alternatively. Specifically, in the $j$th step, one can update $\mathbf{x}$ to be

\begin{equation}
\label{4-6}
\mathbf{x}^{j+1}=\arg \min_{\mathbf{x}} L(\mathbf{x},\pmb{\lambda}^{j},\mathbf{v}^{j})
\end{equation}

Afterwards, one can update $\pmb{\lambda},\mathbf{v}$ to be

\begin{equation}
\pmb{\lambda}^{j+1}=\arg \max_{\pmb{\lambda}} L(\mathbf{x}^{j+1},\pmb{\lambda},\mathbf{v}^{j})
\end{equation}

\begin{equation}
\mathbf{v}^{j+1}=\arg \max_{\mathbf{v}} L(\mathbf{x}^{j+1},\pmb{\lambda}^{j+1},\mathbf{v}).
\end{equation}

For convenience, the gradient ascent method is often used in updating the dual variables, i.e.,

\begin{equation}
\label{4-8a}\pmb{\lambda}^{j+1}=\pmb{\lambda}^{j}+\eta^j\Delta_{\pmb{\lambda}} L(\mathbf{x}^{j+1},\pmb{\lambda},\mathbf{v}^j)
\end{equation}

\begin{equation}
\label{4-8b}\mathbf{v}^{j+1}=\mathbf{v}^{j}+\eta^j\Delta_\mathbf{v}L(\mathbf{x}^{j+1},\pmb{\lambda}^{j+1},\mathbf{v})
\end{equation}

where $\eta^j$ is the step size in the $j$-th iteration, $\Delta_{\pmb{\lambda}}L(\mathbf{x}^{j+1},\pmb{\lambda},\mathbf{v}^j)$ denotes the derivative with respect to $\pmb{\lambda}$, and $\Delta_{\mathbf{v}}L(\mathbf{x}^{j+1},\pmb{\lambda}^{j+1},\mathbf{v})$ denotes the derivative with respect to $\mathbf{v}$. By iteratively updating these parameters using equations (<ref>), (<ref>) and (<ref>), the solutions to the dual problem and the primal problem can be obtained.

§.§ ALDA Algorithm for Solving (P8)

For the inequality constraint in (P8), we equivalently transform it into the following optimization problem by introducing slack variables,

\begin{align}
(\textbf{P9}):~\mathrm{Given}: &~M,N,P,\pmb{\Lambda},\mathbf{z}_0,\epsilon,d\notag\\
\mathrm{Find}:&~\mathbf{z},\mathbf{r},\mathbf{t}\notag\\
\mathrm{Minimize}:&~\mathbf{z}^T\mathbf{z}\\
\mathrm{Subject~to}:&~d-\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}+r_{k,l}=0 , \forall k,l, k\neq l,\notag\\
&~\|\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\|_2^2-\epsilon^2+t_k= 0, \forall k\notag,\\
\end{align}

where $\mathbf{r}\triangleq\{r_{k,l}\},\mathbf{t}\triangleq \{t_k\}$ are the introduced slack variables.
Then, we can write the augmented Lagrangian function for (P9) as

\begin{equation}
\begin{split}
\label{eq34}
&L(\mathbf{z},\mathbf{r},\mathbf{t},\pmb{\lambda},\mathbf{v},\mu)\\&~=\mathbf{z}^T\mathbf{z}-\sum_{l=1}^{M}\sum_{k=1,k\neq l}^{M} \lambda_{l,k}\left(\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}-r_{k,l}-d \right)\\
&~~~~~+\sum_{k=1}^{M}v_k \left(\|\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\|_2^2-\epsilon^2+t_k\right)\\
&~~~~~+\frac{\mu}{2}\left[\sum_{l=1}^{M}\sum_{k=1,k\neq l}^{M} \left(\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}-r_{k,l}-d \right)^2\right.\\
&~~~~~~~~~~~~~+\left.\sum_{k=1}^{M} \left(\|\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\|_2^2+t_k-\epsilon^2\right)^2\right],
\end{split}
\end{equation}

where $\pmb{\lambda}\triangleq\{\lambda_{l,k}\},\mathbf{v}\triangleq\{v_k\}$ are the dual variables, the quadratic part is the penalty item, and $\mu>0$ is the penalty parameter.

According to the principle of the dual ascent method, we will update the parameters alternately. In the $j$-th iteration step, given the dual variables $\pmb{\lambda}^j,\mathbf{v}^j$ and the penalty parameter $\mu^j$, we can update $\mathbf{z}$, $\mathbf{r}$, and $\mathbf{t}$ as follows.

§.§.§ Update $\mathbf{z}$, $\mathbf{r}$, and $\mathbf{t}$

According to (<ref>), we need to solve the following problem to find $\mathbf{z}^{j+1}$, $\mathbf{r}^{j+1}$ and $\mathbf{t}^{j+1}$:

\begin{equation}
\begin{split}
\label{4-11}
\min_{\mathbf{z},\mathbf{r},\mathbf{t}} &~L(\mathbf{z},\mathbf{r},\mathbf{t},\pmb{\lambda}^j,\mathbf{v}^j,\mu^j)\\
\mathrm{subject~to:}&~\mathbf{r}\succeq\mathbf{0},\mathbf{t}\succeq\mathbf{0}.
\end{split}
\end{equation}

The dominant terms in (<ref>) to update slack variables $\mathbf{r}$ and $\mathbf{t}$ can be written as

\begin{equation}
\begin{split}
&\min_{\mathbf{r}\succeq \mathbf{0}} ~\sum_{l=1}^{M}\sum_{k=1,k\neq l}^{M} -\lambda_{l,k}^j\left(\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}-r_{k,l}-d \right)\\
&~~~~~+\frac{\mu^j}{2}\sum_{l=1}^{M}\sum_{k=1,k\neq l}^{M} \left(\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}-r_{k,l}-d \right)^2,
\end{split}
\end{equation}

\begin{equation}
\begin{split}
&\min_{\mathbf{t}\succeq \mathbf{0}}~ \sum_{k=1}^{M}v_k^j \left(\|\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\|_2^2-\epsilon^2+t_k\right)\\
&~~~~~+\frac{\mu^j}{2}\sum_{k=1}^{M} \left(\|\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\|_2^2+t_k-\epsilon^2\right)^2.
\end{split}
\end{equation}

Taking the derivatives of the objective functions with respect to $r_{l,k},v_k$, we can derive the zero point to obtain the updates as

\begin{equation}\label{eq38}
\end{equation}

\begin{equation}\label{eq39}
\end{equation}

\begin{align}
\label{4-14}
&-\lambda_{l,k}^j\left(\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}-r_{k,l}-d \right)+\frac{\mu^j}{2}\left(\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}-r_{k,l}-d \right)^2\notag\\
\begin{array}{l}
-\lambda_{l,k}^j\left(\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}-d \right)+\frac{\mu^j}{2}\left(\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}-d \right)^2,~\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}-d-\frac{\lambda_{l,k}^j}{\mu^j}\leq0,\\-\frac{1}{2\mu^j}(\lambda_{l,k}^j)^2,~\rm{otherwise}.\quad\quad\quad\quad
\end{array}\right.
\\ \label{4-15}
v_k^j\left(\|\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\|_2^2-\epsilon^2\right)+\frac{\mu^j}{2}\left(\|\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\|_2^2-\epsilon^2\right)^2,~-\|\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\|_2^2+\epsilon^2-\frac{v_k^j}{\mu^j}\leq0,  \\-\frac{1}{2\mu^j}(v_k^j)^2,~\rm{otherwise}.
\end{array}\right.
\end{align}

Substituting (<ref>) (<ref>) into (<ref>) yields an equivalent form $L(\mathbf{z},\pmb{\lambda}^j,\mathbf{v}^j,\mu^j)$ without $\mathbf{r},\mathbf{t}$. The substitution terms of (<ref>) can be expressed as (<ref>) and (<ref>), based on which we can simplify $L(\mathbf{z},\pmb{\lambda}^j,\mathbf{v}^j,\mu^j)$ without $\mathbf{r},\mathbf{t}$ as

\begin{equation}
\begin{split}
\label{4-16}
&L(\mathbf{z},\pmb{\lambda}^j,\mathbf{v}^j,\mu^j)\\&~~~=\mathbf{z}^T\mathbf{z}+\sum_{l=1}^{M}\sum_{k=1,k\neq l}^{M} \phi\left(-\mathbf{z}^T\mathbf{Q}_{kl}\mathbf{z}+d,\lambda_{l,k}^j,\mu^j\right)\\&~~~~~~~~~~~~+\sum_{k=1}^{M}\phi\left(\|\mathbf{E}_k(\mathbf{z}-\mathbf{z}_0)\|_2^2-\epsilon^2,v_k^j,\mu^j\right),
\end{split}
\end{equation}

where the function $\phi(a,b,c)$ is defined as

\begin{equation}
\label{4-17}
\phi(a,b,c)\triangleq\left\{\begin{array}{l}
\end{array}\right.
\end{equation}

Then, the problem in (<ref>) becomes finding $\mathbf{z}$ such that $\mathbf{z}^{j+1}=\arg \min_\mathbf{z}L(\mathbf{z},\pmb{\lambda}^j,\mathbf{v}^j,\mu^j)$. Since the closed-form solution is difficult to obtain, we use the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm in <cit.> to solve it.

§.§.§ Update dual variables $\pmb{\lambda}$ and $\mathbf{v}$

Once the approximate solution $\mathbf{z}^{j+1}$ is obtained, according to the Karush–Kuhn–Tucker (KKT) conditions of the augmented Lagrangian <cit.>, the dual variables and the penalty parameters can be updated using the following formulas

\begin{equation}\label{eq44}
\lambda_{l,k}^{j+1}=\left[\lambda_{l,k}^{j}-\mu^j\left((\mathbf{z}^{j+1})^T\mathbf{Q}_{kl}\mathbf{z}^{j+1}-d\right)\right]^+,
\end{equation}

\begin{equation}\label{eq45}
\end{equation}

\begin{equation}\label{eq46}
\mu^{j+1}=\rho\mu^j,
\end{equation}

where $\rho>0$ is the parameter controlling the growth of the penalty parameter.
The detailed steps of the algorithm are summarized in Algorithm 1. It is worth noting that  the result $\mathbf{z}$ from Algorithm 1 will be scaled to meet the total
power constraint at the transmitter.

ALDA Algorithm for Solving (P8).

Initialize $\mathbf{z}^0=\mathbf{1}_{MN} ,  \lambda_{l,k}^0=0.5,\forall l, k, l \neq k, v_k^0=0.5,\forall k,\mu^0=10, \rho=2, $ and $j=0.$
1.Update the primal vector $\mathbf{z}^{j+1}$ by BFGS in <cit.>

\begin{equation}
\mathbf{z}^{j+1}=\arg \min_\mathbf{z}L(\mathbf{z}^j,\pmb{\lambda}^j,\mathbf{v}^j,\mu^j).\notag
\end{equation}

2. Update the dual vector $\lambda_{l,k}^{j+1}$ and $v_k^{k+1}$ through (<ref>) and (<ref>), respectively.
3. Update the penalty parameter through (<ref>).
4. $j=j+1$.

convergence criterion is met.

§.§ Complexity Analysis for ALDA

The main complexity comes from two aspects: 1) the computational complexity of ALDA algorithm at each iteration, and 2) the number of iterations required for the convergence denoted as $K_{ALDA}$. The complexity of each iteration mainly comes from the BFGS solving process and the update of the dual variables. The resulting complexity of the BFGS algorithm for a variable dimension of $MN$ is $O(M^2N^2K_{BFGS})$ with $K_{BFGS}$ being the number of iterations. The resulting complexity of the update of the dual variable is $O(M^4N^2)+O(M^3N^2)$. The total complexity is $O(M^2N^2K_{BFGS}K_{ALDA})+O(M^4N^2K_{ALDA})+O(M^3N^2K_{ALDA})$.

§ LOW-COMPLEXITY SOLUTION WITH BDPS

The computational complexity associated with ADLA is notably high, particularly in the context of adaptive signal shaping that depends on CSIT.  To reduce the computational complexity, we consider decomposing the original problem (P5) by BDPS. In detail, we split $N$ signal dimensions into $G$ groups, i.e., $\mathcal{N}_1,\mathcal{N}_2,\cdots, \mathcal{N}_G$, to carry $m_1,m_2,\cdots,m_G$ bits subjecting to separate power constraints $P_1, P_2,\cdots, P_G$ and  separate sensing constraints  $\epsilon_1, \epsilon_2,\cdots, \epsilon_G$.
The separate constraints are limited by the overall constraints, which can be expressed as

\begin{equation}
P=\sum_{g=1}^G P_g,
\end{equation}

\begin{equation}
\epsilon^2=\sum_{g=1}^G \epsilon_g^2.
\end{equation}

Let $\mathcal{S}_g$ represent the signal set of the $g$th signal dimension group and $M_g$ represent its set sizes  with $M_g=|\mathcal{S}_g|=2^{m_g}$, $\forall g=1,2,\cdots,G$, we have the following relationships

\begin{equation}
\end{equation}

\begin{equation}
\end{equation}

Let $\mathbf{v}_g$ represent the splitting signal vector whose elements are extracted from $\mathbf{s}_0$ corresponding to $g$th signal dimension group. Then, we can formulate small-size problems of finding $\mathcal{S}_g$ given $M_g$, $N_g=|\mathcal{N}_g|$, $\epsilon_g$, $\mathbf{v}_g$, and $\pmb{\Lambda}_g$, where $\pmb{\Lambda}^g\in\mathbb{R}^{2N_r\times N_g}$ is a matrix comprised of the $N_g$ columns corresponding the signal dimension group $\mathcal{N}_g$.  These small-size sub-problems can be transformed similarly to QCQP problems like (P5) and solved by the proposed ALDA algorithm.  The problem size is greatly reduced, leading to computational complexity reduction. Besides,  these subproblems can be solved in parallel, leading to a further reduction of computation time.

Based on the aforementioned BDPS,
the objective of (P5) can be written as

\begin{equation}
\min_{k,l\in [1,M],k\neq l} ||\pmb{\Lambda}(\mathbf{s}_k-\mathbf{s}_l)||_2=\sqrt{\sum_{g=1}^Gd_g^2}
\end{equation}

where $d_g$ denotes the minimum distance of signals in the $g$th signal dimension group and can be obtained by solving the $g$th subproblem. Then, the optimization can be expressed as

\begin{align}
(\textbf{P10}):~\mathrm{Given}: &~M,N,P,\epsilon,\mathbf{s}_0,\pmb{\Lambda}\notag\\
\mathrm{Find}:&~M_1,M_2,\cdots,M_G,\notag \\&~N_1,N_2,\cdots,N_G,\notag\\&~P_1,P_2,\cdots,P_G,\notag\\&~\epsilon_1,\epsilon_2,\cdots,\epsilon_G,\notag\\
\mathrm{Maximize}:&~ \sqrt{\sum_{g=1}^G d_g^2}\\
\mathrm{Subject~to}:&~\prod_{g=1}^G M_g=M,~\sum_{g=1}^{G}N_g=N\notag,\\
&~\sum_{g=1}^G P_g\leq P\notag,~\sum_{g=1}^G \epsilon_g=G\notag.
\end{align}

It is worth mentioning that (P10) is a mixed integer programming and its objective function involves solving a sequence of QCQP problems. To deal with it, we resort to heuristic algorithms such as genetic algorithm (GA), and particle swarm optimization (PSO) algorithm.

Let $N_{iter}$ represent the number of iterations and $N_{pop}$ represent the population size in GA or PSO algorithms, the computational complexity is around $O(N_{iter}N_{pop}\sum_{g=1}^GM_g^4N_g^2)$. Compared with the computational complexity of directly solving (P5), its computational complexity is greatly reduced.

§ SIMULATION RESULTS AND DISCUSSIONS

In this section, we provide thorough comparisons with existing approaches in ISAC signal design. As far as we are aware, our work is the first to address the generalized problem of ISAC signal design. i.e., (P2). For comparison,
we choose existing works on precoding/beamforming for ISAC because all modifications on the precoder/beamformer will affect the signal set as

\begin{equation}
\mathbf{x}_i=\mathbf{F}\mathbf{s}_i,
\end{equation}

where $\mathbf{F}$ can be regarded as a beamforming matrix (in the spatial domain) or a precoding matrix (in the time/frequency domain).

Specifically, we compare the existing signal designs that optimize the beamformer/precoder in MIMO scenarios, as the signal design considering the ICI is the most complicated. The following schemes are chosen for comparison.

* The symbol-level ISAC beamforming with interference reduction (SL-IR) in <cit.>.

* The block-level ISAC beamforming for maximizing minimum Euclidean distance (BL-MMED) with finite alphabet input in <cit.>.

* The block-level ISAC beamforming for Gaussian inputs (BL-Gaussian) in <cit.>.

In our simulation, we configure the carrier frequency at $24$ GHz. The transmitter is equipped with $N/2 = 32$ antennas, while the receiver has $N_r = 8$ antennas. The transmission power is set at $P = 1$. The signal set size is chosen as $M=4$. We model the downlink channel $\mathbf{H}$ as a flat Rayleigh fading channel, where each element of the channel matrix $\mathbf{H}$ follows a standard Complex Gaussian distribution. The CSIT is assumed to be perfectly known.
Furthermore, we utilize a spatial waveform for the sensing optimal waveform, designed to scan within a broad main lobe spanning $[-20^\circ, 20^\circ]$. We ensure that the sidelobe levels (SLL) are maintained below $15$ dB. These simulation parameters will remain constant unless stated otherwise.

CDF of the minimum Euclidean distances among the noise-free received vectors for different schemes, where $d=14$ and $\epsilon=0.3$.

Coverage performance of the proposed BDPS.

§.§ Communications Performance Evaluation

We conducted an evaluation of the minimum Euclidean distance performance across various signal design strategies. Over 1000 channel realizations, we plotted the cumulative distribution function (CDF) of the minimum Euclidean distances among all noise-free received vectors for different schemes and algorithms, as shown in Fig. <ref>, where the minimum distance $d=14$ in (P8) and $\epsilon=0.3$. The results demonstrate that our proposed signal design scheme outperforms others, consistently achieving a minimum Euclidean distance greater than $12.5$. The steep curve corresponding to the ALDA algorithm signifies its stability. In contrast, the minimum Euclidean distances for other schemes like BL-Gaussian, BL-MMD, and SL-IR barely reach $5$, with the probability Pr$\{d>5\}$ being almost zero.
This performance gap is attributed to their constrained modulation schemes, which hinder reaching a globally optimal solution. Our proposed approach, treating ISAC signal as a high-dimensional modulation signal design problem, effectively overcomes this limitation. When compared to traditional SDR-solving methods, it's observable that, despite SDR with randomization providing superior solutions than other schemes, they lack the stability of those achieved by ALDA. This highlights the performance limitations of the SDR algorithm in large-scale QCQP problems <cit.>.
Moreover, as a complexity-reducing algorithm, BDPS closely approaches ALDA’s performance. Fig. <ref> illustrates the convergence curve of the BDPS algorithm, showing rapid convergence within a few iterations.

SER of different signal design schemes.

The Euclidean distance is a critical metric in our study as it signifies the distinguishability between different symbols and is directly correlated with the SER. A larger Euclidean distance typically results in a lower SER. In Fig. <ref>, we present the SER performance of various design schemes. The findings echo our previous conclusions. At $4$ dB, the SER for our proposed design scheme reaches as low as $10^{-4}$, which is notably lower by two orders of magnitude compared to benchmark schemes. The proposed BDPS approach demonstrates SER performance comparable to the ALDA approach. Furthermore, we extended our analysis to simulate the SER curve for a signal set size of $M=8$. It was observed that increasing the signal set size adversely affects the SER. This deterioration is attributed to the reduction in the minimum Euclidean distance under the constant similarity constraints of the sensing waveform, which occurs as the number of ISAC symbols increases. Despite the higher transmission rate, our proposed scheme maintains superior performance over comparative algorithms at an SNR greater than $2$ dB.

§.§ Sensing Performance Evaluation

ISAC symbols designed using the proposed scheme.

Beampattern of ISAC symbols designed using the proposed scheme.

Secondly, we assess the sensing performance through the analysis of waveforms and beampatterns of the ISAC signals, as illustrated in Figs. <ref> and <ref>. The proposed scheme's ISAC symbols, as seen in Fig. <ref>, adhere closely to the reference signal under similarity constraints. This indicates that transmitting any ISAC symbol from the designed set ensures robust communication while satisfying sensing waveform requirements.

In the spatial domain, as observed in Fig. <ref>, the ISAC symbols align well with the main lobe range. By transmitting each ISAC symbol an equal number of times within a communication frame, we plotted the average beampattern over this period, shown as the red line in the figure. This average beampattern stabilizes and closely approximates the reference beampattern after signal accumulation. Fig. <ref> contrasts this with the average beampatterns of other schemes, highlighting the superiority of our proposed approach.

Moreover, wide beams, vital for scanning tasks, are evaluated using detection probability as a performance metric. Assuming a far-field target, we conducted $1000$ Monte Carlo simulations with the target's angle appearing randomly within the main lobe range of $[-20^\circ, 20^\circ]$. The results, depicted in Fig. <ref>, reveal that the proposed scheme's curve closely matches that of the reference waveform, owing to its high waveform similarity. Notably, the ALDA and BDPS schemes under our design achieve a significant gain of at least $6$ dB in detection probability compared to other schemes, affirming their effectiveness at a detection probability of $1$.

Average beampatterns of different schemes.

To further examine the influence of different channel conditions on the sensing performance of ISAC signals, we graphed the waveform similarity across various channel realizations, as depicted in Fig. <ref>. The data shows that the waveform similarity for the benchmark algorithms consistently remained above $0.6$, yet these values exhibited significant fluctuations in response to channel changes. Such variability implies that radar performance is highly sensitive to channel variations, leading to potential instability which is undesirable in radar applications. A similar issue is observed with the SDR-solving approach.
In contrast, the waveform similarity for the proposed ALDA scheme demonstrates remarkable stability, maintaining around $0.3$ with minimal fluctuation across $100$ channel realizations. This consistency underscores the exceptional stability of the sensing performance of our proposed scheme. In the face of dynamically changing communication channels, our ISAC signal design adeptly fulfills sensing requirements, showcasing outstanding robustness and adaptability to varying channel conditions.

Detection probability comparisons among different schemes.

Waveform similarity of different schemes under different channel realizations.

§.§ Performance Trade-off Between Sensing and

Our proposed scheme aims to identify ISAC signals that fulfill a set minimum Euclidean distance within a space constrained by sensing waveform similarity. This approach allows for a balanced performance trade-off between communication and sensing by adjusting the minimum Euclidean distance $d$ or the similarity constraint $\epsilon$. In our simulation, we kept the similarity constraint constant at $\epsilon=0.3$ while varying $d$ from 0 to 14. For comparison, we adjusted the trade-off parameter of the benchmark algorithms from $0.01$ to $0.99$.

The results of this performance trade-off are illustrated in Fig. <ref>. We observe that an increase in the minimum Euclidean distance correlates with a higher tolerance in similarity. As hypothesized, the proposed ALDA and BDPS algorithms substantially outperform other schemes. Our proposed scheme consistently exhibits a significantly lower signal similarity tolerance for any given minimum Euclidean distance. The smoothness of our scheme's curves reflects its effective balance between communication and sensing.

To provide further insight, we plotted additional curves in Fig. <ref>, using beampattern MSE and SER as metrics. The majority of the points on the ALDA and BDPS curves fall within a delineated black ellipse region. In this zone, both a SER below $10^{-4}$ and a beampattern MSE under $-10$ dB are attainable, further evidencing the proposed scheme’s ability to achieve an optimal compromise between reliable communication and accurate sensing.

Trade-off of the similarity and the minimum Euclidean distance.

Trade-off of the beampattern MSE and the SER, SNR = $6$ dB.

§.§ Impact of imperfect CSIT on ISAC Signal Performance

The signal design schemes discussed thus far are based on the premise of perfect CSIT, a condition that is often unattainable in real-world scenarios. Therefore, it is imperative to assess the signal's robustness in the face of channel estimation errors or imperfect CSIT. We consider the channel error model as defined in <cit.>:

\begin{equation}
\mathbf{H}_{im} = \mathbf{H} + \mathbf{H}_e,
\end{equation}

where $\mathbf{H}_e$ represents the matrix of channel errors, with each element following a complex Gaussian distribution with zero mean and variance $\sigma_e^2$. Assuming the application of Least Squares (LS) channel estimation, the channel estimation error $\sigma_e^2$ is proportional to the noise variance $\sigma_n^2$, specifically, $\sigma_e^2 = \eta\sigma_n^2$. In our simulations, we set $\eta = 0.3$ and compared these results with those obtained under perfect CSIT conditions, as depicted in Fig. <ref>. The introduction of channel errors leads to a degradation in the SER across all schemes. However, our proposed scheme still surpasses comparative schemes in performance, demonstrating its robustness in communication even under imperfect CSIT.

Furthermore, it is noteworthy that channel estimation errors can be interpreted as channel variations, transitioning from $\mathbf{H}$ to $\mathbf{H}_{im}$. Consequently, as indicated in Fig. <ref>, our proposed schemes exhibit resilience to such channel variations. This implies that the sensing performance of our proposed scheme remains consistently effective for both $\mathbf{H}$ and $\mathbf{H}_{im}$.

SER of different signal design schemes under imperfect CSIT.

The phase of the time signal obtained by the proposed scheme, where $d=2$ and $\epsilon=0.5$.

The phase of the time signal obtained by the proposed scheme for different similarity tolerance.

Ambiguity function of the LFM signal.

Ambiguity function of the designed ISAC symbols.

§.§ Extension to Time Domain ISAC Signal Set Design

As outlined in Section II, our research adopts a unified approach to ISAC signal design. In this section, we extend our proposed scheme to the design of time domain signals. We assume that each ISAC symbol is transmitted over $\frac{N}{2} = 32$ time slots. For our sensing reference, we employ the classical LFM signal, renowned for its excellent pulse compression and ambiguity properties, as detailed in <cit.>.
The time sequence of the LFM signal $\mathbf{x}_0$ is given by [3, eq.(1.21)],

\begin{equation}
x_0(n)=\sqrt{\frac{2P}{N}}e^{j\pi \frac{2n(n-1)}{N}},~n=1,2,\cdots,\frac{N}{2}.
\end{equation}

With the minimum Euclidean distance set to $d=2$ and the similarity constraint at $\epsilon=0.5$, Fig. <ref> displays the phase of signals generated by our proposed scheme. To assess similarity, we also include the phase of the reference LFM signal. The comparison reveals that each ISAC symbol closely approximates the reference signal, effectively balancing between meeting communication requirements and maintaining similarity.

Subsequently, in Fig. <ref>, we examine the impact of varying similarity constraints. For simplicity, we illustrate this using just one ISAC signal under each constraint. As anticipated, a smaller $\epsilon$ value implies a stronger similarity constraint, resulting in signals more closely resembling the LFM signal.

Finally, the ambiguity functions of both the LFM and the ISAC signals are presented in Figs. <ref> and <ref>. The ambiguity functions of the designed ISAC symbols exhibit a knife-edge shape, which is indicative of their enhanced preservation of the ambiguity properties inherent to the LFM signal.

§ CONCLUSION

This paper presents a unified approach to the optimization of ISAC signal set design. We developed a unified signal design framework applicable both with and without knowledge of CSIT. This framework encompasses both parallel and cross-talk channel scenarios. We transformed this unified challenge into a large-scale QCQP problem. To address its significant computational complexity, we introduced an efficient ALDA algorithm tailored for large-scale QCQP. Additionally, to further reduce complexity, we proposed a low-complexity problem decomposition method utilizing BDPS. Our simulation results highlight the effectiveness of our proposed signal designs, demonstrating their superiority over existing signal design methods in the ISAC domain.