Spreading Code Sequence Design via Mixed-Integer Convex Optimization

Cross-Correlation

For any two i, j ∊ 0, ..., m–1, we define the cross-correlation between the two codes x_i and x_j as a length-n vector (xⁱ * x^j), given by the following:

${\left(x^i\star x^j\right)}_k=\sum _{s=0}^{n-1}{x}_s^i{x}^j{}_{{\left(s+k\right)}_{\text{mod}n}},\hspace{0.17em}\hspace{0.17em}k=0,\dots ,\hspace{0.17em}n-1$1

Here, (xⁱ * x^j)_k is the inner product between xⁱ and a k-circularly shifted version of x^j. The cross-correlation may also be defined in terms of negative shifts, i.e., the cross-correlation of xⁱ with x^j at shift –k is equal to (xⁱ * x^j)_n–k for any k. Note that the periodic cross-correlation in Equation (1) also exhibits the following symmetry: (xⁱ * x^j)_k = (x^j * xⁱ)_n–k.

In this work, we exclusively consider the periodic cross-correlation, although an extension to the aperiodic case is straightforward (see, e.g., the work by Song et al. (2016)). The cross-correlation levels are sometimes referred to as cross-correlation sidelobes.

Autocorrelation

We refer to the cross-correlation of xⁱ with itself as the autocorrelation of xⁱ. Note that for any binary sequence y ∊ {±1}ⁿ, the zero-shift autocorrelation is as follows:

${(y\star y)}_0=\sum _{s=0}^{n-1}{y}_s^2=n$2

which is a constant that does not depend on y.

Spreading Code Design

The goal here is to choose x such that (xⁱ * x^j)_k is suitably close to zero for all i, j, and k, with the exception of zero-shift autocorrelation values. In this work, we achieve this goal by minimizing a convex function of the relevant autocorrelation and cross-correlation values. To formalize this, we consider the pairwise cross-correlation function h : {±1}^n×m → ℝ^n×m×m, which we define as follows:

$h{\left(x\right)}_i^{j,k}={\left(x^i\star x^j\right)}_k,\hspace{1em}i,j=0,\dots ,m-1,\hspace{1em}k=0,\dots ,n-1$3

Here, the first dimension of h(x) is indexed by a subscript, and the last two dimensions are indexed by superscripts. The function h maps the spreading code family x to all pairwise cross-correlations, at every shift. Minimizing every component of h(x) involves a set of competing objectives, as reducing one cross-correlation value may come at the expense of increasing a different cross-correlation value.

For ease of notation, we do not exclude the zero-shift autocorrelation values in the definition of h unless otherwise specified, because they are constant. In addition, we note that by symmetry, the autocorrelation may be represented by just ⎿n/2⏌ unique terms.

Mixed-Integer Convex Optimization Problem

Here, we aim to achieve a trade-off by minimizing a convex scalar function of the autocorrelation and cross-correlation values, given by g : ℝ^n×m×m → ℝ ∪ {∞}. In addition, we may have a convex regularization or constraint function r : {±1}^n×m → ℝ ∪ {∞}. The functions g and r may represent constraints by taking a value of infinity for disallowed values of the cross-correlation and code, respectively. A comprehensive background on convex functions and convex optimization may be found in the book by Boyd and Vandenberghe (2004). We may write the spreading code design problem as the following optimization problem:

$\begin{array}{cc}\text{minimize}& g\left(h\right(x\left)\right)+r\left(x\right)\\ \text{subject}\text{to}& x∊{\{\pm 1\}}^{n\times m}\end{array}$4

In general, this combinatorial optimization problem is challenging. We note that the composition of h and g is nonconvex even though g is convex, because h is nonconvex.

In Section 3.2, we show how Equation (4) may be converted into the minimization of a convex function subject to binary constraints, following the approach of Yang et al. (2023a). This conversion is achieved by introducing auxiliary variables and constraints and replacing h with an equivalent linear function of auxiliary variables. Treating the problem as the minimization of a convex function over the nonconvex set of binary variables enables search heuristics based on convex optimization. Those methods, which are discussed in Section 4, utilize the fact that when the binary constraints are removed, the problem in Equation (4) becomes a convex optimization problem, which can be efficiently solved (Boyd & Vandenberghe, 2004).

Product of Binary Variables

First, let us consider two binary variables: u, V ∊ {±1}. The product uv is not a convex function of u and v. However, suppose that we introduce a new auxiliary variable w. Then, we have w = uv for any of the four possible combinations of u and v, if and only if w satisfies the four linear inequality constraints (Yang et al., 2023a, 2023b):

$w\le -u+v+1$5

$w\le -v+u+1$6

$w\ge -1-u-v$7

$w\ge -1+u+v$8

We refer to these constraints as linking constraints, as they relate the auxiliary variable w to the binary variables u and v. Note that the linking constraints define a convex set, corresponding to the intersection of four halfspaces.

Cross-Correlation of Binary Variables

Recall that the cross-correlation between two binary sequences xⁱ and x^j, given by Equation (1), is a sum of products of binary variables. Thus, we may introduce auxiliary variables z^i,j ∊ {±1}^n×n that each satisfy the linking constraints in Equations (5)–(8). This approach leads to the following linking constraints:

$z_{s,k}^{i,j}\le -x_s^i+x_{{(s+k)}_{\hspace{0.17em}\mathrm{mod}\hspace{0.17em}n}}^j+1$9

$z_{s,k}^{i,j}\le -x_{{(s+k)}_{\hspace{0.17em}\mathrm{mod}\hspace{0.17em}n}}^j+x_s^i+1$10

$z_{s,k}^{i,j}\ge -1-x_s^i-x_{{(s+k)}_{\hspace{0.17em}\mathrm{mod}\hspace{0.17em}n}}^j$11

$z_{s,k}^{i,j}\ge -1+x_s^i+x_{{(s+k)}_{\hspace{0.17em}\mathrm{mod}\hspace{0.17em}n}}^j$12

that is, $z_{s,k}^{i,j}=x_s^i{x^j}_{{\left(s+k\right)}_{\text{mod}n}}$. Then, subject to these constraints, we may represent the cross-correlation in Equation (1) using the following linear function:

${\left(x^i\star x^j\right)}_k=\sum _{s=0}^{n-1}{z}_{s,k}^{i,j}\hspace{1em}k=0,\dots ,n-1$13

We will write the auxiliary variables as a four-dimensional array z ∊ ℝ^n×n×m×m, where the first two components are indexed by subscripts and the last two components are indexed by superscripts. Note that the number of unique auxiliary variables is the number of unique pairs of binary variables, which is $\left(\begin{array}{c}nm\\ 2\end{array}\right)$.

We may now define a linearized pairwise cross-correlation function f:ℝ^n×n×m×m → ℝ^n×m×m as follows:

$f{\left(z\right)}_i^{j,k}=\sum _{s=0}^{n-1}{z}_{s,k}^{i,j},\hspace{1em}i,j=0,\dots ,m-1,\hspace{1em}s,k=0,\dots ,n-1$14

Similar to Equation (3), the first dimension of f(x) is indexed by a subscript, and the last two dimensions are indexed by superscripts.

We may now replace the nonconvex function h in Equation (4) with the linearized function f. Let C denote the convex set such that (x, z) ∊ C implies that all of the linking constraints in Equations (9)–(12) are satisfied for all i, j, k, s. Then, the spreading code optimization problem in Equation (4) may be written as the following equivalent problem:

$\begin{array}{cc}\text{minimize}& g\left(f\right(z\left)\right)+r\left(x\right)\\ \text{subject}\text{to}& x∊{\{\pm 1\}}^{m\times n},\hspace{1em}(x,z)∊C\end{array}$15

The problem in Equation (15) is a MICP, as we have integer-valued constraints on only some of the variables (x) and the problem is convex when the binary constraints are removed. The latter conclusion follows from the fact that the composition of a linear function f with a convex function g is convex (Boyd & Vandenberghe, 2004, Section 3.2.2).

Mean of Squares

We first consider the mean of squared autocorrelation and cross-correlation values, given as follows:

$g\left(c\right)=\frac{1}{n\left(m+\left(\begin{array}{c}m\\ 2\end{array}\right)\right)}\sum _{i=0}^{m-1}\sum _{j=i}^{m-1}\sum _{k=0}^{n-1}{\left(c_k^{i,j}\right)}^2$19

The argument c ∊ ℝ^n×m×m represents the auto- and cross-correlation values, where $c_k^{i,j}={\left(x^i\star x^j\right)}_k$. Note that we sum from j = i to m – 1 because of the symmetry property with periodic cross-correlation, as discussed in Section 3.1, where (xⁱ * x^j)_k = (x^j * xⁱ)_n–k. In some contexts, this objective is referred to as the ISL (He et al., 2012). There are $\left(\begin{array}{c}m\\ 2\end{array}\right)$ cross-correlation and m autocorrelation vectors, each of length n. The mean-of-squares objective g(h(x)) is a nonconvex quartic function of x. However, with the use of the linearized cross-correlation, g(f(x)) becomes a convex quadratic function of the auxiliary variables z, and Equation (15) becomes a MIQP.

Balanced Mean of Squares

In the mean-of-squares objective, there are more cross-correlation terms than autocorrelation terms, which can lead to codes with cross-correlation values that are relatively larger than the autocorrelation values. To address this, Mina and Gao (2022) proposed minimizing the larger of the mean of squared autocorrelations and the mean of squared cross-correlations. We refer to this as the balanced mean of squares, which corresponds to the following choice of g :

$g\left(c\right)=\max \left\{\frac{1}{n\left(\begin{array}{c}m\\ 2\end{array}\right)}\sum _{i=0}^{m-1}\sum _{j=i+1}^{m-1}\sum _{k=0}^{n-1}{\left(c_k^{i,j}\right)}^2,\frac{1}{m(n-1)}\sum _{i=0}^{m-1}\sum _{k=1}^{n-1}{\left(c_k^{i,i}\right)}^2\right\}$20

Because g is the maximum of two convex quadratic functions, it is also convex. Moreover, the optimization problem in Equation (15) becomes a mixed-integer second-order cone program (MISOCP) (Boyd & Vandenberghe, 2004; Gurobi Optimization, LLC, 2022).

Peak Sidelobe Level

Another choice of objective is the PSL, which is the maximum absolute value of the autocorrelation and cross-correlation values (He et al., 2012). This objective corresponds to the function g given as follows:

$g\left(c\right)=\underset{{}_{i,j,k}}{\max }\left|c_k^{i,j}\right|$21

Because g is convex and piecewise linear in this case, the problem in Equation (15) becomes a mixed-integer linear program.

Mean of the p-th Power Absolute Correlation Levels.

Instead of minimizing a mean of squares, we may consider minimizing a mean of absolute autocorrelation and cross-correlation values taken to the power of some positive value p, based on the generalized correlation objective defined by Winkel (2011). Although Winkel (2011) defined a generalized correlation objective for even values of p, by considering the absolute correlation magnitude, we can choose odd values of p as well, leading to the following choice of g:

$g\left(c\right)=\frac{1}{n\left(m+\left(\begin{array}{c}m\\ 2\end{array}\right)\right)}\sum _{i=0}^{m-1}\sum _{j=i}^{m-1}\sum _{k=0}^{n-1}{\left|c_k^{i,j}\right|}^p$22

Similar to the PSL objective, this objective may be used to more heavily penalize large correlation peaks by selecting larger values of p (e.g., p ≥ 4). Furthermore, this objective can be easier to minimize than the PSL objective, because it is smooth. When p = 2, we recover the mean of squares. Here, the optimization problem in Equation (15) may be represented as a MISOCP.

Balance Constraints

For a particular code sequence, the difference between the number of +1 and −1 entries is called the balance of the sequence (Winkel, 2011). In many cases, it is desirable for the number of +1 and −1 entries in each code to be close to the same or to have a balance that is close to 0. For example, for the GPS L1C signal, the length-10,230 spreading codes are designed to have a balance of zero (Rushanan, 2007). Incorporating balance constraints ensures that the statistical properties of each code are similar to those of random noise. Moreover, reducing the mean absolute level of the codes reduces the power needed to transmit the signal.

Balance constraints may be implemented by adding linear inequality constraints to the optimization problem in Equation (15). Formally, we set r to be an indicator function:

$r\left(x\right)=\{\begin{array}{cc}0& a\le \sum _{k=0}^{n-1}{x}_k^i\le b,\hspace{1em}i=0,\dots ,m-1\\ \infty & \text{otherwise}\end{array}$23

where a and b are fixed parameters. Here, a and b represent the maximum allowable excess −1 and +1 entries, respectively.

A scalar penalty term may also be used instead of constraints. For example, we may use a quadratic imbalance penalty:

$r\left(x\right)=\rho \sum _{i=0}^{m-1}{\left(\sum _{s=0}^{n-1}{x}_s^i\right)}^2$24

where ρ > 0 is a penalty parameter. When n is even and ρ is sufficiently large, the optimal solution will have codes with equal numbers of +1 and −1.

Autocorrelation Sidelobe Zero

A binary sequence y ∊ {±1}ⁿ that satisfies the (ASZ) property exhibits minimal autocorrelation at a shift of 1 (and equivalently, at a shift of −1). This corresponds to the condition (y * y)₁ = 0 in the even-length case and |(y * y)₁| = 1 in the odd-length case. The ASZ property can improve the ability of a receiver to accurately identify the originating satellite of a received signal. If a receiver detects similar correlations at shifts of 0 and 1, it may erroneously lose the signal lock. The Galileo constellation employs even-length spreading codes that have been designed to satisfy the ASZ property (Wallner et al., 2007).

The ASZ condition is described by a set of linear equality constraints. In the even-length case, we may take g to be the following indicator function:

$g\left(c\right)=\{\begin{array}{cc}\tilde{g}\left(c\right)& c_1^{i,i}=0,\hspace{1em}i=0,\dots ,m-1\\ \infty & \text{otherwise}\end{array}$25

where $\tilde{g}:{ℝ}^{n\times m\times m}\to ℝ\cup \left\{\infty \right\}$ is a convex function, such as the mean-of-squares objective in Equation (19). In the odd-length case, we may set the right-hand sides of the equality constraints to either +1 or –1. The ASZ constraints may also be converted into a penalty term, similar to the balance constraints. For example, we may take the following:

$g\left(c\right)=\tilde{g}\left(c\right)+\rho \sum _{i=0}^{m-1}{\left(c_1^{i,i}\right)}^2$26

where $\tilde{g}$ is convex and ρ > 0 is a penalty parameter.

Peak Correlation Level Constraint

Minimizing the mean of squares in Equation (19) or its variants, such as the balanced mean of squares in Equation (20), can lead to codes with autocorrelation and cross-correlation values that are small on average. However, the codes may have a relatively large peak absolute correlation, in comparison with codes designed to minimize the PSL objective (Equation (21)).

Conversely, the PSL objective seeks to minimize the worst-case correlation. The Gold codes, for example, have been proven to exhibit favorable PSL properties. However, PSL minimization can be overly conservative, resulting in codes with average correlation values that exceed those of a code family specifically optimized to have a small mean of squares.

One strategy for achieving a trade-off is to minimize the mean of squares, subject to a maximum absolute correlation constraint. That is, we may take the following:

$g\left(c\right)=\{\begin{array}{cc}\tilde{g}\left(c\right)& -c^{\max }\le c_k^{i,j}\le c^{\max },\hspace{1em}(i,j,k)∊J\\ \infty & \text{otherwise}\end{array}$27

where $\tilde{g}$ is convex, c^max is a fixed parameter that represents the maximum allowable PSL value, and J is the index set on which the constraints are enforced. Typically, J would exclude the zero-shift autocorrelation values, as those are constant and equal to n. Thus, we would take the following:

$J=\left\{\right(i,j,k)\mid i\ne j\text{or}(i=j\text{and}k=0\left)\right\}$28

Block Selection Strategy

We consider a random block selection strategy in our experiments. In each iteration, we select B of the nm binary variables at random. However, we limit the selected variables to be members of only K of the codes, in order to reduce the difficulty of solving the block optimization problem. We have applied this limitation because the number of auxiliary variables that must be considered grows quadratically with K (Yang et al., 2023a). In this work, we select blocks by first choosing K codes at random and then randomly choosing B / K variables from each of the K selected codes.

In general, the block size B and number of codes K should be chosen such that the block optimization problem can be solved quickly, for example, using Gurobi. For problems with very long code lengths, for example, with n on the order of 10, 000, it may not be possible to choose blocks such that the block optimization problems can be effectively solved, owing to the complexity of the resulting MICPs. In this case, it may be more efficient to simply solve the block optimization problems using an exhaustive search.

Figure 1 illustrates the block selection scheme for an example with m = 10 codes, each of length n = 50. Each cell represents one of the binary variables. Variables that have been selected as part of the current block are shaded in black, and variables that are held fixed are not shaded. In this example, the block size is B = 15, and the variables to be chosen are restricted by choosing 5 variables at random, from only 3 of the codes.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 1

Illustration of block selection in BCD for the case of m = 10 codes, with code length n = 50, block size B = 15, and block variables selected from K = 3 codes

Black cells represent selected variables, and white cells represent variables that are not selected.

Handling Constraints

The BCD method may be modified to handle the balance and ASZ constraints described in Subsection 3.4. One potential issue is that the block optimization problem may be infeasible, if the initial code family used to initialize BCD does not satisfy the constraints. However, as long as the initial code family satisfies the constraints, all following BCD iterations will also satisfy the constraints. An initial code that satisfies the balance constraints may be generated by taking any code family and flipping the signs of a subset of code entries, such that the sum of each code is zero (in the even-length case) or +1 or −1 (in the odd-length case). The method proposed by Yang et al. (2024) may be used to generate an initial code family that satisfies the ASZ property.

4.2 Theoretical Results

The code design optimization problem in Equation (4) is an NP-hard combinatorial optimization problem that, in the worst case, requires an assessment of 2^nm possible combinations of binary code entries. Because the objective function is non-increasing between iterations, the BCD method is guaranteed to converge. However, it is not guaranteed to converge to the global optimum. Theoretical lower bounds on the optimal objective value are available for certain choices of the objective function. For example, in the absence of constraints, the Welch bound provides a lower bound on the PSL (Welch, 1974).

Comparison Against NES and Gold Codes

We minimized the balanced mean-of-squares objective using the BCD method and compared the results against the Gold codes and the NES method proposed by Mina and Gao (2022). Figure 2 shows the objective value per iteration. The BCD method finds a code family that outperforms the Gold codes after 59 iterations and a code family that outperforms the one found by NES after 296 iterations.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Balanced ISL objective value as a function of BCD iteration

We also note that for the mean-of-squares objective in Equation (19), the BCD method found a code family with an objective value of 909.85, thus outperforming the Gold codes, which have a mean-of-squares objective value of 953.87. The NES method did not consider the mean-of-squares objective value.

Balanced vs. Nominal Mean of Squares

We also compared codes optimized for the balanced mean-of-squares objective (Equation (20)) against codes optimized for the nominal mean-of-squares objective (Equation (19)). Figure 3 presents histograms of the obtained absolute autocorrelation and cross-correlation values for all shifts. The compared codes were both found using the same BCD method described above. We see that the code optimized for the balanced objective has similar distributions of absolute autocorrelation and cross-correlation values, as expected. In contrast, the code optimized for the nominal mean of squares has larger values for cross-correlation than autocorrelation.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Frequency of absolute autocorrelation and cross-correlation values when minimizing the mean of squares (top) vs. balanced mean of squares (bottom) using BCD with a subset size of 15

The problem involved finding m = 31 codes, each of length 1023.

The code optimized for the balanced objective has a slightly lower average cross-correlation value (1012.17 vs. 992.10), but as shown in Figure 4, the distributions of cross-correlation values are nearly identical for the two codes. Whereas the balanced objective can be used to trade off autocorrelation performance for cross-correlation performance when m is small (e.g., 3) (Mina & Gao, 2022), the results suggest that this trade-off may not be achievable for m = 31.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Frequency of absolute cross-correlation values when minimizing the mean of squares vs. balanced mean of squares using the same BCD method

Balance Constraints

Finally, we compared the BCD method with and without balance constraints, for both the mean-of-squares and balanced mean-of-squares objectives. Because the code length is odd, the balance constraints ensure that the sum of each code is either +1 or −1.

For the balanced mean-of-squares objective, BCD found a code with an objective value of 990.27 when no balance constraints were enforced. The resulting code family has a maximum sum of 83, which is a significant imbalance. When balance constraints were enforced, the BCD method found a code family with a balanced mean-of-squares value of 991.18, which is only slightly larger.

For the mean-of-squares objective, BCD found a code family with an objective value of 909.85 without balance constraints and 911.34 with balance constraints. The code family optimized without balance constraints had a maximum sum of 47.

In both cases, imposing the balance constraints lead to optimized codes with slightly larger but comparable objective values.