diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 00000000..deb9cdc6 --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,62 @@ +# Agent instructions + +## Commit messages: Conventional Commits (drive the releases) + +Releases are fully automated: every push to `master` runs semantic-release +(`.github/workflows/Release.yml`), which derives the version bump, the +changelog, and the JuliaRegistrator call from the commit messages. Rebase +merge is enforced, so **every commit in a PR lands verbatim on `master`** and +must follow [Conventional Commits](https://www.conventionalcommits.org): + +``` +type(scope): short imperative summary + +optional body explaining what and why +``` + +- `feat` → minor release, `fix` / `perf` → patch release; `refactor`, `test`, + `docs`, `build`, `chore` → no release. +- Pick the scope from the area touched (e.g. `track`, `dc`, `onebit`, + `bench`); look at `git log --oneline` for precedent. + +### Breaking changes must be marked and explained + +A breaking change **must** carry both: + +1. a `!` after the type/scope — e.g. `feat(track)!: …`, and +2. a `BREAKING CHANGE:` footer paragraph that spells out *what* breaks and + *how to migrate* (removed/renamed exported types, functions, fields, or + kwargs; changed constructor signatures; changed observable behavior of + public API). + +The footer is what triggers the major version bump and becomes the release +notes — an unmarked breaking change silently ships as a minor/patch release. +When in doubt: it is breaking if released user code — including custom +`AbstractDopplerEstimator` implementations and direct +`downconvert_and_correlate(!)` callers — could stop compiling or change +behavior after the release. + +Breaking is always judged **against the latest release**, not against earlier +commits of the same PR: renaming or removing API that was itself introduced in +the same (unmerged) PR is an ordinary `refactor`/`fix` commit, not a breaking +change — either fold it into the introducing commit or say in the body that +the touched API is unreleased. Reserve `!` + `BREAKING CHANGE:` for changes +that affect users of a published version. + +### Do not hand-maintain release artifacts + +The changelog and the version are produced by CI: on `master`, +semantic-release generates `CHANGELOG.md` entries and the release notes from +the commit messages and bumps `Project.toml`'s `version` accordingly. Never +edit `CHANGELOG.md` or bump `version` manually in a PR — write the commit +message well instead, since that text *is* the changelog. + +## Formatting + +CI enforces JuliaFormatter with the exact version pinned in +`.github/workflows/format.yml` (currently 2.8.5) against +`.JuliaFormatter.toml`. Before committing: + +```julia +using JuliaFormatter; format(["src", "test"]) +``` diff --git a/docs/plans/2026-07-21-chunked-doppler-update.md b/docs/plans/2026-07-21-chunked-doppler-update.md new file mode 100644 index 00000000..4ca2eff2 --- /dev/null +++ b/docs/plans/2026-07-21-chunked-doppler-update.md @@ -0,0 +1,151 @@ +# Chunked Doppler-update tracking loop + +**Date:** 2026-07-21 +**Status:** Implemented + +## Goal + +Decouple the Doppler-estimation (NCO-update) rate from the code period. Process +each measurement in fixed-size **time chunks**; within a chunk the NCO Doppler +is held fixed and every correlator output that completes is collected — tagged +with the sample index at which it was taken — into a per-signal array. After the +chunk, the estimator folds over those outputs in order and updates every +satellite's NCO once, at a common epoch. + +## Motivation + +Previously `track!` advanced each satellite to its *next code-block boundary* +and immediately re-estimated and rewrote that satellite's NCO Doppler. Three +limitations: + +- The update rate was welded to the code period — no way to trade + Doppler-estimation compute for a slower update rate. +- Correlator outputs were consumed one at a time and never collected; the + estimator could only ever see the single just-completed correlator. Vector + tracking needs a *batch* of outputs across channels, each tagged with its + sample epoch, on a common time grid. +- NCO updates were staggered — each satellite updated its NCO the instant *it* + finished a code period, so different satellites updated at different times. + +## Design + +### Chunk grid + +The chunk length is the Doppler-update interval, a **time** (`doppler_update_interval` +kwarg on [`track`](@ref) / [`track!`](@ref)). Default `nothing` ⇒ auto = the +smallest primary-code period across all tracked signals (`_smallest_code_period`; +e.g. 1 ms for GPS L1 C/A, or for a mixed L1 C/A + Galileo E1B state). For band +`b` at sampling frequency `fs_b`, chunk `k` (0-based) ends at sample +`min(round((k+1)·Δt·fs_b), num_samples_b)` (`_chunk_last_sample`). The boundary +is re-anchored to the absolute chunk index each call (not accumulated), so +rounding never drifts and different bands stay time-aligned to within a sample. +A one-sample floor is validated up front (`_validate_doppler_update_interval`) so the +loop always makes progress. + +### Correlate phase (records outputs, two passes per chunk) + +`_update_tracked_sat_correlator` runs an inner loop over the chunk. Each +sub-step integrates to the next code-block boundary of any signal (or the chunk +end), reusing `_calc_min_samples_and_completed` with the chunk end substituted +for the buffer end (the true buffer length is still used for replica sizing). +When a signal's integration completes, `update` / `_build_new_signals` snapshots +its raw accumulator into a [`CorrelatorOutput`](@ref) — `(correlator, +integrated_samples, sample_index)` — appended to the reused +per-signal `correlator_outputs` buffer, and resets the accumulator so the next +integration in the chunk starts fresh. A partial integration at the chunk (or +buffer) boundary carries in the accumulator. The NCO Doppler is untouched here. +A chunk yields 0, 1, or several outputs per signal depending on code phase and +Doppler. + +`track!` runs this correlate step **once per chunk** with +`stop_before_partial = true`: integrate every completion inside the chunk, but +stop at the last code-block boundary instead of integrating the trailing +chunk-clamped partial (a sub-step that completes no signal is exactly that +partial). The estimator then writes the new NCO Doppler, and the **next** +chunk's pass picks up from the boundary — so each integration runs boundary → +boundary in a single kernel window, entirely at the just-updated Doppler. +After the last chunk a final pass without the flag drains the buffer's +trailing partial into each satellite's live accumulator (at the final chunk's +Doppler) so it carries into the next `track!` call, followed by one trailing +fold for the rare completion landing exactly on the buffer end. + +This keeps every completed integration on a single NCO Doppler and applies each +Doppler correction right at the boundary where its integration completed — the +same loop timing as the classic per-completion update — while still estimating +once per chunk at a common epoch. (Integrating the residue at the pre-update +Doppler instead measurably tightened the FLL pull-in edge; an interim design +that integrated the residue in a separate per-chunk pass restored the edge but +split each code period into two kernel invocations, costing ~20 % throughput +on long buffers.) The outer loop walks the chunk grid itself (`_chunks_left`), +not per-sat progress: satellites lagging behind a chunk boundary — e.g. an +E1B sat whose 4 ms integration spans several 1 ms chunks — are caught up by +whichever later pass contains their boundary, and finally by the drain. + +### Estimate phase (folds, one NCO write) + +`_apply_correlator_output` applies one record: normalize by its sample count, +post-corr filter, push the filtered prompt, advance CN0 and bit buffer, move the +record to `last_fully_integrated_*` (it no longer resets the accumulator — the +correlate phase already did). The driver processor +`_process_estimator_driver_signal` folds over `correlator_outputs` in order, +threading the loop-filter state and the FLL `previous_prompt` across records +(per-record `integration_time` and `1/N` bandwidth scaling), computing a Doppler +each record and keeping the **last**; the NCO is written once. Passenger signals +apply the per-record advance only. Each signal's `correlator_outputs` is +`empty!`-ed after folding. With no outputs in a chunk the Doppler holds. + +### Phase-snap at sync + +Sync is detected during the estimate pass, *after* the whole chunk was +correlated. The one-time secondary-code phase snap +(`_snap_code_phase_from_synced_signal`) still runs against the end-of-chunk +`code_phase`, preserving the within-primary-block phase (issue #117). Because any +in-flight partial for the sync chunk was accumulated at the *pre-snap* phase (and +pre-sync without a secondary overlay), the snap now also resets every signal's +in-flight accumulator (`_reset_inflight_integration`): the phase bookkeeping is +kept, but the next chunk re-integrates cleanly from the snapped phase. Otherwise +a flipped overlay chip can cancel the partial to zero and feed a 0/0 into the +discriminators. (With the two-pass chunk the estimate usually runs exactly on a +boundary, so the reset is a no-op there; it still matters for signals whose +integration spans multiple chunks and thus carry a partial at estimate time.) + +Records that *follow* the syncing record within the same fold were correlated +with pre-sync replicas (no secondary-code wipe-off), so accumulating them into +the bit buffer as if wiped would feed sign-corrupted prompts into the first +post-sync bits. The folds detect the mid-fold `found` transition and apply +those records with `skip_bit_buffer = true` (`_apply_correlator_output`): +prompt filter, CN0, and loop-filter processing still run, only the bit-buffer +update is skipped. Irrelevant at the default interval (a chunk then holds at +most one record past the sync), it protects enlarged `doppler_update_interval`s. + +### Data model + +- [`CorrelatorOutput`](@ref)`{C}` in `src/correlators/correlator.jl`. +- [`TrackedSignal`](@ref) gains `correlator_outputs::Vector{CorrelatorOutput{C}}`, + preallocated and reused like `filtered_prompts` (`sizehint!`-ed, `empty!`-ed at + each `track` start and after each chunk's fold), so steady-state tracking stays + allocation-free. + +## Backends + +All new logic lives in the shared `_update_tracked_sat_correlator` / `update` / +estimator layer, above the backend boundary `_correlate_signals`. The CPU-fused +and Int16 backends (single + threaded) need no changes; the per-`_correlate_signals` +call count is unchanged from before. + +## Behavior compatibility + +With the two-pass chunk, every completed integration is produced by a single +NCO Doppler and each correction takes effect at its completing boundary — the +same loop timing as the previous per-completion update. The suite's edge +probes match the pre-chunking behavior exactly: the FLL pull-in edge sits at +its historical 240 Hz (converges) / 250 Hz (fails), and the 1800-block L1C-P +overlay sync locks after exactly one overlay cycle. (An earlier single-pass +design integrated the chunk residue at the pre-update Doppler, which tightened +the FLL edge to ~210 Hz and could land the overlay sync one block late.) + +## Deferred follow-ups + +- A global (cross-`track!`-call) sample index for streamed real-time chunks + (currently `sample_index` is relative to the current measurement buffer). +- The vector-tracking consumer of the batched, sample-indexed correlator outputs. diff --git a/docs/src/correlator.md b/docs/src/correlator.md index e76ffd86..241043a2 100644 --- a/docs/src/correlator.md +++ b/docs/src/correlator.md @@ -119,6 +119,17 @@ get_num_accumulators get_num_ants(::AbstractCorrelator) ``` +## Correlator outputs + +Each completed integration within a processing chunk is recorded as a +`CorrelatorOutput`, collected per signal and consumed by the Doppler estimator +after the chunk (see [Chunked Doppler updates](track.md#Chunked-Doppler-updates)). + +```@docs +CorrelatorOutput +get_correlator_outputs +``` + ## Sample shifts The downconvert/correlate inner loop walks the incoming sample buffer once diff --git a/docs/src/custom_doppler_estimator.md b/docs/src/custom_doppler_estimator.md index 13f30977..c2b01386 100644 --- a/docs/src/custom_doppler_estimator.md +++ b/docs/src/custom_doppler_estimator.md @@ -78,12 +78,23 @@ per-sat fields directly and rewraps `doppler_estimator_state` unchanged. 5. **An `estimate_dopplers_and_filter_prompt` method** dispatched on `TrackState{<:Any, <:MyEstimator}`. This is where the actual update - logic runs, once per integration completion. It walks each group in - `track_state.groups`, reads the band's `BandMeasurement` from the - `measurements::BandMeasurements` NamedTuple via `get_band_id(group.band)`, - and produces new `TrackedSat`s with updated + logic runs, once per **chunk** (see [Chunked Doppler updates](track.md#Chunked-Doppler-updates)). + It walks each group in `track_state.groups`, reads the band's + `BandMeasurement` from the `measurements::BandMeasurements` NamedTuple via + `get_band_id(group.band)`, and produces new `TrackedSat`s with updated `carrier_doppler`/`code_doppler` and updated per-sat estimator state. + Each signal's correlator outputs completed during the chunk are in its + `correlator_outputs::Vector{`[`CorrelatorOutput`](@ref)`}` (a chunk may hold + zero, one, or several per signal), each carrying the raw correlator, its + integrated-sample count, and the end sample index. + Fold over them in order — threading whatever filter state you carry — and + write the NCO Doppler once (typically from the last output). Empty each + signal's `correlator_outputs` when done, and remember the NCO Doppler is held + fixed across the whole chunk, so a discriminator that needs the replica's + Doppler should read the satellite's `code_doppler`/`carrier_doppler` (the + value that generated the chunk), not an intermediate per-output estimate. + The matching mutating method `estimate_dopplers_and_filter_prompt!(track_state, measurements, prefer)` is what [`track!`](@ref) calls. To support real-time loops, define diff --git a/docs/src/track.md b/docs/src/track.md index dc70a7a0..81d34953 100644 --- a/docs/src/track.md +++ b/docs/src/track.md @@ -11,12 +11,27 @@ track! - `downconvert_and_correlator` — the downconversion and correlation implementation. Defaults to `CPUThreadedDownconvertAndCorrelator()`. **For real-time loops, hoist this outside the loop** (see below). - `intermediate_frequency` — the IF of the signal. Defaults to `0.0Hz`. Only accepted on the bare-buffer form `track!(buf, state, fs; intermediate_frequency = ...)`; on the [`BandMeasurement`](@ref) and multi-band forms the IF lives on each `BandMeasurement`. +- `doppler_update_interval` — the Doppler-estimation / NCO-update interval, a time (e.g. `1u"ms"`). Defaults to `nothing` ⇒ auto = the smallest primary-code period across all tracked signals (1 ms for GPS L1 C/A). Each measurement is processed in fixed-size chunks of this length: within a chunk the NCO Doppler is held fixed and every correlator output that completes is collected, then the estimator processes them in order and updates **every** satellite's NCO once, at a common epoch (see [Chunked Doppler updates](#Chunked-Doppler-updates)). Pick a longer interval to reduce Doppler-estimation cost at the expense of update rate. + The **coherent-integration length** is not a `track!` argument — it is a per-signal setting on each [`TrackedSignal`](@ref) (its `preferred_num_code_blocks_to_integrate` field), changed with [`set_preferred_num_code_blocks_to_integrate!`](@ref). It defaults to `1`, is capped per integration by the signal's bit/secondary-code period, and only takes effect once bit/secondary-code synchronization has been achieved. For data-bearing signals the length must evenly divide the number of code blocks that form one bit (e.g. a divisor of 20 for GPS L1 C/A, of 10 for GPS L5I) so integrations stay aligned to bit boundaries; other values throw an `ArgumentError`. With the conventional estimator the loop bandwidth auto-scales by `1/N` so longer integration stays stable without re-tuning. ```@docs set_preferred_num_code_blocks_to_integrate! ``` +## Chunked Doppler updates + +`track` / `track!` walk each measurement in fixed-size time chunks of length `doppler_update_interval` (default: the smallest code period across all signals). Each chunk runs one correlate pass and one estimate: + +1. **Correlate to the last completed boundary** — each satellite integrates from wherever it stands up to its last coherent-integration boundary inside the chunk; every completed integration is collected into that signal's `correlator_outputs` buffer, tagged with the sample index at which it ended (a [`CorrelatorOutput`](@ref); the sample index is important for vector tracking). A 1 ms-code signal in a 1 ms chunk yields 0, 1, or 2 outputs; a signal whose coherent integration is longer than the chunk yields outputs only on the chunks where it completes. +2. **Estimate** — the Doppler estimator processes the collected outputs **in order** (threading the loop-filter state across them) and writes the resulting Doppler to the NCO **once per chunk** — all satellites' NCOs update at the same point in the processing, a common epoch. + +The chunk's trailing partial — from each satellite's last completed boundary to the chunk end — is *not* integrated separately: the **next** chunk's pass starts right at that boundary, so each integration runs boundary → boundary in one kernel window, entirely at the freshly updated Doppler. Every completed integration is therefore produced by a single NCO Doppler and each correction takes effect right at the boundary where its integration completed — the same loop timing as a classic per-code-period update. A final pass after the last chunk drains the buffer's trailing partial into each satellite's live accumulator so it carries into the next `track!` call. + +Read the collected outputs for the most recent chunk with [`get_correlator_outputs`](@ref) (they are cleared after each chunk's estimate). + +Choosing a larger `doppler_update_interval` batches more correlator outputs per NCO update, trading Doppler-tracking bandwidth for lower estimation cost. At the default interval the behavior is numerically very close to updating the NCO once per code period. + ## Real-time loops A typical receiver loop builds the `TrackState` once, hoists the correlator outside the loop, and calls `track!` per chunk: diff --git a/src/Tracking.jl b/src/Tracking.jl index b861cc45..b7eb083f 100644 --- a/src/Tracking.jl +++ b/src/Tracking.jl @@ -20,7 +20,7 @@ using Polyester # sync-detection-redesign plan in docs/plans for the comparison). BitIntegers.@define_integers 1800 -using Unitful: upreferred, uconvert, Hz, dBHz, ms, s +using Unitful: upreferred, uconvert, dimension, NoUnits, Hz, dBHz, ms, s import Base.zero, Base.length, Base.resize! export get_early, @@ -40,6 +40,8 @@ export get_early, get_last_fully_integrated_correlator, get_last_fully_integrated_filtered_prompt, get_filtered_prompts, + get_correlator_outputs, + CorrelatorOutput, get_bit_buffer, get_bits, get_num_bits, diff --git a/src/conventional_pll_and_dll.jl b/src/conventional_pll_and_dll.jl index 9486ac97..ee98f0e8 100644 --- a/src/conventional_pll_and_dll.jl +++ b/src/conventional_pll_and_dll.jl @@ -321,34 +321,47 @@ function _update_tracked_sat_doppler(sat::TrackedSat, sampling_frequency) # secondary-chip window. # # This is a *one-time* anchoring applied only on the iteration a - # signal transitions `found == false → true`. The snap aligns - # `code_phase` to a primary-code boundary within the secondary-code - # window, which is only valid when `code_phase` actually sits on - # such a boundary — i.e. just after a completed integration, which - # is exactly when sync is detected. Re-running it on later - # iterations would clobber the within-primary-block phase that - # `update` advances during a partial (chunk-bounded) integration, - # pinning `code_phase` to the boundary and wedging the satellite - # into a state where no chunk can ever complete the block again - # (see issue #117). After sync, `update`'s `mod(…, current_code_wrap)` - # maintains the secondary alignment on its own. + # signal transitions `found == false → true`. It preserves the + # within-primary-block phase (`mod(code_phase, primary)`) so the loop + # keeps the current chunk-bounded position; re-running it on later + # iterations would wedge the satellite (see issue #117). After sync, + # `update`'s `mod(…, current_code_wrap)` maintains the alignment. + # + # Because sync is detected in this estimate pass — *after* the whole + # chunk was correlated — any in-flight partial integration for this + # chunk was accumulated at the pre-snap code phase (and, pre-sync, with + # no secondary-code overlay). Once the snap jumps `code_phase` into the + # secondary window that partial's data is phase-inconsistent with the + # new alignment (for a flipped NH chip it can even cancel to zero and + # feed a 0/0 into the discriminators). So on the snap we also reset every + # signal's in-flight accumulator: the phase bookkeeping is kept, but the + # next chunk re-integrates cleanly from the snapped phase to the next + # boundary. Block-aligned starts have no residue, so this is a no-op there. new_signals = (new_head, new_tail...) - snapped_code_phase = if _any_signal_just_synced(sat.signals, new_signals) - _snap_code_phase_from_synced_signal(new_signals, sat.code_phase) - else + just_synced = _any_signal_just_synced(sat.signals, new_signals) + snapped_code_phase = + just_synced ? _snap_code_phase_from_synced_signal(new_signals, sat.code_phase) : sat.code_phase - end + final_signals = + just_synced ? map(_reset_inflight_integration, new_signals) : new_signals TrackedSat( sat; code_phase = snapped_code_phase, carrier_doppler = new_carrier_doppler, code_doppler = new_code_doppler, - signals = new_signals, + signals = final_signals, doppler_estimator_state = new_doppler_estimator_state, ) end +# Drop an in-flight (partial) integration: zero the accumulator and its sample +# counter, leaving all other per-signal state intact. Used at the sync-transition +# phase snap, where the shared `code_phase` moves and any partial accumulated at +# the old phase must not be carried into the re-anchored window. +@inline _reset_inflight_integration(s::TrackedSignal) = + TrackedSignal(s; correlator = zero(s.correlator), integrated_samples = 0) + # De-rotation applied to a component's bit-buffer prompt so its own energy is # real again, given the loops lock the driver onto the real axis. The rotation # is `cis(driver_carrier_phase − get_carrier_phase(signal))`, where the @@ -362,39 +375,46 @@ end @inline _carrier_phase_derotation(driver_carrier_phase::Real, signal) = cis(driver_carrier_phase - get_carrier_phase_offset(signal)) -# Shared per-signal advance applied after a completed integration — -# identical for the estimator-driver signal and the passenger signals: -# normalize the correlator, update/apply the post-corr filter, record the -# filtered prompt, advance the CN0 estimator and bit buffer, and rebuild -# the `TrackedSignal` with the correlator moved to -# `last_fully_integrated_*`. Returns the rebuilt signal plus the -# intermediate values the driver's loop-filter section needs -# (`filtered_correlator`, `integrated_code_blocks`). Keeping this in one -# place is what guarantees driver and passengers cannot drift (issue #133). +# Apply one completed `CorrelatorOutput` record to a signal — shared by the +# estimator-driver and passenger folds so they cannot drift (issue #133): +# normalize the record's (raw) correlator by its sample count, update/apply the +# post-corr filter, record the filtered prompt, advance the CN0 estimator and +# bit buffer, and rebuild the `TrackedSignal` with the record moved to +# `last_fully_integrated_*`. Returns the rebuilt signal plus the intermediate +# values the driver's loop-filter section needs (`filtered_correlator`, +# `integrated_code_blocks`). +# +# Unlike the old per-integration advance, this does NOT reset the live +# accumulator or `integrated_samples`: the correlate phase already reset them +# when it snapshotted this record and began (or is carrying) the next +# integration. It consumes only the record's stored correlator. # # The bit accumulator is credited with the blocks *actually* integrated -# (`calc_num_code_blocks_for_bit_buffer`), not the intended -# `integrated_code_blocks`: post-sync the first integration is truncated to -# land on the data-bit boundary, so crediting the intended length would -# misalign the decoded bits (issue #125). The intended count is still -# returned for the driver's `1/N` loop-bandwidth scaling. -@inline function _advance_signal_after_integration( +# (`calc_num_code_blocks_for_bit_buffer`), recovered from the record's sample +# count: post-sync the first integration is truncated to land on the data-bit +# boundary, so crediting the intended length would misalign the decoded bits +# (issue #125). The block count returned for the driver's `1/N` loop-bandwidth +# scaling is the same actual count (floored at 1): the bandwidth must pair with +# the record's true integration time, so it only switches when the integration +# actually lengthened — not already on the fold where sync was detected but the +# records were still single-block (see `_process_estimator_driver_signal`). +# `skip_bit_buffer = true` applies everything EXCEPT the bit-buffer update. +# Used for records that follow a bit/secondary sync detected earlier in the +# same fold: those records were correlated with pre-sync replicas (no +# secondary-code wipe-off), so accumulating them as if wiped would feed +# sign-corrupted prompts into the first post-sync bits. They re-enter cleanly +# next chunk, correlated at the snapped phase with the overlay in the replica. +@inline function _apply_correlator_output( tracked_signal::TrackedSignal, + output::CorrelatorOutput, prn::Integer, sampling_frequency, - driver_carrier_phase::Real = 0.0, + driver_carrier_phase::Real = 0.0; + skip_bit_buffer::Bool = false, ) signal = tracked_signal.signal - integrated_code_blocks = calc_num_code_blocks_to_integrate( - signal, - tracked_signal.preferred_num_code_blocks_to_integrate, - has_bit_or_secondary_code_been_found(tracked_signal.bit_buffer), - ) - normalized_correlator = normalize( - tracked_signal.correlator, - tracked_signal.integrated_samples, - get_code_amplitude(signal), - ) + normalized_correlator = + normalize(output.correlator, output.integrated_samples, get_code_amplitude(signal)) post_corr_filter = update(tracked_signal.post_corr_filter, get_prompt(normalized_correlator)) filtered_correlator = apply(post_corr_filter, normalized_correlator) @@ -403,37 +423,41 @@ end cn0_estimator = update(get_cn0_estimator(tracked_signal), prompt) bit_block_count = calc_num_code_blocks_for_bit_buffer( signal, - tracked_signal.integrated_samples, + output.integrated_samples, sampling_frequency, has_bit_or_secondary_code_been_found(tracked_signal.bit_buffer), ) + # Blocks this record actually covered, for the driver's `1/N` bandwidth + # scaling. The floor at 1 covers the fractional-block record right after a + # sync phase-snap accumulator reset, whose rounded block count can be 0. + integrated_code_blocks = max(1, bit_block_count) # De-rotate the prompt onto the driver's (real) phase frame before both the # secondary/bit sync search and the coherent bit accumulation inside # `buffer`, so a quadrature component's data lands on the real axis it is # decided on. No-op for the driver and for co-phased pairs. bit_prompt = prompt * _carrier_phase_derotation(driver_carrier_phase, signal) - bit_buffer = buffer(signal, prn, tracked_signal.bit_buffer, bit_block_count, bit_prompt) + bit_buffer = + skip_bit_buffer ? tracked_signal.bit_buffer : + buffer(signal, prn, tracked_signal.bit_buffer, bit_block_count, bit_prompt) new_signal = TrackedSignal( tracked_signal; - integrated_samples = 0, - is_integration_completed = false, last_fully_integrated_filtered_prompt = prompt, bit_buffer, cn0_estimator, post_corr_filter, - correlator = zero(tracked_signal.correlator), - last_fully_integrated_correlator = tracked_signal.correlator, + last_fully_integrated_correlator = output.correlator, ) return new_signal, filtered_correlator, integrated_code_blocks end -# Process the estimator-driver signal (signals[1]): if its integration -# completed, run the PLL/DLL plus prompt filter / CN0 / bit-buffer update -# and return new values for carrier_doppler, code_doppler, and -# doppler_estimator_state. Otherwise return unchanged values. This is -# where ConventionalPLLAndDLL hard-codes the "signals[1] drives the loop -# filter" rule — a custom AbstractDopplerEstimator may use any/all -# signals' state. +# Process the estimator-driver signal (signals[1]): fold over every +# `CorrelatorOutput` collected during this chunk, in order — running the PLL/DLL +# plus prompt filter / CN0 / bit-buffer update per record and threading the loop +# filters and FLL `previous_prompt` across them — then return the new +# doppler_estimator_state and the *last* record's carrier/code Doppler (the NCO +# is written once per chunk). With no outputs the Doppler holds. This is where +# ConventionalPLLAndDLL hard-codes the "signals[1] drives the loop filter" rule — +# a custom AbstractDopplerEstimator may use any/all signals' state. @inline function _process_estimator_driver_signal( tracked_signal::TrackedSignal, sat::TrackedSat, @@ -441,61 +465,93 @@ end sampling_frequency, driver_carrier_phase::Real = 0.0, ) - if !tracked_signal.is_integration_completed || tracked_signal.integrated_samples == 0 + outputs = tracked_signal.correlator_outputs + if isempty(outputs) return tracked_signal, pll_and_dll_state, sat.carrier_doppler, sat.code_doppler end signal = tracked_signal.signal - # The previous fully-integrated prompt must be read off the *old* - # signal — the shared advance overwrites it with this integration's. - previous_prompt = get_last_fully_integrated_filtered_prompt(tracked_signal) - integration_time = tracked_signal.integrated_samples / sampling_frequency - # The driver de-rotates against itself (offset 0), so this is a no-op for it; - # passed for symmetry with the passenger path. - new_signal, filtered_correlator, integrated_code_blocks = - _advance_signal_after_integration( - tracked_signal, + ts = tracked_signal + carrier_loop_filter = pll_and_dll_state.carrier_loop_filter + code_loop_filter = pll_and_dll_state.code_loop_filter + carrier_doppler = sat.carrier_doppler + code_doppler = sat.code_doppler + found_before_fold = has_bit_or_secondary_code_been_found(ts.bit_buffer) + @inbounds for k in eachindex(outputs) + output = outputs[k] + # FLL needs the previous record's filtered prompt; the first record of + # the chunk chains from the sat's carried-over + # `last_fully_integrated_filtered_prompt` (the previous chunk's last). + # Read it off `ts` BEFORE the advance overwrites it. + previous_prompt = get_last_fully_integrated_filtered_prompt(ts) + # Per-record integration time — the block time, NOT the chunk time. + integration_time = output.integrated_samples / sampling_frequency + # A record that follows a sync detected earlier in THIS fold was + # correlated with pre-sync replicas — keep it out of the bit buffer + # (see `_apply_correlator_output`). + synced_earlier_in_fold = + !found_before_fold && has_bit_or_secondary_code_been_found(ts.bit_buffer) + # The driver de-rotates against itself (offset 0), so the derotation is a + # no-op for it; passed for symmetry with the passenger path. + ts, filtered_correlator, integrated_code_blocks = _apply_correlator_output( + ts, + output, sat.prn, sampling_frequency, - driver_carrier_phase, + driver_carrier_phase; + skip_bit_buffer = synced_earlier_in_fold, ) - # The configured bandwidths are referenced to a one-primary-code-period - # integration. Coherently integrating `integrated_code_blocks` periods - # grows the loop update interval by that factor, so scale the effective - # bandwidth by 1/integrated_code_blocks to hold the loop's BL·Δt stability - # product at its single-period value. For the N=1 path this divides by 1 - # and is bit-identical to before. - carrier_bandwidth = - pll_and_dll_state.carrier_loop_filter_bandwidth / integrated_code_blocks - code_bandwidth = pll_and_dll_state.code_loop_filter_bandwidth / integrated_code_blocks - - carrier_freq_update, carrier_loop_filter = calculate_carrier_frequency_update( - signal, - pll_and_dll_state.carrier_loop_filter, - filtered_correlator, - previous_prompt, - integration_time, - carrier_bandwidth, - ) - code_freq_update, code_loop_filter = calculate_code_frequency_update( - signal, - pll_and_dll_state.code_loop_filter, - filtered_correlator, - sat.code_doppler, - sampling_frequency, - integration_time, - code_bandwidth, - ) - carrier_doppler, code_doppler = aid_dopplers( - signal, - pll_and_dll_state.init_carrier_doppler, - pll_and_dll_state.init_code_doppler, - carrier_freq_update, - code_freq_update, - ) + # The configured bandwidths are referenced to a one-primary-code-period + # integration. Coherently integrating N periods grows the loop update + # interval by that factor, so scale the effective bandwidth by 1/N to + # hold the loop's BL·Δt stability product at its single-period value. + # N (`integrated_code_blocks`) is the number of blocks this record + # ACTUALLY covered — recovered from its sample count — not the intended + # integration length: the bandwidth pairs with the record's true + # `integration_time`, so it only switches once the integrations really + # lengthen. Scaling by the intended length instead would under-gain the + # loop for single-block records folded after a mid-fold sync detection + # (correlated pre-sync, but the live bit buffer already reports the + # post-sync length) and for the first post-sync integration, which is + # truncated to land on the data-bit boundary. For the N=1 path this + # divides by 1 and is bit-identical to before. + carrier_bandwidth = + pll_and_dll_state.carrier_loop_filter_bandwidth / integrated_code_blocks + code_bandwidth = + pll_and_dll_state.code_loop_filter_bandwidth / integrated_code_blocks + + carrier_freq_update, carrier_loop_filter = calculate_carrier_frequency_update( + signal, + carrier_loop_filter, + filtered_correlator, + previous_prompt, + integration_time, + carrier_bandwidth, + ) + # `dll_disc` is fed the chunk-fixed `sat.code_doppler` — the code Doppler + # that actually generated this chunk's replicas — for every record; + # only the loop-filter *state* threads across records. + code_freq_update, code_loop_filter = calculate_code_frequency_update( + signal, + code_loop_filter, + filtered_correlator, + sat.code_doppler, + sampling_frequency, + integration_time, + code_bandwidth, + ) + carrier_doppler, code_doppler = aid_dopplers( + signal, + pll_and_dll_state.init_carrier_doppler, + pll_and_dll_state.init_code_doppler, + carrier_freq_update, + code_freq_update, + ) + end + empty!(outputs) new_doppler_estimator_state = SatConventionalPLLAndDLL(pll_and_dll_state; carrier_loop_filter, code_loop_filter) - return new_signal, new_doppler_estimator_state, carrier_doppler, code_doppler + return ts, new_doppler_estimator_state, carrier_doppler, code_doppler end # Process the non-driver signals (signals[2:end]): the shared per-signal @@ -528,17 +584,28 @@ end sampling_frequency, driver_carrier_phase::Real = 0.0, ) - if !tracked_signal.is_integration_completed || tracked_signal.integrated_samples == 0 - return tracked_signal + outputs = tracked_signal.correlator_outputs + isempty(outputs) && return tracked_signal + ts = tracked_signal + found_before_fold = has_bit_or_secondary_code_been_found(ts.bit_buffer) + @inbounds for k in eachindex(outputs) + # Same rule as the driver fold: records after a sync detected earlier + # in this fold stay out of the bit buffer. + synced_earlier_in_fold = + !found_before_fold && has_bit_or_secondary_code_been_found(ts.bit_buffer) + ts = first( + _apply_correlator_output( + ts, + outputs[k], + prn, + sampling_frequency, + driver_carrier_phase; + skip_bit_buffer = synced_earlier_in_fold, + ), + ) end - first( - _advance_signal_after_integration( - tracked_signal, - prn, - sampling_frequency, - driver_carrier_phase, - ), - ) + empty!(outputs) + ts end """ diff --git a/src/correlators/correlator.jl b/src/correlators/correlator.jl index 45a26219..f310fd7d 100644 --- a/src/correlators/correlator.jl +++ b/src/correlators/correlator.jl @@ -1,6 +1,33 @@ abstract type AbstractCorrelator{M} end abstract type AbstractEarlyPromptLateCorrelator{M} <: AbstractCorrelator{M} end +""" +$(SIGNATURES) + +A single completed correlator output produced within one processing chunk. + +`track!` processes each measurement in fixed-size time chunks. Every time a +signal's coherent integration completes inside a chunk, the (raw, +un-normalized) accumulator is snapshotted into a [`CorrelatorOutput`] and +appended to that signal's `correlator_outputs` array; the Doppler estimator +then folds over those records after the chunk. Storing the raw correlator plus +`integrated_samples` lets the estimator normalize it (dividing by the sample +count) and matches `last_fully_integrated_correlator`. + +Fields: + + - `correlator`: the raw accumulated correlator at completion (not normalized). + - `integrated_samples`: samples integrated into this output (for `normalize`, + the loop-filter `integration_time`, and the bit-buffer block count). + - `sample_index`: buffer-relative sample index at which this integration ended + (the epoch of the measurement; used by vector tracking). +""" +struct CorrelatorOutput{C<:AbstractCorrelator} + correlator::C + integrated_samples::Int + sample_index::Int +end + type_for_num_ants(num_ants::NumAnts{1}) = ComplexF64 type_for_num_ants(num_ants::NumAnts{N}) where {N} = SVector{N,ComplexF64} diff --git a/src/downconvert_and_correlate.jl b/src/downconvert_and_correlate.jl index 62070b5a..bfcb8f48 100644 --- a/src/downconvert_and_correlate.jl +++ b/src/downconvert_and_correlate.jl @@ -1,10 +1,16 @@ -# Per-sat update applied after a single (downconvert + correlate) step over -# `integrated_samples` samples. Advances the shared carrier and code phase, -# each per-signal correlator accumulator, integration counter, and the -# `is_integration_completed` flag (per signal). `signal_start_sample` -# advances on the satellite. Returns the new `TrackedSat`. +# Per-sat update applied after a single (downconvert + correlate) sub-step over +# `integrated_samples` samples. Advances the shared carrier and code phase and +# each per-signal correlator accumulator. `signal_start_sample` advances on the +# satellite. Returns the new `TrackedSat`. # -# `new_signals_data` is a tuple of `(new_correlator, is_integration_completed)` +# When a signal's coherent integration completes on this sub-step, its raw +# accumulator is snapshotted into a `CorrelatorOutput` (appended to the reused +# per-signal `correlator_outputs` buffer) and the accumulator is reset — so the +# next integration within the same processing chunk starts fresh. The Doppler +# estimator later folds over `correlator_outputs`; see +# `estimate_dopplers_and_filter_prompt!`. +# +# `new_signals_data` is a tuple of `(new_correlator, completed)` # pairs, one per element of `sat.signals`. Built via tuple recursion in # `_update_tracked_sat_correlator` so the heterogeneous walk stays # allocation-free and inferable. @@ -38,37 +44,60 @@ function update( code_frequency * integrated_samples / sampling_frequency + sat.code_phase, code_length, ) + new_signal_start_sample = sat.signal_start_sample + integrated_samples - new_signals = _build_new_signals(sat.signals, new_signals_data, integrated_samples) + new_signals = _build_new_signals( + sat.signals, + new_signals_data, + integrated_samples, + new_signal_start_sample, + ) TrackedSat( sat; code_phase, carrier_phase, - signal_start_sample = sat.signal_start_sample + integrated_samples, + signal_start_sample = new_signal_start_sample, signals = new_signals, ) end -# Recursive tuple walk that rebuilds each TrackedSignal with its -# corresponding (new_correlator, is_integration_completed) from -# `new_signals_data`. Allocation-free and inference-friendly. -@inline _build_new_signals(::Tuple{}, ::Tuple{}, ::Int) = () +# Recursive tuple walk that rebuilds each TrackedSignal from its +# `(new_correlator, completed)` in `new_signals_data`. +# Allocation-free and inference-friendly. +# +# On completion, snapshot the (raw) accumulator into the shared +# `correlator_outputs` vector — tagged with `sample_index` (the end sample of +# this integration) — and reset the accumulator. +# The `push!` mutates the same vector the copy-update constructor threads +# through unchanged, so it stays allocation-free after the buffer's capacity is +# seated. On a partial (chunk/buffer-bounded) sub-step, carry the accumulator. +@inline _build_new_signals(::Tuple{}, ::Tuple{}, ::Int, ::Int) = () @inline function _build_new_signals( signals::Tuple, new_data::Tuple, integrated_samples::Int, + signal_start_sample::Int, ) s = first(signals) (corr, completed) = first(new_data) - new_s = TrackedSignal( - s; - integrated_samples = s.integrated_samples + integrated_samples, - is_integration_completed = completed, - correlator = corr, - ) + total_integrated = s.integrated_samples + integrated_samples + if completed + push!( + s.correlator_outputs, + CorrelatorOutput(corr, total_integrated, signal_start_sample - 1), + ) + new_s = TrackedSignal(s; integrated_samples = 0, correlator = zero(corr)) + else + new_s = TrackedSignal(s; integrated_samples = total_integrated, correlator = corr) + end ( new_s, - _build_new_signals(Base.tail(signals), Base.tail(new_data), integrated_samples)..., + _build_new_signals( + Base.tail(signals), + Base.tail(new_data), + integrated_samples, + signal_start_sample, + )..., ) end diff --git a/src/downconvert_and_correlate_cpu.jl b/src/downconvert_and_correlate_cpu.jl index 3aaef294..319e02b5 100644 --- a/src/downconvert_and_correlate_cpu.jl +++ b/src/downconvert_and_correlate_cpu.jl @@ -385,52 +385,76 @@ function _update_tracked_sat_correlator( dc::AbstractDownconvertAndCorrelator, signal, num_samples_signal, + chunk_last_sample, sampling_frequency, intermediate_frequency, + stop_before_partial::Bool = false, ) - # MIN samples-to-next-boundary across all signals on this sat, clamped - # to remaining buffer. Each signal's coherent-integration length comes - # from its own `preferred_num_code_blocks_to_integrate` field; the - # signal-level boundary calc reads it per signal. - samples_to_integrate, per_signal_completed = _calc_min_samples_and_completed( - sat.signals, - sat.signal_start_sample, - sampling_frequency, - sat.code_doppler, - sat.code_phase, - num_samples_signal, - ) - if samples_to_integrate == 0 - return sat + # Integrate this sat forward through the current chunk. Each sub-step + # advances to the next code-block boundary of any signal (or the chunk + # end, whichever is nearer); `update` snapshots a `CorrelatorOutput` and + # resets the accumulator for every signal that just completed, so a chunk + # can yield 0, 1, or several outputs per signal. The NCO Doppler is + # untouched here. A partial integration at the chunk (or buffer) boundary + # carries in the accumulator. + # + # `stop_before_partial = true` stops at the last code-block boundary inside + # the chunk instead of integrating the trailing partial up to the chunk + # end. `track!`'s per-chunk pass uses this: the residue is left for the + # NEXT chunk's pass, which runs boundary → boundary in one kernel window, + # entirely at the Doppler the estimator wrote in between — so every + # completed integration is produced by a single Doppler and the NCO + # correction takes effect right at the completing boundary, like the + # pre-chunking per-completion update. A final pass without the flag drains + # the buffer's trailing partial into the accumulator. + while sat.signal_start_sample <= chunk_last_sample + # MIN samples-to-next-boundary across all signals on this sat, clamped + # to the chunk end. Each signal's coherent-integration length comes from + # its own `preferred_num_code_blocks_to_integrate`; replica sizing still + # uses the true buffer length `num_samples_signal`. + samples_to_integrate, per_signal_completed = _calc_min_samples_and_completed( + sat.signals, + sat.signal_start_sample, + sampling_frequency, + sat.code_doppler, + sat.code_phase, + num_samples_signal, + chunk_last_sample, + ) + samples_to_integrate == 0 && break + # A sub-step that completes no signal is exactly the chunk-clamped + # trailing partial — defer it to the post-estimate pass. + stop_before_partial && !_any_of_tuple(per_signal_completed) && break + carrier_frequency = sat.carrier_doppler + intermediate_frequency + new_signals_data = _correlate_signals( + sat.signals, + per_signal_completed, + dc, + signal, + sat.code_doppler, + sat.code_phase, + carrier_frequency, + sat.carrier_phase, + sampling_frequency, + sat.signal_start_sample, + samples_to_integrate, + sat.prn, + num_samples_signal, + ) + sat = update( + sat, + samples_to_integrate, + intermediate_frequency, + sampling_frequency, + new_signals_data, + ) end - carrier_frequency = sat.carrier_doppler + intermediate_frequency - new_signals_data = _correlate_signals( - sat.signals, - per_signal_completed, - dc, - signal, - sat.code_doppler, - sat.code_phase, - carrier_frequency, - sat.carrier_phase, - sampling_frequency, - sat.signal_start_sample, - samples_to_integrate, - sat.prn, - num_samples_signal, - ) - update( - sat, - samples_to_integrate, - intermediate_frequency, - sampling_frequency, - new_signals_data, - ) + return sat end # Compute (samples_to_integrate, per_signal_completed_tuple) via tuple # recursion. Returns the MIN across all signals' samples-to-next-boundary, -# clamped to remaining buffer samples. Per-signal `completed` flags are +# clamped to the current chunk end. Per-signal `completed` flags are # derived after the MIN is known. @inline function _calc_min_samples_and_completed( signals::Tuple, @@ -439,6 +463,7 @@ end code_doppler, code_phase, num_samples_signal, + chunk_last_sample, ) per_signal_to_boundary = _per_signal_samples_to_boundary( signals, @@ -449,8 +474,10 @@ end num_samples_signal, ) samples_to_integrate = _min_of_tuple(per_signal_to_boundary) - signal_samples_left = num_samples_signal - signal_start_sample + 1 - samples_to_integrate = min(samples_to_integrate, signal_samples_left) + # Clamp the integration to the end of the current chunk (not the buffer): + # the boundary calc / replica sizing above still see the true buffer length. + samples_left = chunk_last_sample - signal_start_sample + 1 + samples_to_integrate = min(samples_to_integrate, samples_left) per_signal_completed = _flag_completed(per_signal_to_boundary, samples_to_integrate) return samples_to_integrate, per_signal_completed end @@ -506,6 +533,9 @@ end @inline _min_of_tuple(t::Tuple{Any}) = first(t) @inline _min_of_tuple(t::Tuple) = min(first(t), _min_of_tuple(Base.tail(t))) +@inline _any_of_tuple(::Tuple{}) = false +@inline _any_of_tuple(t::Tuple) = first(t) || _any_of_tuple(Base.tail(t)) + # For each signal, `is_completed = (chosen_samples == samples_to_boundary)`. @inline _flag_completed(::Tuple{}, _) = () @inline _flag_completed(t::Tuple, chosen) = @@ -799,23 +829,71 @@ backend, different `@batch` iterations write to disjoint slots, so no synchronization is needed. Returns the same `track_state`. Allocation-free in steady state — see [`track!`](@ref). +`chunk_duration` and `chunk_index` restrict the pass to one chunk of a fixed +per-band time grid: each satellite integrates only up to sample +`min(round((chunk_index + 1) * chunk_duration * sampling_frequency), num_samples)`. +The boundary is re-anchored to the absolute `chunk_index` each call (not +accumulated), so rounding never drifts and different bands stay time-aligned. +`chunk_duration = nothing` (the default) disables chunking — the whole buffer +is consumed as one chunk. Every completed integration is snapshotted into its +signal's `correlator_outputs`; only the trailing partial stays in the live +correlator. + +`stop_before_partial = true` additionally stops each satellite at its last +completed code-block boundary inside the chunk instead of integrating the +chunk-clamped trailing partial. `track!`'s per-chunk pass uses this so the +residue is integrated by the *next* chunk's pass, entirely at the Doppler the +estimator writes in between; leave it `false` (the default) to consume the +whole window. + `samples_unchanged = true` promises that every band's sample buffer holds the same content as on the previous call with this `dc`; backends may then reuse sample-derived caches (the bit-wise backends skip re-packing their shared -band sign planes). `track!` passes it on every loop iteration after the first -so the pack happens once per call; leave it `false` (the default) whenever -the buffers may have been refilled. +band sign planes). `track!` passes it on every chunk after the first so the +pack happens once per call; leave it `false` (the default) whenever the +buffers may have been refilled. """ function downconvert_and_correlate!( dc::AbstractDownconvertAndCorrelator, measurements::BandMeasurements, track_state::TrackState; + chunk_index::Int = 0, + chunk_duration = nothing, + stop_before_partial::Bool = false, samples_unchanged::Bool = false, ) - _foreach_group!(_dc_one_group!, track_state.groups, dc, measurements, samples_unchanged) + _foreach_group!( + _dc_one_group!, + track_state.groups, + dc, + measurements, + chunk_index, + chunk_duration, + stop_before_partial, + samples_unchanged, + ) return track_state end +# Last sample index (inclusive) this satellite may integrate up to in the +# current chunk. `chunk_duration === nothing` means "no chunking" — consume the +# whole buffer, i.e. today's behavior and the default for direct callers. When a +# chunk duration is given, the boundary lies on a shared per-band time grid; it +# is re-anchored to the absolute `chunk_index` each call (not accumulated) so +# rounding never drifts and different bands stay time-aligned to within a +# sample. +@inline _chunk_last_sample(::Nothing, chunk_index, sampling_frequency, num_samples) = + num_samples +@inline function _chunk_last_sample( + chunk_duration, + chunk_index, + sampling_frequency, + num_samples, +) + grid = uconvert(NoUnits, (chunk_index + 1) * chunk_duration * sampling_frequency) + min(round(Int, grid), num_samples) +end + # Optional per-backend sample-type check, run once per group before the # per-sat loop. No-op by default; the integer backends override it to reject # non-`Complex{Int16}` sample buffers with a helpful `ArgumentError` (see @@ -834,6 +912,9 @@ end g::SignalGroup, dc::AbstractDownconvertAndCorrelator, measurements::BandMeasurements, + chunk_index::Int, + chunk_duration, + stop_before_partial::Bool, # The float backends derive nothing from the raw samples worth caching; # only the bit-wise backends act on `samples_unchanged`. samples_unchanged::Bool, @@ -842,13 +923,18 @@ end isempty(vals) && return nothing m = measurements[get_band_id(g.band)] _check_sample_type(dc, m) + num_samples = get_num_samples(m) + chunk_last_sample = + _chunk_last_sample(chunk_duration, chunk_index, m.sampling_frequency, num_samples) _dc_group_loop!( dc, vals, m.samples, - get_num_samples(m), + num_samples, + chunk_last_sample, m.sampling_frequency, m.intermediate_frequency, + stop_before_partial, ) end @@ -863,17 +949,17 @@ struct _BatchLoop end @inline _threading(::AbstractDownconvertAndCorrelator) = _SerialLoop() @inline _threading(::CPUThreadedDownconvertAndCorrelator) = _BatchLoop() -@inline _dc_group_loop!(dc::AbstractDownconvertAndCorrelator, vals, args::Vararg{Any,4}) = +@inline _dc_group_loop!(dc::AbstractDownconvertAndCorrelator, vals, args::Vararg{Any,6}) = _dc_group_loop!(_threading(dc), dc, vals, args...) -@inline function _dc_group_loop!(::_SerialLoop, dc, vals, args::Vararg{Any,4}) +@inline function _dc_group_loop!(::_SerialLoop, dc, vals, args::Vararg{Any,6}) @inbounds for i in eachindex(vals) vals[i] = _update_tracked_sat_correlator(vals[i], dc, args...) end return nothing end -@inline function _dc_group_loop!(::_BatchLoop, dc, vals, args::Vararg{Any,4}) +@inline function _dc_group_loop!(::_BatchLoop, dc, vals, args::Vararg{Any,6}) @batch for i = 1:length(vals) @inbounds vals[i] = _update_tracked_sat_correlator(vals[i], dc, args...) end diff --git a/src/downconvert_and_correlate_onebit.jl b/src/downconvert_and_correlate_onebit.jl index b4fc4ad4..43ac689d 100644 --- a/src/downconvert_and_correlate_onebit.jl +++ b/src/downconvert_and_correlate_onebit.jl @@ -1237,52 +1237,18 @@ end map(tuple, new_corrs, per_signal_completed) end -function _update_tracked_sat_correlator( - sat::TrackedSat, - dc::_OneBitDC, - signal, - num_samples_signal, - sampling_frequency, - intermediate_frequency, -) - samples_to_integrate, per_signal_completed = _calc_min_samples_and_completed( - sat.signals, - sat.signal_start_sample, - sampling_frequency, - sat.code_doppler, - sat.code_phase, - num_samples_signal, - ) - samples_to_integrate == 0 && return sat - carrier_frequency = sat.carrier_doppler + intermediate_frequency - new_signals_data = _correlate_signals( - sat.signals, - per_signal_completed, - dc, - signal, - sat.code_doppler, - sat.code_phase, - carrier_frequency, - sat.carrier_phase, - sampling_frequency, - sat.signal_start_sample, - samples_to_integrate, - sat.prn, - num_samples_signal, - ) - update( - sat, - samples_to_integrate, - intermediate_frequency, - sampling_frequency, - new_signals_data, - ) -end +# The per-sat chunk loop is the generic +# `_update_tracked_sat_correlator(sat, dc::AbstractDownconvertAndCorrelator, …)` +# in downconvert_and_correlate_cpu.jl — the backend boundary is +# `_correlate_signals`, which dispatches on `dc::_OneBitDC` above. @inline function _dc_one_group!( g::SignalGroup, dc::_OneBitDC, measurements::BandMeasurements, + chunk_index::Int, + chunk_duration, + stop_before_partial::Bool, samples_unchanged::Bool, ) vals = g.satellites.values @@ -1304,8 +1270,8 @@ end samples = m.samples if samples_unchanged # Same buffer content as the previous call with this dc (track! passes - # this on every loop iteration after the first): the shared band pack — - # or the emptied single-sat state — is still valid, keep it. + # this on every chunk/pass after the first): the shared band pack — or + # the emptied single-sat state — is still valid, keep it. elseif length(vals) > 1 nsamp = get_num_samples(m) M = samples isa AbstractMatrix ? size(samples, 2) : 1 @@ -1322,20 +1288,25 @@ end resize!(dc.band.mrband, 0) resize!(dc.band.miband, 0) end + num_samples = get_num_samples(m) + chunk_last_sample = + _chunk_last_sample(chunk_duration, chunk_index, m.sampling_frequency, num_samples) _dc_group_loop!( dc, vals, m.samples, - get_num_samples(m), + num_samples, + chunk_last_sample, m.sampling_frequency, m.intermediate_frequency, + stop_before_partial, ) end @inline function _dc_group_loop!( dc::OneBitDownconvertAndCorrelator, vals, - args::Vararg{Any,4}, + args::Vararg{Any,6}, ) @inbounds for i in eachindex(vals) vals[i] = _update_tracked_sat_correlator(vals[i], dc, args...) @@ -1346,7 +1317,7 @@ end @inline function _dc_group_loop!( dc::OneBitThreadedDownconvertAndCorrelator, vals, - args::Vararg{Any,4}, + args::Vararg{Any,6}, ) @batch for i = 1:length(vals) @inbounds vals[i] = _update_tracked_sat_correlator(vals[i], dc, args...) @@ -1379,8 +1350,20 @@ function downconvert_and_correlate!( dc::_OneBitDC, measurements::BandMeasurements, track_state::TrackState; + chunk_index::Int = 0, + chunk_duration = nothing, + stop_before_partial::Bool = false, samples_unchanged::Bool = false, ) - _foreach_group!(_dc_one_group!, track_state.groups, dc, measurements, samples_unchanged) + _foreach_group!( + _dc_one_group!, + track_state.groups, + dc, + measurements, + chunk_index, + chunk_duration, + stop_before_partial, + samples_unchanged, + ) return track_state end diff --git a/src/downconvert_and_correlate_twobit.jl b/src/downconvert_and_correlate_twobit.jl index c418a0bf..18caca4d 100644 --- a/src/downconvert_and_correlate_twobit.jl +++ b/src/downconvert_and_correlate_twobit.jl @@ -1722,52 +1722,18 @@ end map(tuple, new_corrs, per_signal_completed) end -function _update_tracked_sat_correlator( - sat::TrackedSat, - dc::_TwoBitDC, - signal, - num_samples_signal, - sampling_frequency, - intermediate_frequency, -) - samples_to_integrate, per_signal_completed = _calc_min_samples_and_completed( - sat.signals, - sat.signal_start_sample, - sampling_frequency, - sat.code_doppler, - sat.code_phase, - num_samples_signal, - ) - samples_to_integrate == 0 && return sat - carrier_frequency = sat.carrier_doppler + intermediate_frequency - new_signals_data = _correlate_signals( - sat.signals, - per_signal_completed, - dc, - signal, - sat.code_doppler, - sat.code_phase, - carrier_frequency, - sat.carrier_phase, - sampling_frequency, - sat.signal_start_sample, - samples_to_integrate, - sat.prn, - num_samples_signal, - ) - update( - sat, - samples_to_integrate, - intermediate_frequency, - sampling_frequency, - new_signals_data, - ) -end +# The per-sat chunk loop is the generic +# `_update_tracked_sat_correlator(sat, dc::AbstractDownconvertAndCorrelator, …)` +# in downconvert_and_correlate_cpu.jl — the backend boundary is +# `_correlate_signals`, which dispatches on `dc::_TwoBitDC` above. @inline function _dc_one_group!( g::SignalGroup, dc::_TwoBitDC, measurements::BandMeasurements, + chunk_index::Int, + chunk_duration, + stop_before_partial::Bool, samples_unchanged::Bool, ) vals = g.satellites.values @@ -1789,8 +1755,8 @@ end samples = m.samples if samples_unchanged # Same buffer content as the previous call with this dc (track! passes - # this on every loop iteration after the first): the shared band pack — - # or the emptied single-sat state — is still valid, keep it. + # this on every chunk/pass after the first): the shared band pack — or + # the emptied single-sat state — is still valid, keep it. elseif length(vals) > 1 nsamp = get_num_samples(m) M = samples isa AbstractMatrix ? size(samples, 2) : 1 @@ -1811,20 +1777,25 @@ end resize!(dc.band.gmrband, 0) resize!(dc.band.gmiband, 0) end + num_samples = get_num_samples(m) + chunk_last_sample = + _chunk_last_sample(chunk_duration, chunk_index, m.sampling_frequency, num_samples) _dc_group_loop!( dc, vals, m.samples, - get_num_samples(m), + num_samples, + chunk_last_sample, m.sampling_frequency, m.intermediate_frequency, + stop_before_partial, ) end @inline function _dc_group_loop!( dc::TwoBitDownconvertAndCorrelator, vals, - args::Vararg{Any,4}, + args::Vararg{Any,6}, ) @inbounds for i in eachindex(vals) vals[i] = _update_tracked_sat_correlator(vals[i], dc, args...) @@ -1835,7 +1806,7 @@ end @inline function _dc_group_loop!( dc::TwoBitThreadedDownconvertAndCorrelator, vals, - args::Vararg{Any,4}, + args::Vararg{Any,6}, ) @batch for i = 1:length(vals) @inbounds vals[i] = _update_tracked_sat_correlator(vals[i], dc, args...) @@ -1868,8 +1839,20 @@ function downconvert_and_correlate!( dc::_TwoBitDC, measurements::BandMeasurements, track_state::TrackState; + chunk_index::Int = 0, + chunk_duration = nothing, + stop_before_partial::Bool = false, samples_unchanged::Bool = false, ) - _foreach_group!(_dc_one_group!, track_state.groups, dc, measurements, samples_unchanged) + _foreach_group!( + _dc_one_group!, + track_state.groups, + dc, + measurements, + chunk_index, + chunk_duration, + stop_before_partial, + samples_unchanged, + ) return track_state end diff --git a/src/sat_state.jl b/src/sat_state.jl index 4df79d08..3e3b22e5 100644 --- a/src/sat_state.jl +++ b/src/sat_state.jl @@ -18,7 +18,6 @@ struct TrackedSignal{ } signal::Sig integrated_samples::Int - is_integration_completed::Bool correlator::C last_fully_integrated_correlator::C last_fully_integrated_filtered_prompt::ComplexF64 @@ -26,6 +25,12 @@ struct TrackedSignal{ bit_buffer::BitBuffer{B} post_corr_filter::PCF filtered_prompts::Vector{ComplexF64} + # Correlator outputs completed within the current processing chunk, each + # tagged with its end sample index. Preallocated and reused (like + # `filtered_prompts`): the correlate phase `push!`es a record per completed + # integration, the Doppler estimator folds over them and `empty!`s the + # array. Empty at chunk boundaries; see [`CorrelatorOutput`]. + correlator_outputs::Vector{CorrelatorOutput{C}} # Preferred coherent-integration length for THIS signal, in primary code # blocks. The actual length is capped per integration by the signal's # bit/secondary-code period and held at 1 until bit/secondary sync (see @@ -98,10 +103,15 @@ function TrackedSignal( # Picking the type here makes `B` concrete in the resulting # `TrackedSignal{Sig, B, C, PCF}`. B = get_code_block_buffer_type(signal) + # Preallocate the per-chunk correlator-output buffer. A default-length + # chunk (= smallest code period) yields at most one record per chunk for + # the driver; size to a small constant so `push!` never grows it after + # warmup even for a moderately enlarged `doppler_update_interval`. + correlator_outputs = CorrelatorOutput{typeof(correlator)}[] + sizehint!(correlator_outputs, 4) TrackedSignal( signal, 0, - false, correlator, correlator, complex(0.0, 0.0), @@ -109,6 +119,7 @@ function TrackedSignal( BitBuffer{B}(), post_corr_filter, ComplexF64[], + correlator_outputs, preferred_num_code_blocks_to_integrate, ) end @@ -122,7 +133,6 @@ function TrackedSignal( t::TrackedSignal{Sig,B,C,PCF}; signal = nothing, integrated_samples = nothing, - is_integration_completed = nothing, correlator::Maybe{C} = nothing, last_fully_integrated_correlator::Maybe{C} = nothing, last_fully_integrated_filtered_prompt = nothing, @@ -130,6 +140,7 @@ function TrackedSignal( bit_buffer::Maybe{BitBuffer{B}} = nothing, post_corr_filter::Maybe{PCF} = nothing, filtered_prompts::Maybe{Vector{ComplexF64}} = nothing, + correlator_outputs::Maybe{Vector{CorrelatorOutput{C}}} = nothing, preferred_num_code_blocks_to_integrate = nothing, ) where { Sig<:AbstractGNSSSignal, @@ -145,8 +156,6 @@ function TrackedSignal( TrackedSignal{Sig,B,C,PCF}( isnothing(signal) ? t.signal : signal, isnothing(integrated_samples) ? t.integrated_samples : integrated_samples, - isnothing(is_integration_completed) ? t.is_integration_completed : - is_integration_completed, isnothing(correlator) ? t.correlator : correlator, isnothing(last_fully_integrated_correlator) ? t.last_fully_integrated_correlator : last_fully_integrated_correlator, @@ -156,6 +165,7 @@ function TrackedSignal( isnothing(bit_buffer) ? t.bit_buffer : bit_buffer, isnothing(post_corr_filter) ? t.post_corr_filter : post_corr_filter, isnothing(filtered_prompts) ? t.filtered_prompts : filtered_prompts, + isnothing(correlator_outputs) ? t.correlator_outputs : correlator_outputs, isnothing(preferred_num_code_blocks_to_integrate) ? t.preferred_num_code_blocks_to_integrate : preferred_num_code_blocks_to_integrate, ) @@ -167,6 +177,17 @@ get_last_fully_integrated_correlator(t::TrackedSignal) = t.last_fully_integrated get_last_fully_integrated_filtered_prompt(t::TrackedSignal) = t.last_fully_integrated_filtered_prompt get_filtered_prompts(t::TrackedSignal) = t.filtered_prompts + +""" +$(SIGNATURES) + +The [`CorrelatorOutput`](@ref)s this signal completed during the most recent +processing chunk, in order. Populated by the correlate phase and consumed + +cleared by the Doppler estimator after each chunk, so it is empty between +`track!` calls; read it inside a custom estimator, or right after a bare +`downconvert_and_correlate!`. See [Chunked Doppler updates](@ref). +""" +get_correlator_outputs(t::TrackedSignal) = t.correlator_outputs get_post_corr_filter(t::TrackedSignal) = t.post_corr_filter get_cn0_estimator(t::TrackedSignal) = t.cn0_estimator get_bit_buffer(t::TrackedSignal) = t.bit_buffer @@ -702,6 +723,7 @@ end @inline function _reset_signal(t::TrackedSignal) empty!(t.filtered_prompts) + empty!(t.correlator_outputs) TrackedSignal(t; bit_buffer = reset(t.bit_buffer)) end diff --git a/src/track.jl b/src/track.jl index b525a105..f2e32d5a 100644 --- a/src/track.jl +++ b/src/track.jl @@ -26,11 +26,14 @@ each group's key set and slot vector are copied, so [`add_satellite!`](@ref) / [`remove_satellite!`](@ref) and tracking itself on either state never affect the other's satellites. The copy is shallow, however — per-satellite scratch vectors (each signal's -`filtered_prompts`, the soft-bit buffer, and the CN0 estimator's -prompt buffer) are shared with the input and are overwritten by the -next `track` call on either state. Treat the input as a stale handle -after the call; `deepcopy` it first if you need to snapshot those -buffers. +`filtered_prompts` and `correlator_outputs`, the soft-bit buffer, and +the CN0 estimator's prompt buffer) are shared with the input and are +overwritten by the next `track` call on either state. Treat the input +as a stale handle after the call; `deepcopy` it first if you need to +snapshot those buffers. The same applies to a bare +`downconvert_and_correlate`: the returned state's +`correlator_outputs` alias the input's, so reuse of one input state +across several calls appends to the same buffers. For real-time loops processing many chunks of signal in sequence, **construct the correlator once outside the loop** and pass it via the @@ -171,25 +174,60 @@ function track!( measurements::BandMeasurements, track_state::TrackState; downconvert_and_correlator::AbstractDownconvertAndCorrelator = CPUThreadedDownconvertAndCorrelator(), + doppler_update_interval = nothing, ) _validate_measurements(track_state, measurements) reset_start_sample_and_bit_buffer!(track_state) - # The measurement buffers are fixed for the whole call, so sample-derived - # backend caches (the bit backends' shared band pack) are built on the - # first iteration and reused ever after (`samples_unchanged`). - first_iteration = true - while true - _all_groups_reached_end(track_state, measurements) && break - + # Resolve the Doppler-update / chunk interval. `nothing` => auto: the + # smallest primary-code period across all tracked signals, so a default + # chunk holds one code period of the shortest signal. The measurement is + # walked chunk by chunk: each chunk correlates (collecting every completed + # correlator output per signal into its `correlator_outputs` buffer) with + # the NCO Doppler held fixed, then the estimator folds over those outputs + # and updates every sat's NCO once — a common epoch across all sats. + chunk_duration = _resolve_doppler_update_interval(doppler_update_interval, track_state) + _validate_doppler_update_interval(chunk_duration, measurements) + # One correlate pass + one estimate per chunk. The pass runs each satellite + # from wherever it stands to its last completed code-block boundary inside + # the chunk (`stop_before_partial` — the chunk-clamped trailing partial is + # NOT integrated); the estimator then folds the collected outputs and + # writes the new NCO Doppler. The residue is picked up by the NEXT chunk's + # pass, which therefore covers boundary → boundary in a single kernel + # window, entirely at the just-updated Doppler. So every completed + # integration is produced by a single NCO Doppler and each correction takes + # effect right at its completing boundary (the classic per-completion loop + # timing), while the estimator still runs once per chunk at a common epoch + # — without splitting each code period into two kernel invocations. + # + # `samples_unchanged`: the measurement buffers are fixed for the whole + # call, so sample-derived backend caches (the bit backends' shared band + # pack) are built on the very first pass and reused ever after. + chunk_index = 0 + while _chunks_left(chunk_duration, chunk_index, measurements) downconvert_and_correlate!( downconvert_and_correlator, measurements, track_state; - samples_unchanged = !first_iteration, + chunk_index, + chunk_duration, + stop_before_partial = true, + samples_unchanged = chunk_index > 0, ) - first_iteration = false estimate_dopplers_and_filter_prompt!(track_state, measurements) + chunk_index += 1 end + # Drain the buffer: consume every satellite's trailing partial — from its + # last completed boundary to the buffer end — into its live accumulator (at + # the final chunk's Doppler), so the integration carries into the next + # `track!` call. A boundary landing exactly on the buffer end completes + # here, so fold once more; a no-op (per-sat early return) otherwise. + downconvert_and_correlate!( + downconvert_and_correlator, + measurements, + track_state; + samples_unchanged = chunk_index > 0, + ) + estimate_dopplers_and_filter_prompt!(track_state, measurements) return track_state end @@ -210,25 +248,76 @@ function track!(measurement::BandMeasurement, track_state::TrackState; kwargs... track!(_single_band_measurements(measurement, track_state), track_state; kwargs...) end -# Loop termination: every group has consumed its band's measurement to the -# end. Each group's `signal_start_sample` advances to `num_samples + 1` -# when that group's measurement is fully integrated; the outer `while` -# exits once every group has reached its own band's measurement end. -@inline function _all_groups_reached_end( - track_state::TrackState, - measurements::BandMeasurements, +# Loop termination: the chunk grid is walked until the previous chunk's end +# already reached the buffer end on every band (`_chunk_last_sample` clamps at +# `num_samples`, so the count is finite and independent of per-sat progress — +# satellites lagging behind a chunk boundary are caught up by later passes and +# by `track!`'s final buffer-draining pass). For `chunk_index == 0` the +# convention `_chunk_last_sample(…, -1, …) == 0` makes this "is the buffer +# non-empty on any band". +@inline function _chunks_left(chunk_duration, chunk_index::Int, measurements) + for m in measurements + n = get_num_samples(m) + _chunk_last_sample(chunk_duration, chunk_index - 1, m.sampling_frequency, n) < n && + return true + end + false +end + +# Resolve the per-chunk update interval to a concrete time. `nothing` => auto: +# the smallest primary-code period across every signal in every group, so the +# default chunk holds exactly one code period of the shortest signal (e.g. 1 ms +# for a GPS L1 C/A + Galileo E1B track state). +@inline _resolve_doppler_update_interval(doppler_update_interval, ::TrackState) = + doppler_update_interval +@inline function _resolve_doppler_update_interval(::Nothing, track_state::TrackState) + _smallest_code_period(track_state) +end + +# Primary-code period (a time) of one signal. +@inline _code_period(signal::AbstractGNSSSignal) = + get_code_length(signal) / get_code_frequency(signal) + +@inline _min_signal_code_period(::Tuple{}, m) = m +@inline _min_signal_code_period(signals::Tuple, m) = + _min_signal_code_period(Base.tail(signals), min(m, _code_period(first(signals)))) + +@inline _min_group_code_period(::Tuple{}, m) = m +@inline _min_group_code_period(groups::Tuple, m) = _min_group_code_period( + Base.tail(groups), + _min_signal_code_period(first(groups).signals, m), ) - _check_all_groups_at_end(Tuple(track_state.groups), measurements) + +function _smallest_code_period(track_state::TrackState) + groups = Tuple(track_state.groups) + # Every TrackState has at least one group and every group at least one + # signal, so seeding from the first signal's period is safe and keeps the + # reduction type-stable across heterogeneous signal tuples. + init = _code_period(first(first(groups).signals)) + _min_group_code_period(groups, init) end -# Recursive tuple-walk: each step has fully concrete types. -@inline _check_all_groups_at_end(::Tuple{}, ::BandMeasurements) = true -@inline function _check_all_groups_at_end(t::Tuple, measurements::BandMeasurements) - g = first(t) - m = measurements[get_band_id(g.band)] - target = get_num_samples(m) + 1 - @inbounds for sat in g.satellites.values - sat.signal_start_sample == target || return false +# A chunk must cover at least one sample on every band, otherwise the chunk +# grid could fail to advance and `track!` would not terminate. The default +# (smallest code period) is always many samples; this only guards against a +# user-supplied `doppler_update_interval` shorter than a sample period. The dimension +# check turns a plain number (e.g. `doppler_update_interval = 1e-3`) into a clear +# ArgumentError instead of a cryptic Unitful conversion error. +function _validate_doppler_update_interval(chunk_duration, measurements::BandMeasurements) + dimension(chunk_duration) == dimension(1.0s) || throw( + ArgumentError( + "doppler_update_interval must be a time quantity, e.g. `1u\"ms\"` or `1e-3u\"s\"` " * + "(with `using Unitful`); got $chunk_duration.", + ), + ) + for m in measurements + samples_per_chunk = uconvert(NoUnits, chunk_duration * m.sampling_frequency) + samples_per_chunk >= 1 || throw( + ArgumentError( + "doppler_update_interval $chunk_duration is shorter than one sample period " * + "at sampling frequency $(m.sampling_frequency); pick a longer interval.", + ), + ) end - _check_all_groups_at_end(Base.tail(t), measurements) + nothing end diff --git a/src/tracking_state.jl b/src/tracking_state.jl index 4934c851..866bdde7 100644 --- a/src/tracking_state.jl +++ b/src/tracking_state.jl @@ -301,8 +301,8 @@ function reset_start_sample_and_bit_buffer!(track_state::TrackState) end # Loop-termination helper for `track`/`track!` lives in `track.jl` as -# `_all_groups_reached_end` — it iterates the per-band measurement -# lengths so each group can terminate against its own band's chunk. +# `_chunks_left` — it iterates the per-band measurement lengths so the +# chunk grid terminates against each band's own buffer end. """ $(SIGNATURES) diff --git a/test/conventional_pll_and_dll.jl b/test/conventional_pll_and_dll.jl index 26b2dc17..cdf53eeb 100644 --- a/test/conventional_pll_and_dll.jl +++ b/test/conventional_pll_and_dll.jl @@ -2,6 +2,7 @@ module ConventionalPLLAndDLLTest using Test: @test, @testset, @inferred, @test_throws using Unitful: Hz +import Tracking using GNSSSignals: GPSL1CA, get_code_center_frequency_ratio using TrackingLoopFilters: ThirdOrderBilinearLF, SecondOrderBilinearLF using StaticArrays: SVector @@ -25,6 +26,7 @@ using Tracking: get_sat_states, update_accumulator, get_default_correlator, + CorrelatorOutput, merge_sats # Build a stub `(L1 = BandMeasurement(...),)` NamedTuple to pass to the @@ -33,6 +35,14 @@ using Tracking: # suffices. _meas_l1(fs) = (L1 = BandMeasurement(ComplexF64[], fs),) +# Build a signal carrying one completed integration as a `CorrelatorOutput` +# record — the estimate phase folds over `correlator_outputs`. `sample_index` +# is metadata here (no vector-tracking consumer). +_completed_signal(sig, correlator, num_samples) = TrackedSignal( + sig; + correlator_outputs = [CorrelatorOutput(correlator, num_samples, num_samples)], +) + @testset "Doppler aiding" begin gpsl1 = GPSL1CA() init_carrier_doppler = 10Hz @@ -140,14 +150,7 @@ end num_samples = 5000 sat_state_after_full_integration = TrackedSat( sat_state; - signals = ( - TrackedSignal( - only(sat_state.signals); - is_integration_completed = true, - integrated_samples = num_samples, - correlator, - ), - ), + signals = (_completed_signal(only(sat_state.signals), correlator, num_samples),), ) track_state = TrackState(gpsl1, sat_state_after_full_integration; doppler_estimator) @@ -171,6 +174,69 @@ end end end +@testset "loop bandwidth follows the record's actual integration length" begin + # The `1/N` bandwidth scaling must pair with the blocks a record ACTUALLY + # covered (recovered from its sample count), not the intended integration + # length: a single-block record folded when the bit buffer already reports + # sync — a mid-fold sync detection with an enlarged + # `doppler_update_interval`, or the truncated first post-sync integration — + # must be filtered at the full single-period bandwidth. The Doppler update + # therefore depends only on the record itself, not on the preferred + # integration length or the sync state. + sampling_frequency = 5e6Hz + gpsl1 = GPSL1CA() + carrier_doppler = 100.0Hz + code_phase = 0.5 + correlator = update_accumulator( + get_default_correlator(gpsl1), + SVector(1000.0 + 10im, 2000.0 + 20im, 750.0 + 10im), + ) + doppler_estimator = ConventionalPLLAndDLL() + + # Same-typed bit buffer with the bit/secondary sync already found + # (secondary phase 0, polarity +1). + synced(::Tracking.BitBuffer{B}) where {B} = Tracking.BitBuffer{B}( + zero(B), + 0, + true, + 0, + Int8(+1), + zero(UInt128), + 0, + complex(0.0, 0.0), + 0, + Float32[], + Tracking.PhaseAccumulators(), + ) + + doppler_after(preferred, found, integrated_samples) = begin + sat = TrackedSat(gpsl1, 1, code_phase, carrier_doppler; doppler_estimator) + base = only(sat.signals) + sig = TrackedSignal( + base; + bit_buffer = found ? synced(base.bit_buffer) : base.bit_buffer, + preferred_num_code_blocks_to_integrate = preferred, + correlator_outputs = [ + CorrelatorOutput(correlator, integrated_samples, integrated_samples), + ], + ) + ts = TrackState(gpsl1, TrackedSat(sat; signals = (sig,)); doppler_estimator) + get_carrier_doppler( + estimate_dopplers_and_filter_prompt(ts, _meas_l1(sampling_frequency)), + ) + end + + one_block = 5000 # one 1 ms L1 C/A code period at 5 MHz + # Synced with a 20-block preferred length, but the record covered a single + # block: full bandwidth — exactly the plain single-block baseline. (Scaling + # by the intended length would divide the bandwidth by 20 here.) + @test doppler_after(20, true, one_block) == doppler_after(1, false, one_block) + + # A record that actually covered 20 blocks is scaled by 1/20 regardless of + # the preferred integration length. + @test doppler_after(20, true, 20 * one_block) == doppler_after(1, true, 20 * one_block) +end + @testset "Per-satellite bandwidths drive the loop filters" begin sampling_frequency = 5e6Hz gpsl1 = GPSL1CA() @@ -189,25 +255,11 @@ end sat2_initial = TrackedSat(gpsl1, 2, code_phase, carrier_doppler; doppler_estimator) sat1 = TrackedSat( sat1_initial; - signals = ( - TrackedSignal( - only(sat1_initial.signals); - is_integration_completed = true, - integrated_samples = num_samples, - correlator, - ), - ), + signals = (_completed_signal(only(sat1_initial.signals), correlator, num_samples),), ) sat2_pre = TrackedSat( sat2_initial; - signals = ( - TrackedSignal( - only(sat2_initial.signals); - is_integration_completed = true, - integrated_samples = num_samples, - correlator, - ), - ), + signals = (_completed_signal(only(sat2_initial.signals), correlator, num_samples),), ) # Bump sat2's bandwidths by replacing its doppler estimator state with a # custom-configured SatConventionalPLLAndDLL. diff --git a/test/doppler_update_interval.jl b/test/doppler_update_interval.jl new file mode 100644 index 00000000..09b9f115 --- /dev/null +++ b/test/doppler_update_interval.jl @@ -0,0 +1,151 @@ +module UpdateIntervalTest + +using Test: @test, @testset, @test_throws +using Unitful: Hz, s, uconvert, NoUnits +import Tracking +using GNSSSignals: + GPSL1CA, + GalileoE1B, + gen_code, + get_code_frequency, + get_code_length, + get_code_center_frequency_ratio +using Tracking: + TrackedSat, + TrackState, + BandMeasurement, + track, + track!, + CPUDownconvertAndCorrelator, + downconvert_and_correlate!, + get_correlator_outputs, + get_sat_state, + get_carrier_doppler, + get_code_doppler, + CorrelatorOutput + +@testset "smallest-code-period default resolution" begin + # `nothing` doppler_update_interval => smallest primary-code period across all + # signals in all groups. GPS L1 C/A: 1023 chips / 1.023 MHz = 1 ms. + ts = TrackState(; signal = GPSL1CA()) + dt = Tracking._smallest_code_period(ts) + @test uconvert(NoUnits, dt * get_code_frequency(GPSL1CA())) ≈ get_code_length(GPSL1CA()) + @test dt ≈ 1e-3 * s + + # Mixed L1 C/A (1 ms) + Galileo E1B (4 ms): the smaller wins. + ts_mixed = TrackState(; signals = (l1ca = (GPSL1CA(),), e1b = (GalileoE1B(),))) + @test Tracking._smallest_code_period(ts_mixed) ≈ 1e-3 * s +end + +@testset "chunk-grid arithmetic (`_chunk_last_sample`)" begin + # Boundaries lie on a shared time grid, re-anchored to the absolute chunk + # index so different bands stay time-aligned without cumulative drift. + @test Tracking._chunk_last_sample(1e-3 * s, 0, 5e6Hz, 100_000) == 5000 + @test Tracking._chunk_last_sample(1e-3 * s, 1, 5e6Hz, 100_000) == 10_000 + @test Tracking._chunk_last_sample(1e-3 * s, 9, 5e6Hz, 100_000) == 50_000 + # A second band at a different sampling rate, same 1 ms time grid. + @test Tracking._chunk_last_sample(1e-3 * s, 0, 13e6Hz, 200_000) == 13_000 + @test Tracking._chunk_last_sample(1e-3 * s, 4, 13e6Hz, 200_000) == 65_000 + # Clamps to the buffer end on the final chunk. + @test Tracking._chunk_last_sample(1e-3 * s, 20, 5e6Hz, 12_345) == 12_345 + # `nothing` => no chunking, whole buffer at once. + @test Tracking._chunk_last_sample(nothing, 0, 5e6Hz, 12_345) == 12_345 +end + +@testset "doppler_update_interval shorter than a sample throws" begin + gpsl1 = GPSL1CA() + ts = TrackState(gpsl1, [TrackedSat(gpsl1, 1, 0.0, 100.0Hz)]) + signal = zeros(ComplexF32, 4000) + @test_throws ArgumentError track(signal, ts, 4e6Hz; doppler_update_interval = 1e-9 * s) +end + +@testset "doppler_update_interval without time units throws a clear error" begin + gpsl1 = GPSL1CA() + ts = TrackState(gpsl1, [TrackedSat(gpsl1, 1, 0.0, 100.0Hz)]) + signal = zeros(ComplexF32, 4000) + @test_throws ArgumentError track(signal, ts, 4e6Hz; doppler_update_interval = 1e-3) + err = try + track(signal, ts, 4e6Hz; doppler_update_interval = 1e-3) + catch e + e + end + @test occursin("time quantity", err.msg) +end + +@testset "a chunk collects every completed correlator output, sample-indexed" begin + gpsl1 = GPSL1CA() + fs = 5e6Hz + prn = 1 + start_code_phase = 0.0 + n_periods = 3 + period = 5000 # 1 ms at 5 MHz = one L1 C/A period + num_samples = n_periods * period + signal = ComplexF32.( + gen_code(num_samples, gpsl1, prn, fs, get_code_frequency(gpsl1), start_code_phase), + ) + ts = TrackState(gpsl1, [TrackedSat(gpsl1, prn, start_code_phase, 0.0Hz)]) + dc = CPUDownconvertAndCorrelator() + measurements = (L1 = BandMeasurement(signal, fs),) + + # One chunk spanning all three periods, NCO held fixed, estimator not run, + # so the collected outputs survive for inspection. + downconvert_and_correlate!( + dc, + measurements, + ts; + chunk_index = 0, + chunk_duration = n_periods * 1e-3 * s, + ) + outputs = get_correlator_outputs(only(get_sat_state(ts, prn).signals)) + + @test length(outputs) == n_periods + @test all(o -> o isa CorrelatorOutput, outputs) + for (i, o) in enumerate(outputs) + # Each integration ends on its code-period boundary; the recorded end + # sample index and integrated-sample count reflect that. + @test isapprox(o.sample_index, i * period; atol = 2) + @test isapprox(o.integrated_samples, period; atol = 2) + end +end + +@testset "estimator empties correlator_outputs after each chunk" begin + gpsl1 = GPSL1CA() + fs = 5e6Hz + prn = 1 + num_samples = 15_000 + signal = + ComplexF32.(gen_code(num_samples, gpsl1, prn, fs, get_code_frequency(gpsl1), 0.0)) + ts = TrackState(gpsl1, [TrackedSat(gpsl1, prn, 0.0, 0.0Hz)]) + track!(signal, ts, fs) + @test isempty(get_correlator_outputs(only(get_sat_state(ts, prn).signals))) +end + +@testset "enlarged doppler_update_interval still converges" begin + # A long continuous buffer with a constant true Doppler; the initial + # estimate is offset. Both the default (per-code-period) and a 4 ms + # (4-periods-per-chunk) update interval must pull the carrier Doppler in. + gpsl1 = GPSL1CA() + fs = 4e6Hz + prn = 1 + true_doppler = 300.0Hz + start_code_phase = 0.0 + code_frequency = + true_doppler * get_code_center_frequency_ratio(gpsl1) + get_code_frequency(gpsl1) + num_samples = 400_000 # 100 ms + range = 0:(num_samples-1) + signal = ComplexF32.( + cis.(2π .* true_doppler .* range ./ fs) .* + gen_code(num_samples, gpsl1, prn, fs, code_frequency, start_code_phase), + ) + + converged(doppler_update_interval) = begin + ts = TrackState(gpsl1, [TrackedSat(gpsl1, prn, start_code_phase, true_doppler - 30Hz)]) + ts = track(signal, ts, fs; doppler_update_interval) + abs(get_carrier_doppler(ts) / Hz - true_doppler / Hz) + end + + @test converged(nothing) < 5 # default = 1 ms + @test converged(4e-3 * s) < 5 # 4 ms => 4 correlator outputs per chunk +end + +end diff --git a/test/downconvert_and_correlate.jl b/test/downconvert_and_correlate.jl index 9bf35808..2340f0da 100644 --- a/test/downconvert_and_correlate.jl +++ b/test/downconvert_and_correlate.jl @@ -17,6 +17,8 @@ using Tracking: downconvert_and_correlate, get_accumulators, get_correlator, + get_correlator_outputs, + get_sat_state, get_correlator_sample_shifts, update_accumulator @@ -38,9 +40,17 @@ end num_samples_signal = 5000 intermediate_frequency = 0.0Hz - sats = [TrackedSat(gpsl1, 1, code_phase, 1000.0Hz), TrackedSat(gpsl1, 2, 11.0, 500.0Hz)] downconvert_and_correlator = DC() - track_state = TrackState(gpsl1, sats) + + # A bare `downconvert_and_correlate` (no `doppler_update_interval`) treats the whole + # buffer as one chunk. Here one code period (~4949 samples) completes, so it + # is snapshotted into `correlator_outputs`; the small residue stays in the + # live accumulator. The completed integration's correlator equals the old + # single-step value — only its location moved (from `get_correlator` to the + # recorded output). Each check uses its own single-sat TrackState so the + # (shared, reused) `correlator_outputs` buffer is not carried between calls. + only_output(track_state, prn) = + only(get_correlator_outputs(only(get_sat_state(track_state, prn).signals))) signal = gen_code( @@ -52,6 +62,7 @@ end code_phase, ) .* cis.(2π * (0:(num_samples_signal-1)) * 1000.0Hz / sampling_frequency) + track_state = TrackState(gpsl1, [TrackedSat(gpsl1, 1, code_phase, 1000.0Hz)]) measurements = (L1 = BandMeasurement(signal, sampling_frequency, intermediate_frequency),) next_track_state = @inferred downconvert_and_correlate( @@ -60,7 +71,8 @@ end track_state, ) - @test real.(get_correlator(next_track_state, 1).accumulators) ≈ [2921, 4949, 2917] rtol=1e-3 + @test real.(only_output(next_track_state, 1).correlator.accumulators) ≈ + [2921, 4949, 2917] rtol = 1e-3 signal = gen_code( @@ -72,6 +84,7 @@ end 11.0, ) .* cis.(2π * (0:(num_samples_signal-1)) * 500.0Hz / sampling_frequency) + track_state = TrackState(gpsl1, [TrackedSat(gpsl1, 2, 11.0, 500.0Hz)]) measurements = (L1 = BandMeasurement(signal, sampling_frequency, intermediate_frequency),) next_track_state = @inferred downconvert_and_correlate( @@ -80,7 +93,8 @@ end track_state, ) - @test real.(get_correlator(next_track_state, 2).accumulators) ≈ [2919, 4947, 2915] rtol=1e-3 + @test real.(only_output(next_track_state, 2).correlator.accumulators) ≈ + [2919, 4947, 2915] rtol = 1e-3 end @testset "Early return when signal_start_sample past end with $DC" for DC in [ diff --git a/test/downconvert_and_correlate_int16.jl b/test/downconvert_and_correlate_int16.jl index 748e857b..e6415c50 100644 --- a/test/downconvert_and_correlate_int16.jl +++ b/test/downconvert_and_correlate_int16.jl @@ -84,7 +84,19 @@ function correlate_once(dc, sig, fs, nsamp, cdopp, cphase; correlator = nothing) TrackedSat(sig, 1, cphase, cdopp; correlator) ts = TrackState(sig, [sat]) ts2 = downconvert_and_correlate(dc, meas, ts) - first(get_sat_state(ts2, 1).signals).correlator + _completed_or_partial_correlator(first(get_sat_state(ts2, 1).signals)) +end + +# A bare `downconvert_and_correlate` treats the whole buffer as one chunk. If a +# code period completed, its (raw) correlator was snapshotted into +# `correlator_outputs` and the live correlator holds only the residue — return +# the first completed integration. If the buffer was shorter than one code +# period (e.g. Galileo E1B's 4 ms period in a 1 ms buffer) nothing completed, so +# the live correlator holds the whole partial integration — return that. Either +# way this matches the value the old single-step call left in the live correlator. +function _completed_or_partial_correlator(sig) + outs = sig.correlator_outputs + isempty(outs) ? sig.correlator : first(outs).correlator end @testset "Int16 downconvert and correlate" begin @@ -180,8 +192,9 @@ end cap = make_capture_mat(sig, fs, nsamp, 200Hz, 100.0, M) meas = (L1 = BandMeasurement(cap, fs, 0.0Hz),) mk() = TrackState(sig, [TrackedSat(sig, 1, 100.0, 200Hz; num_ants = NumAnts(M))]) - run(dc) = - first(get_sat_state(downconvert_and_correlate(dc, meas, mk()), 1).signals).correlator + run(dc) = _completed_or_partial_correlator( + first(get_sat_state(downconvert_and_correlate(dc, meas, mk()), 1).signals), + ) cf = run(CPUThreadedDownconvertAndCorrelator()) ci = run(Int16ThreadedDownconvertAndCorrelator(TEST_MAX_MEAS)) ef, pf, lf = get_early(cf), get_prompt(cf), get_late(cf) # SVector{M,Complex} diff --git a/test/downconvert_and_correlate_onebit.jl b/test/downconvert_and_correlate_onebit.jl index 68117368..2c668d38 100644 --- a/test/downconvert_and_correlate_onebit.jl +++ b/test/downconvert_and_correlate_onebit.jl @@ -89,7 +89,19 @@ function correlate_once( TrackedSat(sig, 1, cphase, cdopp; correlator) ts = TrackState(sig, [sat]) ts2 = downconvert_and_correlate(dc, meas, ts) - first(get_sat_state(ts2, 1).signals).correlator + _completed_or_partial_correlator(first(get_sat_state(ts2, 1).signals)) +end + +# A bare `downconvert_and_correlate` treats the whole buffer as one chunk. If a +# code period completed, its (raw) correlator was snapshotted into +# `correlator_outputs` and the live correlator holds only the residue — return +# the first completed integration. If the buffer was shorter than one code +# period (e.g. Galileo E1B's 4 ms period in a 1 ms buffer) nothing completed, so +# the live correlator holds the whole partial integration — return that. Either +# way this matches the value the old single-step call left in the live correlator. +function _completed_or_partial_correlator(sig) + outs = sig.correlator_outputs + isempty(outs) ? sig.correlator : first(outs).correlator end # Track a noisy GPS L1CA capture at a known C/N0 with backend `dc`, returning the @@ -367,10 +379,11 @@ _std(x) = (m = _mean(x); sqrt(sum(v -> abs2(v - m), x) / (length(x) - 1))) TrackState(sig, satN; doppler_estimator = est), ) for s in get_sat_state(tsN, 1).signals - @test length(get_prompt(s.correlator)) == M - @test get_prompt(s.correlator) == get_prompt(cs) - @test get_early(s.correlator) == get_early(cs) - @test get_late(s.correlator) == get_late(cs) + c = _completed_or_partial_correlator(s) + @test length(get_prompt(c)) == M + @test get_prompt(c) == get_prompt(cs) + @test get_early(c) == get_early(cs) + @test get_late(c) == get_late(cs) end end diff --git a/test/downconvert_and_correlate_twobit.jl b/test/downconvert_and_correlate_twobit.jl index e79fa785..d3ee4699 100644 --- a/test/downconvert_and_correlate_twobit.jl +++ b/test/downconvert_and_correlate_twobit.jl @@ -91,7 +91,19 @@ function correlate_once( TrackedSat(sig, 1, cphase, cdopp; correlator) ts = TrackState(sig, [sat]) ts2 = downconvert_and_correlate(dc, meas, ts) - first(get_sat_state(ts2, 1).signals).correlator + _completed_or_partial_correlator(first(get_sat_state(ts2, 1).signals)) +end + +# A bare `downconvert_and_correlate` treats the whole buffer as one chunk. If a +# code period completed, its (raw) correlator was snapshotted into +# `correlator_outputs` and the live correlator holds only the residue — return +# the first completed integration. If the buffer was shorter than one code +# period (e.g. Galileo E1B's 4 ms period in a 1 ms buffer) nothing completed, so +# the live correlator holds the whole partial integration — return that. Either +# way this matches the value the old single-step call left in the live correlator. +function _completed_or_partial_correlator(sig) + outs = sig.correlator_outputs + isempty(outs) ? sig.correlator : first(outs).correlator end # Direct (per-sample) reference of the two-bit quantised correlation: quantise the @@ -540,10 +552,11 @@ _std(x) = (m = _mean(x); sqrt(sum(v -> abs2(v - m), x) / (length(x) - 1))) TrackState(sig, satN; doppler_estimator = est), ) for s in get_sat_state(tsN, 1).signals - @test length(get_prompt(s.correlator)) == M - @test get_prompt(s.correlator) == get_prompt(cs) - @test get_early(s.correlator) == get_early(cs) - @test get_late(s.correlator) == get_late(cs) + c = _completed_or_partial_correlator(s) + @test length(get_prompt(c)) == M + @test get_prompt(c) == get_prompt(cs) + @test get_early(c) == get_early(cs) + @test get_late(c) == get_late(cs) end end @@ -736,11 +749,18 @@ _std(x) = (m = _mean(x); sqrt(sum(v -> abs2(v - m), x) / (length(x) - 1))) round.(Int16, A .* imag.(s) .+ σ .* randn(N)), ) meas = (L1 = BandMeasurement(cap, fs, 0.0Hz),) - ts = TrackState(sig, [TrackedSat(sig, prn, cphase, cdopp)]) - c1 = - first(get_sat_state(downconvert_and_correlate(d1, meas, ts), 1).signals).correlator - c2 = - first(get_sat_state(downconvert_and_correlate(d2, meas, ts), 1).signals).correlator + # Separate track states per backend: `downconvert_and_correlate` + # shares each signal's (reused) `correlator_outputs` buffer with its + # input, so running both backends on one `ts` would append d2's + # records after d1's and make `first(...)` return d1's for both. + ts1 = TrackState(sig, [TrackedSat(sig, prn, cphase, cdopp)]) + ts2 = TrackState(sig, [TrackedSat(sig, prn, cphase, cdopp)]) + c1 = _completed_or_partial_correlator( + first(get_sat_state(downconvert_and_correlate(d1, meas, ts1), 1).signals), + ) + c2 = _completed_or_partial_correlator( + first(get_sat_state(downconvert_and_correlate(d2, meas, ts2), 1).signals), + ) push!(err1, angle(get_prompt(c1)) - ref) push!(err2, angle(get_prompt(c2)) - ref) end diff --git a/test/runtests.jl b/test/runtests.jl index dce265ad..d1ae4c20 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,6 +25,7 @@ include("cn0_estimation.jl") include("tracking_state.jl") include("track.jl") include("track_in_place.jl") +include("doppler_update_interval.jl") include("multi_signal.jl") include("carrier_phase.jl") include("multi_band.jl") diff --git a/test/sat_state.jl b/test/sat_state.jl index 3361770b..22c774b8 100644 --- a/test/sat_state.jl +++ b/test/sat_state.jl @@ -173,7 +173,6 @@ GNSSSignals.get_data_frequency(::FakeWrapSignal) = 0Hz fake_tracked_signal(code_length, secondary_length, found) = Tracking.TrackedSignal( FakeWrapSignal(code_length, secondary_length), base.integrated_samples, - base.is_integration_completed, base.correlator, base.last_fully_integrated_correlator, base.last_fully_integrated_filtered_prompt, @@ -181,6 +180,7 @@ GNSSSignals.get_data_frequency(::FakeWrapSignal) = 0Hz found ? synced_buffer(base.bit_buffer) : base.bit_buffer, base.post_corr_filter, base.filtered_prompts, + base.correlator_outputs, base.preferred_num_code_blocks_to_integrate, ) diff --git a/test/track.jl b/test/track.jl index 6ed0e462..8c239ee5 100644 --- a/test/track.jl +++ b/test/track.jl @@ -206,7 +206,11 @@ end @test test_convergence(90Hz, false) == true @test test_convergence(100Hz, false) == false - # ConventionalAssistedPLLAndDLL: converges at 240Hz offset, fails at 250Hz + # ConventionalAssistedPLLAndDLL: converges at 240Hz offset, fails at 250Hz. + # The two-pass chunk (completions → NCO update → residue at the updated + # Doppler) keeps every integration on a single Doppler and applies each + # correction at its completing boundary, so the FLL pull-in edge matches + # the classic per-code-period update. @test test_convergence(240Hz, true) == true @test test_convergence(250Hz, true) == false end