Est. read time: 4 minutes | Last updated: March 02, 2026 by John Gentile

Contents

Open In Colab

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from rfproto import plot

import datetime
import dateutil.parser

GPS

# Parameters for GPS L1 C/A signal
code_length = 1023                    # Length of C/A Gold code / PRN (chips)
prn_rate    = 1000                    # PRN repitition rate (Hz)
prn_period  = 1.0 / prn_rate          # Time to correlate PRN over (seconds)
chip_rate   = prn_rate * code_length  # C/A code chip rate (Hz)
bit_rate    = 50                      # Navigation data bit rate (bit/sec)
fc_l1       = 1575.42e6               # GPS L1 Center Frequency (Hz)

sps = 2 # Samples/symbol (oversampling rate) of input I/Q data
fs = sps * chip_rate # input sample rate of I/Q file (Hz)
utc_start_time = dateutil.parser.parse("2023-11-03T18:48:00+00:00")

num_seconds = 1
num_samples = int(num_seconds * fs)
# NumPy complex64 type is a complex number data type composed of two single-precision (32-bit) floating-point numbers, 
# one for the real component and one for the imaginary component
iq_in = np.fromfile("./data/nov_3_time_18_48_st_ives-short.fc32", dtype=np.complex64, count=num_samples, offset=0)

# G2 phase selector taps for each SV (1-37) from GPS ICD-200
# Each SV uses two taps from the G2 register to generate its unique C/A code
G2_TAPS: dict[int, tuple[int, int]] = {
    1: (2, 6),
    2: (3, 7),
    3: (4, 8),
    4: (5, 9),
    5: (1, 9),
    6: (2, 10),
    7: (1, 8),
    8: (2, 9),
    9: (3, 10),
    10: (2, 3),
    11: (3, 4),
    12: (5, 6),
    13: (6, 7),
    14: (7, 8),
    15: (8, 9),
    16: (9, 10),
    17: (1, 4),
    18: (2, 5),
    19: (3, 6),
    20: (4, 7),
    21: (5, 8),
    22: (6, 9),
    23: (1, 3),
    24: (4, 6),
    25: (5, 7),
    26: (6, 8),
    27: (7, 9),
    28: (8, 10),
    29: (1, 6),
    30: (2, 7),
    31: (3, 8),
    32: (4, 9),
    33: (5, 10),
    34: (4, 10),
    35: (1, 7),
    36: (2, 8),
    37: (4, 10),
}

def generate_ca_code(sv: int):
    """Generate the GPS C/A (Coarse/Acquisition) code for a given satellite vehicle.

    The C/A code is a 1023-chip Gold code generated from two 10-bit Linear Feedback
    Shift Registers (G1 and G2). Each satellite has a unique code determined by
    different tap selections on the G2 register.

    G1 polynomial: x^10 + x^3 + 1 (taps at 3 and 10)
    G2 polynomial: x^10 + x^9 + x^8 + x^6 + x^3 + x^2 + 1 (taps at 2, 3, 6, 8, 9, 10)

    Reference: https://www.mathworks.com/matlabcentral/fileexchange/14670-gps-c-a-code-generator
    """
    if sv < 1 or sv > 37:
        raise ValueError(f"SV must be in range 1-37, got {sv}")

    # Get tap selections for this SV (convert to 0-indexed)
    tap1, tap2 = G2_TAPS[sv]
    tap1_idx = tap1 - 1
    tap2_idx = tap2 - 1

    # Initialize both registers with all 1s
    g1 = np.ones(10, dtype=np.int8)
    g2 = np.ones(10, dtype=np.int8)

    # Generate the C/A code
    ca_code = np.zeros(code_length, dtype=np.int8)

    for i in range(code_length):
        # Output is G1[9] XOR G2[tap1] XOR G2[tap2]
        g1_out = g1[9]
        g2_out = g2[tap1_idx] ^ g2[tap2_idx]
        ca_code[i] = g1_out ^ g2_out

        # Compute feedback bits before shifting
        # G1 feedback: positions 3 and 10 (indices 2 and 9)
        g1_feedback = g1[2] ^ g1[9]

        # G2 feedback: positions 2, 3, 6, 8, 9, 10 (indices 1, 2, 5, 7, 8, 9)
        g2_feedback = g2[1] ^ g2[2] ^ g2[5] ^ g2[7] ^ g2[8] ^ g2[9]

        # Shift registers right and insert feedback at position 0
        g1[1:] = g1[:-1]
        g1[0] = g1_feedback

        g2[1:] = g2[:-1]
        g2[0] = g2_feedback

    ## Convert from {0, 1} to {1, -1} (0 -> 1, 1 -> -1)
    #return 1 - 2 * ca_code
    return (2 * ca_code) - 1

# Precompute freq. domain PRN code for each SV at proper SPS
prn_fft_sv = []
for i in range(1, 38):
    prn_i = generate_ca_code(i)
    # Upsample PRN code to match input rate and convert to complex with a zero imaginary part
    prn_up = np.repeat(prn_i, sps).astype(complex)
    # In frequency domain cross-correlation, one part needs to be conjugated, so do that now
    prn_freq = np.conj(np.fft.fft(prn_up))
    
    prn_fft_sv.append(prn_freq)

# Compute number of samples to cross-correlate against
num_corr = int(fs * prn_period)
assert len(prn_fft_sv[0]) == num_corr

First, we show the “traditional” parallel code phase search approach, which efficiently searches across Doppler offsets and correlates in the frequency domain:

r[n]⊛x[−n]=IFFT{FFT{r[n]}(FFT{x[n]})∗}r[n] \circledast x[-n] = \text{IFFT} \bigg\{ \text{FFT}\{r[n]\} (\text{FFT}\{x[n]\})^{*} \bigg\} 
# Define Doppler shift range to search (in Hz)
doppler_shifts = np.linspace(-5000, 5000, 100)  # -5 kHz to 5 kHz, 100 steps

# Array to store correlation results: [doppler shifts, time delays]
correlation_results = np.zeros((len(doppler_shifts), num_corr))

# TODO: show search via:
#for sv_idx, sv in enumerate(prn_fft_sv):
sv_idx = 31
sv = prn_fft_sv[sv_idx]

for dop_idx, doppler_shift in enumerate(doppler_shifts):
    # Time vector for Doppler adjustment
    t = np.arange(num_corr) / fs
    # Carrier wipeoff: apply Doppler shift with complex exponential (phasor)
    doppler_adjustment = np.exp(-1j * 2 * np.pi * doppler_shift * t)
    adjusted_signal = iq_in[:num_corr] * doppler_adjustment
    
    # Frequency domain cross-correlation
    fft_signal = np.fft.fft(adjusted_signal)
    cross_corr_freq = fft_signal * sv
    cross_corr_time = np.fft.ifft(cross_corr_freq)

    # Store magnitude of correlation
    correlation_results[dop_idx, :] = np.abs(cross_corr_time)

# Normalize correlation results to [0, 1] to make peak stand out
Z = correlation_results / np.max(correlation_results)

# Define axes: convert samples to C/A code chips
x = np.arange(num_corr) / sps  # PRN code phase in chips
y = doppler_shifts / 1000       # Doppler shifts in kHz
X, Y = np.meshgrid(x, y)

# Create figure with larger size for better detail
fig = plt.figure(figsize=(14, 9))
ax = fig.add_subplot(111, projection='3d')

# Plot with jet colormap, fine stride, and shading
surf = ax.plot_surface(
    X, Y, Z,
    cmap='jet',
    rstride=1,
    cstride=4,
    antialiased=True,
    linewidth=0,
    alpha=0.95,
)

ax.set_xlabel('Code Phase (chips)', fontsize=12, labelpad=10)
ax.set_ylabel('Doppler Shift (kHz)', fontsize=12, labelpad=10)
ax.set_zlabel('Normalized Correlation', fontsize=12, labelpad=10)
ax.set_title(f'GPS L1 C/A Acquisition — SV {sv_idx + 1}', fontsize=14, fontweight='bold')

# Viewing angle: elevated and rotated to show the peak clearly
ax.view_init(elev=15, azim=225)

# Tighten z-axis to emphasize the peak
ax.set_zlim(0, 1.15)

# Find peak and annotate
peak_dop_idx, peak_phase_idx = np.unravel_index(np.argmax(correlation_results), correlation_results.shape)
peak_phase_chips = peak_phase_idx / sps
peak_doppler_khz = doppler_shifts[peak_dop_idx] / 1000
peak_mag = correlation_results[peak_dop_idx, peak_phase_idx]

annotation = (
    f"Corr: {peak_mag:.1f}\n"
    f"Phase: {peak_phase_chips:.1f} chips\n"
    f"Doppler: {doppler_shifts[peak_dop_idx]:.1f} Hz"
)
ax.text(
    peak_phase_chips, peak_doppler_khz, 1.12,
    annotation,
    fontsize=8, fontfamily='monospace',
    ha='center', va='bottom',
    bbox=dict(boxstyle='round,pad=0.4', facecolor='white', edgecolor='black', alpha=0.85),
    zorder=10,
)
ax.plot(
    [peak_phase_chips, peak_phase_chips],
    [peak_doppler_khz, peak_doppler_khz],
    [1.0, 1.12],
    color='black', linewidth=1, linestyle='--', zorder=9,
)

fig.colorbar(surf, shrink=0.55, aspect=12, pad=0.1, label='Normalized Correlation')
plt.tight_layout()
plt.show()

png

Instead of iterating over every doppler bin, performing carrier wipeoff, then performing frequency domain cross-correlation, we can perform joint frequency and delay correlation by:

  1. FFT the received signal once (length N)
  2. Build a Hankel index matrix where each column is the FFT circularly shifted by a different Doppler offset- this is the key trick. A circular shift of ll DFT samples in the frequency domain, X((k−l))NX((k-l))_{N}, is equivalent to multiplying by a complex exponential (frequency shift) in the time domain, x(n)∗ej2πNlnx(n) * e^{j\frac{2\pi}{N}ln}.
  3. Vectorized correlation where we multiply fft(b) * conj(fft_shifted(a)) for all Dopplers simultaneously via matrix operations, then IFFT column-wise.

However it should be seen that, unlike the above traditional method where we can arbitrarily define the Doppler search space (e.g. more or less carrier wipeoff frequency steps), this joint approach has a fixed Doppler bin spacing dictated by the PRN period (NN samples long) of:

Δf=fs/N\Delta f = f_{s} / N

For coherent correlation, this will almost always equal the PRN rate (1000 Hz in this GPS case). We can calculate coherent integration loss of a given Doppler bin spacing as it’s the same as the DFT (each bin in the DFT is a correlation over the time interval of the sequence to the particular frequency in that bin): when your Doppler hypothesis is offset from truth by Δf\Delta f, the correlation magnitude follows a sinc envelope:

∣R(Δf)∣=∣sinc(π⋅Δf⋅Tcoh)∣| R(\Delta f)| = | \text{sinc}(\pi \cdot \Delta f \cdot T_{coh})|

Where TcohT_{coh} is the coherent integration time, which for GPS of 1ms gives: | Point | Formula | Doppler Bin Spacing | | —– | ——- | ——————- | | First null | 1/Tcoh1/T_{coh} | 1000 Hz | | 3 dB point |  0.44/Tcoh~0.44/T_{coh} | 440 Hz | | 1 dB loss |  0.29/Tcoh~0.29/T_{coh} | 290 Hz |

Essentially longer coherent integration time, and thereby narrower Doppler bins, improves acquisition SNR at the cost of increased search space and computational cost (e.g. acquisition time). Since DSSS processing gain follows 10⋅log⁡10(N)10\cdot \log_{10}(N), we can expect that doubling or halving the coherent integration period results in +3dB and -3dB processing gain, respectively. The tradeoff also needs to factor in Probability of False Alarm (PFA) vs probability of detection.

def fdcorr(a, b, frange=None):
    """Joint frequency and delay correlation using the circular shift property of the DFT.

    Efficiently computes circular cross-correlation between two vectors over all
    possible frequency offsets (Dopplers) by exploiting the fact that a frequency
    shift in time is equivalent to a circular index shift in the frequency domain.
    Instead of looping over each Doppler bin, a Hankel-like index matrix creates
    all Doppler-shifted FFTs simultaneously, enabling a fully vectorized correlation.

    Based on: Dan Boschen, "FDCORR Joint frequency and delay correlation", 2009

    Parameters
    ----------
    a : array_like
        Received signal vector (sampled at fs = 2*pi).
    b : array_like
        Correlation code vector (sampled at fs = 2*pi).
    frange : tuple of (float, float), optional
        Frequency search range as (low, high) in radians, within [-pi, pi].
        Limiting this improves efficiency when Doppler offsets are small
        relative to the Nyquist bandwidth.

    Returns
    -------
    fd_surface : ndarray, shape (N, num_freqs)
        2D correlation surface: rows are delay bins, columns are frequency bins.
    freq_axis : ndarray
        Normalized frequency values (rad/sample) aligned to the frequency axis.
    """
    a = np.asarray(a).ravel()
    b = np.asarray(b).ravel()
    N = max(len(a), len(b))

    fft_a = np.fft.fft(a, N)

    # Frequency axis: centered from -pi to pi (equivalent to fftshift ordering)
    freq_index = np.arange(-np.floor((N - 1) / 2), np.floor(N / 2) + 1, dtype=int)
    freq_axis = freq_index * 2 * np.pi / N

    if frange is not None:
        mask = (freq_axis >= frange[0]) & (freq_axis <= frange[1])
        freq_axis = freq_axis[mask]
        freq_index = freq_index[mask]

    num_freqs = len(freq_index)

    # Build Hankel index matrix: each column is the FFT of 'a' circularly shifted
    # by a different Doppler offset. This is the core trick — a circular shift in
    # the frequency domain corresponds to multiplication by a complex exponential
    # (frequency shift) in the time domain.
    fft_idx = np.arange(N).reshape(-1, 1)   # (N, 1)
    shifts = freq_index.reshape(1, -1)      # (1, num_freqs)
    hankel_idx = (fft_idx + shifts) % N     # (N, num_freqs) circular indices

    fft_dopplers = fft_a[hankel_idx]        # All Doppler-shifted FFTs at once

    # Circular correlation: ifft(fft(b) * conj(fft_shifted(a))) for all Dopplers
    fft_b = np.fft.fft(b, N).reshape(-1, 1)  # (N, 1) broadcast across freq columns
    fd_surface = np.fft.ifft(fft_b * np.conj(fft_dopplers), axis=0)

    return fd_surface, freq_axis


# --- Run joint frequency-delay acquisition for a single SV ---
sv_idx = 31  # SV 32 (0-indexed)
prn_code = generate_ca_code(sv_idx + 1)
prn_up = np.repeat(prn_code, sps).astype(complex)

# Limit Doppler search to +/-5 kHz expressed in normalized frequency (rad/sample)
# Normalized freq = 2*pi * f_hz / fs
max_doppler_hz = 5000
frange_norm = (-2 * np.pi * max_doppler_hz / fs, 2 * np.pi * max_doppler_hz / fs)

fd_surface, freq_axis_norm = fdcorr(iq_in[:num_corr], prn_up, frange=frange_norm)

# Convert axes to physical units
freq_axis_hz = freq_axis_norm * fs / (2 * np.pi)  # Normalized rad/sample -> Hz
delay_chips = np.arange(fd_surface.shape[0]) / sps  # Samples -> chips

Z = np.abs(fd_surface)
Z = Z / np.max(Z)  # Normalize

X, Y = np.meshgrid(freq_axis_hz / 1000, delay_chips)  # kHz and chips

fig = plt.figure(figsize=(14, 9))
ax = fig.add_subplot(111, projection='3d')

surf = ax.plot_surface(
    X, Y, Z,
    cmap='jet',
    rstride=4,
    cstride=1,
    antialiased=True,
    linewidth=0,
    alpha=0.95,
)

ax.set_xlabel('Doppler Shift (kHz)', fontsize=12, labelpad=10)
ax.set_ylabel('Code Phase (chips)', fontsize=12, labelpad=10)
ax.set_zlabel('Normalized Correlation', fontsize=12, labelpad=10)
ax.set_title(f'GPS L1 C/A Acquisition (fdcorr) — SV {sv_idx + 1}', fontsize=14, fontweight='bold')
ax.view_init(elev=20, azim=225)
ax.set_zlim(0, 1.15)

# Find peak location and annotate on the surface
peak_idx = np.unravel_index(np.argmax(np.abs(fd_surface)), fd_surface.shape)
peak_delay_chips = peak_idx[0] / sps
peak_freq_hz = freq_axis_hz[peak_idx[1]]
peak_mag_fdcorr = np.abs(fd_surface[peak_idx])

annotation = (
    f"Corr: {peak_mag_fdcorr:.1f}\n"
    f"Phase: {peak_delay_chips:.1f} chips\n"
    f"Doppler: {peak_freq_hz:.1f} Hz"
)
ax.text(
    peak_freq_hz / 1000, peak_delay_chips, 1.05,
    annotation,
    fontsize=8, fontfamily='monospace',
    ha='center', va='bottom',
    bbox=dict(boxstyle='round,pad=0.4', facecolor='white', edgecolor='black', alpha=0.85),
    zorder=10,
)
ax.plot(
    [peak_freq_hz / 1000, peak_freq_hz / 1000],
    [peak_delay_chips, peak_delay_chips],
    [1.0, 1.0],
    "s",
    color='black', zorder=9,
)

fig.colorbar(surf, shrink=0.55, aspect=12, pad=0.1, label='Normalized Correlation')
plt.tight_layout()
plt.show()

png

NOTE: Doppler resolution is limited to 1000 Hz in this method (versus previous) and the FFT’ed signal is conjugated (versus the PRN code in previous) which leads to opposite, time-reversed Code Phase output (e.x. code length of 1023−788.5=2341023 - 788.5 = 234).

fdcorr_loss = 10.0 * np.log10(peak_mag / peak_mag_fdcorr)
print(f"Processing loss of fdcorr() vs traditional approach due to larger doppler bin: {fdcorr_loss:.2f} dB")

Processing loss of fdcorr() vs traditional approach due to larger doppler bin: 0.62 dB

Further improvements can be made to the circular correlation approach below:

  • Reuse double-buffered approach from optimal SW FIR filtering to not have to actually circularly shift FFT samples, or create Hankel matrix, but just index/offset the already-FFT’ed samples directly!
  • Parallelize each correlation (multi-threading with crates like rayon given completely parallel)
  • Show possible binary correlation method: (quantize r_bin = 1 if r > 0 else 0 -> PN code as 0/1 bin -> XOR to detect sign disagreement -> popcount() number of mismatches -> correlation = N - 2 * popcount() (N = total samples)). Note this still happens after carrier wiping, with ~1.96 dB (2/π2/\pi) theoretical loss Van Vleck for 1-bit quantization
    • Though since this is to improve time-domain correlation (multiply-accumulate), and frequency domain method is simple point-wise multiplication, does this even matter/help?
  • Handing off frequency and phase information to tracking loops
sv_idx = 31
sv_fft = prn_fft_sv[sv_idx] # conj(FFT(PRN[sv]))
assert len(sv_fft) == num_corr

sig_fft = np.fft.fft(iq_in[:num_corr])
assert len(sig_fft) == num_corr

# Similar write FFT(signal) to both halfs of a 2x larger buffer
sig_fft_dbl = np.concatenate((sig_fft, sig_fft))
assert len(sig_fft_dbl) == 2 * num_corr

# Limit Doppler search to +/- 5kHz
max_doppler_hz = 5000

doppler_bin_size = fs / num_corr
print(f"Doppler bin size = {doppler_bin_size:.2f} Hz")

# Get whole indices used for frequency shift (NOTE: floor() + 1 useful for symmetrical indices)
num_sym_bins = np.floor(max_doppler_hz / doppler_bin_size)
freq_idx = np.arange(-num_sym_bins, num_sym_bins + 1, dtype=int)
print(f"Frequency indices (shift in Freq domain): {freq_idx}")

# Instead of going through forming whole Hankel matrix, reassigning values, etc.
max_corr = 0.0
freq_offset = 0.0
code_offset = 0
for freq in freq_idx:
    if freq >= 0:
        x_corr_f = sig_fft_dbl[freq:num_corr + freq] * sv_fft
    else: # negative index value, use second half
        x_corr_f = sig_fft_dbl[freq + num_corr:freq + (2*num_corr)] * sv_fft
    x_corr_t = np.abs(np.fft.ifft(x_corr_f))

    freq_bin = freq*doppler_bin_size/1e3
    plt.plot(x_corr_t, linewidth=0.5, alpha=0.8, label=f'{freq_bin} kHz')

    highest_corr = np.max(x_corr_t)
    if highest_corr > max_corr:
        max_corr = highest_corr
        freq_offset = freq_bin
        code_offset = np.argmax(x_corr_t)

print(f"Max corr = {max_corr:.2f} [{freq_offset:.2f} kHz offset | {int(code_offset/2)} code phase]")
plt.plot(code_offset, max_corr, "s", color='black')
plt.legend(loc='upper right', title="Doppler")
plt.title("Double-buffered Frequency Domain Cross-Correlation")
plt.xlabel("Code Phase (chips)")
plt.ylabel("Correlation Magnitude")
plt.tight_layout()
plt.show()

Doppler bin size = 1000.00 Hz Frequency indices (shift in Freq domain): [-5 -4 -3 -2 -1 0 1 2 3 4 5] Max corr = 8.17 [2.00 kHz offset | 234 code phase]

png

References

Other Methods of PNT

Blind Signal Receiver for PNT

Known Terrestrial Signals

Visual



Maintained by John Gentile • Found a problem on this site? File an Issue
This work is licensed under a Creative Commons Attribution 4.0 International License.