Noise estimation

Abstract supertype for per-signal noise estimators — the source of the noise density N₀ that NoiseRefCN0Estimator divides each record's prompt power by.

One instance is held per signal in Tracking.TrackState's noise_estimators NamedTuple (keyed by GNSSSignals.get_signal_id), never per satellite: every satellite of a signal shares one floor, and averaging it once per signal is what makes the reference's own variance negligible against the per-record prompt statistics.

Why per signal and not per RF band

What a record actually divides by is the post-correlation density. For an interferer of PSD S_I(f) on top of the thermal floor that is

N₀,eff = N₀,thermal + ∫ S_I(f) · |G(f)|² df

— the spectral separation coefficient, weighted by the despreading modulation's own spectrum |G(f)|². Two signals sharing one band, one antenna and one front end therefore see different floors the moment the interference is not white, and the difference is not small: BPSK(1) has its peak at DC and a null at ±1.023 MHz, BOC(1,1) is the reverse, so a CW tone at band centre is rejected by Galileo E1B and lands squarely in GPS L1 C/A, while a tone at ±1.023 MHz does the opposite. Front-end tilt, filter roll-off at the band edge and adjacent-band leakage all colour the floor the same way, more quietly.

Because CorrelatorNoiseEstimator despreads rather than metering power, keying by signal makes that integral measured rather than modelled: the reference runs the consumer's own code, so its spectral weighting is the consumer's by construction. This is the same argument that makes the reference backend-free — it traverses the identical path as the prompt — extended from the quantiser to the interference environment.

It also matches the hardware. A noise channel is a tracking channel with a wrong PRN, and a tracking channel is configured with a code; per-signal is what an FPGA would build anyway.

Why a density and not a power

For a correlator normalised the way TrackingLoops.normalize normalises — by integrated_samples * code_amplitude — white input noise of per-sample variance σ² gives E|P|² = σ²/N = N₀/T. So N₀ = σ²/f_s is independent of the integration time, which is what lets one per-signal figure serve records of any length. It is stored as a Unitful quantity of dimension 1/Hz, so the consumer's ⟨|P|²⟩/N₀ − 1/T is dimension-checked rather than trusted.

The interface

Three required methods, all with a default on this abstract type:

• update_noise! — measure one signal's slice of its band's samples and append the resulting observations. This is the software fill path; it has exactly one call site, inside downconvert_and_correlate!.
• append_noise_observation! — append an observation built elsewhere. This is the hardware fill path (FPGA/ASIC correlator or a front-end power monitor), parallel to Tracking.append_correlator_output!.
• get_noise_density — the signal's current density, or nothing while nothing has been measured yet. A read, not a drain: the window keeps sliding.

A fourth, TrackingLoops.noise_window_looks, is optional and only matters to a source that reports a covariance: an M×M estimate averaged from fewer than M independent looks is rank-deficient by construction, so the fold withholds it until there are enough. A source that leaves the default in place reports no look count and is never withheld.

The two fill paths live on disjoint call graphs — a hardware producer never calls downconvert_and_correlate! — so one concrete type serves both and no type parameter distinguishes them. See CorrelatorNoiseEstimator, the only shipped implementation.

A source that must not measure from samples (a pure power-monitor integration, say) subtypes this and gives update_noise! a no-op method; that is why the interface is exported.

NoiseDensity

The type of a single-antenna noise density N₀, typeof(1.0 / 1.0Hz): what every shipped observation builder produces, and what get_noise_density returns for a single-antenna estimator. A multi-antenna window measures a spatial covariance matrix instead; see noise_density_type.

Type alias for a NamedTuple of AbstractNoiseEstimators keyed by signal id — the shape Tracking.TrackState holds them in. A NamedTuple and not a Dictionary, because a dictionary would need an abstract value type as soon as two signals hold different estimator types, which is type-unstable and would allocate on every chunk.

One noise measurement, the producer-side value parallel to CorrelatorOutput. It carries the already-normalised density plus the two pieces of bookkeeping a sliding window needs, so the internal contract is one scalar and every unit conversion happens at the boundary where the caller knows their hardware.

Build one with noise_observation, noise_observation_from_correlator or noise_observation_from_samples rather than by hand — they are the three shapes a producer actually has, and all three reduce to the same N₀ on the same white input.

Fields:

• noise_density — N₀, of dimension 1/Hz. For an antenna array it is the M×M spatial covariance R̂ instead, of the same dimension elementwise: the diagonal is each antenna's own N₀ and the off-diagonals their noise correlation. A satellite reduces it to its own scalar floor through its beamforming weights, wᴴR̂w (see AbstractPostCorrFilter).
• num_sub_integrations — M, the number of independent looks that went into it, and therefore its statistical weight. Not the sample count: for Σ_m |B_m|² over M sub-integrations of N samples each, the density needs only M·N, but its relative variance is 1/M. A producer handing back one 1 ms coherent despread and one handing back 64 16-chip despreads report the same total samples and wildly different precision, so the window is averaged weighted by M.
• duration — the wall-clock span of samples the observation covers, which is what bounds the sliding window. It is not M · N / f_s in general: the software source's M looks are the correlator's taps, which are simultaneous rather than consecutive.
• prn — the PRN whose code was used to measure it. Legitimate metadata (the observation really was measured with that code), and it is what carries the software source's PRN-rotation position without any scalar state.

append_noise_observation!(estimator, _)

Append one externally built NoiseObservation to estimator's sliding window and return estimator. This is the hardware fill path, parallel to Tracking.append_correlator_output! — see there for how the two differ.

The window is mutated in place and the struct is not rebuilt, so this is allocation-free in steady state and works through the immutable TrackState.

The Tracking.TrackState form selects the signal:

append_noise_observation!(track_state, obs)             # single-signal TrackState
append_noise_observation!(track_state, obs, :GPSL1CA)   # explicit signal id
append_noise_observation!(track_state, obs, GPSL1CA)    # or the signal type

get_noise_density(_)

The signal's current noise density N₀, or nothing while the window holds nothing yet.

A read, not a drain: the window keeps sliding across chunks and across track! calls, so every record of every satellite tracking that signal divides by the same figure. The returned quantity has dimension 1/Hz — see AbstractNoiseEstimator for why a density rather than a power, and why per signal rather than per band.

noise_density_type(_)

The concrete type get_noise_density returns for estimator once its window holds anything. Used to keep the fold's density argument a plain scalar — never a Union{Nothing,…} — so update stays devirtualised and allocation-free even while the window is still empty.

Defaults to typeof(1.0/1.0Hz), which every shipped builder produces.

noise_observation(
    accumulation,
    integrated_samples,
    sampling_frequency;
    code_amplitude,
    prn
)

Build a NoiseObservation from one dump given as a raw complex accumulation Σ x·c over integrated_samples samples — the most faithful form a producer can report, since Tracking does the squaring, the scaling and all the averaging itself.

code_amplitude is the RMS amplitude of the sampled replica (see TrackingLoops.normalize); leave it at 1 for a ±1 code, pass the code's RMS for a multi-level one (CBOC). prn records which code measured it.

Pass an SVector of M per-antenna accumulations for an antenna array, and the observation carries the spatial covariance b·bᴴ instead of the scalar |b|². Report every element the front end has: a satellite's floor is wᴴR̂w under its own beamforming weights, and an array collapsed to one number cannot answer that.

noise_observation_from_correlator(
    accumulated_power,
    num_sub_integrations,
    total_samples,
    sampling_frequency;
    code_amplitude,
    prn,
    duration
)

Build a NoiseObservation from M = num_sub_integrations dumps pre-summed on chip as Σ_m |B_m|², covering total_samples = M · N samples in all.

Pre-summing must be incoherent (Σ|B_m|², never |Σ B_m|²), which is why this builder takes a power rather than a complex value. M > 1 exists only so a high-rate producer can cut host traffic; the natural granularity is one observation per dump with M = 1, which noise_observation covers.

duration defaults to total_samples / sampling_frequency, which is right when the M sub-integrations are consecutive in time. Pass it explicitly when they are simultaneous — the software source's M looks are the reference correlator's taps, which all integrate the same N samples, so their span is N / sampling_frequency and not M times that.

For an antenna array, pass the pre-summed spatial covariance Σ_m B_m·B_mᴴ as an SMatrix. The incoherence requirement is the same one term for term: sum the outer products, never form the outer product of the sum.

noise_observation_from_samples(
    accumulated_power,
    num_samples,
    sampling_frequency;
    prn
)

Build a NoiseObservation from a front-end / AGC power monitor: Σ|x|² over num_samples raw samples.

This is noise_observation_from_correlator at a sub-integration of one sample — the replica degenerates to a single ±1 value, |x·c|² = |x|², and the window mean is exactly σ̂² — so M == num_samples and no code_amplitude applies. It is the minimum-variance, flat-weighting end of the same continuum the despread source sits at the spectral-fidelity end of; a tilted front end is the one thing that tells them apart.

For an antenna array, pass Σ x·xᴴ as an SMatrix — the array's raw spatial covariance, which is what a beamformer's weights are reduced against.

noise_window_looks(_)

How many independent looks the estimator's current density is averaged from, or nothing if the source does not report it.

Only the rank gate below reads this: an M×M covariance averaged from fewer than M looks is rank-deficient by construction, whatever the samples were, so it cannot answer wᴴR̂w for every w. A source that cannot report a look count is taken at its word and not gated.

update_noise!(
    estimator,
    measurement,
    first_sample,
    last_sample,
    context
)

Measure one signal's noise on its band's samples and append the resulting observations to estimator's window, returning estimator.

measurement is the band's Tracking.BandMeasurement — the samples are a band property, one front end feeding every signal on it; only the despreading code, and therefore the measured floor, is per signal. first_sample and last_sample bound the slice of it this call may consume (the current chunk). context is a NoiseUpdateContext, carrying what a source may need but a BandMeasurement does not have — the signal being measured, which PRNs are currently tracked on it, and the chunk index.

This is the software fill path and it has exactly one call site, inside downconvert_and_correlate!. A source whose observations arrive from outside (see append_noise_observation!) never reaches it, so the default here is a no-op: it mutates the window in place and leaves the struct alone, which is what lets TrackState hold it immutably.

The signal's noise reference, measured by despreading an untracked PRN — the only AbstractNoiseEstimator Tracking ships, and the one to configure on a hardware path too (you simply fill it with append_noise_observation! instead of letting update_noise! fill it).

Why a despread and not a power meter

The reference traverses the identical quantise → downconvert → despread → accumulate path as the prompt, reusing the same kernel, the same replica generator and — because the estimator is keyed by signal — the same code. So the measured N₀ already contains every imperfection of that chain — the 1-bit / 2-bit quantisation loss, the quantiser's operating point under load, the input scaling an AGC step moves, the code amplitude of a CBOC replica — with no per-backend model and no closed-form correction. A Σ|x|² power meter would need one line per backend and would still miss the second-order coupling where a strong signal shifts that operating point.

It is also what makes the spectral weighting right rather than assumed. Sharing the consumer's code means the despread evaluates that signal's own ∫ S_I(f)·|G(f)|² df — the coloured-interference case a power meter, or a reference despread with some other signal's code, both get wrong. See AbstractNoiseEstimator.

It is open-loop

The reference has no feedback of any kind: it runs at the band's nominal IF, at a randomly dithered offset from zero Doppler, with a random code phase and a rotating PRN. There is no discriminator, no loop filter and no NCO update ever written to it — and, since the randomisation removed the need to know which PRNs are in use, it reads nothing from satellite state, not even the tracked key set. It is therefore correct with zero satellites tracked, and on an FPGA a noise channel is a strict subset of a tracking channel rather than an addition to one.

Why the code phase and the Doppler are randomised

A reference at a fixed code phase and exactly zero Doppler is a stationary target, and the failure that matters is not an attacker — it is that a hit never goes away. The replica is re-anchored to code phase 0 on a grid running at the nominal chip rate, so the relative phase against any incoming signal at zero Doppler is frozen: whatever it lands on, it stays. A signal that is present but untracked — a spoofed PRN, or simply a visible satellite the receiver has not acquired — has a 3 × 1.5 / 1023 ≈ 0.44 % chance per PRN of sitting inside one of the three taps, and if it does it contributes T·(C/N₀)/3 ≈ 10.5·N₀ at 45 dB-Hz to every observation on that PRN, indefinitely (≈+1.5 dB on N̂₀ after the rotation's dilution). Worse, the geometry is perverse: relative phase drifts only at the signal's own code Doppler (f_d/1540), so low Doppler is both what makes a hit possible and what makes it permanent.

Drawing a fresh code phase and a fresh carrier offset per sub-integration turns that standing bias into an independent per-observation trial. For a full sky at 45 dB-Hz the residual is ≈0.07 dB — a handful of 1 %-sized outliers scattered through a 1000-entry window, rather than a permanent shift of the floor — and, because a hit now needs both draws to land, the two randomisations multiply.

It also makes the untracked-PRN restriction unnecessary, which is why the reference now rotates over all PRNs of the family: a random phase lands within ±1 chip of a tracked satellite's peak with probability ≈0.6 %, worth 10.5/1000 ≈ 0.045 dB for the one observation it touches — an order of magnitude under the whole-sky leakage the estimator already carries uncorrected. A constant 32-PRN pool also stops the rotation from shortening (and the dilution from worsening) exactly as the receiver acquires more satellites.

Neither draw biases the measurement. N̂₀ = |B|²/(N·A_c²·f_s) is unbiased for any starting phase, and for the carrier the reference measures the same noise power wherever it sits; carrier_dither of ±5 kHz smears the spectral weighting ∫ S_I(f)·|G(f)|² df by 0.25 % of a 2 MHz main lobe, which no coloured interferer resolves. On an FPGA an arbitrary code phase is easier than phase 0 — a free-running code generator gives you one, where phase 0 needs a reset.

Fields / configuration

• window_duration — how far back the sliding window reaches (1 s by default). The one genuine tunable, because it trades: longer means lower variance, but a longer smear across AGC changes. At the default and one sub-integration per code period the window holds K_n ≈ 1000 looks per tap, which costs ≤0.08 dB against a variance-free reference at every C/N₀ — below which no further tuning is warranted.

• tap_code_shift — the reference correlator's tap spacing in chips (1.5 by default). Taps are only worth having if they are statistically independent, and at the tracking default of ±0.5 chip they are not: for jointly circular Gaussian accumulations corr(|Bᵢ|²,|Bⱼ|²) = |ρᵢⱼ|², and the sampled-code correlation at half a chip is ≈0.5, so three taps are worth 2.25 independent looks. At ≥1 chip they are worth 2.98. 1.5 chips sits in the autocorrelation null and clear of both the 1- and 2-chip sidelobe values.

• carrier_dither — half-width of the uniform offset added to the band's nominal IF, per sub-integration (5 kHz by default, i.e. the terrestrial GNSS Doppler spread). Zero pins the reference at the IF exactly, which restores half of the stationary target described above — useful to isolate the code-phase draw in a test, and not otherwise.

• rng — the source of the code-phase and carrier draws, seeded by default (Xoshiro(0), one stream per estimator, advanced in place like buffered so the struct is never rebuilt). What the randomisation has to defeat is a stationary reference, not a reader: an attacker cannot observe the chunk grid the phase would be measured against, so scattering the draws is the whole requirement and unpredictability buys nothing on top. A seeded default keeps a run repeatable, which is worth having in a library whose other guarantees are phrased as bit-identical arithmetic. Pass rng = Random.default_rng() for a task-local stream (also the one to use if you ever share an estimator across threads).

  Repeatable on one Julia version, and no further: Xoshiro's stream is not part of Julia's compatibility guarantee and does change across releases. Do not build a tolerance around a particular draw — size the window so the assertion holds for any draw. Every N̂₀ here is a 1/√(3K) estimate over K observations, and that is the number to design against.

• buffered — the sliding window itself, a length-managed FIFO written in place. The Vector's own length is the position, so there is no ring index to write back and the struct is never rebuilt — which is what lets per-signal state live in an immutable Tracking.TrackState.

• totals — the window's running sums (span, M-weighted density, looks), maintained as entries are pushed and dropped. They are what keeps both append_noise_observation! and get_noise_density O(1); recomputing them per call made each an O(K) scan, and since the whole design buys its variance by making K large (≈2500 entries at a 1 s window and a 0.4 ms chunk), that put a per-chunk cost directly proportional to the accuracy asked for. Written in place through Ref cells, exactly as buffered is, so the "never rebuilt" property still holds — and a cache of buffered rather than state of its own, recomputed from it exactly often enough that the incremental arithmetic cannot drift.

The sub-integration length is deliberately not a field: it is the primary code period of the signal it measures, derived rather than configured. Coherent integration buys nothing for noise-power estimation — one dump carries 100 % relative error however long it is, because Var(|B|²) = (E|B|²)² — so all the information lives in the number of looks, and shortening the dump is the only way to buy more of them. Three things say not to: SIMD (64 kernel calls per 1 ms chunk instead of one, and a shorter dump would need a new kernel variant, forfeiting the bit-identical-arithmetic property the whole approach rests on), DC balance (a full-period Gold code is balanced, E[(Σc)²] ≈ 1; any sub-period window behaves like iid signs, ≈16 at 16 chips, so an ADC offset biases a short despread and not a full-period one), and cadence parity with a hardware correlator, which naturally dumps on the code epoch. The variance is bought back with window_duration instead, which is nearly free — a 1 s window is ~1000 Float64 densities per signal.

CorrelatorNoiseEstimator(
;
    window_duration,
    tap_code_shift,
    carrier_dither,
    num_ants,
    rng
)

Construct a CorrelatorNoiseEstimator averaging over the last window_duration of observations, with the reference correlator's taps spaced tap_code_shift chips apart and its carrier dithered by up to carrier_dither either side of the band's nominal IF. See the type's docstring for what the parameters buy and why nothing else is configurable.

rng is drawn from for the per-sub-integration code phase and carrier offset. It defaults to a seeded Xoshiro(0), so a run repeats on a given Julia version; pass Random.default_rng() for a task-local stream. See the type's docstring for why scattering the draws — rather than making them unguessable — is the whole requirement, and for why the seed is not something to calibrate against.

The window is sizehint!-ed to four times the number of code-period observations it expects to hold. The headroom is what makes the FIFO's push!/popfirst! pair measure exactly zero bytes: at 1× or 2× Julia periodically shifts the front offset back and reallocates.

num_ants must match the antenna count of the signal group this estimator is keyed to. At NumAnts(1) the window holds scalar densities; above it, the M×M spatial covariance R̂ of the array's noise, which each satellite reduces to its own scalar floor through its own beamforming weights (wᴴR̂w, see update_noise! and AbstractPostCorrFilter). TrackState provisions the right count automatically when noise_estimators is left at nothing, and rejects a mismatch when the estimators are passed explicitly.

append_noise_observation!(estimator, observation)

Append observation to the signal's sliding window, dropping entries off the front while the remainder still spans window_duration. Returns estimator.

The window is bounded in time, not in observation count, which is what lets one producer report 16-chip accumulations and another a single pre-averaged 0.2 s figure under the same configuration and with no scalar state — the span is a property of the FIFO rather than of a configured count.

O(1) per call, amortised, whatever the window holds — see totals in the type's docstring for why that matters here rather than being a micro-optimisation.

Any NoiseObservation is accepted and retyped onto the window's own field types, which costs nothing for one the builders produced (they already emit the canonical pair) and is what keeps a hand-assembled or Float32 one from matching no method here and falling through to the abstract no-op — where it would be dropped silently and leave the window empty forever.

get_noise_density(estimator)

The window's M-weighted mean density, or nothing while it is empty.

Weighted by num_sub_integrations because that is the number of independent looks each entry represents: the density itself needs only the sample count, but its relative variance is 1/M, so a one-dump entry and a 64-dump entry combine correctly only when weighted this way.

Read straight off the window's running totals, so this is O(1). It is called once per signal per chunk from the Doppler-estimator fold, which is why it may not walk the window: doing so charged every chunk for the window length, i.e. for the very K the estimator's accuracy is bought with.

update_noise!(
    estimator,
    measurement,
    first_sample,
    last_sample,
    context
)

Measure this signal's noise over samples first_sample:last_sample of measurement and append the resulting observations, returning estimator.

The measurement itself — despreading an untracked PRN through the caller's own correlator kernel — is the software receiver's, so it lives with the software backends: this method forwards to despread_noise! on context.downconvert_and_correlator, which Tracking.jl implements for its backends. A loop process fills the window through append_noise_observation! instead and never reaches this method.