Long GPS Memory Code Design via the Cross-Entropy Method

2.1 Binary Codes

In this paper, we design a set of binary spreading codes, where each code is a binary sequence. The set of such binary sequences is also referred to as a family of codes. In the context of GPS and GNSSs, each satellite in the constellation is assigned one of the unique codes in the family, and this code overlays the satellite’s navigation signal.

In this work, we use ±1 binary notation, with each bit represented as either a +1 or a –1. The purpose of this notation is for ease of representing mathematical expressions for auto- and cross-correlation, as defined in Section 2.3. Let L be the length of each binary code, and let K be the size of the code family, or the number of codes in the set. We represent the complete space of binary code families as χ = {−1, +1}^LK. Let the length-KL vector x ∊ χ represent a particular binary code family. We represent the i-th bit, or chip, in the complete code family as x[i] ∊ {−1, +1}. Additionally, we represent the k-th code in the family as x^(k), and the i-th chip of this particular code is represented as x^(k)[i].

2.2 Probability and Proposal Distributions

As discussed in more detail in Section 3, the CEM is a stochastic optimization method, which represents and optimizes a probability distribution over the design space, also called a proposal distribution. We represent the probability of a particular event E as ℙ{E}. Additionally 𝟙 {·} represents the indicator function that outputs a 1 if the inner argument is true and 0 if false.

We denote the proposal distribution over the design space as p_θ(x), where x ∊ χ. The proposal distribution p_θ(x) represents the focus of the search region of the optimization algorithm and is parameterized by the vector θ. During the optimization process, each iteration updates the parameter vector of the proposal distribution in order to shift the search region to better-performing codes in the design space χ. At a given iteration t, we denote the parameter vector at that iteration as θ^(t) and the corresponding proposal distribution as p_θ(t). More details about the proposal distribution representation specific to this work of designing binary spreading codes are provided in Section 3.1.

2.3 Auto- and Cross-Correlation

Given a code family x ∊ χ, we can examine the correlation between any two codes in the family. In particular, if we select code indices k, j ∊ {1, …, K} from the code family x, the periodic correlation between the corresponding respective code sequences x^(k) and x^(j) is defined as follows:

Rk,j(δ):=∑s=0L−1x(k)[s]⁢x(j)[(s+δ)⁢modL]1

where δ ∊ {0, 1, …, L – 1} represents the relative delay between the two codes being correlated.

When k ≠ j, the correlation in Equation (1) is with respect to two different code sequences in the code family x and is often referred to as a periodic cross-correlation. When k = j, the correlation is between a particular code k and a δ-shifted version of itself, which is often referred to as a periodic auto-correlation. For simplicity, we often represent the auto-correlation with respect to code k, or R_k,k(δ), as simply R_k(δ).

When considering an auto-correlation with a non-zero relative delay (i.e., δ ≠ 0), we refer to these correlation values as auto-correlation sidelobes. Additionally, we refer to all cross-correlations for all delays δ as cross-correlation sidelobes. For CDMA applications, including GPS/GNSS applications, we seek to minimize the magnitude of both the auto- and cross-correlation sidelobes.

3.1 Proposal Distribution Representation

A natural representation for the proposal distribution p_θ is a multivariate Bernoulli distribution, with the parameter vector θ (Botev et al., 2013). In particular, each chip in the code set x[i] is represented as a Bernoulli random variable, with parameter θ[i] ∊ [0,1] representing its probability of being sampled as a 1. Mathematically, the binary bit x[i] ~ Bern (θ[i]), where we have the following:

P{x[i]=c}={θ[i],for c=11−θ[i],for c=−12

For a particular code family sample x_j, the complete proposal distribution can be represented as a vector of independently sampled Bernoulli random variables:

pθ(xj)=∏i=0L⁢K−1θ[i]1⁢{xj[i]=1}⁢(1−θ[i])1⁢{xj[i]=−1}3

Correspondingly, for a binary family of K length-L codes, the code family has a total of KL chips, and we can thus represent the proposal distribution parameter vector as a length-(KL) vector of individual probabilities, i.e., θ ∊ [0,1]^KL.

3.2 Initializing and Sampling the Proposal Distribution

When setting up the CEM, we initialize the proposal distribution as follows:

θ(0):=0.5⋅1K⁢L4

where 1^d represents a length-d vector of ones. At each iteration t, we sample a population of N code families from the current proposal distribution:

X(t):={xj:xj∼pθ(t)∀j∈{1,…,N}}5

where each x_j is independently and identically distributed ∀j.

3.3 Evaluating Sampled Candidates

After we initialize the proposal distribution, at each iteration, the CEM evaluates the performance of each sample in the population X^(t) according to a user-defined objective f. In this work, we examine different objectives to evaluate the performance of a particular candidate spreading code family. The choice of code objective affects the weight assigned to different correlation levels in the code family, thereby affecting which correlation properties are preferred in the design optimization process.

In this work, we design families of binary spreading code to minimize the following overall objective:

f=max(fauto,fcross)6

where f_auto and f_cross represent the components of the objective evaluating only auto- and cross-correlation performance, respectively. Given that both the auto- and cross-correlation performances are simultaneously important for satellite navigation, this objective formulation balances between both correlation objectives simultaneously, without compromising one type of correlation performance for the other.

For the individual auto-correlation and cross-correlation performances, similar to code objectives in prior works (Raei et al., 2022; Winkel, 2011), we consider the mean auto- and cross-correlation sidelobes to the power p (MMAC-p):

fauto=(1K⁢(L−1)⁢∑k=1K∑δ=1L−1|Rk(δ)|p)1/p7

fcross=(1(K2)⁢L⁢∑k=1K∑j=k+1K∑δ=0L−1|Rk,j(δ)|p)1/p8

where K represents the size of the code family, L represents the length of each sequence, R_k represents the auto-correlation of the k-th code in the family, and R_k,j represents the cross-correlation of codes k and j. The summation of the auto-correlation terms from relative delays starts at δ = 1, rather than δ = 0; this is done to capture the auto-correlation sidelobe levels in the cost function, given that all binary codes will have a normalized auto-correlation of 1 when perfectly self-aligned (corresponding to δ = 0). The external 1/p power allows for the final correlation value to be less than L and allows for a more meaningful comparison to code correlation performance for different power levels p, as shown in the results in Section 4 and Table 1. While the 1/p power helps with the interpretation of the final correlation values, it is worth noting that the objective in Equation (6) represents an order-isomorphic transformation without the application of the 1/p power. Thus, the CEM-based optimization can be carried out without this power transformation while still achieving the same correlation performance for the final, optimized set of codes.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Correlation histograms for a spreading code family of 1023 chips and 3 codes, optimized with the MMAC-p objective for different power levels p

A corresponding table of correlation performance metrics for these results is also shown in Table 1. The CEM was run for a population size of N = 200, an elite rate of к = 0.5, and a smoothing parameter of α = 0.5. We observe that the higher-p correlation objective results in a greater penalization of larger correlation values, along with a reduced frequency of very low correlation levels as compared with the squared correlation objective (p = 2).

The different choices of power level p in Equations (7) and (8) reflect the weight given to different sizes of absolute correlations, with larger p values assigning greater weights to larger correlation levels. In the experimental results in Section 4, we explore the effect of different choices of p on the shape of the correlation histogram after optimization with the CEM.

Additional desired metrics could be incorporated within the code design objective function, including odd correlation, relative Doppler frequency correlation metrics, and intersystem interference with other GNSS signals (Wallner et al., 2006). These additional metrics could be incorporated as new arguments within the maximization function of Equation (6) with relative scaling factors if desired, or these objectives could be combined via a weighted summation as reported by Soualle et al. (2005). Our prior work examined code design in the presence of relative Doppler effects (Yang et al., 2024). In the present work, we consider code design using the objective in Equation (6), under zero-Doppler conditions. The reason for this shift is two-fold. First, the primary focus of this work is to examine the performance of the CEM, making a more straightforward choice of code design objective a good starting point. Second, from our own empirical observations and those of other prior works, it has been found that for codes with a Gaussian-like histogram of correlation values (as observed with memory codes, but not Gold codes), the peak correlation values of the histogram tend, on average, to shift toward zero under non-zero Doppler conditions (Winkel, 2011).

3.4 Finding the Elite Set

At this point, we have evaluated each candidate code according to the user-defined objective f, discussed in Section 3.3. Now, among all of the evaluated candidates, we sort the objective values from smallest to largest:

f(1)≤f(2)≤⋯≤f(N−1)≤f(N)9

where better performance is indicated by a lower objective value.

For an elite ratio κ, we determine the top N_elite = [κN] candidates in the population, or the elite set. More formally, we can define the elite set with respect to a critical performance level of γ(t):=f(Nelite), where the elite set Xelite(t) consists of the set of all candidates in the population with equivalent or better performance than this critical level:

Xelite(t):={xj∈X(t):f(xj)≤γ(t)}10

3.5 Updating the Proposal Distribution via Maximum Likelihood Estimation

Once the elite set has been found, the CEM updates the proposal distribution by evaluating the maximum likelihood distribution over the elite set, which corresponds to solving the following program:

θ^ML(t):=argmaxθ∑j=1N1⁢{xj∈Xelite(t)}⁢lnpθ(xj)11

The classical algorithm for the CEM directly uses this maximum likelihood solution to update the proposal distribution in the next iteration, i.e., θ(t+1)←θ^ML(t). In this work, we generalize this update step to incorporate smoothing as discussed in Section 3.6.

Solving the maximum likelihood estimate for the proposal distribution parameter in Equation (11) may not always be possible analytically. However, for the rep-resentation of the proposal distribution specified in Section 3.1, as a multivariate Bernoulli distribution with independent variables, the maximum likelihood program in Equation (11) can be found both analytically and in a computationally efficient manner, by simply finding the mean vector among the elite set:

θ^ML(t)=∑j=1N𝟙⁢{xj∈Xelite(t)}⁢xj|Xelite(t)|12

Note that the denominator will generally be equivalent to N_elite, unless the critical performance level γ^(t) is shared by multiple candidate solutions in the population. For continuous objective values, however, this case is rare unless the algorithm has nearly converged owing to the degeneracy of the proposal distribution. As a result, in this work, we take the top N_elite candidate solutions after sorting the population by objective, for simplicity in implementation, thus maintaining a consistent elite set size of N_elite.

The population size N and elite ratio к are chosen by the algorithm designer. A larger population allows for more exploration of the design space at each iteration, but corresponds to more computations per iteration. The elite ratio к ∊ (0,1) represents a balance between exploration and exploitation for the algorithm. Because a larger elite ratio considers more high-performing candidates when fitting the proposal distribution, this case generally results in greater variance in the updated proposal distribution, thereby inducing greater exploration in the next iteration. Similarly, because a smaller elite ratio fits the updated proposal distribution to a fewer number of candidates, this case induces the CEM to narrow its search more quickly and to rely more heavily on information obtained from previous iterations. Generally, increasing the elite ratio results in increased exploration and thus better final performance on the user-defined objective, although at the cost of a greater number of algorithm iterations.

3.6 Incorporating Smoothing in the Proposal Distribution Update

Smoothing can be incorporated within the CEM (Botev et al., 2013) to improve the final optimization solution and avoid premature convergence. With smoothing, the proposal distribution parameter θ^(t) is not directly updated to be the maximum likelihood solution as in Equation (11); rather, a convex weighting is applied between the current distribution parameter and the maximum likelihood solution:

θ(t+1)←α⁢θ^ML(t)+(1−α)⁢θ(t)13

where θ^ML(t) denotes the maximum likelihood estimate as in Equation (11) and α ∊ (0, 1] denotes the smoothing parameter. Generally, a smaller smoothing parameter α reduces the algorithm’s convergence rate owing to smaller updates of the proposal distribution parameter vector at each iteration.

Prior research has examined the impact of various smoothing parameter selections on CEM convergence properties for discrete binary optimization (Costa et al., 2007; Hu & Hu, 2009). In particular, Costa et al. (2007) derived constant and iteration-dependent definitions for the smoothing parameter that ensure that the CEM will find the optimal solution with high probability. However, despite these theoretically guaranteed convergence results, the resulting significantly slower convergence oftentimes renders these smoothing parameter definitions impractical for large-scale combinatorial optimization problems; thus, for the navigation code design problem, we focus on selecting a constant smoothing parameter.

Several multi-objective optimization strategies also exist to search for Pareto-optimal solutions across multiple objectives (Deb, 2011; Deb et al., 2002). However, upon examining these strategies in our prior work (Mina & Gao, 2019), when scaling the problem size for navigation spreading code applications, we experienced challenges in maintaining a broad local Pareto-optimal front. In particular, maintaining an accurate representation of the Pareto-optimal subspace, which encompasses different levels of non-dominated, joint auto- and cross-correlation performance via a spread of candidate solutions, was observed to be increasingly challenging as the problem size expanded to a higher-dimensional space for navigation code design, with at least 1000 chips in length and family sizes of at least 30 codes. Performing large-scale multi-objective optimization is currently an area of active research (Qin et al., 2021; Tian et al., 2020, 2021), and indeed, we believe that a multi-objective navigation code design strategy would be a fruitful avenue for future work.

4.1 Spreading Code Correlation Histogram for Different Objectives

Figure 3 depicts auto- and cross-correlation histograms for CEM-optimized code families, with different powers p of the objective function defined in Section 3.3. In particular, each bin of the histogram corresponds to an absolute correlation value, and the bin height corresponds to the number of occurrences of that correlation value in the set of all possible correlation sidelobes within the code family (for all δ, k, and j in the auto- and cross-correlations in Equations (7) and (8)). We additionally compare the histogram for the CEM-optimized codes with that of a randomly generated code family, with each chip equally likely to be assigned either +1 or −1 in the family. The purpose of this comparison is to demonstrate how the CEM shapes the correlation histogram of the final, optimized code family from the initialized set of codes, via the initial proposal distribution defined in Equation (4). For this run of the CEM, we used a population size, elite rate, and smoothing parameter of N = 200, к = 0.5 (corresponding to N_elite = 100), and α = 0.5, respectively. The same code family size could be run for a larger population size and smaller smoothing parameter, which generally results in better final performance levels.

In Figure 3, we observe that the higher powers of p result in greater penalization of the larger correlation objectives and a corresponding reduction in the maximum correlation level. With different values of p, we did not observe significant changes in the number of iterations required for CEM convergence or in the total computation time. We notice from Figure 3 that, in order to compensate for the reduced number of larger correlations, the optimization for higher powers of p = 4 and p = 8 results in a greater fraction of medium-sized correlation values, in comparison to the optimization with p = 2 and with random codes.

For these same results, Table 1 presents the corresponding correlation metrics. For each power level p, the values in the left auto-correlation columns are evaluated as in Equation (7) whereas the values in the right cross-correlation column are evaluated as in Equation (8). Table 1 further indicates that when the CEM is applied with the MMAC objective for larger values of p, we observe smaller maximum correlation values; meanwhile, when we optimize with the MMAC-2 objective, we observe smaller mean absolute and squared correlation levels. Interestingly, the codes optimized with the MMAC-4 objective had a smaller mean eighth power cross-correlation in comparison with the MMAC-8 objective; however, we also observe that the overall MMAC-8 objective performance (maximum between the mean auto- and cross-correlation values to the eighth power as in Equation (6)) for the MMAC-8-optimized codes is better, with a lower value of 40.48.

Indeed, we observe in both Figure 3 and Table 1 that the choice of power level p reflects a balance between the average correlation performance and the worst-case correlation performance. Table 1 also demonstrates that all of the CEM-optimized code families outperform random codes on all correlation objectives.

4.2 Comparison Against the NES and GA

In our prior work, we performed code optimization using the NES as well as a GA baseline under the MMAC-2 objective. We examined the performance for a range of code families, including those with lengths of up to 1031 bits and family sizes of 3 to 31 codes (Mina & Gao, 2022). For these algorithms, we had a population size of 100 and allowed both algorithms to run for 10000 iterations. The maximization objective specified for the NES and GA in our prior work is equivalent to a squared and normalized version of the objective specified in Equation (6) with p = 2. As a result, given that squaring and normalization operations are order-preserving for positive values, the optimization objective in our prior work is an order-isomorphic transformation of the MMAC-2 objective and is thus equivalent to performing optimization with respect to the MMAC-2 objective. Thus, for simplicity, we present the results in terms of the MMAC objective (with p = 2), used in the remainder of our experimental results.

In Figure 4, we compare the MMAC-2 performance of the CEM with respect to the NES and GA. We further indicate the average squared correlation level for random codes, i.e., L. For the CEM, we utilized an elite rate of к = 0.5 and a smoothing parameter of α = 0.5, while the results for the NES and GA optimization results were obtained from our prior work (Mina & Gao, 2022). All three optimization strategies were provided the same population size of 100 code families. The primary computational cost for all three optimization strategies arises from evaluating the auto- and cross-correlations of each candidate. As a result, by using the same population size, each iteration of the GA, NES, and CEM requires approximately the same computation time.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Comparison of the CEM with the NES and a GA baseline, from Mina and Gao (2022) (a) length-1023 (b) length-1031

The NES and GA were allowed to run for 10000 iterations, whereas the CEM converged in fewer than 1000 iterations. For both code lengths, we observe that the NES demonstrates a significant performance improvement in comparison to the GA, with lower mean-squared correlation values. However, for families of 10 codes and larger, we observe that the CEM implementation consistently obtains better performance than the NES, despite having fewer optimization iterations.

While the NES and GA were allowed to run for 10000 iterations, we observed that our CEM converged in significantly fewer iterations, with fewer than 1000 iterations. For both the length-1023 and length-1031 codes, in Figure 4, we observe that the NES demonstrates a significant performance improvement with respect to random codes, especially in comparison to the GA; however, for families of 10 codes and larger, we observe that the CEM consistently obtains a lower MMAC-2 objective than the NES, without the need for the network architecture of Mina and Gao (2022) and with fewer optimization iterations. While we did not stress-test the choice of the hyperparameters к, α, greater performance improvements could likely be achieved by the CEM by reducing the smoothing parameter and potentially increasing the elite rate as well.

4.3 Comparison Against Gold and Weil Codes

With the CEM, we further optimized binary code families on the MMAC-2 objective for code lengths ranging from 1023 chips to longer codes of 8191 and 10223 chips. For each code length, we designed code family sizes of 3 to 50 codes, comparing against well-chosen families of Gold and Weil codes of the same length. With the CEM, we chose a smaller population size of N = 100, an elite rate of к = 0.5, and a smoothing parameter of α = 0.5; however, better performance can be expected with a larger population size and/or a smaller smoothing parameter (e.g., N = 1000 and/or α = 0.1).

Based on the definitions for Gold (Gold, 1967; Misra & Enge, 2012) and Weil (Rushanan, 2006) code families, the complete set of Gold codes corresponds to a total of K = L + 2, whereas the complete set of Weil codes corresponds to a total of K = (L − 1)/2 codes. To provide a concrete example, these sets correspond to a total of 8193 Gold codes for length L = 8191 and a total of 5111 Weil codes for length L = 10223. Thus, to select code families of size K = 3 to 50 codes, we sampled 10000 random families of Gold and Weil codes at these lengths, with each sampled family consisting of K codes. Then, we selected the best-performing code families to compare against in this section. We further compared the performance of the CEM with the original versions of the Gold and Weil codes at their original lengths, without any truncation or extension, which disturbs their desirable correlation properties.

Figures 5 and 6 demonstrate the performance of the CEM codes optimized on the MMAC-2 objective function, while comparing against the best-performing equal-length Gold and Weil codes from 10000 randomly sampled families. In particular, Figure 5 compares the CEM with length-1023 and length-2047 Gold codes as well as length-1031 and length-2053 Weil codes, while Figure 6 compares the CEM with length-8191 Gold and Weil codes as well as length-10223 Weil codes.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

Comparison of the CEM (N = 100, к = 0.5, α = 0.5) with the best-performing Gold and Weil codes from 10000 random sampled families for lengths of 1023 to 2053 bits (a) length-1023 (b) length-1031 (c) length-2047 (d) length-2053

We note that the CEM finds better-performing codes on the MMAC-2 objective for all code lengths and family sizes, in comparison to the well-chosen families of Gold and Weil codes. Furthermore, the CEM-optimized codes simultaneously maintain both lower mean-squared auto- and cross-correlation sidelobes than the Gold and Weil codes for all code lengths and family sizes.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

Comparison of the CEM (N = 100, к = 0.5, α = 0.5) with the best-performing Gold and Weil codes from 10000 random sampled families for long memory codes of 8191 and 10223 bits (a) length-8191 (b) length-10223

While the Gold and Weil codes maintain nearly consistent auto- and cross-correlation performance, because the CEM is tailored to the particular chosen code family size, the algorithm is able to find better-performing solutions for smaller subsets of code families. We observe that the CEM-optimized codes obtain better squared auto- and cross-correlation performance in comparison to all families of Gold and Weil codes.

Interestingly, in both Figures 5 and 6, we observe that the Weil codes maintain a consistent mean-squared auto-correlation performance for all tested code family sizes. In particular, the length-1031 Weil codes in Figure 5(b) have a distinctively low mean-squared auto-correlation objective; however, in this instance, the corresponding cross-correlation objective is significantly greater, leading to a larger overall MMAC-2 objective. We further notice that the CEM balances the auto- and cross-correlation objectives in order to perform well on the overall maximization objective, by finding codes that have nearly identical mean- squared auto- and cross-correlation performance. This behavior of balancing between the two correlation objectives was similarly observed in our prior work (Mina & Gao, 2022). Overall, in both Figures 5 and 6, we observe that the CEM consistently finds code families that outperform the best-performing Gold and Weil codes on the MMAC-2 objective, for all code lengths and family sizes. Furthermore, for all Gold codes and for all Weil codes of length-2053 or longer, the CEM-optimized codes simultaneously demonstrate lower mean-squared sidelobes for both the auto- and cross-correlation objectives.

The Gold and Weil algorithmic codes are defined for a large number of codes, with a total of 8193 codes for length-8191 Gold codes and a total of 5111 codes for length-10223 Weil codes. However, when tailoring the design to a smaller family size of ~ 50 codes for satellite navigation, we expect to see an improvement in the correlation performance. Indeed, we observe in Figure 6 that the Gold and Weil codes maintain nearly consistent auto- and cross-correlation performance for the different family sizes. However, because the CEM is tailored to the particular chosen code family size, the algorithm is able to find better performing solutions for smaller subsets of code families, and we observe improved squared auto- and cross-correlation performance in comparison to the best Gold and Weil code families.