fix: resolve GPS L1 week-rollover ambiguity via approximate_year

Project.toml

@@ -6,6 +6,7 @@ authors = ["Soeren Schoenbrod <soeren.schoenbrod@nav.rwth-aachen.de>", "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"

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
 

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

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

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")
+include("pvt.jl")
+include("get_week.jl")