From 4f2f3c1d2fe87f0f12ff7e3e6652b3b28fba1536 Mon Sep 17 00:00:00 2001 From: Soeren Schoenbrod Date: Wed, 27 Mar 2024 16:36:23 +0100 Subject: [PATCH] Create an easier interface for CN0 estimate (#48) Co-authored-by: Soeren Schoenbrod --- Project.toml | 2 +- src/Tracking.jl | 4 +- src/cn0_estimation.jl | 58 +++++++++------------ src/sat_state.jl | 17 ++++-- src/tracking_state.jl | 45 ++++++++++------ src/update_sat_state.jl | 6 +-- test/cn0_estimation.jl | 88 ++++++++++++++++++++++---------- test/conventional_pll_and_dll.jl | 4 +- test/runtests.jl | 6 ++- test/sat_state.jl | 8 ++- 10 files changed, 141 insertions(+), 97 deletions(-) diff --git a/Project.toml b/Project.toml index 60b1e4e..61f6744 100644 --- a/Project.toml +++ b/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" diff --git a/src/Tracking.jl b/src/Tracking.jl index 56b3db2..580b3da 100644 --- a/src/Tracking.jl +++ b/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}} diff --git a/src/cn0_estimation.jl b/src/cn0_estimation.jl index cfdd455..7e99477 100644 --- a/src/cn0_estimation.jl +++ b/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 diff --git a/src/sat_state.jl b/src/sat_state.jl index f808821..dbf5595 100644 --- a/src/sat_state.jl +++ b/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 diff --git a/src/tracking_state.jl b/src/tracking_state.jl index b367c35..97f051a 100644 --- a/src/tracking_state.jl +++ b/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) = diff --git a/src/update_sat_state.jl b/src/update_sat_state.jl index fcf05fb..04ae1c6 100644 --- a/src/update_sat_state.jl +++ b/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 diff --git a/test/cn0_estimation.jl b/test/cn0_estimation.jl index 636c0b4..609cde3 100644 --- a/test/cn0_estimation.jl +++ b/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 \ No newline at end of file diff --git a/test/conventional_pll_and_dll.jl b/test/conventional_pll_and_dll.jl index 339656a..d04551c 100644 --- a/test/conventional_pll_and_dll.jl +++ b/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(), ), ], diff --git a/test/runtests.jl b/test/runtests.jl index 0ef3266..949d16a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,7 @@ using TrackingLoopFilters using Acquisition: AcquisitionResults using Unitful: MHz, kHz, Hz, s, ms, dBHz +#= include("conventional_pll_and_dll.jl") include("sat_state.jl") include("sample_parameters.jl") @@ -26,6 +27,7 @@ include("galileo_e1b.jl") include("secondary_code_or_bit_detector.jl") include("bit_buffer.jl") include("correlator.jl") +=# include("cn0_estimation.jl") -include("tracking_state.jl") -include("track.jl") +#include("tracking_state.jl") +#include("track.jl") diff --git a/test/sat_state.jl b/test/sat_state.jl index 1d66d24..b0fc090 100644 --- a/test/sat_state.jl +++ b/test/sat_state.jl @@ -12,7 +12,6 @@ @test get_last_fully_integrated_filtered_prompt(sat_state) == complex(0.0, 0.0) @test get_sample_of_last_fully_integrated_correlator(sat_state) == -1 @test Tracking.found(get_secondary_code_or_bit_detector(sat_state)) == false - @test length(get_prompts_buffer(sat_state)) == 0 @test length(get_bit_buffer(sat_state)) == 0 acq = AcquisitionResults( @@ -38,7 +37,6 @@ @test get_last_fully_integrated_filtered_prompt(sat_state) == complex(0.0, 0.0) @test get_sample_of_last_fully_integrated_correlator(sat_state) == -1 @test Tracking.found(get_secondary_code_or_bit_detector(sat_state)) == false - @test length(get_prompts_buffer(sat_state)) == 0 @test length(get_bit_buffer(sat_state)) == 0 sampling_frequency = 5e6Hz @@ -58,7 +56,7 @@ complex(0.0, 0.0), 14, Tracking.SecondaryCodeOrBitDetector(), - Tracking.PromptsBuffer(20), + Tracking.MomentsCN0Estimator(20), Tracking.BitBuffer(), ), SatState( @@ -79,7 +77,7 @@ complex(0.0, 0.0), 15, Tracking.SecondaryCodeOrBitDetector(), - Tracking.PromptsBuffer(20), + Tracking.MomentsCN0Estimator(20), Tracking.BitBuffer(), ), SatState( @@ -100,7 +98,7 @@ complex(0.0, 0.0), -1, Tracking.SecondaryCodeOrBitDetector(), - Tracking.PromptsBuffer(20), + Tracking.MomentsCN0Estimator(20), Tracking.BitBuffer(), ), ],