Create an easier interface for CN0 estimate (#48)

Co-authored-by: Soeren Schoenbrod <soeren.schoenbrod@juliahub.com>

Project.toml

@@ -20,7 +20,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Acquisition = "0.1.1"
 CUDA = "4.0"
 DocStringExtensions = "0.6, 0.7, 0.8, 0.9"
-GNSSSignals = "0.17"
+GNSSSignals = "0.17.1"
 LoopVectorization = "0.8, 0.9, 0.10, 0.11, 0.12"
 StaticArrays = "0.9, 0.10, 0.11, 0.12, 1.0"
 StructArrays = "0.4, 0.6.5"

src/Tracking.jl

@@ -26,7 +26,6 @@ export get_early,
  get_last_fully_integrated_filtered_prompt,
  get_sample_of_last_fully_integrated_correlator,
  get_secondary_code_or_bit_detector,
- get_prompts_buffer,
  get_bit_buffer,
  get_bits,
  get_accumulators,
@@ -53,7 +52,8 @@ export get_early,
  remove_sats!,
  get_sat_states,
  get_sat_state,
- get_system
+ get_system,
+ estimate_cn0
 
 const StructVecOrMat{T} = Union{StructVector{T},StructArray{T,2}}
 

src/cn0_estimation.jl

@@ -5,57 +5,45 @@ $(SIGNATURES)
 
 MomentsCN0Estimator to estimate the CN0
 """
-struct MomentsCN0Estimator <: AbstractCN0Estimator end
-
-"""
-$(SIGNATURES)
-
-MomentsCN0Estimator to estimate the CN0
-"""
-struct PromptsBuffer
+struct MomentsCN0Estimator <: AbstractCN0Estimator
  prompt_buffer::Vector{ComplexF64}
- current_index::Int
- length::Int
+ buffer_current_index::Int
+ filled_buffer_length::Int
 end
 
-function PromptsBuffer(N)
- PromptsBuffer(zeros(ComplexF64, N), 0, 0)
+function MomentsCN0Estimator(N)
+ MomentsCN0Estimator(zeros(ComplexF64, N), 0, 0)
 end
 
-length(buffer::PromptsBuffer) = buffer.length
-get_prompt_buffer(buffer::PromptsBuffer) = buffer.prompt_buffer
-get_current_index(buffer::PromptsBuffer) = buffer.current_index
+length(estimator::MomentsCN0Estimator) = estimator.filled_buffer_length
+get_prompt_buffer(estimator::MomentsCN0Estimator) = estimator.prompt_buffer
+get_current_index(estimator::MomentsCN0Estimator) = estimator.buffer_current_index
 
 """
 $(SIGNATURES)
 
-Updates the `cn0_estimator` to include the information of the current prompt value.
+Buffers the prompts such that they can be used to estimate the CN0
 """
-function update(buffer::PromptsBuffer, prompt)
- buffer_length = length(buffer.prompt_buffer)
- next_index = mod(get_current_index(buffer), buffer_length) + 1
- buffer.prompt_buffer[next_index] = prompt
- next_length = min(length(buffer) + 1, buffer_length)
- PromptsBuffer(buffer.prompt_buffer, next_index, next_length)
+function update(estimator::MomentsCN0Estimator, prompt)
+ buffer_length = length(estimator.prompt_buffer)
+ next_index = mod(get_current_index(estimator), buffer_length) + 1
+ estimator.prompt_buffer[next_index] = prompt
+ filled_buffer_length = min(length(estimator) + 1, buffer_length)
+ MomentsCN0Estimator(estimator.prompt_buffer, next_index, filled_buffer_length)
 end
 
 """
 $(SIGNATURES)
 
-Estimates the CN0 based on the struct `cn0_estimator`.
+Estimates the CN0 based on the struct `MomentsCN0Estimator`.
 """
-function estimate_cn0(
- buffer::PromptsBuffer,
- cn0_estimator::MomentsCN0Estimator,
- integration_time,
-)
- length(buffer) == 0 && return 0.0dBHz
- abs2_prompt_buffer = abs2.(get_prompt_buffer(buffer))
- M₂ = 1 / length(buffer) * sum(abs2_prompt_buffer)
- M₄ = 1 / length(buffer) * sum(abs2_prompt_buffer .^ 2)
+function estimate_cn0(estimator::MomentsCN0Estimator, integration_time)
+ length(estimator) == 0 && return 0.0dBHz
+ prompt_buffer = get_prompt_buffer(estimator)
+ M₂ = 1 / length(estimator) * sum(abs2, prompt_buffer)
+ M₄ = 1 / length(estimator) * sum(x -> abs2(x)^2, prompt_buffer)
  Pd = sqrt(abs(2 * M₂^2 - M₄))
- noise_power = (M₂ - Pd)
- noise_power_non_neg = noise_power - 2 * (noise_power < 0) * noise_power
- SNR = Pd / noise_power_non_neg
+ noise_power = abs(M₂ - Pd)
+ SNR = Pd / noise_power
  dBHz(SNR / integration_time)
 end

src/sat_state.jl

@@ -10,7 +10,7 @@ struct SatState{C<:AbstractCorrelator}
  last_fully_integrated_filtered_prompt::ComplexF64
  sample_of_last_fully_integrated_correlator::Int
  sc_bit_detector::SecondaryCodeOrBitDetector
- prompts_buffer::PromptsBuffer
+ cn0_estimator::MomentsCN0Estimator
  bit_buffer::BitBuffer
 end
 
@@ -28,7 +28,7 @@ get_last_fully_integrated_filtered_prompt(s::SatState) =
 get_sample_of_last_fully_integrated_correlator(s::SatState) =
  s.sample_of_last_fully_integrated_correlator
 get_secondary_code_or_bit_detector(s::SatState) = s.sc_bit_detector
-get_prompts_buffer(s::SatState) = s.prompts_buffer
+get_cn0_estimator(s::SatState) = s.cn0_estimator
 get_bit_buffer(s::SatState) = s.bit_buffer
 
 function SatState(
@@ -40,7 +40,7 @@ function SatState(
  num_ants::NumAnts = NumAnts(1),
  carrier_phase = 0.0,
  code_doppler = carrier_doppler * get_code_center_frequency_ratio(system),
- num_prompts_buffer::Int = 20,
+ num_prompts_for_cn0_estimation::Int = 100,
 )
  SatState(
  prn,
@@ -54,7 +54,7 @@ function SatState(
  complex(0.0, 0.0),
  -1,
  SecondaryCodeOrBitDetector(),
- PromptsBuffer(num_prompts_buffer),
+ MomentsCN0Estimator(num_prompts_for_cn0_estimation),
  BitBuffer(),
  )
 end
@@ -77,3 +77,12 @@ end
 
 get_system(sss::SystemSatsState) = sss.system
 get_states(sss::SystemSatsState) = sss.states
+get_sat_state(sss::SystemSatsState, sat_idx::Integer) = sss.states[sat_idx]
+
+function estimate_cn0(sss::SystemSatsState, sat_idx::Integer)
+ system = sss.system
+ estimate_cn0(
+ get_cn0_estimator(get_sat_state(sss, sat_idx)),
+ get_code_length(system) / get_code_frequency(system),
+ )
+end

src/tracking_state.jl

@@ -51,18 +51,25 @@ function TrackState(
  )
 end
 
-function get_sat_states(
+function get_sats_states(
  track_state::TrackState{<:TupleLike{<:NTuple{N,SystemSatsState}}},
  system_idx::Union{Symbol,Integer},
 ) where {N}
- track_state.system_sats_states[system_idx].states
+ track_state.system_sats_states[system_idx]
 end
 
-function get_sat_states(
+function get_sats_states(
  track_state::TrackState{<:TupleLike{<:NTuple{N,SystemSatsState}}},
  system_idx::Val{M},
 ) where {N,M}
- track_state.system_sats_states[M].states
+ track_state.system_sats_states[M]
+end
+
+function get_sat_states(
+ track_state::TrackState{<:TupleLike{<:NTuple{N,SystemSatsState}}},
+ system_idx::Union{Symbol,Integer,Val},
+) where {N}
+ get_sats_states(track_state, system_idx).states
 end
 
 function get_sat_states(track_state::TrackState{<:TupleLike{<:NTuple{1,SystemSatsState}}})
@@ -71,16 +78,9 @@ end
 
 function get_system(
  track_state::TrackState{<:TupleLike{<:NTuple{N,SystemSatsState}}},
- system_idx::Union{Symbol,Integer},
+ system_idx::Union{Symbol,Integer,Val},
 ) where {N}
- track_state.system_sats_states[system_idx].system
-end
-
-function get_system(
- track_state::TrackState{<:TupleLike{<:NTuple{N,SystemSatsState}}},
- system_idx::Val{M},
-) where {N,M}
- track_state.system_sats_states[M].system
+ get_sats_states(track_state, system_idx).system
 end
 
 function get_system(track_state::TrackState{<:TupleLike{<:NTuple{1,SystemSatsState}}})
@@ -92,14 +92,29 @@ function get_sat_state(
  system_idx::Union{Symbol,Integer,Val},
  sat_idx::Integer,
 ) where {N}
- get_sat_states(track_state, system_idx)[sat_idx]
+ get_sat_state(get_sats_states(track_state, system_idx), sat_idx)
 end
 
 function get_sat_state(
  track_state::TrackState{<:TupleLike{<:NTuple{1,SystemSatsState}}},
  sat_idx::Integer,
 )
- get_sat_states(track_state, 1)[sat_idx]
+ get_sat_state(track_state, 1, sat_idx)
+end
+
+function estimate_cn0(
+ track_state::TrackState{<:TupleLike{<:NTuple{N,SystemSatsState}}},
+ system_idx::Union{Symbol,Integer,Val},
+ sat_idx::Integer,
+) where {N}
+ estimate_cn0(get_sats_states(track_state, system_idx), sat_idx)
+end
+
+function estimate_cn0(
+ track_state::TrackState{<:TupleLike{<:NTuple{1,SystemSatsState}}},
+ sat_idx::Integer,
+)
+ estimate_cn0(track_state, 1, sat_idx)
 end
 
 create_buffer(buffers::Vector{<:CPUBuffers}, system, sat_states, num_samples) =

src/update_sat_state.jl

@@ -48,7 +48,7 @@ function update(
  sat_state.last_fully_integrated_filtered_prompt,
  sample_of_last_fully_integrated_correlator,
  sat_state.sc_bit_detector,
- sat_state.prompts_buffer,
+ sat_state.cn0_estimator,
  sat_state.bit_buffer,
  )
  end
@@ -77,7 +77,7 @@ function update(
  filtered_prompt = dopplers_and_filtered_prompt.filtered_prompt
  carrier_doppler = dopplers_and_filtered_prompt.carrier_doppler
  code_doppler = dopplers_and_filtered_prompt.code_doppler
- prompts_buffer = update(state.prompts_buffer, filtered_prompt)
+ cn0_estimator = update(get_cn0_estimator(state), filtered_prompt)
  bit_buffer = buffer(
  system_sats.system,
  state.bit_buffer,
@@ -103,7 +103,7 @@ function update(
  sample_params.signal_start_sample +
  sample_params.signal_samples_to_integrate - 1,
  sc_bit_detector,
- prompts_buffer,
+ cn0_estimator,
  bit_buffer,
  )
  new_state

test/cn0_estimation.jl

@@ -1,24 +1,24 @@
-@testset "Prompts buffer" begin
- prompts_buffer = Tracking.PromptsBuffer(20)
- @test @inferred(Tracking.get_prompt_buffer(prompts_buffer)) ==
+@testset "Moments CN0 estimator" begin
+ cn0_estimator = Tracking.MomentsCN0Estimator(20)
+ @test @inferred(Tracking.get_prompt_buffer(cn0_estimator)) ==
  zero(SVector{20,ComplexF64})
- @test @inferred(Tracking.get_current_index(prompts_buffer)) == 0
- @test @inferred(Tracking.length(prompts_buffer)) == 0
-
- next_prompts_buffer = @inferred Tracking.update(prompts_buffer, 1 + 2im)
- @test @inferred(Tracking.get_prompt_buffer(next_prompts_buffer))[1] == 1 + 2im
- @test @inferred(Tracking.get_current_index(next_prompts_buffer)) == 1
- @test @inferred(Tracking.length(next_prompts_buffer)) == 1
-
- prompts_buffer = Tracking.PromptsBuffer(ones(SVector{20,ComplexF64}), 20, 20)
- next_prompts_buffer = @inferred Tracking.update(prompts_buffer, 1 + 2im)
- @test @inferred(Tracking.get_prompt_buffer(next_prompts_buffer))[1] == 1 + 2im
- @test @inferred(Tracking.get_current_index(next_prompts_buffer)) == 1
- @test @inferred(Tracking.length(next_prompts_buffer)) == 20
-
- prompts_buffer = Tracking.PromptsBuffer(ones(SVector{20,ComplexF64}), 19, 20)
- @test @inferred(Tracking.get_current_index(prompts_buffer)) == 19
- @test @inferred(Tracking.length(prompts_buffer)) == 20
+ @test @inferred(Tracking.get_current_index(cn0_estimator)) == 0
+ @test @inferred(Tracking.length(cn0_estimator)) == 0
+
+ next_cn0_estimator = @inferred Tracking.update(cn0_estimator, 1 + 2im)
+ @test @inferred(Tracking.get_prompt_buffer(next_cn0_estimator))[1] == 1 + 2im
+ @test @inferred(Tracking.get_current_index(next_cn0_estimator)) == 1
+ @test @inferred(Tracking.length(next_cn0_estimator)) == 1
+
+ cn0_estimator = Tracking.MomentsCN0Estimator(ones(SVector{20,ComplexF64}), 20, 20)
+ next_cn0_estimator = @inferred Tracking.update(cn0_estimator, 1 + 2im)
+ @test @inferred(Tracking.get_prompt_buffer(next_cn0_estimator))[1] == 1 + 2im
+ @test @inferred(Tracking.get_current_index(next_cn0_estimator)) == 1
+ @test @inferred(Tracking.length(next_cn0_estimator)) == 20
+
+ cn0_estimator = Tracking.MomentsCN0Estimator(ones(SVector{20,ComplexF64}), 19, 20)
+ @test @inferred(Tracking.get_current_index(cn0_estimator)) == 19
+ @test @inferred(Tracking.length(cn0_estimator)) == 20
 end
 
 @testset "CN0 estimation" begin
@@ -30,13 +30,12 @@ end
  prn = 1
  range = 0:3999
  start_carrier_phase = π / 2
- cn0_estimator = MomentsCN0Estimator()
- prompts_buffer = Tracking.PromptsBuffer(20)
+ cn0_estimator = Tracking.MomentsCN0Estimator(100)
  start_sample = 1
  num_samples = 4000
  gpsl1 = GPSL1()
 
- for i = 1:20
+ for i = 1:100
  signal =
  get_code.(
  gpsl1,
@@ -53,11 +52,44 @@ end
  signal_struct = StructArray(signal)
  correlator =
  Tracking.correlate(correlator, signal_struct, code, start_sample, num_samples)
- prompts_buffer = Tracking.update(prompts_buffer, get_prompt(correlator))
+ cn0_estimator = Tracking.update(cn0_estimator, get_prompt(correlator))
  end
- @test @inferred(Tracking.get_current_index(prompts_buffer)) == 20
- @test @inferred(Tracking.length(prompts_buffer)) == 20
- cn0_estimate = @inferred Tracking.estimate_cn0(prompts_buffer, cn0_estimator, 1ms)
+ @test @inferred(Tracking.get_current_index(cn0_estimator)) == 100
+ @test @inferred(Tracking.length(cn0_estimator)) == 100
+ cn0_estimate = @inferred estimate_cn0(cn0_estimator, 1ms)
 
- @test cn0_estimate ≈ 45dBHz atol = 2.05dBHz
+ @test cn0_estimate ≈ 45dBHz atol = 1.0dBHz
 end
+
+@testset "CN0 estimation integration test" begin
+ Random.seed!(1234)
+ carrier_doppler = 0Hz
+ start_code_phase = 0
+ code_frequency = 1023kHz
+ sampling_frequency = 4MHz
+ prn = 1
+ range = 0:3999
+ start_carrier_phase = π / 2
+ cn0_estimator = Tracking.MomentsCN0Estimator(100)
+ start_sample = 1
+ num_samples = 4000
+ gpsl1 = GPSL1()
+
+ track_state = @inferred TrackState(
+ gpsl1,
+ [SatState(gpsl1, prn, sampling_frequency, start_code_phase, carrier_doppler)];
+ num_samples,
+ )
+
+ for i = 1:100
+ signal =
+ get_code.(
+ gpsl1,
+ code_frequency .* range ./ sampling_frequency .+ start_code_phase,
+ prn,
+ ) .* 10^(45 / 20) .+ randn(ComplexF64, length(range)) .* sqrt(4e6)
+ track_state = @inferred track(signal, track_state, sampling_frequency)
+ end
+ cn0_estimate = @inferred estimate_cn0(track_state, 1)
+ @test cn0_estimate ≈ 45dBHz atol = 1.0dBHz
+end

test/conventional_pll_and_dll.jl

@@ -49,7 +49,7 @@ end
  complex(0.0, 0.0),
  14,
  Tracking.SecondaryCodeOrBitDetector(),
- Tracking.PromptsBuffer(20),
+ Tracking.MomentsCN0Estimator(20),
  Tracking.BitBuffer(),
  ),
  ],
@@ -100,7 +100,7 @@ end
  complex(0.0, 0.0),
  14,
  Tracking.SecondaryCodeOrBitDetector(),
- Tracking.PromptsBuffer(20),
+ Tracking.MomentsCN0Estimator(20),
  Tracking.BitBuffer(),
  ),
  ],