From e7e219128346626da70432b618ca4a86ef9914c6 Mon Sep 17 00:00:00 2001 From: Soeren Schoenbrod Date: Thu, 7 May 2026 16:54:29 +0200 Subject: [PATCH] fix: resolve GPS L1 week-rollover ambiguity via approximate_year MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The legacy GPS L1 C/A LNAV message (IS-GPS-200, §20.3.3.3) broadcasts only a 10-bit week number — modulo 1024 — so the receiver cannot determine which 1024-week cycle the recording is in from the broadcast data alone. Each cycle is ~19.6 years, with rollover boundaries at 1999-08-22, 2019-04-07, 2038-11-21, and 2058-07-08. The previous implementation hardcoded `2048 + decoder.data.trans_week`, which works only for recordings inside the post-2019 cycle. Recordings from earlier cycles produced PVT timestamps offset by N×19.6 years (observed: ION RTL-SDR sample data from 2017-09-10 reported as 2037-04-26). This change adds an `approximate_year::Integer` keyword to `calc_pvt` (default `year(now(UTC))`) that selects the cycle whose absolute date is closest to the anchor. Anything within ~9 years of the actual observation date works. Live receivers keep working with the default. Post-processing of archived data passes the rough year of the recording; for example: calc_pvt(states; approximate_year = 2017) Pinning is also necessary on test fixtures that bake in expected TAI timestamps — without it those tests would silently break in 2030+ when the wall-clock-derived default starts picking the next cycle. Co-Authored-By: Claude Opus 4.7 (1M context) --- Project.toml | 2 + src/PositionVelocityTime.jl | 58 ++++++++++++++++++++++++---- test/get_week.jl | 75 +++++++++++++++++++++++++++++++++++++ test/pvt.jl | 14 +++++-- test/runtests.jl | 5 ++- 5 files changed, 141 insertions(+), 13 deletions(-) create mode 100644 test/get_week.jl diff --git a/Project.toml b/Project.toml index 01d2a35..e450a6e 100755 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ authors = ["Soeren Schoenbrod ", "Michael [deps] AstroTime = "c61b5328-d09d-5e37-a9a8-0eb41c39009c" CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" GNSSDecoder = "bf6b1ec8-12ad-4501-972c-3dec0d675c2b" GNSSSignals = "52c80523-2a4e-5c38-8979-05588f836870" @@ -22,6 +23,7 @@ Aqua = "0.8" AstroTime = "0.7" BitIntegers = "0.2.6, 0.3.5" CoordinateTransformations = "0.6" +Dates = "1" DocStringExtensions = "0.6, 0.7, 0.8, 0.9" GNSSDecoder = "0.1.5, 0.2, 1" GNSSSignals = "0.17, 1" diff --git a/src/PositionVelocityTime.jl b/src/PositionVelocityTime.jl index 62b63e7..500fce7 100755 --- a/src/PositionVelocityTime.jl +++ b/src/PositionVelocityTime.jl @@ -10,7 +10,8 @@ using CoordinateTransformations, StaticArrays, Tracking, Unitful, - Statistics + Statistics, + Dates using Unitful: s, Hz @@ -196,7 +197,9 @@ function get_sat_enu(user_pos_ecef::ECEF, sat_pos_ecef::ECEF) end """ - calc_pvt(states::AbstractVector{<:SatelliteState}, prev_pvt::PVTSolution=PVTSolution()) -> PVTSolution + calc_pvt(states::AbstractVector{<:SatelliteState}, + prev_pvt::PVTSolution = PVTSolution(); + approximate_year::Integer = year(now(UTC))) -> PVTSolution Calculate Position, Velocity, and Time (PVT) from GNSS satellite measurements. @@ -208,6 +211,15 @@ carrier Doppler measurements. - `states`: Vector of [`SatelliteState`](@ref) for observed satellites - `prev_pvt`: Previous PVT solution used as initial guess (default: origin) +# Keyword Arguments +- `approximate_year`: Calendar year of the observation, used to resolve the + GPS L1 C/A 1024-week rollover ambiguity (legacy LNAV broadcasts only a + 10-bit week number, so the receiver needs external information to + determine which 1024-week cycle the recording is in). Anything within + ±9 years of the actual observation date works. Defaults to the current + UTC year, which is correct for live signals; for processing archived + recordings, pass the rough year of the recording. + # Returns A [`PVTSolution`](@ref) containing position, velocity, time, DOP values, and satellite information. Returns `prev_pvt` if fewer than 4 healthy satellites are @@ -218,7 +230,8 @@ available or if the GDOP is negative. """ function calc_pvt( states::AbstractVector{<:SatelliteState}, - prev_pvt::PVTSolution = PVTSolution(), + prev_pvt::PVTSolution = PVTSolution(); + approximate_year::Integer = year(now(UTC)), ) length(states) < 4 && throw(ArgumentError("You'll need at least 4 satellites to calculate PVT")) @@ -255,7 +268,7 @@ function calc_pvt( # See https://github.com/JuliaGNSS/PositionVelocityTime.jl/issues/8 corrected_reference_time = reference_time - time_correction / SPEEDOFLIGHT - week = get_week(first(healthy_states).decoder) + week = get_week(first(healthy_states).decoder; approximate_year) start_time = get_system_start_time(first(healthy_states).decoder) time = TAIEpoch( week * 7 * 24 * 60 * 60 + floor(Int, corrected_reference_time) + start_time.second, @@ -308,11 +321,42 @@ function get_system_start_time( TAIEpoch(1999, 8, 22, 0, 0, (32 - 13.0)) # There were 32 leap seconds at 08/22/1999 compared to UTC end -function get_week(decoder::GNSSDecoder.GNSSDecoderState{<:GNSSDecoder.GPSL1Data}) - 2048 + decoder.data.trans_week +""" + get_week(decoder::GNSSDecoderState{<:GPSL1Data}; approximate_year) + +Return the absolute GPS week number for a GPSL1 decoder, resolving the +1024-week rollover ambiguity using `approximate_year` as a calendar +anchor. + +The legacy GPS L1 C/A LNAV message broadcasts only a 10-bit week number +(0–1023) modulo 1024, so the receiver cannot determine which 1024-week +cycle the recording is in from the data alone (IS-GPS-200, §20.3.3.3). +Each cycle is ~19.6 years, so any anchor within ±9 years of the true +observation date selects the correct cycle. + +GPS week 0 is 1980-01-06; cycle boundaries fall on 1999-08-22, +2019-04-07, 2038-11-21, 2058-07-08, … + +For Galileo, the broadcast WN is 12 bits and does not need this +treatment in any practical operational scenario. +""" +function get_week( + decoder::GNSSDecoder.GNSSDecoderState{<:GNSSDecoder.GPSL1Data}; + approximate_year::Integer = year(now(UTC)), +) + # GPS week 0 begins 1980-01-06. Compute the integer week count from + # there to mid-`approximate_year`, then choose the cycle base such + # that `cycle_base + trans_week` is closest to that anchor. + days_at_anchor = Date(approximate_year, 6, 30) - Date(1980, 1, 6) + weeks_at_anchor = Dates.value(days_at_anchor) ÷ 7 + n_cycles = round(Int, (weeks_at_anchor - decoder.data.trans_week) / 1024) + return n_cycles * 1024 + decoder.data.trans_week end -function get_week(decoder::GNSSDecoder.GNSSDecoderState{<:GNSSDecoder.GalileoE1BData}) +function get_week( + decoder::GNSSDecoder.GNSSDecoderState{<:GNSSDecoder.GalileoE1BData}; + approximate_year::Integer = year(now(UTC)), +) decoder.data.WN end diff --git a/test/get_week.jl b/test/get_week.jl new file mode 100644 index 0000000..1cba2ac --- /dev/null +++ b/test/get_week.jl @@ -0,0 +1,75 @@ +using PositionVelocityTime: get_week +using GNSSDecoder: GNSSDecoderState, GPSL1Data +using GNSSSignals: GPSL1, GalileoE1B + +# Build a minimal GNSSDecoderState{GPSL1Data} with the given broadcast +# 10-bit `trans_week`. `get_week` only reads `decoder.data.trans_week` +# so the rest of the state can stay at its default zero values. +function decoder_with_trans_week(trans_week) + base = GNSSDecoderState(GPSL1(), 1) + GNSSDecoderState(base; data = GPSL1Data(base.data; trans_week)) +end + +@testset "get_week resolves GPS L1 week-rollover ambiguity" begin + # GPS week 0 begins 1980-01-06. + # Cycle boundaries: WN=1024 → 1999-08-22, WN=2048 → 2019-04-07, + # WN=3072 → 2038-11-21, WN=4096 → 2058-07-08. + + # Cycle 0 (1980-01-06 .. 1999-08-22): broadcast WN ∈ [0, 1023] + # ION RTL-SDR sample data was recorded 2017-09-10 in cycle 1, + # broadcast trans_week ≈ 941 → absolute week 1024 + 941 = 1965. + @testset "anchor in 2017 picks cycle 1 (post-1999 rollover)" begin + decoder = decoder_with_trans_week(941) + @test get_week(decoder; approximate_year = 2017) == 1024 + 941 + # ±9 years of slop still picks the correct cycle. + @test get_week(decoder; approximate_year = 2009) == 1024 + 941 + @test get_week(decoder; approximate_year = 2018) == 1024 + 941 + end + + # Same trans_week but anchor in cycle 2 (post-2019 rollover): + # absolute week 2048 + 941 = 2989 → ~April 2037. The picked cycle is + # whichever places the absolute date closer to the anchor; with + # trans_week = 941 (~mid-cycle), cycle 1 covers 2017 and cycle 2 + # covers 2037, so anchors near 2037 pick cycle 2 and anchors before + # ~2027 still pick cycle 1. + @testset "anchor in 2037 picks cycle 2 (post-2019 rollover)" begin + decoder = decoder_with_trans_week(941) + @test get_week(decoder; approximate_year = 2037) == 2048 + 941 + @test get_week(decoder; approximate_year = 2030) == 2048 + 941 # closer to 2037 than 2017 + end + + # Anchor in cycle 0 (between epoch and first rollover). + @testset "anchor in 1995 picks cycle 0" begin + decoder = decoder_with_trans_week(500) + @test get_week(decoder; approximate_year = 1995) == 0 + 500 + end + + # Anchor in cycle 3 (after the upcoming 2038 rollover). + @testset "anchor in 2050 picks cycle 3" begin + decoder = decoder_with_trans_week(100) + @test get_week(decoder; approximate_year = 2050) == 3072 + 100 + end + + # Boundary case: trans_week = 0 with anchor near the rollover. + # 2019-04-07 (rollover to cycle 2) — anchor at 2019 should pick the + # cycle whose week 0 is closest. Mid-2019 is closer to start-of-2019 + # (cycle 1, end) than start-of-2020 (cycle 2, beginning), so cycle 1 + # is selected. + @testset "trans_week = 0 near rollover" begin + decoder = decoder_with_trans_week(0) + # Mid-2019 anchor ≈ week 2050, distance to cycle 1 base (week 1024) + # is 1026 weeks; distance to cycle 2 base (week 2048) is 2 weeks. + # So cycle 2 wins. + @test get_week(decoder; approximate_year = 2019) == 2048 + end + + # Sanity-check that the default anchor matches an explicit current-UTC-year + # anchor — both should pick the same cycle (the default reads `now()` at + # call time, and the explicit comparison value reads it on the same line, + # so they agree by construction). + @testset "default anchor uses current UTC year" begin + decoder = decoder_with_trans_week(500) + current_year = Dates.year(Dates.now(Dates.UTC)) + @test get_week(decoder) == get_week(decoder; approximate_year = current_year) + end +end diff --git a/test/pvt.jl b/test/pvt.jl index b95c23e..3c1388c 100644 --- a/test/pvt.jl +++ b/test/pvt.jl @@ -490,14 +490,18 @@ BitIntegers.@define_integers 320 ), ] - pvt = calc_pvt(states) + # Fixture data was recorded on 2021-05-31. Pin `approximate_year` so + # this test stays valid as the wall clock drifts past 2030 (when the + # default `now()` anchor would start picking the wrong rollover cycle + # for these archived fixtures). + pvt = calc_pvt(states; approximate_year = 2021) expected_lla = LLA(; lat = 50.778851672464015, lon = 6.065568885758519, alt = 289.4069805158367) @test pvt.position ≈ ECEFfromLLA(wgs84)(expected_lla) rtol = 1e-8 @test pvt.time ≈ TAIEpoch(2021, 5, 31, 12, 53, 14.1183385390904732) @test pvt.velocity ≈ ECEF(0.0, 0.0, 0.0) atol = 9 @test get_frequency_offset(pvt, get_center_frequency(galileo_e1b)) ≈ -(1675.63Hz + freq_offset) atol = 0.01Hz - warm_pvt = calc_pvt(states, pvt) + warm_pvt = calc_pvt(states, pvt; approximate_year = 2021) @test get_LLA(warm_pvt) ≈ get_LLA(pvt) @test warm_pvt.time ≈ pvt.time @test warm_pvt.velocity ≈ pvt.velocity atol = 1e-6 @@ -1426,14 +1430,16 @@ end ), ] - pvt = calc_pvt(states) + # Fixture data was recorded on 2021-05-31. See the note on the + # Galileo testset above for why we pin `approximate_year`. + pvt = calc_pvt(states; approximate_year = 2021) expected_lla = LLA(; lat = 50.77885249310784, lon = 6.0656199911189175, alt = 291.95658091689086) @test pvt.position ≈ ECEFfromLLA(wgs84)(expected_lla) rtol = 1e-8 @test pvt.time ≈ TAIEpoch(2021, 5, 31, 12, 53, 14.1491024351271335) @test pvt.velocity ≈ ECEF(0.0, 0.0, 0.0) atol = 2.5 @test get_frequency_offset(pvt, get_center_frequency(gpsl1)) ≈ -(1632.59Hz + freq_offset) atol = 0.01Hz - warm_pvt = calc_pvt(states, pvt) + warm_pvt = calc_pvt(states, pvt; approximate_year = 2021) @test get_LLA(warm_pvt) ≈ get_LLA(pvt) @test warm_pvt.time ≈ pvt.time @test warm_pvt.velocity ≈ pvt.velocity atol = 1e-6 diff --git a/test/runtests.jl b/test/runtests.jl index 3708ced..f5e5c36 100755 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,8 @@ -using Test, PositionVelocityTime, GNSSDecoder, AstroTime, BitIntegers, GNSSSignals, Geodesy +using Test, PositionVelocityTime, GNSSDecoder, AstroTime, BitIntegers, GNSSSignals, Geodesy, Dates using Unitful: Hz include("aqua.jl") include("sat_time.jl") -include("pvt.jl") \ No newline at end of file +include("pvt.jl") +include("get_week.jl") \ No newline at end of file