Usage Guide

Basic Acquisition

The simplest way to acquire GNSS signals is with the acquire function:

using Acquisition, GNSSSignals
import Unitful: Hz

system = GPSL1()

# Generate a synthetic GPS L1 signal for PRN 1 (1 ms, 4 MHz)
(; signal, sampling_freq, interm_freq) = generate_test_signal(system, 1;
    num_samples = 4096, sampling_freq = 4e6Hz)

Acquire multiple PRNs at once:

results = acquire(system, signal, sampling_freq, 1:3; interm_freq)

┌─────┬────────────┬──────────────────────┬────────────────────┐
│ PRN │ CN0 (dBHz) │ Carrier Doppler (Hz) │ Code phase (chips) │
├─────┼────────────┼──────────────────────┼────────────────────┤
│ 1.0 │    44.0886 │            1000.0 Hz │            110.484 │
│ 2.0 │    39.8347 │            2000.0 Hz │            727.609 │
│ 3.0 │    40.6004 │            6000.0 Hz │            125.829 │
└─────┴────────────┴──────────────────────┴────────────────────┘

Acquire a single PRN:

result = acquire(system, signal, sampling_freq, 1; interm_freq)

AcquisitionResults: PRN 1, CN0 = 44.09 dB-Hz, Doppler = 1000.0 Hz, Code phase = 110.484 chips

Each result contains:

• carrier_doppler: Estimated Doppler frequency
• code_phase: Code phase in chips
• CN0: Carrier-to-noise density ratio (dB-Hz)
• peak_to_noise_ratio: Peak correlation power divided by noise power
• power_bins: Correlation power matrix (only populated when store_power_bins=true)

Detecting Satellites

Use is_detected to decide whether a satellite signal is present:

# Filter to detected satellites (1% false alarm probability)
detected = filter(is_detected, results)

┌─────┬────────────┬──────────────────────┬────────────────────┐
│ PRN │ CN0 (dBHz) │ Carrier Doppler (Hz) │ Code phase (chips) │
├─────┼────────────┼──────────────────────┼────────────────────┤
│ 1.0 │    44.0886 │            1000.0 Hz │            110.484 │
└─────┴────────────┴──────────────────────┴────────────────────┘

# Custom false alarm probability
detected = filter(r -> is_detected(r; pfa = 0.001), results)

┌─────┬────────────┬──────────────────────┬────────────────────┐
│ PRN │ CN0 (dBHz) │ Carrier Doppler (Hz) │ Code phase (chips) │
├─────┼────────────┼──────────────────────┼────────────────────┤
│ 1.0 │    44.0886 │            1000.0 Hz │            110.484 │
└─────┴────────────┴──────────────────────┴────────────────────┘

Under the hood, this compares each result's peak_to_noise_ratio against a CFAR (Constant False Alarm Rate) threshold computed by cfar_threshold.

CFAR Threshold — How It Works

The detector uses the statistic

peak_to_noise_ratio = peak_power / noise_power

Under the null hypothesis (noise only), the peak across num_cells = num_doppler_bins × num_code_phases search cells follows a scaled chi-squared distribution with 2M degrees of freedom, where M is the number of non-coherent accumulations stored in the result. cfar_threshold returns the quantile of that distribution at the requested per-grid false alarm probability, using a Bonferroni-like correction across num_cells:

threshold = cfar_threshold(0.01, get_num_cells(result); num_noncoherent_integrations = 1)

15.666787667699174

is_detected(result; pfa) takes care of both arguments for you — it reads num_cells and num_noncoherent_integrations directly from the result — so you only need cfar_threshold when you want to inspect the threshold separately (e.g. for plotting or custom logic).

Choosing pfa:

• pfa = 0.01 (the default) is a reasonable starting point for a cold acquire.
• Lower values (1e-4, 1e-6) are typical when a false track is expensive — the threshold rises modestly because the chi-squared tail is steep.
• Only the ratio matters: the threshold is independent of the absolute noise power, so it works unchanged at any CN0 and any sampling frequency.

Using Acquisition Plans

For repeated acquisitions — e.g. tracking many epochs or processing a file — pre-compute a plan once to avoid repeated FFT planning and memory allocation:

plan = plan_acquire(system, sampling_freq, collect(1:3))

results = acquire!(plan, signal, 1:3; interm_freq)

┌─────┬────────────┬──────────────────────┬────────────────────┐
│ PRN │ CN0 (dBHz) │ Carrier Doppler (Hz) │ Code phase (chips) │
├─────┼────────────┼──────────────────────┼────────────────────┤
│ 1.0 │    44.0886 │            1000.0 Hz │            110.484 │
│ 2.0 │    39.8347 │            2000.0 Hz │            727.609 │
│ 3.0 │    40.6004 │            6000.0 Hz │            125.829 │
└─────┴────────────┴──────────────────────┴────────────────────┘

Multi-threaded Acquisition

When processing many PRNs, start Julia with multiple threads and the PRN loop runs in parallel automatically — no code changes required:

julia -t 4

The plan allocates per-thread scratch buffers at construction time based on the number of threads available when plan_acquire is called. If you create the plan and later run with more threads, restart Julia with the desired thread count before calling plan_acquire.

Non-coherent Integration

At low CN0, accumulate power across multiple successive signal segments:

plan_ni = plan_acquire(system, sampling_freq, [1];
    num_coherently_integrated_code_periods = 10,
    num_noncoherent_accumulations = 8)

long_signal = generate_test_signal(system, 1;
    num_samples = 8 * 10 * 4096,
    sampling_freq = sampling_freq, CN0 = 30).signal

result_ni = acquire!(plan_ni, long_signal, [1])

The signal must contain at least num_noncoherent_accumulations × num_coherently_integrated_code_periods × samples_per_code samples.

Plotting Results

Acquisition results can be plotted directly with Plots.jl. Pass store_power_bins = true to retain the correlation power surface.

The three examples below use the same PRN 1 signal at default CN0 (45 dB-Hz) and show how coherent and non-coherent integration affect the correlation surface. A single 2 ms signal is generated once; the 1 ms baseline uses only the first half.

1 ms coherent integration (baseline)

1 ms gives ~1000 Hz Doppler bin spacing — coarse but fast:

using Plots
plotlyjs()

fs = 4e6Hz
prn = 1
doppler = 1500Hz

# Generate one 2 ms signal (8000 samples at 4 MHz) — reused across all three plots
signal_2ms = generate_test_signal(system, prn;
    num_samples = 8000, sampling_freq = fs, doppler = doppler).signal

plan_1ms = plan_acquire(system, fs, [prn])
result_1ms = acquire!(plan_1ms, signal_2ms[1:4000], [prn]; store_power_bins = true)
plot(result_1ms[1])

2 ms coherent integration

2× longer integration → 2× finer Doppler bins (~500 Hz spacing). The correlation spike narrows visibly in the Doppler dimension:

plan_2ms = plan_acquire(system, fs, [prn];
    num_coherently_integrated_code_periods = 2)
result_2ms = acquire!(plan_2ms, signal_2ms, [prn]; store_power_bins = true)
plot(result_2ms[1])

1 ms coherent + 2 non-coherent accumulations

Non-coherent integration adds power from 2 successive 1 ms segments, improving sensitivity without requiring longer phase-coherent integration. The Doppler resolution stays at ~1000 Hz but the peak-to-noise ratio improves:

plan_ni = plan_acquire(system, fs, [prn]; num_noncoherent_accumulations = 2)
result_ni = acquire!(plan_ni, signal_2ms, [prn]; store_power_bins = true)
plot(result_ni[1])

Algorithm Constraints and Trade-offs

The FM-DBZP algorithm (Heckler & Garrison 2009) has different constraints from a classical serial Doppler search. Understanding them is important for choosing acquisition parameters.

Doppler Resolution and Coverage

The coherent integration time T_coh determines both the Doppler resolution and the Doppler bin spacing:

T_coh = num_coherently_integrated_code_periods × samples_per_code / sampling_freq

Doppler bin spacing = 1 / T_coh

For GPS L1 C/A at 4 MHz (samples_per_code = 4092):

Unlike a classical search where you can choose any Doppler step independently of integration length, FM-DBZP fixes the Doppler step at 1 / T_coh. To get finer Doppler resolution you must integrate longer.

The total Doppler coverage is:

Doppler coverage = num_blocks / T_coh = num_blocks × Doppler bin spacing

where num_blocks is chosen automatically to cover at least min_doppler_coverage on each side (default ±7000 Hz). You can widen the search with:

plan = plan_acquire(system, sampling_freq, prns; min_doppler_coverage = 10_000Hz)

What min_doppler_coverage actually guarantees

The Doppler grid stored in plan.doppler_freqs is

range(-coverage/2, step = bin_spacing, length = num_doppler_bins)

where coverage = num_blocks × (sampling_freq / samples_per_code) and num_doppler_bins = num_coherently_integrated_code_periods × num_blocks. This is a half-open interval — the grid spans [-coverage/2, +coverage/2), so the highest searched bin is +coverage/2 - bin_spacing, not +coverage/2.

min_doppler_coverage is the minimum guaranteed reach on both ends: plan_acquire chooses num_blocks such that

last(plan.doppler_freqs)  ≥ +min_doppler_coverage
first(plan.doppler_freqs) ≤ -min_doppler_coverage

Concretely, with the default min_doppler_coverage = 7000Hz:

The asymmetry is a consequence of the FFT bin layout: a length-N DFT covers exactly N bins worth of bandwidth, and centering those N bins on 0 leaves the upper edge open. It is not a bug — the highest searched Doppler is the last bin, and that bin is guaranteed to be ≥ +min_doppler_coverage.

The num_blocks Divisibility Constraint

num_blocks must divide samples_per_code exactly so that each block has an integer number of samples (block_size = samples_per_code ÷ num_blocks). plan_acquire finds the smallest valid divisor automatically, but this means not all sampling frequencies support all Doppler coverages.

For example, at 2.048 MHz (samples_per_code = 2048) the valid divisors are 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 — so num_blocks will always be a power of two and Doppler coverage is always a multiple of the bin spacing. At 5 MHz (samples_per_code = 5000) the divisors include 5, 10, 20, 25, … — giving more choices but potentially a larger jump to the next valid num_blocks.

If plan_acquire cannot find a valid num_blocks for your sampling frequency and min_doppler_coverage, it throws an ArgumentError. Try a slightly different sampling frequency or reduce min_doppler_coverage.

Why the smallest valid divisor?

Once min_doppler_coverage fixes a lower bound, several divisors of samples_per_code are usually admissible (e.g. at 5 MHz with ±9 kHz requested, both num_blocks = 20 and num_blocks = 25 are valid). It is tempting to pick a larger divisor: that shrinks block_size and thus the inner double-block FFT of size 2 × block_size. In practice the opposite is the right default — plan_acquire picks the smallest valid divisor, which corresponds to the narrowest Doppler coverage that still meets your request.

The reason is that acquisition has two FFT stages with opposite scaling in num_blocks:

• Inner double-block stage. num_blocks forward + inverse FFTs of size 2 × block_size = 2 × samples_per_code / num_blocks. Total work scales as 2 × samples_per_code × log₂(2 × samples_per_code / num_blocks) — it shrinks logarithmically as num_blocks grows.
• Column FFT stage. samples_per_code FFTs of size num_doppler_bins = num_coherently_integrated_code_periods × num_blocks. Total work scales as samples_per_code × num_doppler_bins × log₂(num_doppler_bins) — it grows roughly linearly (× log) in num_blocks. The non-coherent accumulation matrix is num_doppler_bins × samples_per_code, so its memory traffic scales with num_blocks as well.

A concrete comparison at 4 MHz, num_coherently_integrated_code_periods = 1:

The inner stage barely changes; the column stage grows ~35 %, and overall cost rises by ~26 %. The gap widens with longer coherent integration, because num_doppler_bins then equals num_coherently_integrated_code_periods × num_blocks, so every extra divisor unit costs num_coherently_integrated_code_periods extra column-FFT bins.

There is one second-order effect that can invert this trade-off: the inner FFT size is 2 × samples_per_code / num_blocks, so a smaller divisor can land on a block_size with a large prime factor (11, 31, 257, …) and a slow FFTW kernel (see Sampling Frequency and FFT Performance). In that situation a slightly larger divisor with a smoother block_size may win — but this is sampling-frequency-specific and rare in practice. If you suspect you have hit such a case, benchmark both min_doppler_coverage settings explicitly.

Sampling Frequency and FFT Performance

Acquisition runs many in-place FFTs at size 2 × block_size, where block_size = samples_per_code ÷ num_blocks. FFTW is very fast when those sizes factor into small primes (2, 3, 5, 7) and noticeably slower when they contain a large prime factor such as 11, 13, 31, or 257. Since samples_per_code itself is ceil(code_length × fs / code_freq), its prime factorization is decided entirely by your sampling frequency — a 0.1 % change in fs can make the inner FFT 3-5× faster.

As a rule of thumb for GPS L1 C/A:

• Prefer sampling frequencies where samples_per_code (≈ fs / 1000 in Hz, rounded up) has only small prime factors.
• Powers of two (2.048, 4.096, 8.192, 16.384 MHz) and 2^a · 5^b rates (2, 2.5, 5, 10.24 MHz) are always fast.
• Avoid rates whose samples_per_code contains a prime ≥ 11. The commonly recommended 16.368 MHz (= 16 × 1.023 MHz) is a notable offender: 16368 = 2^4 · 3 · 11 · 31, and the radix-31 inner FFT is ~5× slower than nearby smooth sizes.

Measured per-PRN acquire! times on a Ryzen 7 PRO 5850U (16 threads), GPS L1 C/A, min_doppler_coverage = 10 000 Hz:

Notice the near-identical-rate pairs: 8.192 MHz runs at 50 ms vs. 8.184 MHz at 82 ms (1.6× faster); 16.384 MHz runs at 112 ms vs. 16.368 MHz at 167 ms (1.5× faster). If you control the RF front-end clock, picking a smooth rate is usually the single biggest performance lever available.

Picking a Good Sampling Frequency

When you have flexibility in the front-end clock, recommend_sampling_freqs sweeps a frequency range and returns the candidates with the smoothest FFT factorizations and lowest estimated cost. It accounts for the num_blocks divisibility constraint and both FFT stages of the algorithm (inner double-block + column FFT).

using Acquisition, GNSSSignals
import Unitful: Hz

recommend_sampling_freqs(GPSL1();
    fs_min = 2e6Hz,
    fs_max = 5e6Hz,
    min_doppler_coverage = 7000Hz,
    num_coherently_integrated_code_periods = 1,
    num_alternatives = 5,
)

┌─────────────────────┬──────────────────┬────────────┬────────────┬────────────
│ Sampling freq (MHz) │ Samples per code │ Num blocks │ Block size │ Inner FFT ⋯
├─────────────────────┼──────────────────┼────────────┼────────────┼────────────
│                 2.0 │             2000 │         16 │        125 │           ⋯
│               2.016 │             2016 │         16 │        126 │           ⋯
│               2.048 │             2048 │         16 │        128 │           ⋯
│                2.16 │             2160 │         16 │        135 │           ⋯
│                2.24 │             2240 │         16 │        140 │           ⋯
└─────────────────────┴──────────────────┴────────────┴────────────┴────────────
                                                               4 columns omitted

For a custom code geometry (e.g. a non-standard system) pass code_length and code_freq directly:

recommend_sampling_freqs(10_000, 36e6Hz;
    fs_min = 36e6Hz,
    fs_max = 40e6Hz,
    min_doppler_coverage = 10_000Hz,
    num_coherently_integrated_code_periods = 36,
    num_alternatives = 5,
)

┌─────────────────────┬──────────────────┬────────────┬────────────┬────────────
│ Sampling freq (MHz) │ Samples per code │ Num blocks │ Block size │ Inner FFT ⋯
├─────────────────────┼──────────────────┼────────────┼────────────┼────────────
│              36.285 │            10080 │          6 │       1680 │           ⋯
│              36.739 │            10206 │          6 │       1701 │           ⋯
│              37.041 │            10290 │          6 │       1715 │           ⋯
│              37.322 │            10368 │          6 │       1728 │           ⋯
│              37.797 │            10500 │          6 │       1750 │           ⋯
└─────────────────────┴──────────────────┴────────────┴────────────┴────────────
                                                               4 columns omitted

By default candidates are ranked by estimated FFT cost; pass sort_by = :smoothness to rank by largest prime factor of the inner FFT size instead. Use max_prime to tighten or relax the smoothness budget (default 7, FFTW's "fast" regime).

Filtering against an SDR's hardware constraints

Most SDRs can only reach a discrete subset of sampling frequencies — the output of a master clock, divider tree, and PLL. Pass an AbstractSDRClockPlan via sdr_clock_plan to skip rates the hardware can't produce:

recommend_sampling_freqs(GPSL1();
    fs_min = 2e6Hz,
    fs_max = 8e6Hz,
    min_doppler_coverage = 7000Hz,
    num_alternatives = 5,
    sdr_clock_plan = AD9361ClockPlan(),
)

┌─────────────────────┬──────────────────┬────────────┬────────────┬────────────
│ Sampling freq (MHz) │ Samples per code │ Num blocks │ Block size │ Inner FFT ⋯
├─────────────────────┼──────────────────┼────────────┼────────────┼────────────
│                2.16 │             2160 │         16 │        135 │           ⋯
│                2.24 │             2240 │         16 │        140 │           ⋯
│               2.304 │             2304 │         16 │        144 │           ⋯
│               2.352 │             2352 │         16 │        147 │           ⋯
│                 2.4 │             2400 │         16 │        150 │           ⋯
└─────────────────────┴──────────────────┴────────────┴────────────┴────────────
                                                               4 columns omitted

AD9361ClockPlan reproduces the validation done by the AD9361 driver shipped with litex_m2sdr (550 kHz – 61.44 MSPS, divider search through {12, 8, 6, 4, 3, 2, 1} against an ADC clock window of 25–640 MHz, and a reachable BBPLL). Other SDRs can be supported by subtyping AbstractSDRClockPlan and implementing is_valid_sample_rate (and optionally sample_rate_range).

Coherent Integration with Data Bits (GPS L1 C/A)

GPS L1 C/A data bits flip every 20 ms. Integrating across a bit transition cancels signal energy. The library handles this automatically:

Sub-bit integration (num_coherently_integrated_code_periods < 20): At most one bit transition can fall in the window. Use bit_edge_search_steps > 1 to search over candidate alignment positions:

plan = plan_acquire(system, sampling_freq, prns;
    num_coherently_integrated_code_periods = 10,
    bit_edge_search_steps = 10)

Multi-bit integration (num_coherently_integrated_code_periods ≥ 20): The window spans whole data bit periods. Two constraints apply:

1. num_coherently_integrated_code_periods must be divisible by bit_period_codes (= 20 for GPS L1)
2. bit_period_codes must be divisible by bit_edge_search_steps

The algorithm searches all 2^(num_data_bits - 1) sign-flip combinations across the data bits in the window.

Pilot channels (e.g. GPS L5-Q, Galileo E1-C): no data bit constraint; set num_coherently_integrated_code_periods freely.

Sub-sample Interpolation

By default, code phase and Doppler estimates are quantised to the search grid:

• Code-phase step = 1 / sampling_freq in seconds (0.25 µs at 4 MHz ≈ 0.25 chips for GPS L1 C/A).
• Doppler step = 1 / T_coh (see the Doppler Resolution and Coverage table).

Pass subsample_interpolation = true to refine both below the grid spacing using a parabolic fit across the peak bin and its two neighbours:

result_interp = acquire(system, signal, sampling_freq, 1;
    interm_freq, subsample_interpolation = true)

AcquisitionResults: PRN 1, CN0 = 44.09 dB-Hz, Doppler = 1034.576952457428 Hz, Code phase = 110.608 chips

The fit is only applied when the neighbouring bins exceed √noise_power (so it doesn't chase noise on non-detections), and it costs four extra array reads per PRN — negligible compared to the acquisition itself.

Storing the Power Surface

Pass store_power_bins = true to keep the full Doppler × code-phase correlation matrix in the returned result. It's required for plotting (see Plotting Results) and useful for post-hoc analysis; without it the power_bins field is nothing and the plan's internal buffer is reused across calls to keep acquisitions allocation-free.

result = acquire!(plan, signal, [1]; store_power_bins = true)