From 2c5daa0befb26d2bbe74a4bebcf5f3b21361784b Mon Sep 17 00:00:00 2001 From: Soeren Schoenbrod Date: Fri, 19 May 2023 21:18:45 +0200 Subject: [PATCH] Restructure: Allow to track multiple satellites (#43) --- .JuliaFormatter.toml | 8 + Project.toml | 13 +- src/Tracking.jl | 106 ++-- src/bit_buffer.jl | 37 +- src/cn0_estimation.jl | 47 +- src/conventional_pll_and_dll.jl | 130 +++++ src/correlator.jl | 374 ++++++++++---- src/discriminators.jl | 35 +- src/downconvert.jl | 2 +- src/downconvert_and_correlate.jl | 238 +++++---- src/galileo_e1b.jl | 4 +- src/gpsl1.jl | 4 +- src/gpsl5.jl | 4 +- src/gpu_downconvert_and_correlate.jl | 202 ++++++++ src/post_corr_filter.jl | 15 +- src/sample_parameters.jl | 148 ++++++ src/sat_state.jl | 84 ++++ src/track.jl | 60 +++ src/tracking_loop.jl | 515 -------------------- src/tracking_results.jl | 177 ------- src/tracking_state.jl | 376 +++++--------- src/update_sat_state.jl | 107 ++++ test/bit_buffer.jl | 73 +-- test/cn0_estimation.jl | 55 +-- test/code_replica.jl | 6 +- test/conventional_pll_and_dll.jl | 118 +++++ test/correlator.jl | 217 ++++++--- test/cuda/bit_buffer.jl | 93 ---- test/cuda/cn0_estimation.jl | 50 -- test/cuda/discriminators.jl | 61 --- test/cuda/galileo_e1b.jl | 15 - test/cuda/gps_l1.jl | 11 - test/cuda/gps_l5.jl | 11 - test/cuda/secondary_code_or_bit_detector.jl | 40 -- test/cuda/tracking_loop.jl | 267 ---------- test/cuda/tracking_results.jl | 52 -- test/cuda/tracking_state.jl | 60 --- test/discriminators.jl | 43 +- test/downconvert_and_correlate.jl | 151 ++++++ test/galileo_e1b.jl | 8 +- test/gps_l1.jl | 8 +- test/gps_l5.jl | 8 +- test/post_corr_filter.jl | 2 +- test/runtests.jl | 22 +- test/sample_parameters.jl | 146 ++++++ test/sat_state.jl | 196 ++++++++ test/track.jl | 491 +++++++++++++++++++ test/tracking_loop.jl | 393 --------------- test/tracking_results.jl | 52 -- test/tracking_state.jl | 287 ++++++++--- 50 files changed, 3014 insertions(+), 2608 deletions(-) create mode 100644 .JuliaFormatter.toml create mode 100644 src/conventional_pll_and_dll.jl create mode 100644 src/gpu_downconvert_and_correlate.jl create mode 100644 src/sample_parameters.jl create mode 100644 src/sat_state.jl create mode 100644 src/track.jl delete mode 100644 src/tracking_loop.jl delete mode 100644 src/tracking_results.jl create mode 100644 src/update_sat_state.jl create mode 100644 test/conventional_pll_and_dll.jl delete mode 100644 test/cuda/bit_buffer.jl delete mode 100644 test/cuda/cn0_estimation.jl delete mode 100644 test/cuda/discriminators.jl delete mode 100644 test/cuda/galileo_e1b.jl delete mode 100644 test/cuda/gps_l1.jl delete mode 100644 test/cuda/gps_l5.jl delete mode 100644 test/cuda/secondary_code_or_bit_detector.jl delete mode 100644 test/cuda/tracking_loop.jl delete mode 100644 test/cuda/tracking_results.jl delete mode 100644 test/cuda/tracking_state.jl create mode 100644 test/downconvert_and_correlate.jl create mode 100644 test/sample_parameters.jl create mode 100644 test/sat_state.jl create mode 100644 test/track.jl delete mode 100644 test/tracking_loop.jl delete mode 100644 test/tracking_results.jl diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..b32e65e --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,8 @@ +# Configuration file for JuliaFormatter.jl +# For more information, see: https://domluna.github.io/JuliaFormatter.jl/stable/config/ + +margin = 92 +remove_extra_newlines = true +separate_kwargs_with_semicolon = true +format_docstrings = true +format_markdown = true \ No newline at end of file diff --git a/Project.toml b/Project.toml index b34d85e..60d1c0d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,29 +1,30 @@ name = "Tracking" uuid = "10b2438b-ffd4-5096-aa58-44041d5c8f3b" -authors = ["Soeren Schoenbrod ", - "Michael Niestroj "] +authors = ["Soeren Schoenbrod ", "Michael Niestroj "] version = "0.15.2" [deps] Acquisition = "d4bbf60b-102b-5ffb-8f97-a7ea5817e69f" +Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" GNSSSignals = "52c80523-2a4e-5c38-8979-05588f836870" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" TrackingLoopFilters = "0814aff9-93cb-554c-9fff-9bf1cfdb5efa" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] -Acquisition = "0.1.0" -CUDA = "3.5" +Acquisition = "0.1.1" +CUDA = "4.0" DocStringExtensions = "0.6, 0.7, 0.8, 0.9" -GNSSSignals = "0.16" +GNSSSignals = "0.17" 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" TrackingLoopFilters = "0.1" Unitful = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 1.0" -julia = "1.6" \ No newline at end of file +julia = "1.6" diff --git a/src/Tracking.jl b/src/Tracking.jl index 3cb968c..49c3a3e 100644 --- a/src/Tracking.jl +++ b/src/Tracking.jl @@ -12,37 +12,46 @@ module Tracking using Unitful: upreferred, Hz, dBHz, ms import Base.zero, Base.length, Base.resize!, LinearAlgebra.dot - export get_early, get_prompt, get_late - export get_early_index,get_prompt_index, get_late_index - export get_correlator - export get_accumulator, get_accumulators - export get_num_accumulators - export get_carrier_doppler, get_carrier_phase - export get_code_doppler, get_code_phase - export get_correlator_sample_shifts - export get_early_late_sample_spacing - export get_early_late_index_shift - export get_secondary_code_or_bit_found - export get_correlator_carrier_phase, get_correlator_carrier_frequency - export get_state - export get_system - export get_cn0 - export track - export get_bits - export get_num_bits - export get_filtered_prompt - export get_post_corr_filter - - export TrackingState - export NumAnts - export NumAccumulators - export MomentsCN0Estimator - export AbstractCN0Estimator - export EarlyPromptLateCorrelator - #export VeryEarlyPromptLateCorrelator - export SecondaryCodeOrBitDetector - export GainControlledSignal - export AbstractPostCorrFilter + export get_early, get_prompt, get_late, + get_prn, + get_code_phase, + get_code_doppler, + get_carrier_phase, + get_carrier_doppler, + get_integrated_samples, + get_correlator, + get_last_fully_integrated_correlator, + 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, + get_early_late_sample_spacing, + track, + NumAnts, + NumAccumulators, + MomentsCN0Estimator, + AbstractCN0Estimator, + EarlyPromptLateCorrelator, + VeryEarlyPromptLateCorrelator, + SecondaryCodeOrBitDetector, + GainControlledSignal, + AbstractPostCorrFilter, + SatState, + SystemSatsState, + CPUDownconvertAndCorrelator, + GPUDownconvertAndCorrelator, + ConventionalPLLAndDLL, + ConventionalPLLsAndDLLs, + DefaultPostCorrFilter, + TrackState, + add_sats!, + remove_sats!, + get_sat_states, + get_sat_state, + get_system struct NumAnts{x} end @@ -54,20 +63,45 @@ module Tracking NumAccumulators(x) = NumAccumulators{x}() - include("post_corr_filter.jl") + TupleLike{T <: Tuple} = Union{T, NamedTuple{<:Any, T}} + + struct DopplersAndFilteredPrompt + carrier_doppler::typeof(1.0Hz) + code_doppler::typeof(1.0Hz) + filtered_prompt::ComplexF64 + end + + """ + $(SIGNATURES) + + Get the number of samples in the signal. + """ + @inline function get_num_samples(signal) + length(signal) + end + + @inline function get_num_samples(signal::AbstractMatrix) + size(signal, 1) + end + include("code_replica.jl") include("carrier_replica.jl") include("downconvert.jl") include("cn0_estimation.jl") - include("discriminators.jl") include("bit_buffer.jl") include("correlator.jl") + include("discriminators.jl") + include("post_corr_filter.jl") include("secondary_code_or_bit_detector.jl") - include("tracking_state.jl") - include("tracking_results.jl") - include("tracking_loop.jl") include("gpsl1.jl") include("gpsl5.jl") include("galileo_e1b.jl") + include("sat_state.jl") + include("sample_parameters.jl") + include("update_sat_state.jl") include("downconvert_and_correlate.jl") + include("gpu_downconvert_and_correlate.jl") + include("conventional_pll_and_dll.jl") + include("tracking_state.jl") + include("track.jl") end diff --git a/src/bit_buffer.jl b/src/bit_buffer.jl index 8a7cb13..81c21cf 100644 --- a/src/bit_buffer.jl +++ b/src/bit_buffer.jl @@ -6,10 +6,12 @@ BitBuffer to buffer bits struct BitBuffer buffer::UInt64 length::Int + prompt_accumulator::ComplexF64 + prompt_accumulator_integrated_code_blocks::Int end function BitBuffer() - BitBuffer(0, 0) + BitBuffer(0, 0, complex(0.0, 0.0), 0) end @inline get_bits(bit_buffer::BitBuffer) = bit_buffer.buffer @@ -23,24 +25,31 @@ Buffer data bits based on the prompt accumulation and the current prompt value. function buffer( system::AbstractGNSS, bit_buffer, - prompt_accumulator, + integrated_code_blocks, secondary_code_or_bit_found, - prev_code_phase, - code_phase, - integration_time, - prompt_correlator + prompt ) - prompt_accumulator = prompt_accumulator + secondary_code_or_bit_found * - prompt_correlator + prompt_accumulator = bit_buffer.prompt_accumulator + secondary_code_or_bit_found * prompt + prompt_accumulator_integrated_code_blocks = bit_buffer.prompt_accumulator_integrated_code_blocks + + secondary_code_or_bit_found * integrated_code_blocks + + num_code_blocks_that_form_a_bit = Int(get_code_frequency(system) / (get_code_length(system) * get_data_frequency(system))) - if secondary_code_or_bit_found && - (code_phase - prev_code_phase < 0 || integration_time == 1 / get_data_frequency(system)) + return if secondary_code_or_bit_found && + prompt_accumulator_integrated_code_blocks == num_code_blocks_that_form_a_bit bit = real(prompt_accumulator) > 0 - bit_buffer = BitBuffer( + BitBuffer( get_bits(bit_buffer) << 1 + UInt64(bit), - length(bit_buffer) + 1 + length(bit_buffer) + 1, + zero(prompt_accumulator), + 0 + ) + else + BitBuffer( + bit_buffer.buffer, + bit_buffer.length, + prompt_accumulator, + prompt_accumulator_integrated_code_blocks ) - prompt_accumulator = zero(prompt_accumulator) end - bit_buffer, prompt_accumulator end diff --git a/src/cn0_estimation.jl b/src/cn0_estimation.jl index 4913772..f4db1b2 100644 --- a/src/cn0_estimation.jl +++ b/src/cn0_estimation.jl @@ -5,30 +5,39 @@ $(SIGNATURES) MomentsCN0Estimator to estimate the CN0 """ -struct MomentsCN0Estimator{N} <: AbstractCN0Estimator - prompt_buffer::SVector{N, ComplexF64} +struct MomentsCN0Estimator <: AbstractCN0Estimator +end + +""" +$(SIGNATURES) + +MomentsCN0Estimator to estimate the CN0 +""" +struct PromptsBuffer + prompt_buffer::Vector{ComplexF64} current_index::Int length::Int end -function MomentsCN0Estimator(N) - MomentsCN0Estimator{N}(zero(SVector{N, ComplexF64}), 0, 0) +function PromptsBuffer(N) + PromptsBuffer(zeros(ComplexF64, N), 0, 0) end -length(cn0_state::MomentsCN0Estimator) = cn0_state.length -get_prompt_buffer(cn0_state::MomentsCN0Estimator) = cn0_state.prompt_buffer -get_current_index(cn0_state::MomentsCN0Estimator) = cn0_state.current_index +length(buffer::PromptsBuffer) = buffer.length +get_prompt_buffer(buffer::PromptsBuffer) = buffer.prompt_buffer +get_current_index(buffer::PromptsBuffer) = buffer.current_index """ $(SIGNATURES) Updates the `cn0_estimator` to include the information of the current prompt value. """ -function update(cn0_estimator::MomentsCN0Estimator{N}, prompt) where N - next_index = mod(get_current_index(cn0_estimator), N) + 1 - next_prompt_buffer = setindex(get_prompt_buffer(cn0_estimator), prompt, next_index) - next_length = min(length(cn0_estimator) + 1, N) - MomentsCN0Estimator{N}(next_prompt_buffer, next_index, next_length) +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) end """ @@ -36,11 +45,15 @@ $(SIGNATURES) Estimates the CN0 based on the struct `cn0_estimator`. """ -function estimate_cn0(cn0_estimator::MomentsCN0Estimator, integration_time) - length(cn0_estimator) == 0 && return 0.0dBHz - abs2_prompt_buffer = abs2.(get_prompt_buffer(cn0_estimator)) - M₂ = 1 / length(cn0_estimator) * sum(abs2_prompt_buffer) - M₄ = 1 / length(cn0_estimator) * sum(abs2_prompt_buffer .^ 2) +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) Pd = sqrt(abs(2 * M₂^2 - M₄)) noise_power = (M₂ - Pd) noise_power_non_neg = noise_power - 2 * (noise_power < 0) * noise_power diff --git a/src/conventional_pll_and_dll.jl b/src/conventional_pll_and_dll.jl new file mode 100644 index 0000000..c2e1765 --- /dev/null +++ b/src/conventional_pll_and_dll.jl @@ -0,0 +1,130 @@ +abstract type AbstractDopplerEstimator end + +struct ConventionalPLLAndDLL{ + CA <: AbstractLoopFilter, + CO <: AbstractLoopFilter, + F <: AbstractPostCorrFilter +} + init_carrier_doppler::typeof(1.0Hz) + init_code_doppler::typeof(1.0Hz) + carrier_loop_filter::CA + code_loop_filter::CO + carrier_loop_filter_bandwidth::typeof(1.0Hz) + code_loop_filter_bandwidth::typeof(1.0Hz) + post_corr_filter::F +end + +function ConventionalPLLAndDLL( + sat_state::SatState, + carrier_loop_filter = ThirdOrderBilinearLF(), + code_loop_filter = SecondOrderBilinearLF(), + carrier_loop_filter_bandwidth = 18Hz, + code_loop_filter_bandwidth = 1Hz, + post_corr_filter = DefaultPostCorrFilter() +) + ConventionalPLLAndDLL( + sat_state.carrier_doppler, + sat_state.code_doppler, + carrier_loop_filter, + code_loop_filter, + float(carrier_loop_filter_bandwidth), + float(code_loop_filter_bandwidth), + post_corr_filter + ) +end + +struct ConventionalPLLsAndDLLs{ + N, + PD <: TupleLike{<:NTuple{N, Vector{<:ConventionalPLLAndDLL}}}, +} <: AbstractDopplerEstimator + plls_and_dlls::PD +end + +""" +$(SIGNATURES) + +Aid dopplers. That is velocity aiding for the carrier doppler and carrier aiding +for the code doppler. +""" +function aid_dopplers( + system::AbstractGNSS, + init_carrier_doppler, + init_code_doppler, + carrier_freq_update, + code_freq_update +) + carrier_doppler = carrier_freq_update + code_doppler = code_freq_update + carrier_doppler * get_code_center_frequency_ratio(system) + init_carrier_doppler + carrier_doppler, init_code_doppler + code_doppler +end + +function estimate_dopplers_and_filter_prompt( + doppler_estimator::ConventionalPLLsAndDLLs, + system_sats_state::TupleLike{<:NTuple{N, SystemSatsState}}, + sat_sample_params::TupleLike{<:NTuple{N, Vector{SampleParams}}}, + sampling_frequency, + min_integration_time +) where N + results = map(doppler_estimator.plls_and_dlls, system_sats_state, sat_sample_params) do estimator_per_system, system_sats, sat_sample_params_per_system + map(estimator_per_system, system_sats.states, sat_sample_params_per_system) do estimator, state, sample_params + integration_time = state.integrated_samples / sampling_frequency + return if sample_params.signal_samples_to_integrate == sample_params.samples_to_integrate_until_code_end && integration_time >= min_integration_time + correlator = normalize(state.correlator, state.integrated_samples) + post_corr_filter = update(estimator.post_corr_filter, get_prompt(correlator)) + filtered_correlator = apply(post_corr_filter, correlator) + pll_discriminator = pll_disc( + system_sats.system, + filtered_correlator, + ) + dll_discriminator = dll_disc( + system_sats.system, + filtered_correlator, + (state.code_doppler + get_code_frequency(system_sats.system)) / sampling_frequency + ) + carrier_freq_update, carrier_loop_filter = filter_loop( + estimator.carrier_loop_filter, + pll_discriminator, + integration_time, + estimator.carrier_loop_filter_bandwidth + ) + code_freq_update, code_loop_filter = filter_loop( + estimator.code_loop_filter, + dll_discriminator, + integration_time, + estimator.code_loop_filter_bandwidth + ) + carrier_doppler, code_doppler = aid_dopplers( + system_sats.system, + estimator.init_carrier_doppler, + estimator.init_code_doppler, + carrier_freq_update, + code_freq_update, + ) + dopplers_and_filtered_prompt = DopplersAndFilteredPrompt( + carrier_doppler, + code_doppler, + get_prompt(filtered_correlator) + ) + new_estimator = ConventionalPLLAndDLL( + estimator.init_carrier_doppler, + estimator.init_code_doppler, + carrier_loop_filter, + code_loop_filter, + estimator.carrier_loop_filter_bandwidth, + estimator.code_loop_filter_bandwidth, + post_corr_filter + ) + (estimator = new_estimator, dopplers_and_filtered_prompt = dopplers_and_filtered_prompt) + else + (estimator = estimator, dopplers_and_filtered_prompt = DopplersAndFilteredPrompt( + 0.0Hz, + 0.0Hz, + complex(0.0,0.0) + )) + end + end + end + new_doppler_estimator = ConventionalPLLsAndDLLs(map(x -> map(y -> y.estimator, x), results)) + new_dopplers_and_filtered_prompts = map(x -> map(y -> y.dopplers_and_filtered_prompt, x), results) + new_doppler_estimator, new_dopplers_and_filtered_prompts +end \ No newline at end of file diff --git a/src/correlator.jl b/src/correlator.jl index 2889b6d..4ab6a9e 100644 --- a/src/correlator.jl +++ b/src/correlator.jl @@ -1,48 +1,206 @@ -abstract type AbstractCorrelator{T, TA} end +abstract type AbstractCorrelator{M} end +abstract type AbstractEarlyPromptLateCorrelator{M} <: AbstractCorrelator{M} end """ $(SIGNATURES) EarlyPromptLateCorrelator holding a user defined number of correlation values. """ -struct EarlyPromptLateCorrelator{T, TA <: AbstractVector{T}} <: AbstractCorrelator{T, TA} +struct EarlyPromptLateCorrelator{ + M, + T, + TA <: AbstractVector{T}, + S <: AbstractVector{Int} +} <: AbstractEarlyPromptLateCorrelator{M} accumulators::TA - EarlyPromptLateCorrelator{T, TA}(accumulators) where {T, TA <: AbstractVector{T}} = length(accumulators) < 3 ? - error("Early-Prompt-Late-Correlator needs at least 3 accumulators") : new{T, TA}(accumulators) -end - -function EarlyPromptLateCorrelator(accumulators::AbstractVector) - EarlyPromptLateCorrelator{typeof(first(accumulators)), typeof(accumulators)}(accumulators) + shifts::S + early_shift::Int + prompt_shift::Int + late_shift::Int + function EarlyPromptLateCorrelator( + accumulators::AbstractVector{SVector{M, T}}, + shifts::S, + early_shift::Int, + prompt_shift::Int, + late_shift::Int, + ) where {M, T <: Complex, S} + idxs = SVector(early_shift, prompt_shift, late_shift) + length(accumulators) < 3 ? + error("Early-Prompt-Late-Correlator needs at least 3 distinct accumulators") : + length(shifts) != length(accumulators) ? + error("Number of accumulators and shifts must be the same.") : + !allunique(shifts[idxs]) ? + error("Early, Prompt and Late must not have the same sample shift. Either increase sample rate or adjust the preferred code shift.") : + new{M, SVector{M, T}, typeof(accumulators), S}(accumulators, shifts, early_shift, prompt_shift, late_shift) + end + function EarlyPromptLateCorrelator( + accumulators::AbstractVector{T}, + shifts::S, + early_shift::Int, + prompt_shift::Int, + late_shift::Int, + ) where {T <: Complex, S} + idxs = SVector(early_shift, prompt_shift, late_shift) + length(accumulators) < 3 ? + error("Early-Prompt-Late-Correlator needs at least 3 distinct accumulators") : + length(shifts) != length(accumulators) ? + error("Number of accumulators and shifts must be the same.") : + !allunique(shifts[idxs]) ? + error("Early, Prompt and Late must not have the same sample shift. Either increase sample rate or adjust the preferred code shift.") : + new{1, T, typeof(accumulators), S}(accumulators, shifts, early_shift, prompt_shift, late_shift) + end end -type_for_num_ants(num_ants::NumAnts{1}) = ComplexF64 -type_for_num_ants(num_ants::NumAnts{N}) where N = SVector{N, ComplexF64} - """ $(SIGNATURES) EarlyPromptLateCorrelator constructor without parameters and some default parameters. """ function EarlyPromptLateCorrelator( + system::AbstractGNSS, + sampling_frequency; + preferred_code_shift = 0.5, + preferred_early_late_to_prompt_code_shift = 0.5, # TODO: At the momemnt this must be similar to preferred_code_shift num_ants::NumAnts = NumAnts(1), - num_accumulators::NumAccumulators{M} = NumAccumulators(3), -) where M + num_accumulators = NumAccumulators(3), +) + shifts = get_correlator_sample_shifts( + system, + get_num_accumulators(num_accumulators), + sampling_frequency, + preferred_code_shift + ) + prompt_shift = length(shifts) >> 1 + 1 + + sample_shift = preferred_code_shift_to_sample_shift( + preferred_early_late_to_prompt_code_shift, + sampling_frequency, + system + ) + _, late_shift = findmin(abs, shifts .+ sample_shift) + _, early_shift = findmin(abs, shifts .- sample_shift) + EarlyPromptLateCorrelator( - zero(SVector{M, type_for_num_ants(num_ants)}) + get_initial_accumulator(num_ants, num_accumulators), + shifts, + early_shift, + prompt_shift, + late_shift ) end """ $(SIGNATURES) -EarlyPromptLateCorrelator constructor for large number of accumulators. +EarlyPromptLateCorrelator holding a user defined number of correlation values. """ -function EarlyPromptLateCorrelator( +struct VeryEarlyPromptLateCorrelator{ + M, + T, + TA <: AbstractVector{T}, + S <: AbstractVector{Int} +} <: AbstractEarlyPromptLateCorrelator{M} + accumulators::TA + shifts::S + very_early_shift::Int + early_shift::Int + prompt_shift::Int + late_shift::Int + very_late_shift::Int + function VeryEarlyPromptLateCorrelator( + accumulators::AbstractVector{SVector{M, T}}, + shifts::S, + very_early_shift::Int, + early_shift::Int, + prompt_shift::Int, + late_shift::Int, + very_late_shift::Int + ) where {M, T <: Complex, S} + idxs = SVector(very_early_shift, early_shift, prompt_shift, late_shift, very_late_shift) + length(accumulators) < 5 ? + error("Very-Early-Prompt-Late-Correlator needs at least 5 distinct accumulators") : + length(shifts) != length(accumulators) ? + error("Number of accumulators and shifts must be the same.") : + !allunique(shifts[idxs]) ? + error("Early, Prompt and Late must not have the same sample shift. Either increase sample rate or adjust the preferred code shift.") : + new{M, SVector{M, T}, typeof(accumulators), S}(accumulators, shifts, very_early_shift, early_shift, prompt_shift, late_shift, very_late_shift) + end + function VeryEarlyPromptLateCorrelator( + accumulators::AbstractVector{T}, + shifts::S, + very_early_shift::Int, + early_shift::Int, + prompt_shift::Int, + late_shift::Int, + very_late_shift::Int + ) where {T <: Complex, S} + idxs = SVector(very_early_shift, early_shift, prompt_shift, late_shift, very_late_shift) + length(accumulators) < 5 ? + error("Very-Early-Prompt-Late-Correlator needs at least 5 distinct accumulators") : + length(shifts) != length(accumulators) ? + error("Number of accumulators and shifts must be the same.") : + !allunique(shifts[idxs]) ? + error("Early, Prompt and Late must not have the same sample shift. Either increase sample rate or adjust the preferred code shift.") : + new{1, T, typeof(accumulators), S}(accumulators, shifts, very_early_shift, early_shift, prompt_shift, late_shift, very_late_shift) + end +end + +""" +$(SIGNATURES) + +VeryEarlyPromptLateCorrelator constructor without parameters and some default parameters. +Default parameters take from https://gnss-sdr.org/docs/sp-blocks/tracking/#implementation-galileo_e1_dll_pll_veml_tracking +""" +function VeryEarlyPromptLateCorrelator( + system::AbstractGNSS, + sampling_frequency; + preferred_early_late_to_prompt_code_shift = 0.15, + preferred_very_early_late_to_prompt_code_shift = 0.6, + num_ants::NumAnts = NumAnts(1), +) + early_late_sample_shift = preferred_code_shift_to_sample_shift( + preferred_early_late_to_prompt_code_shift, + sampling_frequency, + system + ) + very_early_late_sample_shift = preferred_code_shift_to_sample_shift( + preferred_very_early_late_to_prompt_code_shift, + sampling_frequency, + system + ) + shifts = SVector( + -very_early_late_sample_shift, + -early_late_sample_shift, + 0, + early_late_sample_shift, + very_early_late_sample_shift + ) + + VeryEarlyPromptLateCorrelator( + get_initial_accumulator(num_ants, NumAccumulators(5)), + shifts, + 5, + 4, + 3, + 2, + 1 + ) +end + +type_for_num_ants(num_ants::NumAnts{1}) = ComplexF64 +type_for_num_ants(num_ants::NumAnts{N}) where N = SVector{N, ComplexF64} + +function get_initial_accumulator( + num_ants::NumAnts, + num_accumulators::NumAccumulators{M} +) where M + zero(SVector{M, type_for_num_ants(num_ants)}) +end + +function get_initial_accumulator( num_ants::NumAnts, num_accumulators::Integer ) - EarlyPromptLateCorrelator( - [zero(type_for_num_ants(num_ants)) for i = 1:num_accumulators] - ) + [zero(type_for_num_ants(num_ants)) for i = 1:num_accumulators] end @@ -51,8 +209,7 @@ $(SIGNATURES) Get number of antennas from correlator """ -get_num_ants(correlator::AbstractCorrelator{Complex{T}}) where {T} = 1 -get_num_ants(correlator::AbstractCorrelator{SVector{N,T}}) where {N,T} = N +get_num_ants(correlator::AbstractCorrelator{M}) where {M} = M """ $(SIGNATURES) @@ -60,29 +217,43 @@ $(SIGNATURES) Get number of accumulators """ get_num_accumulators(correlator::AbstractCorrelator) = size(correlator.accumulators, 1) +get_num_accumulators(num_accumulators::Int) = num_accumulators +get_num_accumulators(num_accumulators::NumAccumulators{M}) where M = M + +""" +$(SIGNATURES) + +Get early correlator index +""" +@inline get_very_early_index(corr::VeryEarlyPromptLateCorrelator) = corr.very_early_shift """ $(SIGNATURES) Get early correlator index """ -@inline get_early_index(correlator_sample_shifts, early_late_index_shift) = - get_prompt_index(correlator_sample_shifts) + early_late_index_shift +@inline get_early_index(corr::AbstractEarlyPromptLateCorrelator) = corr.early_shift """ $(SIGNATURES) Get prompt correlator index """ -@inline get_prompt_index(correlator_sample_shifts) = findfirst(iszero, correlator_sample_shifts) +@inline get_prompt_index(corr::AbstractEarlyPromptLateCorrelator) = corr.prompt_shift + +""" +$(SIGNATURES) + +Get late correlator index +""" +@inline get_late_index(corr::AbstractEarlyPromptLateCorrelator) = corr.late_shift """ $(SIGNATURES) Get late correlator index """ -@inline get_late_index(correlator_sample_shifts, early_late_index_shift) = - get_prompt_index(correlator_sample_shifts) - early_late_index_shift +@inline get_very_late_index(corr::VeryEarlyPromptLateCorrelator) = corr.very_late_shift """ $(SIGNATURES) @@ -97,8 +268,17 @@ $(SIGNATURES) Get a specific accumulator with `index` counted negative for late and positive for early accumulators. """ -function get_accumulator(correlator::AbstractCorrelator, correlator_sample_shifts, index::Integer) - correlator.accumulators[index + get_prompt_index(correlator_sample_shifts)] +function get_accumulator(correlator::AbstractCorrelator, index::Integer) + correlator.accumulators[index + get_prompt_index(correlator)] +end + +""" +$(SIGNATURES) + +Get very early correlator +""" +function get_very_early(correlator::VeryEarlyPromptLateCorrelator) + correlator.accumulators[get_very_early_index(correlator)] end """ @@ -106,14 +286,8 @@ $(SIGNATURES) Get early correlator """ -function get_early( - correlator::AbstractCorrelator, - correlator_sample_shifts, - early_late_index_shift -) - correlator.accumulators[ - get_early_index(correlator_sample_shifts, early_late_index_shift) - ] +function get_early(correlator::AbstractEarlyPromptLateCorrelator) + correlator.accumulators[get_early_index(correlator)] end """ @@ -121,8 +295,8 @@ $(SIGNATURES) Get prompt correlator """ -function get_prompt(correlator::AbstractCorrelator, correlator_sample_shifts) - correlator.accumulators[get_prompt_index(correlator_sample_shifts)] +function get_prompt(correlator::AbstractCorrelator) + correlator.accumulators[get_prompt_index(correlator)] end """ @@ -130,14 +304,43 @@ $(SIGNATURES) Get late correlator """ -function get_late( - correlator::AbstractCorrelator, - correlator_sample_shifts, - early_late_index_shift -) - correlator.accumulators[ - get_late_index(correlator_sample_shifts, early_late_index_shift) - ] +function get_late(correlator::AbstractEarlyPromptLateCorrelator) + correlator.accumulators[get_late_index(correlator)] +end + +""" +$(SIGNATURES) + +Get very late correlator +""" +function get_very_late(correlator::VeryEarlyPromptLateCorrelator) + correlator.accumulators[get_very_late_index(correlator)] +end + +""" +$(SIGNATURES) + +Update the Correlator +""" +function update_accumulator(correlator::EarlyPromptLateCorrelator, accumulators) + EarlyPromptLateCorrelator( + accumulators, + correlator.shifts, + correlator.early_shift, + correlator.prompt_shift, + correlator.late_shift + ) +end +function update_accumulator(correlator::VeryEarlyPromptLateCorrelator, accumulators) + VeryEarlyPromptLateCorrelator( + accumulators, + correlator.shifts, + correlator.very_early_shift, + correlator.early_shift, + correlator.prompt_shift, + correlator.late_shift, + correlator.very_late_shift + ) end """ @@ -145,8 +348,8 @@ $(SIGNATURES) Zero the Correlator """ -function zero(correlator::T) where T <: AbstractCorrelator - T(zero(correlator.accumulators)) +function zero(correlator::AbstractCorrelator) + update_accumulator(correlator, zero(correlator.accumulators)) end """ @@ -154,8 +357,8 @@ $(SIGNATURES) Filter the Correlator by the function `post_corr_filter` """ -function apply(post_corr_filter, correlator::T) where T <: AbstractCorrelator - (T.name.wrapper)(map(post_corr_filter, get_accumulators(correlator))) +function apply(post_corr_filter, correlator::AbstractCorrelator) + update_accumulator(correlator, map(post_corr_filter, get_accumulators(correlator))) end """ @@ -167,32 +370,19 @@ ordered from latest to earliest replica. """ function get_correlator_sample_shifts( system::AbstractGNSS, - correlator::AbstractCorrelator{T, <:SVector{M}}, - sampling_frequency, - preferred_code_shift -) where {T,M} - num_corrs = floor(Int, M / 2) - sample_shift = preferred_code_shift_to_sample_shift( - preferred_code_shift, - sampling_frequency, - system - ) - SVector{M}(-num_corrs:num_corrs) .* sample_shift -end -function get_correlator_sample_shifts( - system::AbstractGNSS, - correlator::AbstractCorrelator, + num_correlators::Int, sampling_frequency, preferred_code_shift ) - num_corrs = floor(Int, get_num_accumulators(correlator) / 2) + half_num_correlators = num_correlators >> 1 sample_shift = preferred_code_shift_to_sample_shift( preferred_code_shift, sampling_frequency, system ) - (-num_corrs:num_corrs) .* sample_shift + (-half_num_correlators:half_num_correlators) .* sample_shift end + function preferred_code_shift_to_sample_shift( preferred_code_shift, sampling_frequency, @@ -210,21 +400,16 @@ Calculates the index for the early and late samples function get_early_late_index_shift( system, correlator_sample_shifts, - correlator, sampling_frequency, preferred_code_shift ) - preferred_code_phase_distance = Inf - early_late_index_shift = 1 - for i = get_prompt_index(correlator_sample_shifts) + 1:length(correlator_sample_shifts) - code_phase_shift = correlator_sample_shifts[i] * get_code_frequency(system) / sampling_frequency - preferred_code_phase_distance_temp = abs(code_phase_shift - preferred_code_shift) - if preferred_code_phase_distance_temp < preferred_code_phase_distance - early_late_index_shift = i - get_prompt_index(correlator_sample_shifts) - preferred_code_phase_distance = preferred_code_phase_distance_temp - end - end - early_late_index_shift + sample_shift = preferred_code_shift_to_sample_shift( + preferred_code_shift, + sampling_frequency, + system + ) + min_val, min_idx = findmin(abs, correlator_sample_shifts .- sample_shift) + min_idx end """ @@ -233,11 +418,10 @@ $(SIGNATURES) Calculate the total spacing between early and late correlator in samples. """ function get_early_late_sample_spacing( - correlator_sample_shifts, - early_late_index_shift + correlator::AbstractEarlyPromptLateCorrelator ) - correlator_sample_shifts[get_early_index(correlator_sample_shifts, early_late_index_shift)] - - correlator_sample_shifts[get_late_index(correlator_sample_shifts, early_late_index_shift)] + correlator.shifts[get_early_index(correlator)] - + correlator.shifts[get_late_index(correlator)] end """ @@ -254,26 +438,25 @@ $(SIGNATURES) Perform a correlation for multi antenna systems """ function correlate( - correlator::T, + correlator::AbstractCorrelator{1}, downconverted_signal::AbstractVector, code, - correlator_sample_shifts, start_sample, num_samples -) where {T <: AbstractCorrelator} +) a_re = zero_accumulators(get_accumulators(correlator), downconverted_signal) a_im = zero_accumulators(get_accumulators(correlator), downconverted_signal) d_re = downconverted_signal.re d_im = downconverted_signal.im @avx for i = start_sample:num_samples + start_sample - 1 for j = 1:length(a_re) - sample_shift = correlator_sample_shifts[j] - correlator_sample_shifts[1] + sample_shift = correlator.shifts[j] - correlator.shifts[1] a_re[j] += d_re[i] * code[i + sample_shift] a_im[j] += d_im[i] * code[i + sample_shift] end end accumulators_result = complex.(a_re, a_im) - T(map(+, get_accumulators(correlator), accumulators_result)) + update_accumulator(correlator, map(+, get_accumulators(correlator), accumulators_result)) end function zero_accumulators(accumulators::SVector, signal) @@ -289,13 +472,12 @@ $(SIGNATURES) Perform a correlation for multi antenna systems """ function correlate( - correlator::AbstractCorrelator{<: SVector{N}}, + correlator::AbstractCorrelator{M}, downconverted_signal::AbstractMatrix, code, - correlator_sample_shifts, start_sample, num_samples, -) where {N} +) where {M} a_re = zero_accumulators(get_accumulators(correlator), downconverted_signal) a_im = zero_accumulators(get_accumulators(correlator), downconverted_signal) d_re = downconverted_signal.re @@ -303,14 +485,20 @@ function correlate( @avx for i = start_sample:num_samples + start_sample - 1 for k = 1:size(a_re, 2) for j = 1:size(a_re, 1) - shift = correlator_sample_shifts[k] - correlator_sample_shifts[1] + shift = correlator.shifts[k] - correlator.shifts[1] a_re[j,k] += d_re[i,j] * code[shift + i] a_im[j,k] += d_im[i,j] * code[shift + i] end end end + update_accumulator(correlator, add_to_previous(get_accumulators(correlator), a_re, a_im)) +end - typeof(correlator)(map(+, get_accumulators(correlator), eachcol(complex.(a_re, a_im)))) +function add_to_previous(accumulators::SVector{NC, <:SVector}, a_re, a_im) where NC + SVector{NC, eltype(accumulators)}(map(+, accumulators, eachcol(complex.(a_re, a_im)))) +end +function add_to_previous(accumulators::Vector{<:SVector}, a_re, a_im) + map(+, accumulators, eachcol(complex.(a_re, a_im))) end function zero_accumulators(accumulators::SVector{NC, <:SVector{NA}}, signal) where {NC, NA} diff --git a/src/discriminators.jl b/src/discriminators.jl index 02285ad..e0112f5 100644 --- a/src/discriminators.jl +++ b/src/discriminators.jl @@ -5,25 +5,44 @@ Calculates the code phase error in chips. """ function dll_disc( system::AbstractGNSS, - correlator, - correlator_sample_shifts, - early_late_index_shift, + correlator::EarlyPromptLateCorrelator, code_phase_delta ) - E = abs(get_early(correlator, correlator_sample_shifts, early_late_index_shift)) - L = abs(get_late(correlator, correlator_sample_shifts, early_late_index_shift)) + E = abs(get_early(correlator)) + L = abs(get_late(correlator)) distance_between_early_and_late = - get_early_late_sample_spacing(correlator_sample_shifts, early_late_index_shift) * + get_early_late_sample_spacing(correlator) * code_phase_delta (E - L) / (E + L) / (2 * (2 - distance_between_early_and_late)) end +""" +$(SIGNATURES) +DLL discriminator for very early prompt late correlator. +Implementation from: +https://gnss-sdr.org/docs/sp-blocks/tracking/#implementation-galileo_e1_dll_pll_veml_tracking +""" +function dll_disc( + system::AbstractGNSS, + correlator::VeryEarlyPromptLateCorrelator, + code_phase_delta +) + VE = abs(get_very_early(correlator)) + E = abs(get_early(correlator)) + L = abs(get_late(correlator)) + VL = abs(get_very_late(correlator)) + distance_between_early_and_late = + get_early_late_sample_spacing(correlator) * + code_phase_delta + (VE + E - VL - L) / (VE + E + VL + L) / (2 * (2 - distance_between_early_and_late)) +end + """ $(SIGNATURES) Calculates the carrier phase error in radians. """ -function pll_disc(system::AbstractGNSS, correlator, correlator_sample_shifts) - p = get_prompt(correlator, correlator_sample_shifts) +function pll_disc(system::AbstractGNSS, correlator) + p = get_prompt(correlator) atan(imag(p) / real(p)) end diff --git a/src/downconvert.jl b/src/downconvert.jl index 04ae6b7..e70853c 100644 --- a/src/downconvert.jl +++ b/src/downconvert.jl @@ -51,7 +51,7 @@ function downconvert!( start_phase, start_sample, num_samples -) where {N, T, TS} +) where {T, TS} ds_re = downconverted_signal.re; ds_im = downconverted_signal.im s_re = signal.re; s_im = signal.im carrier_freq = upreferred(carrier_frequency / Hz) diff --git a/src/downconvert_and_correlate.jl b/src/downconvert_and_correlate.jl index b19727b..e949d3f 100644 --- a/src/downconvert_and_correlate.jl +++ b/src/downconvert_and_correlate.jl @@ -31,121 +31,155 @@ function downconvert_and_correlate( end =# -# CUDA Kernel -function downconvert_and_correlate_kernel( - res_re, - res_im, - signal_re, - signal_im, - codes, +function downconvert_and_correlate!( + system, + signal, + correlator, + code_replica, + code_phase, + carrier_replica, + carrier_phase, + downconverted_signal, code_frequency, - correlator_sample_shifts, carrier_frequency, sampling_frequency, - start_code_phase, - carrier_phase, - code_length, - prn, - num_samples, - num_ants, - num_corrs -) - cache = @cuDynamicSharedMem(Float32, (2 * blockDim().x, num_ants, num_corrs)) - sample_idx = 1 + ((blockIdx().x - 1) * blockDim().x + (threadIdx().x - 1)) - antenna_idx = 1 + ((blockIdx().y - 1) * blockDim().y + (threadIdx().y - 1)) - corr_idx = 1 + ((blockIdx().z - 1) * blockDim().z + (threadIdx().z - 1)) - iq_offset = blockDim().x - cache_index = threadIdx().x - 1 - - code_phase = accum_re = accum_im = dw_re = dw_im = carrier_re = carrier_im = 0.0f0 - mod_floor_code_phase = Int(0) + signal_start_sample, + num_samples_left, + prn +) + gen_code_replica!( + code_replica, + system, + code_frequency, + sampling_frequency, + code_phase, + signal_start_sample, + num_samples_left, + correlator.shifts, + prn + ) + gen_carrier_replica!( + carrier_replica, + carrier_frequency, + sampling_frequency, + carrier_phase, + signal_start_sample, + num_samples_left + ) + downconvert!( + downconverted_signal, + signal, + carrier_replica, + signal_start_sample, + num_samples_left + ) + correlate( + correlator, + downconverted_signal, + code_replica, + signal_start_sample, + num_samples_left + ) +end - if sample_idx <= num_samples && antenna_idx <= num_ants && corr_idx <= num_corrs - # generate carrier - carrier_im, carrier_re = CUDA.sincos(2π * ((sample_idx - 1) * carrier_frequency / sampling_frequency + carrier_phase)) - - # downconvert with the conjugate of the carrier - dw_re = signal_re[sample_idx, antenna_idx] * carrier_re + signal_im[sample_idx, antenna_idx] * carrier_im - dw_im = signal_im[sample_idx, antenna_idx] * carrier_re - signal_re[sample_idx, antenna_idx] * carrier_im +abstract type AbstractDownconvertAndCorrelator end - # calculate the code phase - code_phase = code_frequency / sampling_frequency * ((sample_idx - 1) + correlator_sample_shifts[corr_idx]) + start_code_phase +const StructVecOrMat{T} = Union{StructVector{T}, StructArray{T, 2}} - # wrap the code phase around the code length e.g. phase = 1024 -> modfloorphase = 1 - mod_floor_code_phase = 1 + mod(floor(Int32, code_phase), code_length) +struct CPUBuffers{T, CT, DS <: StructVecOrMat{Complex{T}}} + code_replica_buffer::Vector{CT} + carrier_replica_buffer::StructVector{Complex{T}} + downconvert_signal_buffer::DS +end - # multiply elementwise with the code - accum_re += codes[mod_floor_code_phase, prn] * dw_re - accum_im += codes[mod_floor_code_phase, prn] * dw_im - end +function CPUBuffers( + ::Type{T}, + system::AbstractGNSS, + state::SatState, + num_samples, +) where T + code_shifts = get_correlator(state).shifts + CPUBuffers( + Vector{get_code_type(system)}(undef, num_samples + maximum(code_shifts) - minimum(code_shifts)), + StructVector{Complex{T}}(undef, num_samples), + get_downconvert_signal_buffer(T, num_samples, state) + ) +end - cache[1 + cache_index + 0 * iq_offset, antenna_idx, corr_idx] = accum_re - cache[1 + cache_index + 1 * iq_offset, antenna_idx, corr_idx] = accum_im +function CPUBuffers( + system::AbstractGNSS, + state::SatState, + num_samples, +) + CPUBuffers(Float32, system, state, num_samples) +end - ## Reduction - # wait until all the accumulators have done writing the results to the cache - sync_threads() +struct CPUDownconvertAndCorrelator{N, B <: TupleLike{<:NTuple{N, Vector{<:CPUBuffers}}}} <: AbstractDownconvertAndCorrelator + buffers::B + function CPUDownconvertAndCorrelator(buffers::TupleLike{<:NTuple{N, Vector{<:CPUBuffers}}}) where N + new{N, typeof(buffers)}(buffers) + end +end - i::Int = blockDim().x ÷ 2 - @inbounds while i != 0 - if cache_index < i - cache[1 + cache_index + 0 * iq_offset, antenna_idx, corr_idx] += cache[1 + cache_index + 0 * iq_offset + i, antenna_idx, corr_idx] - cache[1 + cache_index + 1 * iq_offset, antenna_idx, corr_idx] += cache[1 + cache_index + 1 * iq_offset + i, antenna_idx, corr_idx] +function CPUDownconvertAndCorrelator( + ::Type{T}, + system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}, + num_samples::Integer, +) where {T, N} + CPUDownconvertAndCorrelator( + map(system_sats_states) do system_sats_state + system = system_sats_state.system + map(system_sats_state.states) do state + CPUBuffers(T, system, state, num_samples) + end end - sync_threads() - i ÷= 2 - end - - if (threadIdx().x - 1) == 0 - res_re[blockIdx().x, antenna_idx, corr_idx] += cache[1 + 0 * iq_offset, antenna_idx, corr_idx] - res_im[blockIdx().x, antenna_idx, corr_idx] += cache[1 + 1 * iq_offset, antenna_idx, corr_idx] - end - return nothing + ) +end +function CPUDownconvertAndCorrelator( + system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}, + num_samples::Integer, +) where N + CPUDownconvertAndCorrelator(Float32, system_sats_states, num_samples) end -function downconvert_and_correlate_kernel_wrapper( - system, +function get_downconvert_signal_buffer(::Type{T}, num_samples::Int, state::SatState{<: AbstractCorrelator{1}}) where T + StructVector{Complex{T}}(undef, num_samples) +end +function get_downconvert_signal_buffer(::Type{T}, num_samples::Int, state::SatState{<: AbstractCorrelator{M}}) where {T, M} + StructArray{Complex{T}}(undef, num_samples, M) +end + +function downconvert_and_correlate( + downconvert_and_correlator::CPUDownconvertAndCorrelator, signal, - correlator, - code_phase, - carrier_phase, - code_frequency, - correlator_sample_shifts, - carrier_frequency, sampling_frequency, - signal_start_sample, - num_samples_left, - prn -) - num_corrs = length(correlator_sample_shifts) - num_ants = size(signal, 2) - num_samples = size(signal, 1) - block_dim_z = num_corrs - block_dim_y = num_ants - # keep num_corrs and num_ants in seperate dimensions, truncate num_samples accordingly to fit - block_dim_x = prevpow(2, 1024 ÷ block_dim_y ÷ block_dim_z) - threads = (block_dim_x, block_dim_y, block_dim_z) - blocks = cld(size(signal, 1), block_dim_x) - res_re = CUDA.zeros(Float32, blocks, block_dim_y, block_dim_z) - res_im = CUDA.zeros(Float32, blocks, block_dim_y, block_dim_z) - shmem_size = sizeof(ComplexF32)*block_dim_x*block_dim_y*block_dim_z - @cuda threads=threads blocks=blocks shmem=shmem_size downconvert_and_correlate_kernel( - res_re, - res_im, - signal.re, - signal.im, - system.codes, - Float32(code_frequency), - correlator_sample_shifts, - Float32(carrier_frequency), - Float32(sampling_frequency), - Float32(code_phase), - Float32(carrier_phase), - size(system.codes, 1), - prn, - num_samples, - num_ants, - num_corrs - ) - return sum(res_re .+ 1im*res_im, dims=1) + intermediate_frequency, + system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}, + params::TupleLike{<:NTuple{N, Vector{SampleParams}}}, +) where N + map(params, system_sats_states, downconvert_and_correlator.buffers) do system_params, system_sats, buffers + map(system_params, system_sats.states, buffers) do sat_params, sat_state, buffer + if sat_params.signal_samples_to_integrate == 0 + return sat_state.correlator + end + carrier_frequency = sat_state.carrier_doppler + intermediate_frequency + code_frequency = sat_state.code_doppler + get_code_frequency(system_sats.system) + downconvert_and_correlate!( + system_sats.system, + signal, + sat_state.correlator, + buffer.code_replica_buffer, + sat_state.code_phase, + buffer.carrier_replica_buffer, + sat_state.carrier_phase, + buffer.downconvert_signal_buffer, + code_frequency, + carrier_frequency, + sampling_frequency, + sat_params.signal_start_sample, + sat_params.signal_samples_to_integrate, + sat_state.prn + )::typeof(sat_state.correlator) + end + end end \ No newline at end of file diff --git a/src/galileo_e1b.jl b/src/galileo_e1b.jl index e9cc74f..05d25a7 100644 --- a/src/galileo_e1b.jl +++ b/src/galileo_e1b.jl @@ -12,6 +12,6 @@ function is_upcoming_integration_new_bit(galileo_e1b::GalileoE1B, prns, num_prns end # TODO: Very early very late correlator? -function get_default_correlator(galileo_e1b::GalileoE1B, num_ants::NumAnts{N}) where N - EarlyPromptLateCorrelator(num_ants) +function get_default_correlator(galileo_e1b::GalileoE1B, sampling_frequency, num_ants::NumAnts) + VeryEarlyPromptLateCorrelator(galileo_e1b, sampling_frequency; num_ants) end diff --git a/src/gpsl1.jl b/src/gpsl1.jl index 065ca03..476ba04 100644 --- a/src/gpsl1.jl +++ b/src/gpsl1.jl @@ -11,6 +11,6 @@ function is_upcoming_integration_new_bit(gpsl1::GPSL1, prns, num_prns) masked_bit_synchronizer == 0xfffff || masked_bit_synchronizer == 0xfffff00000 end -function get_default_correlator(gpsl1::GPSL1, num_ants::NumAnts{N}) where N - EarlyPromptLateCorrelator(num_ants) +function get_default_correlator(gpsl1::GPSL1, sampling_frequency, num_ants::NumAnts) + EarlyPromptLateCorrelator(gpsl1, sampling_frequency; num_ants) end diff --git a/src/gpsl5.jl b/src/gpsl5.jl index d9d9918..d1cf912 100644 --- a/src/gpsl5.jl +++ b/src/gpsl5.jl @@ -11,6 +11,6 @@ function is_upcoming_integration_new_bit(gpsl5::GPSL5, prns, num_prns) xored_bit_synchronizer == 0 || xored_bit_synchronizer == 0x3ff end -function get_default_correlator(gpsl5::GPSL5, num_ants::NumAnts{N}) where N - EarlyPromptLateCorrelator(num_ants) +function get_default_correlator(gpsl5::GPSL5, sampling_frequency, num_ants::NumAnts) + EarlyPromptLateCorrelator(gpsl5, sampling_frequency; num_ants) end diff --git a/src/gpu_downconvert_and_correlate.jl b/src/gpu_downconvert_and_correlate.jl new file mode 100644 index 0000000..49e51d7 --- /dev/null +++ b/src/gpu_downconvert_and_correlate.jl @@ -0,0 +1,202 @@ +struct GPUBuffers{T, DS <: CuArray{Complex{T}, 3}} + downconverted_and_decoded_signal::DS +end + +function GPUBuffers( + ::Type{T}, + system::AbstractGNSS, + state::SatState, + num_samples, +) where T + GPUBuffers( + CuArray{ComplexF32}(undef, (num_samples, get_num_ants(state), get_num_accumulators(get_correlator(state)))) + ) +end + +function GPUBuffers( + system::AbstractGNSS, + state::SatState, + num_samples, +) + GPUBuffers(Float32, system, state, num_samples) +end + +struct GPUDownconvertAndCorrelator{N, CB <: TupleLike{<:NTuple{N, <:AbstractGNSS}}, B <: TupleLike{<:NTuple{N, Vector{<:GPUBuffers}}}} <: AbstractDownconvertAndCorrelator + code_buffers::CB + buffers::B + function GPUDownconvertAndCorrelator(code_buffers::TupleLike{<:NTuple{N, <:AbstractGNSS}}, buffers::TupleLike{<:NTuple{N, Vector{<:GPUBuffers}}}) where N + new{N, typeof(code_buffers), typeof(buffers)}(code_buffers, buffers) + end +end + +function GPUDownconvertAndCorrelator(::Type{T}, system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}, num_samples) where {T, N} + GPUDownconvertAndCorrelator( + map(convert_code_to_texture_memory, map(get_system, system_sats_states)), + map(system_sats_states) do system_sats_state + system = system_sats_state.system + map(system_sats_state.states) do state + GPUBuffers(T, system, state, num_samples) + end + end + ) +end + +function GPUDownconvertAndCorrelator( + system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}, + num_samples::Integer, +) where N + GPUDownconvertAndCorrelator(Float32, system_sats_states, num_samples) +end + +import Adapt + +Adapt.@adapt_structure GPSL1 +Adapt.@adapt_structure GPSL5 +Adapt.@adapt_structure GalileoE1B + +recreate_system_with_texture(system::GPSL1, texture) = GPSL1(texture) +recreate_system_with_texture(system::GPSL5, texture) = GPSL5(texture) +recreate_system_with_texture(system::GalileoE1B, texture) = GalileoE1B(texture) + +function convert_code_to_texture_memory(system::S) where S <: AbstractGNSS + # Get only base code without secondary code, since otherwise code might be too + # large for texture memory. Texture memory has a max size of 65536 in each + # 2D dimension. GPSL5 would have a size of 102300 with secondary code. + # Without secondary code GPSL5 has a code size of 10230. + # The secondary code is multiplied in the kernel instead. + # The same goes for any subcarrier code. + codes = get_codes(system)[1:get_code_length(system), :] + recreate_system_with_texture(system, CuTexture( + CuTextureArray(CuArray(Float32.(codes))), + address_mode = CUDA.ADDRESS_MODE_WRAP, + interpolation = CUDA.NearestNeighbour(), + )) +end + +function downconvert_and_correlate( + downconvert_and_correlator::GPUDownconvertAndCorrelator, + signal, + sampling_frequency, + intermediate_frequency, + system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}, + params::TupleLike{<:NTuple{N, Vector{SampleParams}}}, +) where N + map(params, system_sats_states, downconvert_and_correlator.code_buffers, downconvert_and_correlator.buffers) do system_params, system_sats, system, buffers + map(system_params, system_sats.states, buffers) do sat_params, sat_state, buffer + if sat_params.signal_samples_to_integrate == 0 + return sat_state.correlator + end + carrier_frequency = sat_state.carrier_doppler + intermediate_frequency + code_frequency = sat_state.code_doppler + get_code_frequency(system) + downconvert_and_correlate!( + system, + signal, + sat_state.correlator, + sat_state.code_phase, + sat_state.carrier_phase, + code_frequency, + carrier_frequency, + sampling_frequency, + sat_params.signal_start_sample, + sat_params.signal_samples_to_integrate, + sat_state.prn, + buffer.downconverted_and_decoded_signal, + )::typeof(sat_state.correlator) + end + end +end + +function get_code(system::AbstractGNSS, phase, prn) + get_code(system, get_modulation(system), phase, prn) +end + +function get_code(system::AbstractGNSS, modulation::GNSSSignals.LOC, phase, prn) + # Must add 0.5 because CUDA uses nearest neighbour instead of floor. + system.codes[phase + 0.5f0, prn] * get_secondary_code(system, phase) +end + +function get_code(system::AbstractGNSS, modulation::GNSSSignals.BOC, phase, prn) + # Must add 0.5 because CUDA uses nearest neighbour instead of floor. + system.codes[phase + 0.5f0, prn] * get_secondary_code(system, phase) * + GNSSSignals.get_subcarrier_code(modulation, phase) +end + +# This kernel currently assumes that we have more threads than number of samples +# to precess +# TODO: handle more samples than number of threads available +function downconvert_and_decode_prn_kernel!( + downconverted_and_decoded_signal, + signal, + system::AbstractGNSS, + prn::Int32, + correlator_sample_shifts, + num_samples::Int32, + code_frequency, + carrier_frequency, + sampling_frequency, + start_code_phase::Float32, + start_carrier_phase::Float32, + start_sample::Int32, + num_ants::NumAnts{N}, +) where {N} + sample = ((blockIdx().x - 0x1) * blockDim().x + (threadIdx().x - 0x1)) + index = sample + 0x1 + if sample < num_samples + carrier_wipe_off = cis(-Float32(2π) * (sample * carrier_frequency / sampling_frequency + start_carrier_phase)) + for sample_shift_index in eachindex(correlator_sample_shifts) + sample_shift = correlator_sample_shifts[sample_shift_index] + code = get_code(system, (sample + sample_shift) * code_frequency / sampling_frequency + start_code_phase, prn) + for antenna_index = 0x1:N + @inbounds downconverted_and_decoded_signal[index, antenna_index, sample_shift_index] = signal[sample + start_sample, antenna_index] * carrier_wipe_off * code + end + end + end + return +end + +function downconvert_and_correlate!( + code_buffer, + signal, + correlator::AbstractCorrelator{M}, + code_phase, + carrier_phase, + code_frequency, + carrier_frequency, + sampling_frequency, + signal_start_sample, + num_samples_left, + prn, + downconverted_and_decoded_signal, +) where M + # Assume 1024 to be the max number of threads + # TODO: Evaluate this at run time + threads = min(size(signal, 1), 1024) + blocks = cld(size(signal, 1), threads) + num_correlators = size(downconverted_and_decoded_signal, 3) + @cuda threads=threads blocks=blocks downconvert_and_decode_prn_kernel!( + downconverted_and_decoded_signal, + signal, + code_buffer, + Int32(prn), + correlator.shifts, + Int32(num_samples_left), + Float32(code_frequency / Hz), + Float32(carrier_frequency / Hz), + Float32(sampling_frequency / Hz), + Float32(code_phase), + Float32(carrier_phase), + Int32(signal_start_sample), + NumAnts{M}() + ) + correlated_signal = sum(view(downconverted_and_decoded_signal, 1:num_samples_left, :, :), dims = 1) + result = reshape(Array(correlated_signal), M, num_correlators) + gpu_add_to_accumulators(correlator, result) +end + +function gpu_add_to_accumulators(correlator::AbstractCorrelator{1}, result) + update_accumulator(correlator, SVector(map((a, b) -> a + b[1], get_accumulators(correlator), eachcol(result)))) +end + +function gpu_add_to_accumulators(correlator::AbstractCorrelator{M}, result) where M + update_accumulator(correlator, SVector(map(+, get_accumulators(correlator), eachcol(result)))) +end \ No newline at end of file diff --git a/src/post_corr_filter.jl b/src/post_corr_filter.jl index 8ee92ae..7fbf97f 100644 --- a/src/post_corr_filter.jl +++ b/src/post_corr_filter.jl @@ -5,4 +5,17 @@ struct DefaultPostCorrFilter <: AbstractPostCorrFilter end update(filter::DefaultPostCorrFilter, prompt) = filter (filter::DefaultPostCorrFilter)(x) = x -(filter::DefaultPostCorrFilter)(x::AbstractVector) = last(x) \ No newline at end of file +(filter::DefaultPostCorrFilter)(x::AbstractVector) = last(x) + +""" +$(SIGNATURES) + +Default post corr filter +""" +function get_default_post_corr_filter(correlator::AbstractCorrelator{T}) where T <: SVector + x -> x[1] +end + +function get_default_post_corr_filter(correlator) + x -> x +end \ No newline at end of file diff --git a/src/sample_parameters.jl b/src/sample_parameters.jl new file mode 100644 index 0000000..0e2f58e --- /dev/null +++ b/src/sample_parameters.jl @@ -0,0 +1,148 @@ +""" +$(SIGNATURES) + +Returns the appropiate number of code blocks to integrate. It will be just a single +code block when the secondary code or bit isn't found. The maximum number of code +blocks to integrate is limited by the bit shift. +""" +function calc_num_code_blocks_to_integrate( + system::AbstractGNSS, + preferred_num_code_blocks::Int, + secondary_code_or_bit_found::Bool +) + ifelse( + secondary_code_or_bit_found, + min( + Int(get_code_frequency(system) / (get_code_length(system) * get_data_frequency(system))), + preferred_num_code_blocks + ), + 1 + ) +end + +""" +$(SIGNATURES) + +Calculates the number of chips to integrate. +""" +function calc_num_chips_to_integrate( + system::AbstractGNSS, + num_code_blocks::Int, + current_code_phase +) + max_phase = num_code_blocks * get_code_length(system) + current_phase_mod_max_phase = mod(current_code_phase, max_phase) + return max_phase - current_phase_mod_max_phase +end + +""" +$(SIGNATURES) + +Calculates the number of samples to integrate. +""" +function calc_num_samples_left_to_integrate( + system::AbstractGNSS, + num_code_blocks::Int, + sampling_frequency, + current_code_doppler, + current_code_phase, +) + chips_to_integrate = calc_num_chips_to_integrate( + system, + num_code_blocks, + current_code_phase, + ) + code_frequency = get_code_frequency(system) + current_code_doppler + ceil(Int, chips_to_integrate * sampling_frequency / code_frequency) +end + +struct SampleParams + signal_samples_to_integrate::Int + signal_start_sample::Int + samples_to_integrate_until_code_end::Int + num_code_blocks_to_integrate::Int +end + +function SampleParams( + system::AbstractGNSS, + preferred_num_code_blocks_to_integrate::Int, + secondary_code_or_bit_found::Bool +) + num_code_blocks_to_integrate = calc_num_code_blocks_to_integrate( + system, + preferred_num_code_blocks_to_integrate, + secondary_code_or_bit_found + ) + SampleParams(0, 1, 0, num_code_blocks_to_integrate) +end + +function SampleParams( + system::AbstractGNSS, + prev_sample_params::SampleParams, + num_samples_signal::Int, + sampling_frequency, + code_doppler, + code_phase, + preferred_num_code_blocks_to_integrate::Int, + secondary_code_or_bit_found::Bool +) + num_code_blocks_to_integrate = calc_num_code_blocks_to_integrate( + system, + preferred_num_code_blocks_to_integrate, + secondary_code_or_bit_found + ) + samples_to_integrate_until_code_end = calc_num_samples_left_to_integrate( + system, + num_code_blocks_to_integrate, + sampling_frequency, + code_doppler, + code_phase + ) + signal_start_sample = prev_sample_params.signal_start_sample + prev_sample_params.signal_samples_to_integrate + signal_samples_left = num_samples_signal - signal_start_sample + 1 + signal_samples_to_integrate = min(samples_to_integrate_until_code_end, signal_samples_left) + SampleParams( + signal_samples_to_integrate, + signal_start_sample, + samples_to_integrate_until_code_end, + num_code_blocks_to_integrate + ) +end + +function init_sample_params( + system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}, + preferred_num_code_blocks_to_integrate::Int, +) where N + map(system_sats_states) do states + map(states.states) do state + SampleParams( + states.system, + preferred_num_code_blocks_to_integrate, + found(state.sc_bit_detector) + ) + end + end +end + +function calc_sample_params( + system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}, + prev_system_sats_sample_params, + num_samples_signal, + sampling_frequency, + preferred_num_code_blocks_to_integrate::Int, +) where N + map(system_sats_states, prev_system_sats_sample_params) do states, prev_sats_sample_params + map(states.states, prev_sats_sample_params) do state, prev_sample_params + SampleParams( + states.system, + prev_sample_params, + num_samples_signal, + sampling_frequency, + state.code_doppler, + state.code_phase, + preferred_num_code_blocks_to_integrate, + found(state.sc_bit_detector) + ) + end + end +end \ No newline at end of file diff --git a/src/sat_state.jl b/src/sat_state.jl new file mode 100644 index 0000000..b5eabb5 --- /dev/null +++ b/src/sat_state.jl @@ -0,0 +1,84 @@ +struct SatState{C <: AbstractCorrelator} + prn::Int + code_phase::Float64 + code_doppler::typeof(1.0Hz) + carrier_phase::Float64 + carrier_doppler::typeof(1.0Hz) + integrated_samples::Int + correlator::C + last_fully_integrated_correlator::C + last_fully_integrated_filtered_prompt::ComplexF64 + sample_of_last_fully_integrated_correlator::Int + sc_bit_detector::SecondaryCodeOrBitDetector + prompts_buffer::PromptsBuffer + bit_buffer::BitBuffer +end + +get_prn(s::SatState) = s.prn +get_num_ants(s::SatState{<: AbstractCorrelator{M}}) where M = M +get_code_phase(s::SatState) = s.code_phase +get_code_doppler(s::SatState) = s.code_doppler +get_carrier_phase(s::SatState) = s.carrier_phase * 2π +get_carrier_doppler(s::SatState) = s.carrier_doppler +get_integrated_samples(s::SatState) = s.integrated_samples +get_correlator(s::SatState) = s.correlator +get_last_fully_integrated_correlator(s::SatState) = s.last_fully_integrated_correlator +get_last_fully_integrated_filtered_prompt(s::SatState) = s.last_fully_integrated_filtered_prompt +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_bit_buffer(s::SatState) = s.bit_buffer + + +function SatState( + system::AbstractGNSS, + prn::Int, + sampling_frequency, + code_phase, + carrier_doppler; + num_ants::NumAnts = NumAnts(1), + carrier_phase = 0.0, + code_doppler = carrier_doppler * get_code_center_frequency_ratio(system), + num_prompts_buffer::Int = 20, +) + SatState( + prn, + float(code_phase), + float(code_doppler), + float(carrier_phase) / 2π, + float(carrier_doppler), + 0, + get_default_correlator(system, sampling_frequency, num_ants), + get_default_correlator(system, sampling_frequency, num_ants), + complex(0.0, 0.0), + -1, + SecondaryCodeOrBitDetector(), + PromptsBuffer(num_prompts_buffer), + BitBuffer() + ) +end + +function SatState( + acq::AcquisitionResults; + args... +) + SatState( + acq.system, + acq.prn, + acq.sampling_frequency, + acq.code_phase, + acq.carrier_doppler; + args... + ) +end + +struct SystemSatsState{ + S <: AbstractGNSS, + SS <: SatState +} + system::S + states::Vector{SS} +end + +get_system(sss::SystemSatsState) = sss.system +get_states(sss::SystemSatsState) = sss.states \ No newline at end of file diff --git a/src/track.jl b/src/track.jl new file mode 100644 index 0000000..6c1b462 --- /dev/null +++ b/src/track.jl @@ -0,0 +1,60 @@ +function track( + signal::AbstractVecOrMat, + track_state::TrackState, + sampling_frequency; + intermediate_frequency = 0.0Hz, + preferred_num_code_blocks_to_integrate = 1, + min_integration_time = 0.75ms +) + system_sats_states = track_state.system_sats_states + doppler_estimator = track_state.doppler_estimator + sat_sample_params = init_sample_params( + system_sats_states, + preferred_num_code_blocks_to_integrate, + ) + num_samples_signal = get_num_samples(signal) + while true + sat_sample_params = calc_sample_params( + system_sats_states, + sat_sample_params, + num_samples_signal, + sampling_frequency, + preferred_num_code_blocks_to_integrate, + ) + all(x -> all(y -> y.signal_samples_to_integrate == 0, x), sat_sample_params) && break + correlators = downconvert_and_correlate( + track_state.downconvert_and_correlator, + signal, + sampling_frequency, + intermediate_frequency, + system_sats_states, + sat_sample_params, + ) + system_sats_states = update( + system_sats_states, + sampling_frequency, + intermediate_frequency, + correlators, + sat_sample_params, + num_samples_signal + ) + doppler_estimator, dopplers_and_prompts = estimate_dopplers_and_filter_prompt( + doppler_estimator, + system_sats_states, + sat_sample_params, + sampling_frequency, + min_integration_time + ) + system_sats_states = update( + system_sats_states, + dopplers_and_prompts, + sat_sample_params, + ) + end + TrackState( + system_sats_states, + doppler_estimator, + track_state.downconvert_and_correlator, + track_state.num_samples + ) +end \ No newline at end of file diff --git a/src/tracking_loop.jl b/src/tracking_loop.jl deleted file mode 100644 index 0357746..0000000 --- a/src/tracking_loop.jl +++ /dev/null @@ -1,515 +0,0 @@ -""" -$(SIGNATURES) - -Track the signal `signal` based on the current tracking `state`, the sampling frequency -`sampling_frequency` and PRN `prn`. Optional configurations are: -- Post correlation filter `post_corr_filter` defaults to `get_default_post_corr_filter(...)` -- Intermediate frequency `intermediate_frequency` defaults to 0Hz -- Maximal integration time `max_integration_time` defaults to 1ms. The actual integration - time can be less in order to find the secondary code or the bit shift. Once this is found - the integration time is identical to the maximal integration time. -- Minimal integration time `min_integration_time` defaults to 0.75ms. It's the minimal - integration time, that leads to a valid correlation result. It is only used for the first - integration period. -- Sample shift of the correlator replica `correlator_sample_shifts` defaults to - `get_correlator_sample_shifts(...)` -- Bandwidth of the carrier loop `carrier_loop_filter_bandwidth` defaults to 18Hz -- Bandwidth of the code loop `code_loop_filter_bandwidth` defaults to 1Hz -- Velocity aiding `velocity_aiding` defaults to 0Hz -""" -function track( - signal, - state::TrackingState{S, C, CALF, COLF, CN, DS, CAR, COR, PCF}, - sampling_frequency; - intermediate_frequency = 0.0Hz, - max_integration_time::typeof(1ms) = 1ms, - min_integration_time::typeof(1.0ms) = 0.75ms, - correlator_sample_shifts = get_correlator_sample_shifts(get_system(state), - get_correlator(state), sampling_frequency, 0.5), - early_late_index_shift = get_early_late_index_shift(get_system(state), - correlator_sample_shifts, get_correlator(state), sampling_frequency, 0.5), - carrier_loop_filter_bandwidth = 18Hz, - code_loop_filter_bandwidth = 1Hz, - velocity_aiding = 0Hz, -) where { - S <: AbstractGNSS, - C <: AbstractCorrelator, - CALF <: AbstractLoopFilter, - COLF <: AbstractLoopFilter, - CN <: AbstractCN0Estimator, - DS, - CAR, - COR, - PCF <: AbstractPostCorrFilter, -} - prn = get_prn(state) - system = get_system(state) - if get_data_frequency(system) != 0Hz - @assert rem(1 / get_data_frequency(system), max_integration_time) == 0ms - end - correlator = get_correlator(state) - num_ants = get_num_ants(correlator) - size(signal, 2) == num_ants || throw(ArgumentError("The second dimension of the signal should be equal to the number of antennas specified by num_ants = NumAnts(N) in the TrackingState.")) - if typeof(get_codes(system)) <: CuMatrix - typeof(signal) <: StructArray || throw(ArgumentError("Signal is not a StructArray, initialize the signal properly and try again.")) - typeof(signal.re) <: CuArray && typeof(get_codes(system)) <: CuArray || throw(ArgumentError("Signal and GNSS codes are not of the same type. Please check if CPU or GPU is used.")) - end - downconverted_signal_temp = get_downconverted_signal(state) - downconverted_signal = resize!(downconverted_signal_temp, size(signal, 1), signal) - carrier_replica = get_carrier(state) - resize!(choose(carrier_replica, signal), size(signal, 1)) - code_replica = get_code(state) - resize!(code_replica, size(signal, 1) + correlator_sample_shifts[end]-correlator_sample_shifts[1]) - init_carrier_doppler = get_init_carrier_doppler(state) - init_code_doppler = get_init_code_doppler(state) - carrier_doppler = get_carrier_doppler(state) - code_doppler = get_code_doppler(state) - carrier_phase = get_carrier_phase(state) - code_phase = get_code_phase(state) - sc_bit_detector = get_sc_bit_detector(state) - carrier_loop_filter = get_carrier_loop_filter(state) - code_loop_filter = get_code_loop_filter(state) - prompt_accumulator = get_prompt_accumulator(state) - integrated_samples = get_integrated_samples(state) - cn0_estimator = get_cn0_estimator(state) - post_corr_filter = get_post_corr_filter(state) - signal_start_sample = 1 - bit_buffer = BitBuffer() - filtered_prompt = zero(ComplexF64) - valid_correlator = zero(correlator) - valid_correlator_carrier_phase = 0.0 - valid_correlator_carrier_frequency = 0.0Hz - got_correlator = false - while true - num_samples_left_to_integrate = get_num_samples_left_to_integrate( - system, - max_integration_time, - sampling_frequency, - code_doppler, - code_phase, - found(sc_bit_detector) - ) - signal_samples_left = get_num_samples(signal) - signal_start_sample + 1 - num_samples_left = min(num_samples_left_to_integrate, signal_samples_left) - carrier_frequency = get_current_carrier_frequency( - intermediate_frequency, - carrier_doppler - ) - code_frequency = get_current_code_frequency(system, code_doppler) - correlator = downconvert_and_correlate!( - system, - signal, - correlator, - code_replica, - code_phase, - carrier_replica, - carrier_phase, - downconverted_signal, - code_frequency, - correlator_sample_shifts, - carrier_frequency, - sampling_frequency, - signal_start_sample, - num_samples_left, - prn - ) - integrated_samples += num_samples_left - carrier_phase = update_carrier_phase( - num_samples_left, - carrier_frequency, - sampling_frequency, - carrier_phase, - ) - prev_code_phase = code_phase - code_phase = update_code_phase( - system, - num_samples_left, - code_frequency, - sampling_frequency, - code_phase, - found(sc_bit_detector) - ) - integration_time = integrated_samples / sampling_frequency - if num_samples_left == num_samples_left_to_integrate && - integration_time >= min_integration_time - got_correlator = true - correlator = normalize(correlator, integrated_samples) - valid_correlator = correlator - valid_correlator_carrier_phase = carrier_phase - valid_correlator_carrier_frequency = carrier_frequency - post_corr_filter = update(post_corr_filter, get_prompt(correlator, correlator_sample_shifts)) - filtered_correlator = apply(post_corr_filter, correlator) - filtered_prompt = get_prompt(filtered_correlator, correlator_sample_shifts) - pll_discriminator = pll_disc( - system, - filtered_correlator, - correlator_sample_shifts - ) - dll_discriminator = dll_disc( - system, - filtered_correlator, - correlator_sample_shifts, - early_late_index_shift, - code_frequency / sampling_frequency - ) - carrier_freq_update, carrier_loop_filter = filter_loop( - carrier_loop_filter, - pll_discriminator, - integration_time, - carrier_loop_filter_bandwidth - ) - code_freq_update, code_loop_filter = filter_loop( - code_loop_filter, - dll_discriminator, - integration_time, - code_loop_filter_bandwidth - ) - carrier_doppler, code_doppler = aid_dopplers( - system, - init_carrier_doppler, - init_code_doppler, - carrier_freq_update, - code_freq_update, - velocity_aiding - ) - cn0_estimator = update( - cn0_estimator, - get_prompt(filtered_correlator, correlator_sample_shifts) - ) - bit_buffer, prompt_accumulator = buffer( - system, - bit_buffer, - prompt_accumulator, - found(sc_bit_detector), - prev_code_phase, - code_phase, - max_integration_time, - get_prompt(filtered_correlator, correlator_sample_shifts) - ) - sc_bit_detector = find( - system, - sc_bit_detector, - get_prompt(filtered_correlator, correlator_sample_shifts) - ) - correlator = zero(correlator) - integrated_samples = 0 - end - - num_samples_left == signal_samples_left && break - signal_start_sample += num_samples_left - end - next_state = TrackingState{S, C, CALF, COLF, CN, DS, CAR, COR, PCF}( - prn, - system, - init_carrier_doppler, - init_code_doppler, - carrier_doppler, - code_doppler, - carrier_phase, - code_phase, - correlator, - carrier_loop_filter, - code_loop_filter, - sc_bit_detector, - integrated_samples, - prompt_accumulator, - cn0_estimator, - downconverted_signal, - carrier_replica, - code_replica, - post_corr_filter, - ) - estimated_cn0 = estimate_cn0(cn0_estimator, max_integration_time) - TrackingResults( - next_state, - valid_correlator, - filtered_prompt, - correlator_sample_shifts, - early_late_index_shift, - valid_correlator_carrier_frequency, - valid_correlator_carrier_phase, - got_correlator, - bit_buffer, - estimated_cn0 - ) -end - -function downconvert_and_correlate!( - system, - signal, - correlator, - code_replica, - code_phase, - carrier_replica, - carrier_phase, - downconverted_signal, - code_frequency, - correlator_sample_shifts, - carrier_frequency, - sampling_frequency, - signal_start_sample, - num_samples_left, - prn -) - gen_code_replica!( - code_replica, - system, - code_frequency, - sampling_frequency, - code_phase, - signal_start_sample, - num_samples_left, - correlator_sample_shifts, - prn - ) - gen_carrier_replica!( - choose(carrier_replica, signal), - carrier_frequency, - sampling_frequency, - carrier_phase, - signal_start_sample, - num_samples_left - ) - downconvert!( - choose(downconverted_signal, signal), - signal, - choose(carrier_replica, signal), - signal_start_sample, - num_samples_left - ) - correlate( - correlator, - choose(downconverted_signal, signal), - code_replica, - correlator_sample_shifts, - signal_start_sample, - num_samples_left - ) -end - -# CUDA downconvert_and_correlate for num_ants > 1 -function downconvert_and_correlate!( - system::AbstractGNSS{C}, - signal::AbstractMatrix, - correlator::T, - code_replica, - code_phase, - carrier_replica, - carrier_phase, - downconverted_signal, - code_frequency, - correlator_sample_shifts, - carrier_frequency, - sampling_frequency, - signal_start_sample, - num_samples_left, - prn -) where {C <: CuMatrix, T <: AbstractCorrelator} - accumulator_result = downconvert_and_correlate_kernel_wrapper( - system, - view(signal, signal_start_sample:signal_start_sample - 1 + num_samples_left,:), - correlator, - code_phase, - carrier_phase, - code_frequency, - correlator_sample_shifts, - carrier_frequency, - sampling_frequency, - signal_start_sample, - num_samples_left, - prn - ) - return T(map(+, get_accumulators(correlator), eachcol(Array(accumulator_result[1,:,:])))) -end - -# CUDA downconvert_and_correlate for num_ants = 1 -function downconvert_and_correlate!( - system::AbstractGNSS{C}, - signal::AbstractVector, - correlator::T, - code_replica, - code_phase, - carrier_replica, - carrier_phase, - downconverted_signal, - code_frequency, - correlator_sample_shifts, - carrier_frequency, - sampling_frequency, - signal_start_sample, - num_samples_left, - prn -) where {C <: CuMatrix, T <: AbstractCorrelator} - accumulator_result = downconvert_and_correlate_kernel_wrapper( - system, - view(signal, signal_start_sample:signal_start_sample - 1 + num_samples_left), - correlator, - code_phase, - carrier_phase, - code_frequency, - correlator_sample_shifts, - carrier_frequency, - sampling_frequency, - signal_start_sample, - num_samples_left, - prn - ) - addition(a,b) = a + first(b) - return T(map(addition, get_accumulators(correlator), eachcol(Array(accumulator_result[1,:,:])))) -end - -function choose(replica::CarrierReplicaCPU, signal::AbstractArray{Complex{Float64}}) - replica.carrier_f64 -end -function choose(replica::CarrierReplicaCPU, signal::AbstractArray{Complex{T}}) where T <: Number - replica.carrier_f32 -end -function choose(replica::DownconvertedSignalCPU, signal::AbstractArray{Complex{Float64}}) - replica.downconverted_signal_f64 -end -function choose(replica::DownconvertedSignalCPU, signal::AbstractArray{Complex{T}}) where T <: Number - replica.downconverted_signal_f32 -end -function choose(replica::Nothing, signal::AbstractArray) - nothing -end - -""" -$(SIGNATURES) - -Default post corr filter -""" -function get_default_post_corr_filter(correlator::AbstractCorrelator{T}) where T <: SVector - x -> x[1] -end - -function get_default_post_corr_filter(correlator) - x -> x -end - -""" -$(SIGNATURES) - -Returns the appropiate integration time. It will be the maximum integration time once the -secondary code or the bit shift has been found. -""" -function get_integration_time( - system::AbstractGNSS, - max_integration_time, - secondary_code_or_bit_found::Bool -) - ifelse( - secondary_code_or_bit_found, - max_integration_time, - min( - ceil(typeof(1ms), get_code_length(system) / get_code_frequency(system)), - max_integration_time - ) - ) -end - -""" -$(SIGNATURES) - -Calculates the number of chips to integrate. -""" -function get_num_chips_to_integrate( - system::AbstractGNSS, - max_integration_time, - current_code_phase, - secondary_code_or_bit_found -) - max_phase = Int(upreferred(get_code_frequency(system) * - get_integration_time(system, max_integration_time, secondary_code_or_bit_found))) - current_phase_mod_max_phase = mod(current_code_phase, max_phase) - return max_phase - current_phase_mod_max_phase -end - -""" -$(SIGNATURES) - -Calculates the number of samples to integrate. -""" -function get_num_samples_left_to_integrate( - system::AbstractGNSS, - max_integration_time, - sampling_frequency, - current_code_doppler, - current_code_phase, - secondary_code_or_bit_found -) - phase_to_integrate = get_num_chips_to_integrate( - system, - max_integration_time, - current_code_phase, - secondary_code_or_bit_found - ) - code_frequency = get_code_frequency(system) + current_code_doppler - ceil(Int, phase_to_integrate * sampling_frequency / code_frequency) -end - -""" -$(SIGNATURES) - -Aid dopplers. That is velocity aiding for the carrier doppler and carrier aiding -for the code doppler. -""" -function aid_dopplers( - system::AbstractGNSS, - init_carrier_doppler, - init_code_doppler, - carrier_freq_update, - code_freq_update, - velocity_aiding -) - carrier_doppler = carrier_freq_update + velocity_aiding - code_doppler = code_freq_update + carrier_doppler * get_code_center_frequency_ratio(system) - init_carrier_doppler + carrier_doppler, init_code_doppler + code_doppler -end - -""" -$(SIGNATURES) - -Get the number of samples in the signal. -""" -@inline function get_num_samples(signal) - length(signal) -end - -@inline function get_num_samples(signal::AbstractMatrix) - size(signal, 1) -end - -function resize!(ds::DownconvertedSignalCPU, b::Integer, signal::AbstractVector) - resize!(choose(ds, signal), b) - ds -end - -function resize!(ds::DownconvertedSignalCPU, b::Integer, signal::AbstractMatrix{Complex{Float64}}) - num_ants = size(signal, 2) - DownconvertedSignalCPU( - ds.downconverted_signal_f32, - size(ds.downconverted_signal_f64, 1) == b ? - ds.downconverted_signal_f64 : - StructArray{Complex{Float64}}((Matrix{Float64}(undef, b, num_ants), Matrix{Float64}(undef, b, num_ants))) - ) -end - -function resize!(ds::DownconvertedSignalCPU, b::Integer, signal::AbstractMatrix{Complex{T}}) where T <: Number - num_ants = size(signal, 2) - DownconvertedSignalCPU( - size(ds.downconverted_signal_f32, 1) == b ? - ds.downconverted_signal_f32 : - StructArray{Complex{Float32}}((Matrix{Float32}(undef, b, num_ants), Matrix{Float32}(undef, b, num_ants))), - ds.downconverted_signal_f64 - ) -end - -# No need for resizing when dealing with GPU signals -function resize!(ds::Nothing, b::Integer, signal::AbstractArray) - return ds -end -# No need for resizing the GPU GNSS codes -function resize!(codes::Nothing, b::Integer) - return codes -end \ No newline at end of file diff --git a/src/tracking_results.jl b/src/tracking_results.jl deleted file mode 100644 index 53bcd9b..0000000 --- a/src/tracking_results.jl +++ /dev/null @@ -1,177 +0,0 @@ -""" -$(SIGNATURES) - -TrackingResults that hold the correlation results. -""" -struct TrackingResults{ - TS <: TrackingState, - C <: AbstractCorrelator, - CS, - ELI - } - state::TS - correlator::C - filtered_prompt::ComplexF64 - correlator_sample_shifts::CS - early_late_index_shift::ELI - correlator_carrier_frequency::typeof(1.0Hz) - correlator_carrier_phase::Float64 - got_correlator::Bool - bit_buffer::BitBuffer - cn0::typeof(1.0dBHz) -end - -""" -$(SIGNATURES) - -Get state of the tracking result. -""" -@inline get_state(results::TrackingResults) = results.state - -""" -$(SIGNATURES) - -Get the carrier doppler of the tracking result. -""" -@inline get_carrier_doppler(results::TrackingResults) = get_carrier_doppler(results.state) - -""" -$(SIGNATURES) - -Get the carrier phase of the tracking result. -""" -@inline get_carrier_phase(results::TrackingResults) = get_carrier_phase(results.state) * 2π - -""" -$(SIGNATURES) - -Get the code doppler of the tracking result. -""" -@inline get_code_doppler(results::TrackingResults) = get_code_doppler(results.state) - -""" -$(SIGNATURES) - -Get the code phase of the tracking result. -""" -@inline get_code_phase(results::TrackingResults) = get_code_phase(results.state) - -""" -$(SIGNATURES) - -Get the correlator of the tracking result. -""" -@inline get_correlator(results::TrackingResults) = results.correlator - -""" -$(SIGNATURES) - -Get the correlator sample shifts of the tracking result. -""" -@inline get_correlator_sample_shifts(results::TrackingResults) = results.correlator_sample_shifts - -""" -$(SIGNATURES) - -Get the correlator sample shifts of the tracking result. -""" -@inline get_early_late_index_shift(results::TrackingResults) = results.early_late_index_shift - -""" -$(SIGNATURES) - -Get carrier phase at the same timestamp when the correlator is created of the tracking -result. -""" -@inline function get_correlator_carrier_phase(results::TrackingResults) - results.correlator_carrier_phase * 2π -end - -""" -$(SIGNATURES) - -Get carrier frequency at the same timestamp when the correlator is created of the tracking -result. -""" -@inline function get_correlator_carrier_frequency(results::TrackingResults) - results.correlator_carrier_frequency -end - -""" -$(SIGNATURES) - -Get the early of the tracking result. -""" -@inline get_early(results::TrackingResults) = get_early( - get_correlator(results), - get_correlator_sample_shifts(results), - get_early_late_index_shift(results) -) - -""" -$(SIGNATURES) - -Get the prompt of the tracking result. -""" -@inline get_prompt(results::TrackingResults) = get_prompt(get_correlator(results), get_correlator_sample_shifts(results)) - -""" -$(SIGNATURES) - -Get the late of the tracking result. -""" -@inline get_late(results::TrackingResults) = get_late( - get_correlator(results), - get_correlator_sample_shifts(results), - get_early_late_index_shift(results) -) - -""" -$(SIGNATURES) - -Get the bits of the tracking result. -""" -@inline get_bits(results::TrackingResults) = get_bits(results.bit_buffer) - -""" -$(SIGNATURES) - -Get the number of bits of the tracking result. -""" -@inline get_num_bits(results::TrackingResults) = length(results.bit_buffer) - -""" -$(SIGNATURES) - -Get the carrier to noise density ratio of the tracking result. -""" -@inline get_cn0(results::TrackingResults) = results.cn0 - -""" -$(SIGNATURES) - -Check if the secondary code or bit has been found. -""" -@inline get_secondary_code_or_bit_found(results::TrackingResults) = - found(get_sc_bit_detector(results.state)) - -""" -$(SIGNATURES) - -Get used system. -""" -@inline get_system(results::TrackingResults) = get_system(get_state(results)) - -""" -$(SIGNATURES) - -Get filtered prompt. -""" -@inline get_filtered_prompt(results::TrackingResults) = results.filtered_prompt - -""" -$(SIGNATURES) - -Get Post correlation filter. -""" -@inline get_post_corr_filter(results::TrackingResults) = get_post_corr_filter(get_state(results)) \ No newline at end of file diff --git a/src/tracking_state.jl b/src/tracking_state.jl index 44735d4..b444048 100644 --- a/src/tracking_state.jl +++ b/src/tracking_state.jl @@ -1,280 +1,152 @@ -""" -$(SIGNATURES) +struct TrackState{ + S <: TupleLike{<:NTuple{N, SystemSatsState}} where N, + DE <: AbstractDopplerEstimator, + DC <: AbstractDownconvertAndCorrelator, +} + system_sats_states::S + doppler_estimator::DE + downconvert_and_correlator::DC + num_samples::Int +end -CarrierReplicaCPU that holds possible carrier representations. -The type is unknown, when TrackingState is initialized. -It is determined based on the signal, which is given in the track -function. -""" -struct CarrierReplicaCPU{ - CF32 <: StructArray{ComplexF32}, - CF64 <: StructArray{ComplexF64} - } - carrier_f32::CF32 - carrier_f64::CF64 +function TrackState( + system::AbstractGNSS, + sat_states::Vector{<: SatState}; + num_samples, + doppler_estimator = map(sat_state -> ConventionalPLLAndDLL(sat_state), sat_states), + downconvert_and_correlator = CPUDownconvertAndCorrelator((SystemSatsState(system, sat_states),), num_samples) +) + system_sats_stats = (SystemSatsState(system, sat_states),) + TrackState( + system_sats_stats, + ConventionalPLLsAndDLLs((doppler_estimator,)), + downconvert_and_correlator, + num_samples + ) end -function CarrierReplicaCPU() - CarrierReplicaCPU( - StructArray{Complex{Float32}}(undef, 0), - StructArray{Complex{Float64}}(undef, 0) +function TrackState( + system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}; + num_samples, + doppler_estimators = ConventionalPLLsAndDLLs(map(sat_states -> map(sat_state -> ConventionalPLLAndDLL(sat_state), sat_states.states), system_sats_states)), + downconvert_and_correlator = CPUDownconvertAndCorrelator(system_sats_states, num_samples), +) where N + TrackState( + system_sats_states, + doppler_estimators, + downconvert_and_correlator, + num_samples ) end -""" -$(SIGNATURES) -DownconvertedSignalCPU that holds possible downconverted signal representations. -The type is unknown, when TrackingState is initialized. -It is determined based on the signal, which is given in the track -function. -""" -struct DownconvertedSignalCPU{ - DSF32 <: StructArray{ComplexF32}, - DSF64 <: StructArray{ComplexF64} - } - downconverted_signal_f32::DSF32 - downconverted_signal_f64::DSF64 +function get_sat_states( + track_state::TrackState{<:TupleLike{<:NTuple{N, SystemSatsState}}}, + system_idx::Union{Symbol, Integer} +) where N + track_state.system_sats_states[system_idx].states end -function DownconvertedSignalCPU(num_ants::NumAnts{1}) - DownconvertedSignalCPU( - StructArray{Complex{Float32}}(undef, 0), - StructArray{Complex{Float64}}(undef, 0) - ) +function get_sat_states( + track_state::TrackState{<:TupleLike{<:NTuple{N, SystemSatsState}}}, + system_idx::Val{M} +) where {N, M} + track_state.system_sats_states[M].states end -function DownconvertedSignalCPU(num_ants::NumAnts{N}) where N - DownconvertedSignalCPU( - StructArray{Complex{Float32}}(undef, 0, N), - StructArray{Complex{Float64}}(undef, 0, N) - ) +function get_sat_states( + track_state::TrackState{<:TupleLike{<:NTuple{1, SystemSatsState}}}, +) + get_sat_states(track_state, 1) end -""" -$(SIGNATURES) +function get_system( + track_state::TrackState{<:TupleLike{<:NTuple{N, SystemSatsState}}}, + system_idx::Union{Symbol, Integer} +) where N + track_state.system_sats_states[system_idx].system +end -TrackingState that holds the tracking state. -""" -struct TrackingState{ - S <: AbstractGNSS, - C <: AbstractCorrelator, - CALF <: AbstractLoopFilter, - COLF <: AbstractLoopFilter, - CN <: AbstractCN0Estimator, - DS <: Union{DownconvertedSignalCPU, Nothing}, - CAR <: Union{CarrierReplicaCPU, Nothing}, - COR <: Union{Vector{<:Real}, Nothing}, - PCF <: AbstractPostCorrFilter, - } - prn::Int - system::S - init_carrier_doppler::typeof(1.0Hz) - init_code_doppler::typeof(1.0Hz) - carrier_doppler::typeof(1.0Hz) - code_doppler::typeof(1.0Hz) - carrier_phase::Float64 - code_phase::Float64 - correlator::C - carrier_loop_filter::CALF - code_loop_filter::COLF - sc_bit_detector::SecondaryCodeOrBitDetector - integrated_samples::Int - prompt_accumulator::ComplexF64 - cn0_estimator::CN - downconverted_signal::DS - carrier::CAR - code::COR - post_corr_filter::PCF +function get_system( + track_state::TrackState{<:TupleLike{<:NTuple{N, SystemSatsState}}}, + system_idx::Val{M} +) where {N, M} + track_state.system_sats_states[M].system end -function TrackingState( - track_state::TrackingState; - post_corr_filter, +function get_system( + track_state::TrackState{<:TupleLike{<:NTuple{1, SystemSatsState}}}, ) - TrackingState( - track_state.prn, - track_state.system, - track_state.init_carrier_doppler, - track_state.init_code_doppler, - track_state.carrier_doppler, - track_state.code_doppler, - track_state.carrier_phase, - track_state.code_phase, - track_state.correlator, - track_state.carrier_loop_filter, - track_state.code_loop_filter, - track_state.sc_bit_detector, - track_state.integrated_samples, - track_state.prompt_accumulator, - track_state.cn0_estimator, - track_state.downconverted_signal, - track_state.carrier, - track_state.code, - post_corr_filter - ) + get_system(track_state, 1) end -""" -$(SIGNATURES) - -Convenient TrackingState constructor. Mandatory parameters are the GNSS system, the -carrier doppler `carrier_doppler` and the code phase `code_phase`. Optional parameters are -- Code doppler `code_doppler`, that defaults to carrier doppler multiplied with the code / - center frequency ratio -- Carrier phase `carrier_phase`, that defaults to 0 -- Carrier loop filter `carrier_loop_filter`, that defaults to the `ThirdOrderBilinarLF()` -- Code loop filter `code_loop_filter`, that defaults to the `SecondOrderBilinearLF()` -- Number of antenna elemants `num_ants`, that defaults to `NumAnts(1)` -- Correlator `correlator`, that defaults to `get_default_correlator(...)` -- Integrated samples `integrated_samples`, that defaults to 0 -- Prompt accumulator for bit detection `prompt_accumulator`, that defaults to 0 -- CN0 estimator `cn0_estimator`, that defaults to `MomentsCN0Estimator(20)` -""" -function TrackingState( - prn::Integer, - system::S, - carrier_doppler, - code_phase; - code_doppler = carrier_doppler * get_code_center_frequency_ratio(system), - carrier_phase = 0.0, - carrier_loop_filter::CALF = ThirdOrderBilinearLF(), - code_loop_filter::COLF = SecondOrderBilinearLF(), - sc_bit_detector = SecondaryCodeOrBitDetector(), - num_ants = NumAnts(1), - correlator::C = get_default_correlator(system, num_ants), - integrated_samples = 0, - prompt_accumulator = zero(ComplexF64), - cn0_estimator::CN = MomentsCN0Estimator(20), - post_corr_filter::PCF = DefaultPostCorrFilter(), - num_samples = 0 -) where { - T <: AbstractMatrix, - S <: AbstractGNSS{T}, - C <: AbstractCorrelator, - CALF <: AbstractLoopFilter, - COLF <: AbstractLoopFilter, - CN <: AbstractCN0Estimator, - PCF <: AbstractPostCorrFilter -} - if found(sc_bit_detector) - code_phase = mod(code_phase, get_code_length(system) * - get_secondary_code_length(system)) - else - code_phase = mod(code_phase, get_code_length(system)) - end - downconverted_signal = DownconvertedSignalCPU(num_ants) - carrier = CarrierReplicaCPU() - code = Vector{get_code_type(system)}(undef, 0) +function get_sat_state( + track_state::TrackState{<:TupleLike{<:NTuple{N, SystemSatsState}}}, + system_idx::Union{Symbol, Integer, Val}, + sat_idx::Integer +) where N + get_sat_states(track_state, system_idx)[sat_idx] +end - TrackingState{S, C, CALF, COLF, CN, typeof(downconverted_signal), typeof(carrier), typeof(code), PCF}( - prn, - system, - carrier_doppler, - code_doppler, - carrier_doppler, - code_doppler, - carrier_phase / 2π, - code_phase, - correlator, - carrier_loop_filter, - code_loop_filter, - sc_bit_detector, - integrated_samples, - prompt_accumulator, - cn0_estimator, - downconverted_signal, - carrier, - code, - post_corr_filter, - ) +function get_sat_state( + track_state::TrackState{<:TupleLike{<:NTuple{1, SystemSatsState}}}, + sat_idx::Integer +) + get_sat_states(track_state, 1)[sat_idx] end -# CUDA dispatch -function TrackingState( - prn::Integer, - system::S, - carrier_doppler, - code_phase; - num_samples, - code_doppler = carrier_doppler * get_code_center_frequency_ratio(system), - carrier_phase = 0.0, - carrier_loop_filter::CALF = ThirdOrderBilinearLF(), - code_loop_filter::COLF = SecondOrderBilinearLF(), - sc_bit_detector = SecondaryCodeOrBitDetector(), - num_ants = NumAnts(1), - correlator::C = get_default_correlator(system, num_ants), - integrated_samples = 0, - prompt_accumulator = zero(ComplexF64), - cn0_estimator::CN = MomentsCN0Estimator(20), - post_corr_filter::PCF = DefaultPostCorrFilter(), -) where { - T <: CuMatrix, - S <: AbstractGNSS{T}, - C <: AbstractCorrelator, - CALF <: AbstractLoopFilter, - COLF <: AbstractLoopFilter, - CN <: AbstractCN0Estimator, - PCF <: AbstractPostCorrFilter -} - if found(sc_bit_detector) - code_phase = mod(code_phase, get_code_length(system) * - get_secondary_code_length(system)) +create_buffer(buffers::Vector{<:CPUBuffers}, system, sat_states, num_samples) = + CPUBuffers(system, sat_states, num_samples) + +create_buffer(buffers::Vector{<:GPUBuffers}, system, sat_states, num_samples) = + GPUBuffers(system, sat_states, num_samples) + +function add_sats!( + track_state::TrackState{<:TupleLike{<:NTuple{N, SystemSatsState}}, <: ConventionalPLLsAndDLLs}, + system_idx::Union{Symbol, Integer}, + system::AbstractGNSS, + sat_states::Union{SatState, Vector{<:SatState}}, +) where N + system_sats_state = track_state.system_sats_states[system_idx] + pll_and_dlls = track_state.doppler_estimator.plls_and_dlls[system_idx] + buffers = track_state.downconvert_and_correlator.buffers[system_idx] + if sat_states isa Vector + append!(system_sats_state.states, sat_states) + append!(pll_and_dlls, map(sat_state -> ConventionalPLLAndDLL(sat_state), sat_states)) + append!(buffers, map(sat_state -> create_buffer(buffers, system, sat_state, track_state.num_samples), sat_states)) else - code_phase = mod(code_phase, get_code_length(system)) + push!(system_sats_state.states, sat_states) + push!(pll_and_dlls, ConventionalPLLAndDLL(sat_states)) + push!(buffers, create_buffer(buffers, system, sat_states, track_state.num_samples)) end - downconverted_signal = nothing - carrier = nothing - code = nothing + track_state +end - TrackingState{S, C, CALF, COLF, CN, Nothing, Nothing, Nothing, PCF}( - prn, - system, - carrier_doppler, - code_doppler, - carrier_doppler, - code_doppler, - carrier_phase / 2π, - code_phase, - correlator, - carrier_loop_filter, - code_loop_filter, - sc_bit_detector, - integrated_samples, - prompt_accumulator, - cn0_estimator, - downconverted_signal, - carrier, - code, - post_corr_filter, - ) +function add_sats!( + track_state::TrackState{<:TupleLike{<:NTuple{1, SystemSatsState}}, <: ConventionalPLLsAndDLLs}, + system::AbstractGNSS, + sat_states::Union{SatState, Vector{<:SatState}}, +) + add_sats!(track_state, 1, system, sat_states) end -function TrackingState(acq::AcquisitionResults; args...) - TrackingState( - acq.prn, - acq.system, - acq.carrier_doppler, - acq.code_phase; - args... - ) +function remove_sats!( + track_state::TrackState{<:TupleLike{<:NTuple{N, SystemSatsState}}, <: ConventionalPLLsAndDLLs}, + system_idx::Union{Symbol, Integer}, + prns::Union{Int, Vector{Int}}, +) where N + system_sats_state = track_state.system_sats_states[system_idx] + pll_and_dlls = track_state.doppler_estimator.plls_and_dlls[system_idx] + buffers = track_state.downconvert_and_correlator.buffers[system_idx] + idxs_to_delete = findall(state -> state.prn in prns, system_sats_state.states) + deleteat!(system_sats_state.states, idxs_to_delete) + deleteat!(pll_and_dlls, idxs_to_delete) + deleteat!(buffers, idxs_to_delete) + track_state end -@inline get_system(state::TrackingState) = state.system -@inline get_code_phase(state::TrackingState) = state.code_phase -@inline get_carrier_phase(state::TrackingState) = state.carrier_phase -@inline get_init_code_doppler(state::TrackingState) = state.init_code_doppler -@inline get_init_carrier_doppler(state::TrackingState) = state.init_carrier_doppler -@inline get_code_doppler(state::TrackingState) = state.code_doppler -@inline get_carrier_doppler(state::TrackingState) = state.carrier_doppler -@inline get_correlator(state::TrackingState) = state.correlator -@inline get_sc_bit_detector(state::TrackingState) = state.sc_bit_detector -@inline get_carrier_loop_filter(state::TrackingState) = state.carrier_loop_filter -@inline get_code_loop_filter(state::TrackingState) = state.code_loop_filter -@inline get_prompt_accumulator(state::TrackingState) = state.prompt_accumulator -@inline get_integrated_samples(state::TrackingState) = state.integrated_samples -@inline get_cn0_estimator(state::TrackingState) = state.cn0_estimator -@inline get_downconverted_signal(state::TrackingState) = state.downconverted_signal -@inline get_carrier(state::TrackingState) = state.carrier -@inline get_code(state::TrackingState) = state.code -@inline get_prn(state::TrackingState) = state.prn -@inline get_post_corr_filter(state::TrackingState) = state.post_corr_filter +function remove_sats!( + track_state::TrackState{<:TupleLike{<:NTuple{1, SystemSatsState}}, <: ConventionalPLLsAndDLLs}, + prns::Union{Int, Vector{Int}}, +) + remove_sats!(track_state, 1, prns) +end \ No newline at end of file diff --git a/src/update_sat_state.jl b/src/update_sat_state.jl new file mode 100644 index 0000000..9f7b7db --- /dev/null +++ b/src/update_sat_state.jl @@ -0,0 +1,107 @@ +function update( + system_sats_states::TupleLike{<:NTuple{N, SystemSatsState}}, + sampling_frequency, + intermediate_frequency, + new_correlators::TupleLike{<:NTuple{N, Vector{<:AbstractCorrelator}}}, + params::TupleLike{<:NTuple{N, Vector{SampleParams}}}, + num_samples_signal +) where N + map(new_correlators, params, system_sats_states) do system_correlators, system_params, system_sats + sats_state = map(system_correlators, system_params, system_sats.states) do new_correlator, sat_params, sat_state + carrier_frequency = sat_state.carrier_doppler + intermediate_frequency + code_frequency = sat_state.code_doppler + get_code_frequency(system_sats.system) + carrier_phase = update_carrier_phase( + sat_params.signal_samples_to_integrate, + carrier_frequency, + sampling_frequency, + sat_state.carrier_phase, + ) + code_phase = update_code_phase( + system_sats.system, + sat_params.signal_samples_to_integrate, + code_frequency, + sampling_frequency, + sat_state.code_phase, + found(sat_state.sc_bit_detector) + ) + sample_of_last_fully_integrated_correlator = sat_state.sample_of_last_fully_integrated_correlator - + (sat_params.signal_start_sample == 1 ? num_samples_signal : 0) + SatState( + sat_state.prn, + code_phase, + sat_state.code_doppler, + carrier_phase, + sat_state.carrier_doppler, + sat_state.integrated_samples + sat_params.signal_samples_to_integrate, + new_correlator, + sat_state.last_fully_integrated_correlator, + sat_state.last_fully_integrated_filtered_prompt, + sample_of_last_fully_integrated_correlator, + sat_state.sc_bit_detector, + sat_state.prompts_buffer, + sat_state.bit_buffer + ) + end + SystemSatsState( + system_sats.system, + sats_state + ) + end +end + +function update( + system_sats_state::TupleLike{<:NTuple{N, SystemSatsState}}, + dopplers_and_filtered_prompts::TupleLike{<:NTuple{N, Vector{<:Union{Nothing, DopplersAndFilteredPrompt}}}}, + sat_sample_params::TupleLike{<:NTuple{N, Vector{SampleParams}}}, +) where N + map(system_sats_state, dopplers_and_filtered_prompts, sat_sample_params) do system_sats, dopplers_and_filtered_prompts_per_system, sat_sample_params_per_system + sats_state = map(system_sats.states, dopplers_and_filtered_prompts_per_system, sat_sample_params_per_system) do state, dopplers_and_filtered_prompt, sample_params + return if dopplers_and_filtered_prompt.filtered_prompt != complex(0.0, 0.0) + 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 + ) + bit_buffer = buffer( + system_sats.system, + state.bit_buffer, + sample_params.num_code_blocks_to_integrate, + found(state.sc_bit_detector), + filtered_prompt + ) + sc_bit_detector = find( + system_sats.system, + state.sc_bit_detector, + filtered_prompt + ) + correlator = zero(state.correlator) + integrated_samples = 0 + + new_state = SatState( + state.prn, + state.code_phase, + code_doppler, + state.carrier_phase, + carrier_doppler, + integrated_samples, + correlator, + state.correlator, + dopplers_and_filtered_prompt.filtered_prompt, + sample_params.signal_start_sample + sample_params.signal_samples_to_integrate - 1, + sc_bit_detector, + prompts_buffer, + bit_buffer + ) + new_state + else + state + end + end + SystemSatsState( + system_sats.system, + sats_state + ) + end +end \ No newline at end of file diff --git a/test/bit_buffer.jl b/test/bit_buffer.jl index 0bd2f56..3514411 100644 --- a/test/bit_buffer.jl +++ b/test/bit_buffer.jl @@ -3,91 +3,52 @@ @test @inferred(get_bits(bit_buffer)) == 0 @test @inferred(Tracking.length(bit_buffer)) == 0 @test @allocated(Tracking.BitBuffer()) == 0 + @test bit_buffer.prompt_accumulator == 0.0 + 0.0im + @test bit_buffer.prompt_accumulator_integrated_code_blocks == 0 - prompt_accumulator = 0.0 + 0.0im secondary_code_or_bit_found = false - prev_code_phase = 0.0 - code_phase = 1023 prompt_correlator = 1.0 + 0.0im - integration_time = 5ms + integrated_code_blocks = 1 gpsl1 = GPSL1() - next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( + next_bit_buffer = @inferred Tracking.buffer( gpsl1, bit_buffer, - prompt_accumulator, + integrated_code_blocks, secondary_code_or_bit_found, - prev_code_phase, - code_phase, - integration_time, prompt_correlator ) - @test next_prompt_accumulator == 0.0 + 0.0im + @test next_bit_buffer.prompt_accumulator == 0.0 + 0.0im + @test next_bit_buffer.prompt_accumulator_integrated_code_blocks == 0 @test @inferred(get_bits(next_bit_buffer)) == 0 @test @inferred(Tracking.length(next_bit_buffer)) == 0 - prompt_accumulator = 20.0 + 0.0im + bit_buffer = Tracking.BitBuffer(0, 0, 20.0 + 0.0im, 1) secondary_code_or_bit_found = true - prev_code_phase = 0.0 - code_phase = 1023 - prompt_correlator = 1.0 + 0.0im - integration_time = 5ms - next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( + next_bit_buffer = @inferred Tracking.buffer( gpsl1, bit_buffer, - prompt_accumulator, + integrated_code_blocks, secondary_code_or_bit_found, - prev_code_phase, - code_phase, - integration_time, prompt_correlator ) - @test next_prompt_accumulator == 21.0 + 0.0im + @test next_bit_buffer.prompt_accumulator == 21.0 + 0.0im + @test next_bit_buffer.prompt_accumulator_integrated_code_blocks == 2 @test @inferred(get_bits(next_bit_buffer)) == 0 @test @inferred(Tracking.length(next_bit_buffer)) == 0 - prompt_accumulator = 20.0 + 0.0im - secondary_code_or_bit_found = true - prev_code_phase = 9000 - code_phase = 1023 - prompt_correlator = 1.0 + 0.0im - integration_time = 20ms - - next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( - gpsl1, - bit_buffer, - prompt_accumulator, - secondary_code_or_bit_found, - prev_code_phase, - code_phase, - integration_time, - prompt_correlator - ) - @test next_prompt_accumulator == 0.0 + 0.0im - @test @inferred(get_bits(next_bit_buffer)) == 1 - @test @inferred(Tracking.length(next_bit_buffer)) == 1 - - prompt_accumulator = 20.0 + 0.0im - secondary_code_or_bit_found = true - prev_code_phase = 0 - code_phase = 0 - prompt_correlator = 1.0 + 0.0im - integration_time = 20ms + integrated_code_blocks = 19 - next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( + next_bit_buffer = @inferred Tracking.buffer( gpsl1, bit_buffer, - prompt_accumulator, + integrated_code_blocks, secondary_code_or_bit_found, - prev_code_phase, - code_phase, - integration_time, prompt_correlator ) - @test next_prompt_accumulator == 0.0 + 0.0im + @test next_bit_buffer.prompt_accumulator == 0.0 + 0.0im + @test next_bit_buffer.prompt_accumulator_integrated_code_blocks == 0 @test @inferred(get_bits(next_bit_buffer)) == 1 @test @inferred(Tracking.length(next_bit_buffer)) == 1 - - end diff --git a/test/cn0_estimation.jl b/test/cn0_estimation.jl index 61f1b27..4610d27 100644 --- a/test/cn0_estimation.jl +++ b/test/cn0_estimation.jl @@ -1,26 +1,24 @@ -@testset "Update CN0 Estimator" begin +@testset "Prompts buffer" begin - cn0_estimator = MomentsCN0Estimator(20) - @test @inferred(Tracking.get_prompt_buffer(cn0_estimator)) == zero(SVector{20, ComplexF64}) - @test @inferred(Tracking.get_current_index(cn0_estimator)) == 0 - @test @inferred(Tracking.length(cn0_estimator)) == 0 + prompts_buffer = Tracking.PromptsBuffer(20) + @test @inferred(Tracking.get_prompt_buffer(prompts_buffer)) == zero(SVector{20, ComplexF64}) + @test @inferred(Tracking.get_current_index(prompts_buffer)) == 0 + @test @inferred(Tracking.length(prompts_buffer)) == 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 + 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 - cn0_estimator = 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 + 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 - cn0_estimator = MomentsCN0Estimator(ones(SVector{20, ComplexF64}), 19, 20) - @test @inferred(Tracking.get_current_index(cn0_estimator)) == 19 - @test @inferred(Tracking.length(cn0_estimator)) == 20 - - @test @allocated(Tracking.update(MomentsCN0Estimator(20), 1 + 2im)) == 0 + 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 end @testset "CN0 estimation" begin @@ -33,8 +31,8 @@ end prn = 1 range = 0:3999 start_carrier_phase = π / 2 - cn0_estimator = MomentsCN0Estimator(20) - correlator_sample_shifts = SVector(-2, 0, 2) + cn0_estimator = MomentsCN0Estimator() + prompts_buffer = Tracking.PromptsBuffer(20) start_sample = 1 num_samples = 4000 gpsl1 = GPSL1() @@ -51,24 +49,23 @@ end code_frequency .* (-2:4001) ./ sampling_frequency .+ start_code_phase, prn ) - correlator = EarlyPromptLateCorrelator() + correlator = EarlyPromptLateCorrelator(gpsl1, sampling_frequency) signal_struct = StructArray(signal) correlator = Tracking.correlate( correlator, signal_struct, code, - correlator_sample_shifts, start_sample, num_samples, ) - cn0_estimator = Tracking.update( - cn0_estimator, - get_prompt(correlator, correlator_sample_shifts) + prompts_buffer = Tracking.update( + prompts_buffer, + get_prompt(correlator) ) end - @test @inferred(Tracking.get_current_index(cn0_estimator)) == 20 - @test @inferred(Tracking.length(cn0_estimator)) == 20 - cn0_estimate = @inferred Tracking.estimate_cn0(cn0_estimator, 1ms) + @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 cn0_estimate ≈ 45dBHz atol = 2.05dBHz diff --git a/test/code_replica.jl b/test/code_replica.jl index eadbca7..b317b15 100644 --- a/test/code_replica.jl +++ b/test/code_replica.jl @@ -10,7 +10,7 @@ 2.0, 11, 2480, - SVector(-1, 0, 1), + -1:1, 1 ) @@ -27,7 +27,7 @@ 2.0, 11, 6480, - SVector(-1, 0, 1), + -1:1, 1 ) @@ -45,7 +45,7 @@ 2.0, 11, 2480, - SVector(-1, 0, 1), + -1:1, 1 ) diff --git a/test/conventional_pll_and_dll.jl b/test/conventional_pll_and_dll.jl new file mode 100644 index 0000000..edf06b2 --- /dev/null +++ b/test/conventional_pll_and_dll.jl @@ -0,0 +1,118 @@ +@testset "Doppler aiding" begin + gpsl1 = GPSL1() + init_carrier_doppler = 10Hz + init_code_doppler = 1Hz + carrier_freq_update = 2Hz + code_freq_update = -0.5Hz + + carrier_freq, code_freq = @inferred Tracking.aid_dopplers( + gpsl1, + init_carrier_doppler, + init_code_doppler, + carrier_freq_update, + code_freq_update, + ) + + @test carrier_freq == 10Hz + 2Hz + @test code_freq == 1Hz + 2Hz / 1540 - 0.5Hz +end + +@testset "Conventional PLL and DLL" begin + gpsl1 = GPSL1() + sat_state = SatState(gpsl1, 1, 5e6Hz, 10.0, 500.0Hz) + + pll_and_dll = @inferred ConventionalPLLAndDLL( + sat_state, + ) + + @test pll_and_dll.init_carrier_doppler == 500.0Hz + @test pll_and_dll.init_code_doppler == 500.0Hz * get_code_center_frequency_ratio(gpsl1) + @test pll_and_dll.carrier_loop_filter == ThirdOrderBilinearLF() + @test pll_and_dll.code_loop_filter == SecondOrderBilinearLF() + @test pll_and_dll.carrier_loop_filter_bandwidth == 18.0Hz + @test pll_and_dll.code_loop_filter_bandwidth == 1.0Hz + @test pll_and_dll.post_corr_filter == DefaultPostCorrFilter() + + sampling_frequency = 5e6Hz + + gpsl1 = GPSL1() + system_sats_state = SystemSatsState( + gpsl1, + [ + SatState( + 1, + 0.5, + 1.0Hz, + 0.01, + 100.0Hz, + 150, + EarlyPromptLateCorrelator(complex.([1000, 2000, 1000], 0.0), -1:1, 3, 2, 1), + EarlyPromptLateCorrelator(complex.(zeros(3), zeros(3)), -1:1, 3, 2, 1), + complex(0.0, 0.0), + 14, + Tracking.SecondaryCodeOrBitDetector(), + Tracking.PromptsBuffer(20), + Tracking.BitBuffer() + ), + ] + ) + + doppler_estimators = ConventionalPLLsAndDLLs(map(sat_states -> map(sat_state -> ConventionalPLLAndDLL(sat_state), sat_states.states), (system_sats_state,))) + + system_sats_sample_params = Tracking.init_sample_params( + (system_sats_state,), + 1, + ) + next_system_sats_sample_params = Tracking.calc_sample_params( + (system_sats_state,), + system_sats_sample_params, + 5000, + sampling_frequency, + 1, + ) + + next_doppler_estimators, dopplers_and_filtered_prompts = Tracking.estimate_dopplers_and_filter_prompt( + doppler_estimators, + (system_sats_state,), + next_system_sats_sample_params, + sampling_frequency, + 0.75ms + ) + + @test dopplers_and_filtered_prompts[1][1].carrier_doppler == 0.0Hz + @test dopplers_and_filtered_prompts[1][1].code_doppler == 0.0Hz + @test dopplers_and_filtered_prompts[1][1].filtered_prompt == 0.0 + + system_sats_state = SystemSatsState( + gpsl1, + [ + SatState( + 1, + 0.5, + 1.0Hz, + 0.01, + 100.0Hz, + 4500, + EarlyPromptLateCorrelator(complex.([1000, 2000, 1000], 0.0), -1:1, 3, 2, 1), + EarlyPromptLateCorrelator(complex.(zeros(3), zeros(3)), -1:1, 3, 2, 1), + complex(0.0, 0.0), + 14, + Tracking.SecondaryCodeOrBitDetector(), + Tracking.PromptsBuffer(20), + Tracking.BitBuffer() + ), + ] + ) + + next_doppler_estimators, dopplers_and_filtered_prompts = Tracking.estimate_dopplers_and_filter_prompt( + doppler_estimators, + (system_sats_state,), + next_system_sats_sample_params, + sampling_frequency, + 0.75ms + ) + + @test dopplers_and_filtered_prompts[1][1].carrier_doppler == 100.0Hz + @test dopplers_and_filtered_prompts[1][1].code_doppler == 1.0Hz + @test dopplers_and_filtered_prompts[1][1].filtered_prompt == 2000 / 4500 +end \ No newline at end of file diff --git a/test/correlator.jl b/test/correlator.jl index 09fdb3e..2aa3cf0 100644 --- a/test/correlator.jl +++ b/test/correlator.jl @@ -2,39 +2,63 @@ @testset "Correlator constructor" begin gpsl1 = GPSL1() - correlator = @inferred EarlyPromptLateCorrelator() - correlator_sample_shifts = SVector(-1, 0, 1) - early_late_index_shift = 1 - @test @inferred(get_early(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0 - @test @inferred(get_prompt(correlator, correlator_sample_shifts)) == 0.0 - @test @inferred(get_late(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0 - - correlator = @inferred EarlyPromptLateCorrelator(NumAnts(1), 3) + sampling_frequency = 5e6Hz + correlator = @inferred EarlyPromptLateCorrelator(gpsl1, sampling_frequency) + + @test @inferred(get_early(correlator)) == 0.0 + @test @inferred(get_prompt(correlator)) == 0.0 + @test @inferred(get_late(correlator)) == 0.0 + @test correlator.shifts == -2:2:2 + @test correlator.early_shift == 3 + @test correlator.prompt_shift == 2 + @test correlator.late_shift == 1 + + correlator = @inferred EarlyPromptLateCorrelator(gpsl1, sampling_frequency, num_ants = NumAnts(1), num_accumulators = 3) @test isa(get_accumulators(correlator), Vector) + @test correlator.shifts == -2:2:2 + @test correlator.early_shift == 3 + @test correlator.prompt_shift == 2 + @test correlator.late_shift == 1 - correlator = @inferred EarlyPromptLateCorrelator(NumAnts(1)) + correlator = @inferred EarlyPromptLateCorrelator(gpsl1, sampling_frequency, num_ants = NumAnts(1)) - @test @inferred(get_early(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0 - @test @inferred(get_prompt(correlator, correlator_sample_shifts)) == 0.0 - @test @inferred(get_late(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0 + @test @inferred(get_early(correlator)) == 0.0 + @test @inferred(get_prompt(correlator)) == 0.0 + @test @inferred(get_late(correlator)) == 0.0 + @test correlator.shifts == -2:2:2 + @test correlator.early_shift == 3 + @test correlator.prompt_shift == 2 + @test correlator.late_shift == 1 - correlator = @inferred EarlyPromptLateCorrelator(NumAnts(1), 3) + correlator = @inferred EarlyPromptLateCorrelator(gpsl1, sampling_frequency, num_ants = NumAnts(1), num_accumulators = 3) - @test @inferred(get_early(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0 - @test @inferred(get_prompt(correlator, correlator_sample_shifts)) == 0.0 - @test @inferred(get_late(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0 + @test @inferred(get_early(correlator)) == 0.0 + @test @inferred(get_prompt(correlator)) == 0.0 + @test @inferred(get_late(correlator)) == 0.0 + @test correlator.shifts == -2:2:2 + @test correlator.early_shift == 3 + @test correlator.prompt_shift == 2 + @test correlator.late_shift == 1 - correlator = @inferred EarlyPromptLateCorrelator(NumAnts(2)) + correlator = @inferred EarlyPromptLateCorrelator(gpsl1, sampling_frequency, num_ants = NumAnts(2)) - @test @inferred(get_early(correlator, correlator_sample_shifts, early_late_index_shift)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) - @test @inferred(get_prompt(correlator, correlator_sample_shifts)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) - @test @inferred(get_late(correlator, correlator_sample_shifts, early_late_index_shift)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + @test @inferred(get_early(correlator)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + @test @inferred(get_prompt(correlator)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + @test @inferred(get_late(correlator)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + @test correlator.shifts == -2:2:2 + @test correlator.early_shift == 3 + @test correlator.prompt_shift == 2 + @test correlator.late_shift == 1 - correlator = @inferred EarlyPromptLateCorrelator(NumAnts(2), 3) + correlator = @inferred EarlyPromptLateCorrelator(gpsl1, sampling_frequency, num_ants = NumAnts(2), num_accumulators = 3) - @test @inferred(get_early(correlator, correlator_sample_shifts, early_late_index_shift)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) - @test @inferred(get_prompt(correlator, correlator_sample_shifts)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) - @test @inferred(get_late(correlator, correlator_sample_shifts, early_late_index_shift)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + @test @inferred(get_early(correlator)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + @test @inferred(get_prompt(correlator)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + @test @inferred(get_late(correlator)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + @test correlator.shifts == -2:2:2 + @test correlator.early_shift == 3 + @test correlator.prompt_shift == 2 + @test correlator.late_shift == 1 end @testset "Zeroing correlator" begin @@ -42,8 +66,10 @@ SVector( SVector(1.0 + 0.0im, 1.0 + 0.0im), SVector(1.0 + 0.0im, 1.0 + 0.0im), - SVector(1.0 + 0.0im, 1.0 + 0.0im) - ) + SVector(1.0 + 0.0im, 1.0 + 0.0im), + ), + -1:1, + 3, 2, 1 ) @test @inferred(zero(correlator)) == EarlyPromptLateCorrelator( @@ -51,7 +77,9 @@ SVector(0.0 + 0.0im, 0.0 + 0.0im), SVector(0.0 + 0.0im, 0.0 + 0.0im), SVector(0.0 + 0.0im, 0.0 + 0.0im) - ) + ), + -1:1, + 3, 2, 1 ) correlator = @inferred EarlyPromptLateCorrelator( @@ -59,7 +87,9 @@ 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im - ) + ), + -1:1, + 3, 2, 1 ) @test @inferred(zero(correlator)) == EarlyPromptLateCorrelator( @@ -67,7 +97,9 @@ 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im - ) + ), + -1:1, + 3, 2, 1 ) end @testset "Filter correlator" begin @@ -76,7 +108,9 @@ SVector(1.0 + 0.0im, 1.0 + 1.0im), SVector(1.0 + 0.0im, 2.0 + 0.0im), SVector(1.0 + 1.0im, 1.0 + 3.0im) - ) + ), + -1:1, + 3, 2, 1 ) default_post_corr_filter = Tracking.DefaultPostCorrFilter() filtered_correlator = @inferred Tracking.apply(default_post_corr_filter, correlator) @@ -85,7 +119,9 @@ 1.0 + 1.0im, 2.0 + 0.0im, 1.0 + 3.0im - ) + ), + -1:1, + 3, 2, 1 ) correlator = @inferred EarlyPromptLateCorrelator( @@ -93,7 +129,9 @@ 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im - ) + ), + -1:1, + 3, 2, 1 ) filtered_correlator = @inferred Tracking.apply(x -> 2 * x, correlator) @test filtered_correlator == EarlyPromptLateCorrelator( @@ -101,42 +139,46 @@ 2.0 + 0.0im, 2.0 + 0.0im, 2.0 + 0.0im - ) + ), + -1:1, + 3, 2, 1 ) end - @testset "Correlator sample shifts" begin - gpsl1 = GPSL1() - correlator = EarlyPromptLateCorrelator() - correlator_sample_shifts = @inferred get_correlator_sample_shifts( - gpsl1, - correlator, - 4e6Hz, - 0.5 - ) - @test correlator_sample_shifts ≈ round(0.5 * 4e6 / 1023e3) * SVector(-1, 0, 1) - @test typeof(correlator_sample_shifts) <: SVector{3,<:Integer} - end - @testset "Early late index shift" begin - gpsl1 = GPSL1() - correlator = EarlyPromptLateCorrelator() - correlator_sample_shifts = get_correlator_sample_shifts( - gpsl1, - correlator, - 4e6Hz, - 0.5 + @testset "Early late sample spacing" begin + correlator = @inferred EarlyPromptLateCorrelator( + SVector( + SVector(1.0 + 0.0im, 1.0 + 0.0im), + SVector(1.0 + 0.0im, 1.0 + 0.0im), + SVector(1.0 + 0.0im, 1.0 + 0.0im), + ), + -1:1, + 3, 2, 1 ) - early_late_index_shift = @inferred get_early_late_index_shift( - gpsl1, - correlator_sample_shifts, - correlator, - 4e6Hz, - 0.5 + @test get_early_late_sample_spacing(correlator) == 2 + + correlator = @inferred EarlyPromptLateCorrelator( + SVector( + SVector(1.0 + 0.0im, 1.0 + 0.0im), + SVector(1.0 + 0.0im, 1.0 + 0.0im), + SVector(1.0 + 0.0im, 1.0 + 0.0im), + ), + -2:2:2, + 3, 2, 1 ) + @test get_early_late_sample_spacing(correlator) == 4 - @test early_late_index_shift == 1 - end - @testset "Early late sample spacing" begin - @test get_early_late_sample_spacing(SVector(-2,0,2), 1) == 4 + correlator = @inferred EarlyPromptLateCorrelator( + SVector( + SVector(1.0 + 0.0im, 1.0 + 0.0im), + SVector(1.0 + 0.0im, 1.0 + 0.0im), + SVector(1.0 + 0.0im, 1.0 + 0.0im), + SVector(1.0 + 0.0im, 1.0 + 0.0im), + SVector(1.0 + 0.0im, 1.0 + 0.0im), + ), + -2:2, + 5, 3, 1 + ) + @test get_early_late_sample_spacing(correlator) == 4 end @testset "Normalize correlator" begin correlator = @inferred EarlyPromptLateCorrelator( @@ -144,7 +186,9 @@ 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im - ) + ), + -1:1, + 3, 2, 1 ) normalized_correlator = @inferred Tracking.normalize(correlator, 10) @test normalized_correlator == EarlyPromptLateCorrelator( @@ -152,7 +196,9 @@ 0.1 + 0.0im, 0.1 + 0.0im, 0.1 + 0.0im - ) + ), + -1:1, + 3, 2, 1 ) correlator = EarlyPromptLateCorrelator( @@ -160,7 +206,9 @@ SVector(1.0 + 0.0im, 1.0 + 0.0im), SVector(1.0 + 0.0im, 1.0 + 0.0im), SVector(1.0 + 0.0im, 1.0 + 0.0im) - ) + ), + -1:1, + 3, 2, 1 ) normalized_correlator = @inferred Tracking.normalize(correlator, 10) @test normalized_correlator == EarlyPromptLateCorrelator( @@ -168,45 +216,54 @@ SVector(0.1 + 0.0im, 0.1 + 0.0im), SVector(0.1 + 0.0im, 0.1 + 0.0im), SVector(0.1 + 0.0im, 0.1 + 0.0im) - ) + ), + -1:1, + 3, 2, 1 ) end - @testset "Correlate" begin gpsl1 = GPSL1() + sampling_frequency = 2.5e6Hz signal = StructArray{Complex{Float32}}( - (Float32.(get_code.(gpsl1, (1:2500) * 1023e3 / 2.5e6, 1)), + (Float32.(get_code.(gpsl1, (1:2500) * 1023e3 / (sampling_frequency / Hz), 1)), zeros(Float32, 2500)) ) - correlator_sample_shifts = SVector(-1, 0, 1) code = get_code.( gpsl1, - (1 + correlator_sample_shifts[1]:2500 + correlator_sample_shifts[end]) * 1023e3 / 2.5e6, + (1 + -1:2500 + 1) * 1023e3 / 2.5e6, 1 ) - @testset "↳ $na antennas" for na ∈ [1, 4, 15] - signal_mat = na == 1 ? signal : repeat(signal, outer=(1,na)) + @testset "↳ $na antennas" for na ∈ [1, 4] + signal_mat = na == 1 ? signal : repeat(signal, outer=(1, na)) @testset "↳ $(string(type)) accumulator" for type ∈ [:SVector, :Vector] if type == :SVector - correlator = EarlyPromptLateCorrelator(NumAnts(na)) + correlator = EarlyPromptLateCorrelator( + gpsl1, + sampling_frequency, + num_ants = NumAnts(na) + ) else - correlator = EarlyPromptLateCorrelator(NumAnts(na), 3) + correlator = EarlyPromptLateCorrelator( + gpsl1, + sampling_frequency, + num_ants = NumAnts(na), + num_accumulators = 3 + ) end - correlator_result = Tracking.correlate( + correlator_result = @inferred Tracking.correlate( correlator, signal_mat, code, - correlator_sample_shifts, 1, 2500, ) - early = get_early(correlator_result, correlator_sample_shifts, 1) - prompt = get_prompt(correlator_result, correlator_sample_shifts) - late = get_late(correlator_result, correlator_sample_shifts, 1) + early = get_early(correlator_result) + prompt = get_prompt(correlator_result) + late = get_late(correlator_result) @test early == late @test all(early .== 1476) @test all(prompt .== 2500) diff --git a/test/cuda/bit_buffer.jl b/test/cuda/bit_buffer.jl deleted file mode 100644 index 37bf940..0000000 --- a/test/cuda/bit_buffer.jl +++ /dev/null @@ -1,93 +0,0 @@ -@testset "CUDA: Bit buffer" begin - bit_buffer = Tracking.BitBuffer() - @test @inferred(get_bits(bit_buffer)) == 0 - @test @inferred(Tracking.length(bit_buffer)) == 0 - @test @allocated(Tracking.BitBuffer()) == 0 - - prompt_accumulator = 0.0 + 0.0im - secondary_code_or_bit_found = false - prev_code_phase = 0.0 - code_phase = 1023 - prompt_correlator = 1.0 + 0.0im - integration_time = 5ms - gpsl1 = GPSL1(use_gpu = Val(true)) - next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( - gpsl1, - bit_buffer, - prompt_accumulator, - secondary_code_or_bit_found, - prev_code_phase, - code_phase, - integration_time, - prompt_correlator - ) - - @test next_prompt_accumulator == 0.0 + 0.0im - @test @inferred(get_bits(next_bit_buffer)) == 0 - @test @inferred(Tracking.length(next_bit_buffer)) == 0 - - prompt_accumulator = 20.0 + 0.0im - secondary_code_or_bit_found = true - prev_code_phase = 0.0 - code_phase = 1023 - prompt_correlator = 1.0 + 0.0im - integration_time = 5ms - - next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( - gpsl1, - bit_buffer, - prompt_accumulator, - secondary_code_or_bit_found, - prev_code_phase, - code_phase, - integration_time, - prompt_correlator - ) - @test next_prompt_accumulator == 21.0 + 0.0im - @test @inferred(get_bits(next_bit_buffer)) == 0 - @test @inferred(Tracking.length(next_bit_buffer)) == 0 - - prompt_accumulator = 20.0 + 0.0im - secondary_code_or_bit_found = true - prev_code_phase = 9000 - code_phase = 1023 - prompt_correlator = 1.0 + 0.0im - integration_time = 20ms - - next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( - gpsl1, - bit_buffer, - prompt_accumulator, - secondary_code_or_bit_found, - prev_code_phase, - code_phase, - integration_time, - prompt_correlator - ) - @test next_prompt_accumulator == 0.0 + 0.0im - @test @inferred(get_bits(next_bit_buffer)) == 1 - @test @inferred(Tracking.length(next_bit_buffer)) == 1 - - prompt_accumulator = 20.0 + 0.0im - secondary_code_or_bit_found = true - prev_code_phase = 0 - code_phase = 0 - prompt_correlator = 1.0 + 0.0im - integration_time = 20ms - - next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( - gpsl1, - bit_buffer, - prompt_accumulator, - secondary_code_or_bit_found, - prev_code_phase, - code_phase, - integration_time, - prompt_correlator - ) - @test next_prompt_accumulator == 0.0 + 0.0im - @test @inferred(get_bits(next_bit_buffer)) == 1 - @test @inferred(Tracking.length(next_bit_buffer)) == 1 - - -end diff --git a/test/cuda/cn0_estimation.jl b/test/cuda/cn0_estimation.jl deleted file mode 100644 index 6182cc7..0000000 --- a/test/cuda/cn0_estimation.jl +++ /dev/null @@ -1,50 +0,0 @@ -@testset "CUDA: CN0 estimation" begin - CUDA.allowscalar() do - 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 = MomentsCN0Estimator(20) - correlator_sample_shifts = SVector(-2, 0, 2) - start_sample = 1 - num_samples = 4000 - gpsl1 = GPSL1(use_gpu = Val(true)) - - for i = 1:20 - signal = get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ start_code_phase, - prn - ) .* 10^(45 / 20) .+ - randn(ComplexF64, length(range)) .* sqrt(4e6) - code = get_code.( - gpsl1, - code_frequency .* (-2:4001) ./ sampling_frequency .+ start_code_phase, - prn - ) - correlator = EarlyPromptLateCorrelator() - signal_struct = StructArray(signal) - correlator = Tracking.correlate( - correlator, - signal_struct, - code, - correlator_sample_shifts, - start_sample, - num_samples, - ) - cn0_estimator = Tracking.update( - cn0_estimator, - get_prompt(correlator, correlator_sample_shifts) - ) - end - @test @inferred(Tracking.get_current_index(cn0_estimator)) == 20 - @test @inferred(Tracking.length(cn0_estimator)) == 20 - cn0_estimate = @inferred Tracking.estimate_cn0(cn0_estimator, 1ms) - - @test cn0_estimate ≈ 45dBHz atol = 1.05dBHz -end -end diff --git a/test/cuda/discriminators.jl b/test/cuda/discriminators.jl deleted file mode 100644 index 97c2678..0000000 --- a/test/cuda/discriminators.jl +++ /dev/null @@ -1,61 +0,0 @@ -@testset "CUDA: PLL discriminator" begin - - correlator_minus60off = EarlyPromptLateCorrelator( - SVector( - -0.5 + sqrt(3) / 2im, - -1 + sqrt(3) * 1im, - -0.5 + sqrt(3) / 2im - ) - ) - correlator_0off = EarlyPromptLateCorrelator( - SVector( - 0.5 + 0.0im, - 1 + 0.0im, - 0.5 + 0.0im - ) - ) - correlator_plus60off = EarlyPromptLateCorrelator( - SVector( - 0.5 + sqrt(3) / 2im, - 1 + sqrt(3) * 1im, - 0.5 + sqrt(3) / 2im - ) - ) - gpsl1 = GPSL1(use_gpu = Val(true)) - correlator_sample_shifts = SVector(-1, 0, 1) - @test @inferred(Tracking.pll_disc(gpsl1, correlator_minus60off, correlator_sample_shifts)) == -π / 3 #-60° - @test @inferred(Tracking.pll_disc(gpsl1, correlator_0off, correlator_sample_shifts)) == 0 - @test @inferred(Tracking.pll_disc(gpsl1, correlator_plus60off, correlator_sample_shifts)) == π / 3 #+60° -end - - -@testset "CUDA: DLL discriminator" begin - - sample_shifts = SVector(-2, 0, 2) - index_shift = 1 - delta = 0.25 - @test @inferred(get_early_late_sample_spacing(sample_shifts, index_shift)) == 4 - - gpsl1 = GPSL1(use_gpu = Val(true)) - very_early_correlator = EarlyPromptLateCorrelator(SVector(1.0 + 0.0im, 0.5 + 0.0im, 0.0 + 0.0im)) - early_correlator = EarlyPromptLateCorrelator(SVector(0.75 + 0.0im, 0.75 + 0.0im, 0.25 + 0.0im)) - prompt_correlator = EarlyPromptLateCorrelator(SVector(0.5 + 0.0im, 1 + 0.0im, 0.5 + 0.0im)) - late_correlator = EarlyPromptLateCorrelator(SVector(0.25 + 0.0im, 0.75 + 0.0im, 0.75 + 0.0im)) - very_late_correlator = EarlyPromptLateCorrelator(SVector(0.0 + 0.0im, 0.5 + 0.0im, 1.0 + 0.0im)) - - @test @inferred( - Tracking.dll_disc(gpsl1, very_early_correlator, sample_shifts, index_shift, delta) - ) == -0.5 - @test @inferred( - Tracking.dll_disc(gpsl1, early_correlator, sample_shifts, index_shift, delta) - ) == -0.25 - @test @inferred( - Tracking.dll_disc(gpsl1, prompt_correlator, sample_shifts, index_shift, delta) - ) == 0 - @test @inferred( - Tracking.dll_disc(gpsl1, late_correlator, sample_shifts, index_shift, delta) - ) == 0.25 - @test @inferred( - Tracking.dll_disc(gpsl1, very_late_correlator, sample_shifts, index_shift, delta) - ) == 0.5 -end diff --git a/test/cuda/galileo_e1b.jl b/test/cuda/galileo_e1b.jl deleted file mode 100644 index ffe9d73..0000000 --- a/test/cuda/galileo_e1b.jl +++ /dev/null @@ -1,15 +0,0 @@ -@testset "CUDA: Galileo E1B" begin - galileo_e1b = GalileoE1B(use_gpu = Val(true)) - @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf, 29)) == true - - @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf, 6)) == false - - @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xc, 40)) == false - - @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf, 8)) == true - - @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf0, 8)) == true - - @test @inferred(Tracking.get_default_correlator(galileo_e1b, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_default_correlator(galileo_e1b, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3)) -end diff --git a/test/cuda/gps_l1.jl b/test/cuda/gps_l1.jl deleted file mode 100644 index 45e9684..0000000 --- a/test/cuda/gps_l1.jl +++ /dev/null @@ -1,11 +0,0 @@ -@testset "CUDA: GPS L1" begin - gpsl1 = GPSL1(use_gpu = Val(true)) - @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl1, 0xfffff00000, 50)) == true - - @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl1, 0xfffff, 10)) == false - - @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl1, 0xfffff, 40)) == true - - @test @inferred(Tracking.get_default_correlator(gpsl1, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_default_correlator(gpsl1, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3)) -end diff --git a/test/cuda/gps_l5.jl b/test/cuda/gps_l5.jl deleted file mode 100644 index 4810a54..0000000 --- a/test/cuda/gps_l5.jl +++ /dev/null @@ -1,11 +0,0 @@ -@testset "CUDA: GPS L5" begin - gpsl5 = GPSL5(use_gpu = Val(true)) - @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl5, 0x35, 50)) == true - - @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl5, 0x35, 5)) == false - - @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl5, 0x3ca, 10)) == true # 0x3ca == 1111001010 - - @test @inferred(Tracking.get_default_correlator(gpsl5, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_default_correlator(gpsl5, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3)) -end diff --git a/test/cuda/secondary_code_or_bit_detector.jl b/test/cuda/secondary_code_or_bit_detector.jl deleted file mode 100644 index c128f51..0000000 --- a/test/cuda/secondary_code_or_bit_detector.jl +++ /dev/null @@ -1,40 +0,0 @@ -@testset "CUDA: Secondary code or bit detector" begin - detector = SecondaryCodeOrBitDetector() - gpsl1 = GPSL1(use_gpu = Val(true)) - - @test Tracking.get_buffer(detector) == 0 - @test Tracking.length(detector) == 0 - @test Tracking.found(detector) == false - - next_detector = Tracking.find(gpsl1, detector, -1.0 + 0.0im) - @test Tracking.get_buffer(next_detector) == 0 - @test Tracking.length(next_detector) == 1 - @test Tracking.found(next_detector) == false - - next2_detector = Tracking.find(gpsl1, next_detector, 1.0 + 0.0im) - @test Tracking.get_buffer(next2_detector) == 1 - @test Tracking.length(next2_detector) == 2 - @test Tracking.found(next2_detector) == false - - next3_detector = Tracking.find(gpsl1, next2_detector, 1.0 + 0.0im) - @test Tracking.get_buffer(next3_detector) == 3 - @test Tracking.length(next3_detector) == 3 - @test Tracking.found(next3_detector) == false - - detector = SecondaryCodeOrBitDetector(1, 1, true) - - @test Tracking.get_buffer(detector) == 1 - @test Tracking.length(detector) == 1 - @test Tracking.found(detector) == true - - next_detector = Tracking.find(gpsl1, detector, 1.0 + 0.0im) - @test Tracking.get_buffer(next_detector) == 1 - @test Tracking.length(next_detector) == 1 - @test Tracking.found(next_detector) == true - - detector = SecondaryCodeOrBitDetector(0x7ffff80000, 39, false) - next_detector = Tracking.find(gpsl1, detector, -1.0 + 0.0im) - @test Tracking.get_buffer(next_detector) == 0x000000fffff00000 - @test Tracking.length(next_detector) == 40 - @test Tracking.found(next_detector) == true -end diff --git a/test/cuda/tracking_loop.jl b/test/cuda/tracking_loop.jl deleted file mode 100644 index e2ec163..0000000 --- a/test/cuda/tracking_loop.jl +++ /dev/null @@ -1,267 +0,0 @@ -@testset "CUDA: Integration time" begin - gpsl1 = GPSL1(use_gpu = Val(true)) - galileo_e1b = GalileoE1B(use_gpu = Val(true)) - max_integration_time = 2ms - bit_found = false - time = @inferred Tracking.get_integration_time(gpsl1, max_integration_time, bit_found) - @test time == 1ms - - max_integration_time = 2ms - bit_found = true - time = @inferred Tracking.get_integration_time(gpsl1, max_integration_time, bit_found) - @test time == 2ms - - max_integration_time = 2ms - bit_found = false - time = @inferred Tracking.get_integration_time( - galileo_e1b, - max_integration_time, - bit_found - ) - @test time == 2ms - - max_integration_time = 10ms - bit_found = false - time = @inferred Tracking.get_integration_time( - galileo_e1b, - max_integration_time, - bit_found - ) - @test time == 4ms -end - -@testset "CUDA: Number of chips to integrate GPU" begin - gpsl1 = GPSL1(use_gpu = Val(true)) - max_integration_time = 1ms - code_phase = 200 - bit_found = true - chips = @inferred Tracking.get_num_chips_to_integrate( - gpsl1, - max_integration_time, - code_phase, - bit_found - ) - @test chips == 1023 - 200 - - code_phase = 1200 - chips = @inferred Tracking.get_num_chips_to_integrate( - gpsl1, - max_integration_time, - code_phase, - bit_found - ) - @test chips == 2046 - 1200 -end - -@testset "CUDA: Number of samples to integrate" begin - gpsl1 = GPSL1(use_gpu = Val(true)) - max_integration_time = 1ms - sampling_frequency = 4e6Hz - current_code_doppler = 0Hz - current_code_phase = 0 - bit_found = true - samples = @inferred Tracking.get_num_samples_left_to_integrate( - gpsl1, - max_integration_time, - sampling_frequency, - current_code_doppler, - current_code_phase, - bit_found - ) - @test samples == 4000 -end - -@testset "CUDA: Code phase delta" begin - gpsl1 = GPSL1(use_gpu = Val(true)) - code_doppler = 1Hz - frequency = @inferred Tracking.get_current_code_frequency(gpsl1, code_doppler) - @test frequency == 1023e3Hz + 1Hz -end - -@testset "CUDA: Doppler aiding" begin - gpsl1 = GPSL1(use_gpu = Val(true)) - init_carrier_doppler = 10Hz - init_code_doppler = 1Hz - carrier_freq_update = 2Hz - code_freq_update = -0.5Hz - velocity_aiding = 3Hz - - carrier_freq, code_freq = @inferred Tracking.aid_dopplers( - gpsl1, - init_carrier_doppler, - init_code_doppler, - carrier_freq_update, - code_freq_update, - velocity_aiding - ) - - @test carrier_freq == 10Hz + 2Hz + 3Hz - @test code_freq == 1Hz + (2Hz + 3Hz) / 1540 - 0.5Hz -end - -@testset "CUDA: Tracking with signal of type StructArray of CuArrays ComplexF32" begin - CUDA.allowscalar() do - gpsl1 = GPSL1(use_gpu = Val(true)) - carrier_doppler = 200Hz - start_code_phase = 100 - code_frequency = carrier_doppler / 1540 + 1023e3Hz - sampling_frequency = 4e6Hz - prn = 1 - range = 0:3999 - start_carrier_phase = π / 2 - state = TrackingState(prn, gpsl1, carrier_doppler - 20Hz, start_code_phase, num_samples = 4000) - - carrier_phases = cu(2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase) - code_phases = get_code_frequency(gpsl1) / sampling_frequency .* range .+ start_code_phase - upsampled_code = gpsl1.codes[1 .+ mod.(floor.(Int, code_phases), 1023), prn] - signal = StructArray{ComplexF32}((cos.(carrier_phases), sin.(carrier_phases))) - @. signal *= upsampled_code - - track_result = @inferred track(signal, state, sampling_frequency) - - iterations = 2000 - code_phases = zeros(iterations) - carrier_phases = zeros(iterations) - tracked_code_phases = zeros(iterations) - tracked_carrier_phases = zeros(iterations) - tracked_code_dopplers = zeros(iterations) - tracked_carrier_dopplers = zeros(iterations) - tracked_prompts = zeros(ComplexF64, iterations) - for i = 1:iterations - carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency + - start_carrier_phase + π) - π - code_phase = mod( - code_frequency * 4000 * i / sampling_frequency + start_code_phase, - 1023 - ) - signal_temp = CuArray{ComplexF32}(cis.( - 2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase - ) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ code_phase, - prn - ) - ) - signal = StructArray{ComplexF32}((real.(signal_temp), imag.(signal_temp))) - track_result = @inferred track( - signal, - get_state(track_result), - sampling_frequency - ) - comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) / - sampling_frequency + start_carrier_phase + π) - π - comp_code_phase = mod( - code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, - 1023 - ) - tracked_code_phases[i] = get_code_phase(track_result) - tracked_carrier_phases[i] = get_carrier_phase(track_result) - tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz - tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz - tracked_prompts[i] = get_prompt(track_result) - code_phases[i] = comp_code_phase - carrier_phases[i] = comp_carrier_phase - end - - @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4 - @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-2 - end -# using PyPlot -# pygui(true) -# figure("carrier_phases") -# plot(tracked_carrier_phases) -# plot(carrier_phases) -# grid(true) -# figure("Code Phases") -# plot(300 * (tracked_code_phases .- code_phases)) -# figure("Carrier Doppler") -# plot(tracked_carrier_dopplers) -# figure("Code Doppler") -# plot(tracked_code_dopplers) -# figure("Prompt") -# plot(real.(tracked_prompts)) -# plot(imag.(tracked_prompts)) -end - -@testset "CUDA: Track multiple signals with signal of type StructArray of CuArrays ComplexF32" begin - CUDA.allowscalar() do - gpsl1 = GPSL1(use_gpu = Val(true)) - carrier_doppler = 200Hz - start_code_phase = 100 - code_frequency = carrier_doppler / 1540 + 1023e3Hz - sampling_frequency = 4e6Hz - prn = 1 - range = 0:3999 - start_carrier_phase = π / 2 - state = TrackingState(prn, gpsl1, carrier_doppler - 20Hz, start_code_phase, num_ants = NumAnts(3), num_samples = 4000) - - carrier_phases = cu(2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase) - code_phases = get_code_frequency(gpsl1) / sampling_frequency .* range .+ start_code_phase - upsampled_code = gpsl1.codes[1 .+ mod.(floor.(Int, code_phases), 1023), prn] - carrier_phases_mat = carrier_phases * CUDA.ones(3)' - signal_mat = StructArray{ComplexF32}((cos.(carrier_phases_mat), sin.(carrier_phases_mat))) - @. signal_mat *= upsampled_code - - track_result = @inferred track(signal_mat, state, sampling_frequency) - - iterations = 2000 - code_phases = zeros(iterations) - carrier_phases = zeros(iterations) - tracked_code_phases = zeros(iterations) - tracked_carrier_phases = zeros(iterations) - tracked_code_dopplers = zeros(iterations) - tracked_carrier_dopplers = zeros(iterations) - for i = 1:iterations - carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency + - start_carrier_phase + π) - π - code_phase = mod( - code_frequency * 4000 * i / sampling_frequency + start_code_phase, - 1023 - ) - signal_temp = CuArray{ComplexF32}(cis.( - 2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase - ) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ code_phase, - prn - ) - ) - signal_temp_mat = repeat(signal_temp, outer = (1,3)) - signal_mat = StructArray(signal_temp_mat) - track_result = @inferred track( - signal_mat, - get_state(track_result), - sampling_frequency - ) - comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) / - sampling_frequency + start_carrier_phase + π) - π - comp_code_phase = mod( - code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, - 1023 - ) - tracked_code_phases[i] = get_code_phase(track_result) - tracked_carrier_phases[i] = get_carrier_phase(track_result) - tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz - tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz - code_phases[i] = comp_code_phase - carrier_phases[i] = comp_carrier_phase - end - - @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4 - @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-2 - end -# using PyPlot -# pygui(true) -# figure("carrier_phases") -# plot(tracked_carrier_phases) -# plot(carrier_phases) -# grid(true) -# figure("Code Phases") -# plot(300 * (tracked_code_phases .- code_phases)) -# figure("Carrier Doppler") -# plot(tracked_carrier_dopplers) -# figure("Code Doppler") -# plot(tracked_code_dopplers) -end \ No newline at end of file diff --git a/test/cuda/tracking_results.jl b/test/cuda/tracking_results.jl deleted file mode 100644 index 800cccf..0000000 --- a/test/cuda/tracking_results.jl +++ /dev/null @@ -1,52 +0,0 @@ -@testset "CUDA: Tracking results" begin - gpsl1 = GPSL1(use_gpu = Val(true)) - results = Tracking.TrackingResults( - TrackingState(1, gpsl1, 100Hz, 100, num_samples = 2500), - EarlyPromptLateCorrelator(NumAnts(2)), - complex(1.2,1.3), - SVector(-1, 0, 1), - 1, - 1.0Hz, - 1.0, - true, - Tracking.BitBuffer(), - 45dBHz - ) - - @test @inferred(get_carrier_doppler(results)) == 100Hz - @test @inferred(get_carrier_phase(results)) == 0.0 - @test @inferred(get_code_doppler(results)) == 100Hz / 1540 - @test @inferred(get_code_phase(results)) == 100 - @test @inferred(get_correlator_carrier_phase(results)) ≈ 1.0 * 2π - @test @inferred(get_correlator_carrier_frequency(results)) == 1.0Hz - correlator = @inferred get_correlator(results) - @test correlator == EarlyPromptLateCorrelator(NumAnts(2)) - @test @inferred(get_correlator_sample_shifts(results)) == SVector(-1, 0, 1) - @test @inferred(get_early_late_index_shift(results)) == 1 - @test @inferred(get_filtered_prompt(results)) == complex(1.2,1.3) - @test @inferred(get_early(results)) == [0.0, 0.0] - @test @inferred(get_prompt(results)) == [0.0, 0.0] - @test @inferred(get_late(results)) == [0.0, 0.0] - @test @inferred(get_bits(results)) == 0 - @test @inferred(get_num_bits(results)) == 0 - @test @inferred(get_secondary_code_or_bit_found(results)) == false - @test @inferred(get_cn0(results)) == 45dBHz - - results = Tracking.TrackingResults( - TrackingState(1, gpsl1, 100Hz, 100, num_samples = 2500), - EarlyPromptLateCorrelator(NumAnts(2)), - complex(1.2,1.3), - SVector(-1, 0, 1), - 1, - 1.0Hz, - 1.0, - false, - Tracking.BitBuffer(), - 45dBHz - ) - correlator = @inferred get_correlator(results) - @test @inferred(get_filtered_prompt(results)) == complex(1.2,1.3) - @test @inferred(get_early(results)) == [0, 0] - @test @inferred(get_prompt(results)) == [0, 0] - @test @inferred(get_late(results)) == [0, 0] -end diff --git a/test/cuda/tracking_state.jl b/test/cuda/tracking_state.jl deleted file mode 100644 index 28519a3..0000000 --- a/test/cuda/tracking_state.jl +++ /dev/null @@ -1,60 +0,0 @@ -@testset "CUDA: Tracking state" begin - - carrier_doppler = 100Hz - code_phase = 100 - gpsl1 = GPSL1(use_gpu = Val(true)) - num_samples = 2500 - state = TrackingState(1, gpsl1, carrier_doppler, code_phase, num_samples = num_samples) - - @test_throws UndefKeywordError TrackingState(1, gpsl1, carrier_doppler, code_phase) - - @test @inferred(Tracking.get_prn(state)) == 1 - @test @inferred(Tracking.get_code_phase(state)) == 100 - @test @inferred(Tracking.get_carrier_phase(state)) == 0.0 - @test @inferred(Tracking.get_init_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_init_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_correlator(state)) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_sc_bit_detector(state)) == SecondaryCodeOrBitDetector() - @test @inferred(Tracking.get_carrier_loop_filter(state)) == ThirdOrderBilinearLF() - @test @inferred(Tracking.get_code_loop_filter(state)) == SecondOrderBilinearLF() - @test @inferred(Tracking.get_prompt_accumulator(state)) == 0.0 - @test @inferred(Tracking.get_integrated_samples(state)) == 0 - - state = TrackingState( - 1, - gpsl1, - carrier_doppler, - code_phase; - num_samples = 2500, - code_doppler = carrier_doppler * get_code_center_frequency_ratio(gpsl1), - carrier_phase = 0.0, - carrier_loop_filter = ThirdOrderBilinearLF(), - code_loop_filter = SecondOrderBilinearLF(), - sc_bit_detector = SecondaryCodeOrBitDetector(), - correlator = Tracking.get_default_correlator(gpsl1, NumAnts(1)), - integrated_samples = 0, - prompt_accumulator = zero(ComplexF64) - ) - - @test @inferred(Tracking.get_prn(state)) == 1 - @test @inferred(Tracking.get_code_phase(state)) == 100 - @test @inferred(Tracking.get_carrier_phase(state)) == 0.0 - @test @inferred(Tracking.get_init_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_init_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_correlator(state)) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_sc_bit_detector(state)) == SecondaryCodeOrBitDetector() - @test @inferred(Tracking.get_carrier_loop_filter(state)) == ThirdOrderBilinearLF() - @test @inferred(Tracking.get_code_loop_filter(state)) == SecondOrderBilinearLF() - @test @inferred(Tracking.get_prompt_accumulator(state)) == 0.0 - @test @inferred(Tracking.get_integrated_samples(state)) == 0 - - state = TrackingState(1, gpsl1, carrier_doppler, code_phase, num_samples = num_samples, num_ants = NumAnts(2)) - @test @inferred(Tracking.get_correlator(state)) == EarlyPromptLateCorrelator(NumAnts(2)) - - changed_state = @inferred(TrackingState(state, post_corr_filter = Tracking.DefaultPostCorrFilter())) - @test get_post_corr_filter(changed_state) == Tracking.DefaultPostCorrFilter() -end diff --git a/test/discriminators.jl b/test/discriminators.jl index 9ded338..c4dc88d 100644 --- a/test/discriminators.jl +++ b/test/discriminators.jl @@ -5,57 +5,62 @@ -0.5 + sqrt(3) / 2im, -1 + sqrt(3) * 1im, -0.5 + sqrt(3) / 2im - ) + ), + -1:1, + 3, 2, 1 ) correlator_0off = EarlyPromptLateCorrelator( SVector( 0.5 + 0.0im, 1 + 0.0im, 0.5 + 0.0im - ) + ), + -1:1, + 3, 2, 1 ) correlator_plus60off = EarlyPromptLateCorrelator( SVector( 0.5 + sqrt(3) / 2im, 1 + sqrt(3) * 1im, 0.5 + sqrt(3) / 2im - ) + ), + -1:1, + 3, 2, 1 ) gpsl1 = GPSL1() - correlator_sample_shifts = SVector(-1, 0, 1) - @test @inferred(Tracking.pll_disc(gpsl1, correlator_minus60off, correlator_sample_shifts)) == -π / 3 #-60° - @test @inferred(Tracking.pll_disc(gpsl1, correlator_0off, correlator_sample_shifts)) == 0 - @test @inferred(Tracking.pll_disc(gpsl1, correlator_plus60off, correlator_sample_shifts)) == π / 3 #+60° + @test @inferred(Tracking.pll_disc(gpsl1, correlator_minus60off)) == -π / 3 #-60° + @test @inferred(Tracking.pll_disc(gpsl1, correlator_0off)) == 0 + @test @inferred(Tracking.pll_disc(gpsl1, correlator_plus60off)) == π / 3 #+60° end @testset "DLL discriminator" begin sample_shifts = SVector(-2, 0, 2) - index_shift = 1 delta = 0.25 - @test @inferred(get_early_late_sample_spacing(sample_shifts, index_shift)) == 4 gpsl1 = GPSL1() - very_early_correlator = EarlyPromptLateCorrelator(SVector(1.0 + 0.0im, 0.5 + 0.0im, 0.0 + 0.0im)) - early_correlator = EarlyPromptLateCorrelator(SVector(0.75 + 0.0im, 0.75 + 0.0im, 0.25 + 0.0im)) - prompt_correlator = EarlyPromptLateCorrelator(SVector(0.5 + 0.0im, 1 + 0.0im, 0.5 + 0.0im)) - late_correlator = EarlyPromptLateCorrelator(SVector(0.25 + 0.0im, 0.75 + 0.0im, 0.75 + 0.0im)) - very_late_correlator = EarlyPromptLateCorrelator(SVector(0.0 + 0.0im, 0.5 + 0.0im, 1.0 + 0.0im)) + very_early_correlator = EarlyPromptLateCorrelator(SVector(1.0 + 0.0im, 0.5 + 0.0im, 0.0 + 0.0im), sample_shifts, 3, 2, 1) + early_correlator = EarlyPromptLateCorrelator(SVector(0.75 + 0.0im, 0.75 + 0.0im, 0.25 + 0.0im), sample_shifts, 3, 2, 1) + prompt_correlator = EarlyPromptLateCorrelator(SVector(0.5 + 0.0im, 1 + 0.0im, 0.5 + 0.0im), sample_shifts, 3, 2, 1) + late_correlator = EarlyPromptLateCorrelator(SVector(0.25 + 0.0im, 0.75 + 0.0im, 0.75 + 0.0im), sample_shifts, 3, 2, 1) + very_late_correlator = EarlyPromptLateCorrelator(SVector(0.0 + 0.0im, 0.5 + 0.0im, 1.0 + 0.0im), sample_shifts, 3, 2, 1) + + @test @inferred(Tracking.get_early_late_sample_spacing(prompt_correlator)) == 4 @test @inferred( - Tracking.dll_disc(gpsl1, very_early_correlator, sample_shifts, index_shift, delta) + Tracking.dll_disc(gpsl1, very_early_correlator, delta) ) == -0.5 @test @inferred( - Tracking.dll_disc(gpsl1, early_correlator, sample_shifts, index_shift, delta) + Tracking.dll_disc(gpsl1, early_correlator, delta) ) == -0.25 @test @inferred( - Tracking.dll_disc(gpsl1, prompt_correlator, sample_shifts, index_shift, delta) + Tracking.dll_disc(gpsl1, prompt_correlator, delta) ) == 0 @test @inferred( - Tracking.dll_disc(gpsl1, late_correlator, sample_shifts, index_shift, delta) + Tracking.dll_disc(gpsl1, late_correlator, delta) ) == 0.25 @test @inferred( - Tracking.dll_disc(gpsl1, very_late_correlator, sample_shifts, index_shift, delta) + Tracking.dll_disc(gpsl1, very_late_correlator, delta) ) == 0.5 end diff --git a/test/downconvert_and_correlate.jl b/test/downconvert_and_correlate.jl new file mode 100644 index 0000000..b0a5af5 --- /dev/null +++ b/test/downconvert_and_correlate.jl @@ -0,0 +1,151 @@ +@testset "Downconvert and Correlator" begin + gpsl1 = GPSL1() + + system_sats_states = (SystemSatsState( + gpsl1, + [ + SatState(gpsl1, 1, 5e6Hz, 10.0, 500.0Hz), + ], + ),) + + downconvert_and_correlator = @inferred CPUDownconvertAndCorrelator(system_sats_states, 5000) + @test size(downconvert_and_correlator.buffers[1][1].code_replica_buffer) == (5004,) + @test eltype(downconvert_and_correlator.buffers[1][1].code_replica_buffer) == Int16 + @test size(downconvert_and_correlator.buffers[1][1].carrier_replica_buffer) == (5000,) + @test eltype(downconvert_and_correlator.buffers[1][1].carrier_replica_buffer) == ComplexF32 + @test size(downconvert_and_correlator.buffers[1][1].downconvert_signal_buffer) == (5000,) + @test eltype(downconvert_and_correlator.buffers[1][1].downconvert_signal_buffer) == ComplexF32 + + downconvert_and_correlator_f64 = @inferred CPUDownconvertAndCorrelator(Float64, system_sats_states, 5000) + @test size(downconvert_and_correlator_f64.buffers[1][1].code_replica_buffer) == (5004,) + @test eltype(downconvert_and_correlator_f64.buffers[1][1].code_replica_buffer) == Int16 + @test size(downconvert_and_correlator_f64.buffers[1][1].carrier_replica_buffer) == (5000,) + @test eltype(downconvert_and_correlator_f64.buffers[1][1].carrier_replica_buffer) == ComplexF64 + @test size(downconvert_and_correlator_f64.buffers[1][1].downconvert_signal_buffer) == (5000,) + @test eltype(downconvert_and_correlator_f64.buffers[1][1].downconvert_signal_buffer) == ComplexF64 + + system_sats_states_numants2 = (SystemSatsState( + gpsl1, + [ + SatState(gpsl1, 1, 5e6Hz, 10.0, 500.0Hz, num_ants = NumAnts(2)), + ], + ),) + + downconvert_and_correlator_numants2 = @inferred CPUDownconvertAndCorrelator(system_sats_states_numants2, 5000) + @test size(downconvert_and_correlator_numants2.buffers[1][1].code_replica_buffer) == (5004,) + @test eltype(downconvert_and_correlator_numants2.buffers[1][1].code_replica_buffer) == Int16 + @test size(downconvert_and_correlator_numants2.buffers[1][1].carrier_replica_buffer) == (5000,) + @test eltype(downconvert_and_correlator_numants2.buffers[1][1].carrier_replica_buffer) == ComplexF32 + @test size(downconvert_and_correlator_numants2.buffers[1][1].downconvert_signal_buffer) == (5000, 2) + @test eltype(downconvert_and_correlator_numants2.buffers[1][1].downconvert_signal_buffer) == ComplexF32 + + galileo_e1b = GalileoE1B() + system_sats_states = ( + SystemSatsState( + gpsl1, + [ + SatState(gpsl1, 1, 5e6Hz, 10.0, 500.0Hz), + ], + ), + SystemSatsState( + galileo_e1b, + [ + SatState(galileo_e1b, 1, 5e6Hz, 10.0, 500.0Hz), + ], + ), + ) + + downconvert_and_correlator2 = @inferred CPUDownconvertAndCorrelator(system_sats_states, 5000) + @test size(downconvert_and_correlator2.buffers[1][1].code_replica_buffer) == (5004,) + @test eltype(downconvert_and_correlator2.buffers[1][1].code_replica_buffer) == Int16 + @test size(downconvert_and_correlator2.buffers[1][1].carrier_replica_buffer) == (5000,) + @test eltype(downconvert_and_correlator2.buffers[1][1].carrier_replica_buffer) == ComplexF32 + @test size(downconvert_and_correlator2.buffers[1][1].downconvert_signal_buffer) == (5000,) + @test eltype(downconvert_and_correlator2.buffers[1][1].downconvert_signal_buffer) == ComplexF32 + + @test size(downconvert_and_correlator2.buffers[2][1].code_replica_buffer) == (5006,) + @test eltype(downconvert_and_correlator2.buffers[2][1].code_replica_buffer) == Float32 + @test size(downconvert_and_correlator2.buffers[2][1].carrier_replica_buffer) == (5000,) + @test eltype(downconvert_and_correlator2.buffers[2][1].carrier_replica_buffer) == ComplexF32 + @test size(downconvert_and_correlator2.buffers[2][1].downconvert_signal_buffer) == (5000,) + @test eltype(downconvert_and_correlator2.buffers[2][1].downconvert_signal_buffer) == ComplexF32 +end + +@testset "Downconvert and correlate with $type" for type in (:CPU, :GPU) + arch_function = (CPU = (downconvert_and_correlator = CPUDownconvertAndCorrelator, array_transform = Array), GPU = (downconvert_and_correlator = GPUDownconvertAndCorrelator, array_transform = cu)) + type == :GPU && !CUDA.functional() && return + gpsl1 = GPSL1() + sampling_frequency = 5e6Hz + code_phase = 10.5 + num_samples_signal = 5000 + intermediate_frequency = 0.0Hz + + system_sats_state = SystemSatsState( + gpsl1, + [ + SatState(gpsl1, 1, sampling_frequency, code_phase, 1000.0Hz), + SatState(gpsl1, 2, sampling_frequency, 11.0, 500.0Hz), + ], + ) + system_sats_states = (system_sats_state,) + + downconvert_and_correlator = @inferred arch_function[type].downconvert_and_correlator(system_sats_states, num_samples_signal) + + system_sats_sample_params = Tracking.init_sample_params( + system_sats_states, + 1, + ) + next_system_sats_sample_params = Tracking.calc_sample_params( + system_sats_states, + system_sats_sample_params, + num_samples_signal, + sampling_frequency, + 1, + ) + + signal = arch_function[type].array_transform( + gen_code( + num_samples_signal, + gpsl1, + 1, + sampling_frequency, + get_code_frequency(gpsl1) + 1000Hz * get_code_center_frequency_ratio(gpsl1), + code_phase, + ) .* + cis.(2π * (0:num_samples_signal-1) * 1000.0Hz / sampling_frequency), + ) + + correlators = @inferred Tracking.downconvert_and_correlate( + downconvert_and_correlator, + signal, + sampling_frequency, + intermediate_frequency, + system_sats_states, + next_system_sats_sample_params, + ) + + # GPU uses floating point arithmetic and might differ a little with the fixed point arithmetic + @test real.(correlators[1][1].accumulators) ≈ [2921, 4949, 2917] || real.(correlators[1][1].accumulators) ≈ [2921, 4949, 2921] + + signal = arch_function[type].array_transform(gen_code( + num_samples_signal, + gpsl1, + 2, + sampling_frequency, + get_code_frequency(gpsl1) + 500Hz * get_code_center_frequency_ratio(gpsl1), + 11.0, + ) .* + cis.(2π * (0:num_samples_signal-1) * 500.0Hz / sampling_frequency)) + + correlators = @inferred Tracking.downconvert_and_correlate( + downconvert_and_correlator, + signal, + sampling_frequency, + intermediate_frequency, + system_sats_states, + next_system_sats_sample_params, + ) + + # GPU uses floating point arithmetic and might differ a little with the fixed point arithmetic + @test real.(correlators[1][2].accumulators) ≈ [2919, 4947, 2915] || real.(correlators[1][2].accumulators) ≈ [2919, 4947, 2919] +end diff --git a/test/galileo_e1b.jl b/test/galileo_e1b.jl index 46193b7..86a66ca 100644 --- a/test/galileo_e1b.jl +++ b/test/galileo_e1b.jl @@ -10,6 +10,10 @@ @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf0, 8)) == true - @test @inferred(Tracking.get_default_correlator(galileo_e1b, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_default_correlator(galileo_e1b, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3)) + sampling_frequency = 5e6Hz + + @test @inferred(Tracking.get_default_correlator(galileo_e1b, sampling_frequency, NumAnts(1))) == + VeryEarlyPromptLateCorrelator(galileo_e1b, sampling_frequency, num_ants = NumAnts(1)) + @test @inferred(Tracking.get_default_correlator(galileo_e1b, sampling_frequency, NumAnts(3))) == + VeryEarlyPromptLateCorrelator(galileo_e1b, sampling_frequency, num_ants = NumAnts(3)) end diff --git a/test/gps_l1.jl b/test/gps_l1.jl index 79eacf0..bf98233 100644 --- a/test/gps_l1.jl +++ b/test/gps_l1.jl @@ -6,6 +6,10 @@ @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl1, 0xfffff, 40)) == true - @test @inferred(Tracking.get_default_correlator(gpsl1, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_default_correlator(gpsl1, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3)) + sampling_frequency = 5e6Hz + + @test @inferred(Tracking.get_default_correlator(gpsl1, sampling_frequency, NumAnts(1))) == + EarlyPromptLateCorrelator(gpsl1, sampling_frequency, num_ants = NumAnts(1)) + @test @inferred(Tracking.get_default_correlator(gpsl1, sampling_frequency, NumAnts(3))) == + EarlyPromptLateCorrelator(gpsl1, sampling_frequency, num_ants = NumAnts(3)) end diff --git a/test/gps_l5.jl b/test/gps_l5.jl index 755511e..bffff0c 100644 --- a/test/gps_l5.jl +++ b/test/gps_l5.jl @@ -6,6 +6,10 @@ @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl5, 0x3ca, 10)) == true # 0x3ca == 1111001010 - @test @inferred(Tracking.get_default_correlator(gpsl5, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_default_correlator(gpsl5, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3)) + sampling_frequency = 5e6Hz + + @test @inferred(Tracking.get_default_correlator(gpsl5, sampling_frequency, NumAnts(1))) == + EarlyPromptLateCorrelator(gpsl5, sampling_frequency, num_ants = NumAnts(1)) + @test @inferred(Tracking.get_default_correlator(gpsl5, sampling_frequency, NumAnts(3))) == + EarlyPromptLateCorrelator(gpsl5, sampling_frequency, num_ants = NumAnts(3)) end diff --git a/test/post_corr_filter.jl b/test/post_corr_filter.jl index e9d8f03..3cf6269 100644 --- a/test/post_corr_filter.jl +++ b/test/post_corr_filter.jl @@ -1,5 +1,5 @@ @testset "Post correlation filter" begin - default_post_corr_filter = Tracking.DefaultPostCorrFilter() + default_post_corr_filter = DefaultPostCorrFilter() @test Tracking.update(default_post_corr_filter, randn(ComplexF64)) == default_post_corr_filter diff --git a/test/runtests.jl b/test/runtests.jl index 6e39051..94ffcc6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,10 @@ 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") +include("downconvert_and_correlate.jl") include("post_corr_filter.jl") include("code_replica.jl") include("carrier_replica.jl") @@ -22,20 +26,6 @@ include("galileo_e1b.jl") include("secondary_code_or_bit_detector.jl") include("bit_buffer.jl") include("correlator.jl") -include("tracking_state.jl") -include("tracking_results.jl") -include("tracking_loop.jl") include("cn0_estimation.jl") - -if CUDA.functional() - include("cuda/bit_buffer.jl") - include("cuda/cn0_estimation.jl") - include("cuda/discriminators.jl") - include("cuda/galileo_e1b.jl") - include("cuda/gps_l1.jl") - include("cuda/gps_l5.jl") - include("cuda/secondary_code_or_bit_detector.jl") - include("cuda/tracking_loop.jl") - include("cuda/tracking_results.jl") - include("cuda/tracking_state.jl") -end +include("tracking_state.jl") +include("track.jl") \ No newline at end of file diff --git a/test/sample_parameters.jl b/test/sample_parameters.jl new file mode 100644 index 0000000..229126e --- /dev/null +++ b/test/sample_parameters.jl @@ -0,0 +1,146 @@ +@testset "Sample parameters" begin + + gpsl1 = GPSL1() + + @test @inferred(Tracking.calc_num_code_blocks_to_integrate(gpsl1, 1, false)) == 1 + @test @inferred(Tracking.calc_num_code_blocks_to_integrate(gpsl1, 2, false)) == 1 + @test @inferred(Tracking.calc_num_code_blocks_to_integrate(gpsl1, 2, true)) == 2 + @test @inferred(Tracking.calc_num_code_blocks_to_integrate(gpsl1, 20, true)) == 20 + @test @inferred(Tracking.calc_num_code_blocks_to_integrate(gpsl1, 21, true)) == 20 + + @test @inferred(Tracking.calc_num_chips_to_integrate(gpsl1, 1, 10.5)) == 1023 - 10.5 + @test @inferred(Tracking.calc_num_chips_to_integrate(gpsl1, 1, 1023 + 10.5)) == 1023 - 10.5 + @test @inferred(Tracking.calc_num_chips_to_integrate(gpsl1, 2, 10.5)) == 2 * 1023 - 10.5 + + code_freq = get_code_frequency(gpsl1) + @test @inferred(Tracking.calc_num_samples_left_to_integrate(gpsl1, 1, 5e6Hz, 0.0Hz, 10.5)) == + ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + + sample_params = @inferred Tracking.SampleParams( + gpsl1, + 1, + false + ) + + @test sample_params.signal_samples_to_integrate == 0 + @test sample_params.signal_start_sample == 1 + @test sample_params.samples_to_integrate_until_code_end == 0 + @test sample_params.num_code_blocks_to_integrate == 1 + + num_samples_signal = 5000 + sampling_frequency = 5e6Hz + code_doppler = 0.0Hz + code_phase = 10.5 + preferred_num_code_blocks_to_integrate = 1 + secondary_code_or_bit_found = false + + next_sample_params = @inferred Tracking.SampleParams( + gpsl1, + sample_params, + num_samples_signal, + sampling_frequency, + code_doppler, + code_phase, + preferred_num_code_blocks_to_integrate, + secondary_code_or_bit_found + ) + + @test next_sample_params.signal_samples_to_integrate == ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + @test next_sample_params.signal_start_sample == 1 + @test next_sample_params.samples_to_integrate_until_code_end == ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + @test next_sample_params.num_code_blocks_to_integrate == 1 + + # Bit not found -> It should default to one code block to integrate + next_sample_params = @inferred Tracking.SampleParams( + gpsl1, + sample_params, + num_samples_signal, + sampling_frequency, + code_doppler, + code_phase, + 2, + secondary_code_or_bit_found + ) + + @test next_sample_params.signal_samples_to_integrate == ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + @test next_sample_params.signal_start_sample == 1 + @test next_sample_params.samples_to_integrate_until_code_end == ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + @test next_sample_params.num_code_blocks_to_integrate == 1 + + # Bit found + next_sample_params = @inferred Tracking.SampleParams( + gpsl1, + sample_params, + num_samples_signal, + sampling_frequency, + code_doppler, + code_phase, + 2, + true + ) + + @test next_sample_params.signal_samples_to_integrate == 5000 + @test next_sample_params.signal_start_sample == 1 + @test next_sample_params.samples_to_integrate_until_code_end == ceil(Int, (2 * 1023 - 10.5) * 5e6Hz / code_freq) + @test next_sample_params.num_code_blocks_to_integrate == 2 + + system_sats_state = SystemSatsState( + gpsl1, + [ + SatState(gpsl1, 1, sampling_frequency, code_phase, 0.0Hz), + SatState(gpsl1, 2, sampling_frequency, 11.5, 0.0Hz) + ] + ) + + system_sats_sample_params = @inferred Tracking.init_sample_params( + (system_sats_state,), + 1, + ) + + @test system_sats_sample_params[1][1].signal_samples_to_integrate == 0 + @test system_sats_sample_params[1][1].signal_start_sample == 1 + @test system_sats_sample_params[1][1].samples_to_integrate_until_code_end == 0 + @test system_sats_sample_params[1][1].num_code_blocks_to_integrate == 1 + + @test system_sats_sample_params[1][2].signal_samples_to_integrate == 0 + @test system_sats_sample_params[1][2].signal_start_sample == 1 + @test system_sats_sample_params[1][2].samples_to_integrate_until_code_end == 0 + @test system_sats_sample_params[1][2].num_code_blocks_to_integrate == 1 + + + next_system_sats_sample_params = @inferred Tracking.calc_sample_params( + (system_sats_state,), + system_sats_sample_params, + num_samples_signal, + sampling_frequency, + preferred_num_code_blocks_to_integrate, + ) + + @test next_system_sats_sample_params[1][1].signal_samples_to_integrate == ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + @test next_system_sats_sample_params[1][1].signal_start_sample == 1 + @test next_system_sats_sample_params[1][1].samples_to_integrate_until_code_end == ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + @test next_system_sats_sample_params[1][1].num_code_blocks_to_integrate == 1 + + @test next_system_sats_sample_params[1][2].signal_samples_to_integrate == ceil(Int, (1023 - 11.5) * 5e6Hz / code_freq) + @test next_system_sats_sample_params[1][2].signal_start_sample == 1 + @test next_system_sats_sample_params[1][2].samples_to_integrate_until_code_end == ceil(Int, (1023 - 11.5) * 5e6Hz / code_freq) + @test next_system_sats_sample_params[1][2].num_code_blocks_to_integrate == 1 + + next_next_system_sats_sample_params = @inferred Tracking.calc_sample_params( + (system_sats_state,), + next_system_sats_sample_params, + num_samples_signal, + sampling_frequency, + preferred_num_code_blocks_to_integrate, + ) + + @test next_next_system_sats_sample_params[1][1].signal_samples_to_integrate == 5000 - ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + @test next_next_system_sats_sample_params[1][1].signal_start_sample == ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + 1 + @test next_next_system_sats_sample_params[1][1].samples_to_integrate_until_code_end == ceil(Int, (1023 - 10.5) * 5e6Hz / code_freq) + @test next_next_system_sats_sample_params[1][1].num_code_blocks_to_integrate == 1 + + @test next_next_system_sats_sample_params[1][2].signal_samples_to_integrate == 5000 - ceil(Int, (1023 - 11.5) * 5e6Hz / code_freq) + @test next_next_system_sats_sample_params[1][2].signal_start_sample == ceil(Int, (1023 - 11.5) * 5e6Hz / code_freq) + 1 + @test next_next_system_sats_sample_params[1][2].samples_to_integrate_until_code_end == ceil(Int, (1023 - 11.5) * 5e6Hz / code_freq) + @test next_next_system_sats_sample_params[1][2].num_code_blocks_to_integrate == 1 +end \ No newline at end of file diff --git a/test/sat_state.jl b/test/sat_state.jl new file mode 100644 index 0000000..44e4480 --- /dev/null +++ b/test/sat_state.jl @@ -0,0 +1,196 @@ +@testset "Satellite state" begin + + gpsl1 = GPSL1() + sat_state = @inferred SatState(gpsl1, 1, 5e6Hz, 10.0, 500.0Hz) + @test get_prn(sat_state) == 1 + @test get_code_phase(sat_state) == 10.0 + @test get_code_doppler(sat_state) == 500.0Hz * get_code_center_frequency_ratio(gpsl1) + @test get_carrier_phase(sat_state) == 0.0 + @test get_carrier_doppler(sat_state) == 500.0Hz + @test get_integrated_samples(sat_state) == 0.0 + @test get_correlator(sat_state).accumulators == zeros(3) + @test get_last_fully_integrated_correlator(sat_state).accumulators == zeros(3) + @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(gpsl1, 5, 5e6Hz, 100.0Hz, 524.6, 45.0, 1.0, randn(100,100), -500:100.0:500) + sat_state = @inferred SatState(acq) + @test get_prn(sat_state) == 5 + @test get_code_phase(sat_state) == 524.6 + @test get_code_doppler(sat_state) == 100.0Hz * get_code_center_frequency_ratio(gpsl1) + @test get_carrier_phase(sat_state) == 0.0 + @test get_carrier_doppler(sat_state) == 100.0Hz + @test get_integrated_samples(sat_state) == 0.0 + @test get_correlator(sat_state).accumulators == zeros(3) + @test get_last_fully_integrated_correlator(sat_state).accumulators == zeros(3) + @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 + + system_sats_state = @inferred SystemSatsState( + gpsl1, + [ + SatState( + 1, + 10.5, + 0.0Hz, + 0.49, + 0.0Hz, + 150, + EarlyPromptLateCorrelator(complex.(ones(3), ones(3)), -1:1, 3, 2, 1), + EarlyPromptLateCorrelator(complex.(zeros(3), zeros(3)), -1:1, 3, 2, 1), + complex(0.0, 0.0), + 14, + Tracking.SecondaryCodeOrBitDetector(), + Tracking.PromptsBuffer(20), + Tracking.BitBuffer() + ), + SatState( + 2, + 11.5, + 1.0Hz, + 0.8, + 10.0Hz, + 100, + EarlyPromptLateCorrelator(complex.(ones(3) * 2, ones(3) * 3), -1:1, 3, 2, 1), + EarlyPromptLateCorrelator(complex.(zeros(3), zeros(3)), -1:1, 3, 2, 1), + complex(0.0, 0.0), + 15, + Tracking.SecondaryCodeOrBitDetector(), + Tracking.PromptsBuffer(20), + Tracking.BitBuffer() + ), + SatState( + 3, + 12.5, + 2.0Hz, + -0.5, + 21.0Hz, + 10, + EarlyPromptLateCorrelator(complex.(ones(3) * 3, ones(3) * 4), -1:1, 3, 2, 1), + EarlyPromptLateCorrelator(complex.(ones(3), zeros(3)), -1:1, 3, 2, 1), + complex(0.0, 0.0), + -1, + Tracking.SecondaryCodeOrBitDetector(), + Tracking.PromptsBuffer(20), + Tracking.BitBuffer() + ) + ] + ) + + dopplers_and_filtered_prompts = [ + Tracking.DopplersAndFilteredPrompt(200.0Hz, 2.0Hz, complex(1.0, 2.0)), + Tracking.DopplersAndFilteredPrompt(150.0Hz, 1.5Hz, complex(2.0, 5.0)), + Tracking.DopplersAndFilteredPrompt(0.0Hz, 0.0Hz, complex(0.0, 0.0)), + ] + + system_sats_sample_params = Tracking.init_sample_params( + (system_sats_state,), + 1, + ) + next_system_sats_sample_params = Tracking.calc_sample_params( + (system_sats_state,), + system_sats_sample_params, + 5000, + sampling_frequency, + 1, + ) + + next_system_sats_states = @inferred Tracking.update( + (system_sats_state,), + (dopplers_and_filtered_prompts,), + next_system_sats_sample_params, + ) + + state1 = next_system_sats_states[1].states[1] + @test get_carrier_doppler(state1) == 200.0Hz + @test get_code_doppler(state1) == 2.0Hz + @test get_last_fully_integrated_filtered_prompt(state1) == complex(1.0, 2.0) + @test get_integrated_samples(state1) == 0 + @test get_correlator(state1).accumulators == zeros(3) + @test get_last_fully_integrated_correlator(state1).accumulators == complex.(ones(3), ones(3)) + @test get_sample_of_last_fully_integrated_correlator(state1) == 4949 + + state2 = next_system_sats_states[1].states[2] + @test get_carrier_doppler(state2) == 150.0Hz + @test get_code_doppler(state2) == 1.5Hz + @test get_last_fully_integrated_filtered_prompt(state2) == complex(2.0, 5.0) + @test get_integrated_samples(state2) == 0 + @test get_correlator(state2).accumulators == zeros(3) + @test get_last_fully_integrated_correlator(state2).accumulators == complex.(ones(3) * 2, ones(3) * 3) + @test get_sample_of_last_fully_integrated_correlator(state2) == 4944 + + state3 = next_system_sats_states[1].states[3] + @test get_carrier_doppler(state3) == 21.0Hz + @test get_code_doppler(state3) == 2.0Hz + @test get_last_fully_integrated_filtered_prompt(state3) == complex(0.0, 0.0) + @test get_integrated_samples(state3) == 10 + @test get_correlator(state3).accumulators == complex.(ones(3) * 3, ones(3) * 4) + @test get_last_fully_integrated_correlator(state3).accumulators == complex.(ones(3), zeros(3)) + @test get_sample_of_last_fully_integrated_correlator(state3) == -1 + + correlators = [ + EarlyPromptLateCorrelator(complex.(ones(3) * 3, ones(3) * 4), -1:1, 3, 2, 1) + EarlyPromptLateCorrelator(complex.(ones(3) * 2, ones(3) * 3), -1:1, 3, 2, 1) + EarlyPromptLateCorrelator(complex.(ones(3) * 1, ones(3) * 2), -1:1, 3, 2, 1) + ] + + next_system_sats_states = @inferred Tracking.update( + (system_sats_state,), + sampling_frequency, + 0.0Hz, + (correlators,), + next_system_sats_sample_params, + 5000 + ) + + state1 = next_system_sats_states[1].states[1] + state2 = next_system_sats_states[1].states[2] + state3 = next_system_sats_states[1].states[3] + @test get_code_phase(state1) == mod(10.5 + next_system_sats_sample_params[1][1].signal_samples_to_integrate * (get_code_frequency(gpsl1) + 0.0Hz) / sampling_frequency, 1023) + @test get_code_phase(state2) == mod(11.5 + next_system_sats_sample_params[1][2].signal_samples_to_integrate * (get_code_frequency(gpsl1) + 1.0Hz) / sampling_frequency, 1023) + @test get_code_phase(state3) == mod(12.5 + next_system_sats_sample_params[1][3].signal_samples_to_integrate * (get_code_frequency(gpsl1) + 2.0Hz) / sampling_frequency, 1023) + + @test get_carrier_phase(state1) ≈ mod2pi(2π * (0.49 + next_system_sats_sample_params[1][1].signal_samples_to_integrate * 0.0Hz / sampling_frequency) + π) - π + @test get_carrier_phase(state2) ≈ mod2pi(2π * (0.8 + next_system_sats_sample_params[1][2].signal_samples_to_integrate * 10.0Hz / sampling_frequency) + π) - π + @test get_carrier_phase(state3) ≈ mod2pi(2π * (-0.5 + next_system_sats_sample_params[1][3].signal_samples_to_integrate * 21.0Hz / sampling_frequency) + π) - π + + @test get_correlator(state1) == correlators[1] + @test get_correlator(state2) == correlators[2] + @test get_correlator(state3) == correlators[3] + + @test get_sample_of_last_fully_integrated_correlator(state1) == 14 - 5000 + @test get_sample_of_last_fully_integrated_correlator(state2) == 15 - 5000 + @test get_sample_of_last_fully_integrated_correlator(state3) == -1 - 5000 + + next_next_system_sats_sample_params = @inferred Tracking.calc_sample_params( + (system_sats_state,), + next_system_sats_sample_params, + 5000, + sampling_frequency, + 1, + ) + + next_system_sats_states = @inferred Tracking.update( + (system_sats_state,), + sampling_frequency, + 0.0Hz, + (correlators,), + next_next_system_sats_sample_params, + 5000 + ) + + state1 = next_system_sats_states[1].states[1] + state2 = next_system_sats_states[1].states[2] + state3 = next_system_sats_states[1].states[3] + @test get_sample_of_last_fully_integrated_correlator(state1) == 14 + @test get_sample_of_last_fully_integrated_correlator(state2) == 15 + @test get_sample_of_last_fully_integrated_correlator(state3) == -1 +end \ No newline at end of file diff --git a/test/track.jl b/test/track.jl new file mode 100644 index 0000000..c3967fe --- /dev/null +++ b/test/track.jl @@ -0,0 +1,491 @@ +@testset "Tracking with signal of type $type" for type in (Int16, Int32, Int64, Float32, Float64) + gpsl1 = GPSL1() + carrier_doppler = 200Hz + start_code_phase = 100 + code_frequency = carrier_doppler * get_code_center_frequency_ratio(gpsl1) + get_code_frequency(gpsl1) + sampling_frequency = 4e6Hz + prn = 1 + range = 0:3999 + start_carrier_phase = π / 2 + + num_samples = 4000 + + track_state = @inferred TrackState( + gpsl1, + [ + SatState(gpsl1, 1, sampling_frequency, start_code_phase, carrier_doppler - 20Hz), + ]; + num_samples, + ) + + signal_temp = cis.( + 2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase + ) .* + gen_code(4000, gpsl1, prn, sampling_frequency, code_frequency, start_code_phase) + scaling = 512 + signal = type <: Integer ? complex.(round.(type, real.(signal_temp) * scaling), round.(type, imag.(signal_temp) * scaling)) : Complex{type}.(signal_temp) + + track_state = @inferred track( + signal, + track_state, + sampling_frequency + ) + + iterations = 20000 + code_phases = zeros(iterations) + carrier_phases = zeros(iterations) + tracked_code_phases = zeros(iterations) + tracked_carrier_phases = zeros(iterations) + tracked_code_dopplers = zeros(iterations) + tracked_carrier_dopplers = zeros(iterations) + tracked_prompts = zeros(ComplexF64, iterations) + for i = 1:iterations + carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency + + start_carrier_phase + π) - π + code_phase = mod( + code_frequency * 4000 * i / sampling_frequency + start_code_phase, + 1023 + ) + signal_temp = cis.( + 2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase + ) .* + gen_code(4000, gpsl1, prn, sampling_frequency, code_frequency, code_phase) + signal = type <: Integer ? complex.(round.(type, real.(signal_temp) * scaling), round.(type, imag.(signal_temp) * scaling)) : Complex{type}.(signal_temp) + track_state = @inferred track( + signal, + track_state, + sampling_frequency + ) + comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) / + sampling_frequency + start_carrier_phase + π) - π + comp_code_phase = mod( + code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, + 1023 + ) + tracked_code_phases[i] = get_code_phase(track_state.system_sats_states[1].states[1]) + tracked_carrier_phases[i] = get_carrier_phase(track_state.system_sats_states[1].states[1]) + tracked_carrier_dopplers[i] = get_carrier_doppler(track_state.system_sats_states[1].states[1])/Hz + tracked_code_dopplers[i] = get_code_doppler(track_state.system_sats_states[1].states[1])/Hz + tracked_prompts[i] = get_last_fully_integrated_filtered_prompt(track_state.system_sats_states[1].states[1]) + code_phases[i] = comp_code_phase + carrier_phases[i] = comp_carrier_phase + end + @test tracked_code_phases[end] ≈ code_phases[end] atol = 5e-5 + @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 1e-3 + +# using PyPlot +# pygui(true) +# figure("carrier_phases") +# plot(tracked_carrier_phases) +# plot(carrier_phases) +# grid(true) +# figure("Code Phases") +# plot(300 * (tracked_code_phases .- code_phases)) +# figure("Carrier Doppler") +# plot(tracked_carrier_dopplers) +# figure("Code Doppler") +# plot(tracked_code_dopplers) +# figure("Prompt") +# plot(real.(tracked_prompts)) +# plot(imag.(tracked_prompts)) +end + +@testset "Track multiple systems of type $type" for type in (Int16, Int32, Int64, Float32, Float64) + gpsl1 = GPSL1() + galileo_e1b = GalileoE1B() + carrier_doppler_gps = 200Hz + carrier_doppler_gal = 1200Hz + start_code_phase = 100 + code_frequency_gps = carrier_doppler_gps * get_code_center_frequency_ratio(gpsl1) + get_code_frequency(gpsl1) + code_frequency_gal = carrier_doppler_gal * get_code_center_frequency_ratio(galileo_e1b) + get_code_frequency(galileo_e1b) + sampling_frequency = 4e6Hz + prn = 1 + range = 0:3999 + start_carrier_phase = π / 2 + + num_samples = 4000 + + track_state = @inferred TrackState( + ( + gps = SystemSatsState( + gpsl1, + [ + SatState(gpsl1, prn, sampling_frequency, start_code_phase, carrier_doppler_gps), + ] + ), + gal = SystemSatsState( + galileo_e1b, + [ + SatState(galileo_e1b, prn, sampling_frequency, start_code_phase, carrier_doppler_gal), + ] + ) + ); + num_samples, + ) + + signal_temp = cis.( + 2π .* carrier_doppler_gps .* range ./ sampling_frequency .+ start_carrier_phase + ) .* + gen_code(4000, gpsl1, prn, sampling_frequency, code_frequency_gps, start_code_phase) .+ + cis.( + 2π .* carrier_doppler_gal .* range ./ sampling_frequency .+ start_carrier_phase + ) .* + gen_code(4000, galileo_e1b, prn, sampling_frequency, code_frequency_gal, start_code_phase) + scaling = 5000 + signal = type <: Integer ? complex.(round.(type, real.(signal_temp) * scaling), round.(type, imag.(signal_temp) * scaling)) : Complex{type}.(signal_temp) + + track_state = @inferred track( + signal, + track_state, + sampling_frequency + ) + + iterations = 2000 + for i = 1:iterations + carrier_phase_gps = mod2pi(2π * carrier_doppler_gps * 4000 * i / sampling_frequency + + start_carrier_phase + π) - π + code_phase_gps = mod( + code_frequency_gps * 4000 * i / sampling_frequency + start_code_phase, + get_code_length(gpsl1) + ) + carrier_phase_gal = mod2pi(2π * carrier_doppler_gal * 4000 * i / sampling_frequency + + start_carrier_phase + π) - π + code_phase_gal = mod( + code_frequency_gal * 4000 * i / sampling_frequency + start_code_phase, + get_code_length(galileo_e1b) + ) + signal_temp = cis.( + 2π .* carrier_doppler_gps .* range ./ sampling_frequency .+ carrier_phase_gps + ) .* + gen_code(4000, gpsl1, prn, sampling_frequency, code_frequency_gps, code_phase_gps) .+ + cis.( + 2π .* carrier_doppler_gal .* range ./ sampling_frequency .+ carrier_phase_gal + ) .* + gen_code(4000, galileo_e1b, prn, sampling_frequency, code_frequency_gal, code_phase_gal) + signal = type <: Integer ? complex.(round.(type, real.(signal_temp) * scaling), round.(type, imag.(signal_temp) * scaling)) : Complex{type}.(signal_temp) + track_state = @inferred track( + signal, + track_state, + sampling_frequency + ) + end + comp_carrier_phase_gps = mod2pi(2π * carrier_doppler_gps * 4000 * (iterations + 1) / + sampling_frequency + start_carrier_phase + π) - π + comp_code_phase_gps = mod( + code_frequency_gps * 4000 * (iterations + 1) / sampling_frequency + start_code_phase, + get_code_length(gpsl1) + ) + comp_carrier_phase_gal = mod2pi(2π * carrier_doppler_gal * 4000 * (iterations + 1) / + sampling_frequency + start_carrier_phase + π) - π + comp_code_phase_gal = mod( + code_frequency_gal * 4000 * (iterations + 1) / sampling_frequency + start_code_phase, + get_code_length(galileo_e1b) + ) + @test get_code_phase(track_state.system_sats_states.gps.states[1]) ≈ comp_code_phase_gps atol = 5e-3 + @test mod(get_carrier_phase(track_state.system_sats_states.gps.states[1]), π) ≈ comp_carrier_phase_gps atol = 3e-3 + @test get_code_phase(track_state.system_sats_states.gal.states[1]) ≈ comp_code_phase_gal atol = 5e-3 + @test mod(get_carrier_phase(track_state.system_sats_states.gal.states[1]), π) ≈ comp_carrier_phase_gal atol = 3e-3 +end + +@testset "Tracking with intermediate frequency of $intermediate_frequency" for intermediate_frequency in (0.0Hz, -10000.0Hz, 10000.0Hz, -30000.0Hz, 30000.0Hz) + gpsl1 = GPSL1() + carrier_doppler = 200Hz + start_code_phase = 100 + code_frequency = carrier_doppler / 1540 + get_code_frequency(gpsl1) + sampling_frequency = 4e6Hz + prn = 1 + range = 0:3999 + start_carrier_phase = π / 2 + + num_samples = 4000 + + track_state = @inferred TrackState( + gpsl1, + [ + SatState(gpsl1, 1, sampling_frequency, start_code_phase, carrier_doppler - 20Hz), + ]; + num_samples, + ) + + signal = cis.( + 2π .* (carrier_doppler + intermediate_frequency) .* range ./ sampling_frequency .+ start_carrier_phase + ) .* + get_code.( + gpsl1, + code_frequency .* range ./ sampling_frequency .+ start_code_phase, + prn + ) + + track_state = @inferred track( + signal, + track_state, + sampling_frequency; + intermediate_frequency + ) + + iterations = 2000 + code_phases = zeros(iterations) + carrier_phases = zeros(iterations) + tracked_code_phases = zeros(iterations) + tracked_carrier_phases = zeros(iterations) + tracked_code_dopplers = zeros(iterations) + tracked_carrier_dopplers = zeros(iterations) + tracked_prompts = zeros(ComplexF64, iterations) + for i = 1:iterations + carrier_phase = mod2pi(2π * (carrier_doppler + intermediate_frequency) * 4000 * i / sampling_frequency + + start_carrier_phase + π) - π + code_phase = mod( + code_frequency * 4000 * i / sampling_frequency + start_code_phase, + 1023 + ) + signal = cis.( + 2π .* (carrier_doppler + intermediate_frequency) .* range ./ sampling_frequency .+ carrier_phase + ) .* + get_code.( + gpsl1, + code_frequency .* range ./ sampling_frequency .+ code_phase, + prn + ) + track_state = @inferred track( + signal, + track_state, + sampling_frequency; + intermediate_frequency + ) + comp_carrier_phase = mod2pi(2π * (carrier_doppler + intermediate_frequency) * 4000 * (i + 1) / + sampling_frequency + start_carrier_phase + π) - π + comp_code_phase = mod( + code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, + 1023 + ) + tracked_code_phases[i] = get_code_phase(track_state.system_sats_states[1].states[1]) + tracked_carrier_phases[i] = get_carrier_phase(track_state.system_sats_states[1].states[1]) + tracked_carrier_dopplers[i] = get_carrier_doppler(track_state.system_sats_states[1].states[1])/Hz + tracked_code_dopplers[i] = get_code_doppler(track_state.system_sats_states[1].states[1])/Hz + tracked_prompts[i] = get_last_fully_integrated_filtered_prompt(track_state.system_sats_states[1].states[1]) + code_phases[i] = comp_code_phase + carrier_phases[i] = comp_carrier_phase + end + @test tracked_code_phases[end] ≈ code_phases[end] atol = 5e-5 + @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-5 + +# using PyPlot +# pygui(true) +# figure("carrier_phases") +# plot(tracked_carrier_phases) +# plot(carrier_phases) +# grid(true) +# figure("Code Phases") +# plot(300 * (tracked_code_phases .- code_phases)) +# figure("Carrier Doppler") +# plot(tracked_carrier_dopplers) +# figure("Code Doppler") +# plot(tracked_code_dopplers) +# figure("Prompt") +# plot(real.(tracked_prompts)) +# plot(imag.(tracked_prompts)) +end + +@testset "Track multiple signals with signal of type $type" for type in (Int16, Int32, Int64, Float32, Float64) + gpsl1 = GPSL1() + carrier_doppler = 200Hz + start_code_phase = 100 + code_frequency = carrier_doppler / 1540 + get_code_frequency(gpsl1) + sampling_frequency = 4e6Hz + prn = 1 + range = 0:3999 + start_carrier_phase = π / 2 + + num_samples = 4000 + + track_state = @inferred TrackState( + gpsl1, + [ + SatState(gpsl1, 1, sampling_frequency, start_code_phase, carrier_doppler - 20Hz, num_ants = NumAnts(3)), + ]; + num_samples, + ) + + signal = cis.( + 2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase + ) .* + get_code.( + gpsl1, + code_frequency .* range ./ sampling_frequency .+ start_code_phase, + prn + ) + signal_mat_temp = repeat(signal, outer = (1,3)) + scaling = 512 + signal_mat = type <: Integer ? complex.(floor.(type, real.(signal_mat_temp) * scaling), floor.(type, imag.(signal_mat_temp) * scaling)) : Complex{type}.(signal_mat_temp) + + track_state = @inferred track( + signal_mat, + track_state, + sampling_frequency + ) + + iterations = 2000 + code_phases = zeros(iterations) + carrier_phases = zeros(iterations) + tracked_code_phases = zeros(iterations) + tracked_carrier_phases = zeros(iterations) + tracked_code_dopplers = zeros(iterations) + tracked_carrier_dopplers = zeros(iterations) + tracked_prompts = zeros(ComplexF64, iterations) + for i = 1:iterations + carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency + + start_carrier_phase + π) - π + code_phase = mod( + code_frequency * 4000 * i / sampling_frequency + start_code_phase, + 1023 + ) + signal = cis.( + 2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase + ) .* + get_code.( + gpsl1, + code_frequency .* range ./ sampling_frequency .+ code_phase, + prn + ) + signal_mat = repeat(signal, outer = (1,3)) + track_state = @inferred track( + signal_mat, + track_state, + sampling_frequency + ) + comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) / + sampling_frequency + start_carrier_phase + π) - π + comp_code_phase = mod( + code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, + 1023 + ) + tracked_code_phases[i] = get_code_phase(track_state.system_sats_states[1].states[1]) + tracked_carrier_phases[i] = get_carrier_phase(track_state.system_sats_states[1].states[1]) + tracked_carrier_dopplers[i] = get_carrier_doppler(track_state.system_sats_states[1].states[1])/Hz + tracked_code_dopplers[i] = get_code_doppler(track_state.system_sats_states[1].states[1])/Hz + tracked_prompts[i] = get_last_fully_integrated_filtered_prompt(track_state.system_sats_states[1].states[1]) + code_phases[i] = comp_code_phase + carrier_phases[i] = comp_carrier_phase + end + @test tracked_code_phases[end] ≈ code_phases[end] atol = 5e-5 + @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-5 + +# using PyPlot +# pygui(true) +# figure("carrier_phases") +# plot(tracked_carrier_phases) +# plot(carrier_phases) +# grid(true) +# figure("Code Phases") +# plot(300 * (tracked_code_phases .- code_phases)) +# figure("Carrier Doppler") +# plot(tracked_carrier_dopplers) +# figure("Code Doppler") +# plot(tracked_code_dopplers) +# figure("Prompt") +# plot(real.(tracked_prompts)) +# plot(imag.(tracked_prompts)) +end + +@testset "Track multiple signals with GPU" begin + !CUDA.functional() && return + gpsl1 = GPSL1() + carrier_doppler = 200Hz + start_code_phase = 100 + code_frequency = carrier_doppler / 1540 + get_code_frequency(gpsl1) + sampling_frequency = 4e6Hz + prn = 1 + range = 0:3999 + start_carrier_phase = π / 2 + + num_samples = 4000 + + sat_states = [ + SatState(gpsl1, 1, sampling_frequency, start_code_phase, carrier_doppler - 20Hz, num_ants = NumAnts(3)), + ] + + track_state = @inferred TrackState( + gpsl1, + sat_states; + num_samples, + downconvert_and_correlator = GPUDownconvertAndCorrelator((SystemSatsState(gpsl1, sat_states),), num_samples) + ) + + signal = cis.( + 2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase + ) .* + get_code.( + gpsl1, + code_frequency .* range ./ sampling_frequency .+ start_code_phase, + prn + ) + signal_mat = cu(repeat(signal, outer = (1,3))) + + track_state = @inferred track( + signal_mat, + track_state, + sampling_frequency + ) + + iterations = 2000 + code_phases = zeros(iterations) + carrier_phases = zeros(iterations) + tracked_code_phases = zeros(iterations) + tracked_carrier_phases = zeros(iterations) + tracked_code_dopplers = zeros(iterations) + tracked_carrier_dopplers = zeros(iterations) + tracked_prompts = zeros(ComplexF64, iterations) + for i = 1:iterations + carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency + + start_carrier_phase + π) - π + code_phase = mod( + code_frequency * 4000 * i / sampling_frequency + start_code_phase, + 1023 + ) + signal = cis.( + 2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase + ) .* + get_code.( + gpsl1, + code_frequency .* range ./ sampling_frequency .+ code_phase, + prn + ) + signal_mat = cu(repeat(signal, outer = (1,3))) + track_state = @inferred track( + signal_mat, + track_state, + sampling_frequency + ) + comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) / + sampling_frequency + start_carrier_phase + π) - π + comp_code_phase = mod( + code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, + 1023 + ) + tracked_code_phases[i] = get_code_phase(track_state.system_sats_states[1].states[1]) + tracked_carrier_phases[i] = get_carrier_phase(track_state.system_sats_states[1].states[1]) + tracked_carrier_dopplers[i] = get_carrier_doppler(track_state.system_sats_states[1].states[1])/Hz + tracked_code_dopplers[i] = get_code_doppler(track_state.system_sats_states[1].states[1])/Hz + tracked_prompts[i] = get_last_fully_integrated_filtered_prompt(track_state.system_sats_states[1].states[1]) + code_phases[i] = comp_code_phase + carrier_phases[i] = comp_carrier_phase + end + @test tracked_code_phases[end] ≈ code_phases[end] atol = 1e-2 + @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-5 + +# using PyPlot +# pygui(true) +# figure("carrier_phases") +# plot(tracked_carrier_phases) +# plot(carrier_phases) +# grid(true) +# figure("Code Phases") +# plot(300 * (tracked_code_phases .- code_phases)) +# figure("Carrier Doppler") +# plot(tracked_carrier_dopplers) +# figure("Code Doppler") +# plot(tracked_code_dopplers) +# figure("Prompt") +# plot(real.(tracked_prompts)) +# plot(imag.(tracked_prompts)) +end \ No newline at end of file diff --git a/test/tracking_loop.jl b/test/tracking_loop.jl deleted file mode 100644 index 6c7a117..0000000 --- a/test/tracking_loop.jl +++ /dev/null @@ -1,393 +0,0 @@ -@testset "Resize downconverted signal ($type) for multiple ants" for type in (Int16, Int32, Int64, Float32, Float64) - downconverted_signal_temp = Tracking.DownconvertedSignalCPU(NumAnts(4)) - signal = StructArray(ones(Complex{type}, 2500, 4)) - downconverted_signal = Tracking.resize!(downconverted_signal_temp, 2500, signal) - if type == Float64 - @test size(downconverted_signal.downconverted_signal_f32) == (0, 4) - @test size(downconverted_signal.downconverted_signal_f64) == (2500, 4) - else - @test size(downconverted_signal.downconverted_signal_f32) == (2500, 4) - @test size(downconverted_signal.downconverted_signal_f64) == (0, 4) - end -end - -@testset "Resize downconverted signal ($type) for single ant" for type in (Int16, Int32, Int64, Float32, Float64) - downconverted_signal_temp = Tracking.DownconvertedSignalCPU(NumAnts(1)) - signal = StructArray(ones(Complex{type}, 2500)) - downconverted_signal = Tracking.resize!(downconverted_signal_temp, 2500, signal) - if type == Float64 - @test size(downconverted_signal.downconverted_signal_f32) == (0,) - @test size(downconverted_signal.downconverted_signal_f64) == (2500,) - else - @test size(downconverted_signal.downconverted_signal_f32) == (2500,) - @test size(downconverted_signal.downconverted_signal_f64) == (0,) - - end -end - -@testset "Integration time" begin - gpsl1 = GPSL1() - galileo_e1b = GalileoE1B() - max_integration_time = 2ms - bit_found = false - time = @inferred Tracking.get_integration_time(gpsl1, max_integration_time, bit_found) - @test time == 1ms - - max_integration_time = 2ms - bit_found = true - time = @inferred Tracking.get_integration_time(gpsl1, max_integration_time, bit_found) - @test time == 2ms - - max_integration_time = 2ms - bit_found = false - time = @inferred Tracking.get_integration_time( - galileo_e1b, - max_integration_time, - bit_found - ) - @test time == 2ms - - max_integration_time = 10ms - bit_found = false - time = @inferred Tracking.get_integration_time( - galileo_e1b, - max_integration_time, - bit_found - ) - @test time == 4ms -end - -@testset "Number of chips to integrate" begin - gpsl1 = GPSL1() - max_integration_time = 1ms - code_phase = 200 - bit_found = true - chips = @inferred Tracking.get_num_chips_to_integrate( - gpsl1, - max_integration_time, - code_phase, - bit_found - ) - @test chips == 1023 - 200 - - code_phase = 1200 - chips = @inferred Tracking.get_num_chips_to_integrate( - gpsl1, - max_integration_time, - code_phase, - bit_found - ) - @test chips == 2046 - 1200 -end - -@testset "Number of samples to integrate" begin - gpsl1 = GPSL1() - max_integration_time = 1ms - sampling_frequency = 4e6Hz - current_code_doppler = 0Hz - current_code_phase = 0 - bit_found = true - samples = @inferred Tracking.get_num_samples_left_to_integrate( - gpsl1, - max_integration_time, - sampling_frequency, - current_code_doppler, - current_code_phase, - bit_found - ) - @test samples == 4000 -end - -@testset "Carrier phase delta" begin - intermediate_frequency = 100Hz - carrier_doppler = 50Hz - frequency = @inferred Tracking.get_current_carrier_frequency( - intermediate_frequency, - carrier_doppler - ) - @test frequency == intermediate_frequency + carrier_doppler -end - -@testset "Code phase delta" begin - gpsl1 = GPSL1() - code_doppler = 1Hz - frequency = @inferred Tracking.get_current_code_frequency(gpsl1, code_doppler) - @test frequency == 1023e3Hz + 1Hz -end - -@testset "Doppler aiding" begin - gpsl1 = GPSL1() - init_carrier_doppler = 10Hz - init_code_doppler = 1Hz - carrier_freq_update = 2Hz - code_freq_update = -0.5Hz - velocity_aiding = 3Hz - - carrier_freq, code_freq = @inferred Tracking.aid_dopplers( - gpsl1, - init_carrier_doppler, - init_code_doppler, - carrier_freq_update, - code_freq_update, - velocity_aiding - ) - - @test carrier_freq == 10Hz + 2Hz + 3Hz - @test code_freq == 1Hz + (2Hz + 3Hz) / 1540 - 0.5Hz -end - -@testset "Tracking with signal of type $type" for type in (Int16, Int32, Int64, Float32, Float64) - gpsl1 = GPSL1() - carrier_doppler = 200Hz - start_code_phase = 100 - code_frequency = carrier_doppler / 1540 + 1023e3Hz - sampling_frequency = 4e6Hz - prn = 1 - range = 0:3999 - start_carrier_phase = π / 2 - state = TrackingState(prn, gpsl1, carrier_doppler - 20Hz, start_code_phase) - - signal_temp = cis.( - 2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase - ) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ start_code_phase, - prn - ) - scaling = 512 - signal = type <: Integer ? complex.(round.(type, real.(signal_temp) * scaling), round.(type, imag.(signal_temp) * scaling)) : Complex{type}.(signal_temp) - - track_result = @inferred track(signal, state, sampling_frequency) - - iterations = 2000 - code_phases = zeros(iterations) - carrier_phases = zeros(iterations) - tracked_code_phases = zeros(iterations) - tracked_carrier_phases = zeros(iterations) - tracked_code_dopplers = zeros(iterations) - tracked_carrier_dopplers = zeros(iterations) - tracked_prompts = zeros(ComplexF64, iterations) - for i = 1:iterations - carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency + - start_carrier_phase + π) - π - code_phase = mod( - code_frequency * 4000 * i / sampling_frequency + start_code_phase, - 1023 - ) - signal = cis.( - 2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase - ) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ code_phase, - prn - ) - track_result = @inferred track( - signal, - get_state(track_result), - sampling_frequency - ) - comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) / - sampling_frequency + start_carrier_phase + π) - π - comp_code_phase = mod( - code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, - 1023 - ) - tracked_code_phases[i] = get_code_phase(track_result) - tracked_carrier_phases[i] = get_carrier_phase(track_result) - tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz - tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz - tracked_prompts[i] = get_prompt(track_result) - code_phases[i] = comp_code_phase - carrier_phases[i] = comp_carrier_phase - end - @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4 - @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-2 - -# using PyPlot -# pygui(true) -# figure("carrier_phases") -# plot(tracked_carrier_phases) -# plot(carrier_phases) -# grid(true) -# figure("Code Phases") -# plot(300 * (tracked_code_phases .- code_phases)) -# figure("Carrier Doppler") -# plot(tracked_carrier_dopplers) -# figure("Code Doppler") -# plot(tracked_code_dopplers) -# figure("Prompt") -# plot(real.(tracked_prompts)) -# plot(imag.(tracked_prompts)) -end - -@testset "Tracking with intermediate frequency of $intermediate_frequency" for intermediate_frequency in (0.0Hz, -10000.0Hz, 10000.0Hz, -30000.0Hz, 30000.0Hz) - gpsl1 = GPSL1() - carrier_doppler = 200Hz - start_code_phase = 100 - code_frequency = carrier_doppler / 1540 + 1023e3Hz - sampling_frequency = 4e6Hz - prn = 1 - range = 0:3999 - start_carrier_phase = π / 2 - state = TrackingState(prn, gpsl1, carrier_doppler - 20Hz, start_code_phase) - - signal = cis.( - 2π .* (carrier_doppler + intermediate_frequency) .* range ./ sampling_frequency .+ start_carrier_phase - ) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ start_code_phase, - prn - ) - - track_result = @inferred track(signal, state, sampling_frequency; intermediate_frequency) - - iterations = 2000 - code_phases = zeros(iterations) - carrier_phases = zeros(iterations) - tracked_code_phases = zeros(iterations) - tracked_carrier_phases = zeros(iterations) - tracked_code_dopplers = zeros(iterations) - tracked_carrier_dopplers = zeros(iterations) - tracked_prompts = zeros(ComplexF64, iterations) - for i = 1:iterations - carrier_phase = mod2pi(2π * (carrier_doppler + intermediate_frequency) * 4000 * i / sampling_frequency + - start_carrier_phase + π) - π - code_phase = mod( - code_frequency * 4000 * i / sampling_frequency + start_code_phase, - 1023 - ) - signal = cis.( - 2π .* (carrier_doppler + intermediate_frequency) .* range ./ sampling_frequency .+ carrier_phase - ) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ code_phase, - prn - ) - track_result = @inferred track( - signal, - get_state(track_result), - sampling_frequency; - intermediate_frequency - ) - comp_carrier_phase = mod2pi(2π * (carrier_doppler + intermediate_frequency) * 4000 * (i + 1) / - sampling_frequency + start_carrier_phase + π) - π - comp_code_phase = mod( - code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, - 1023 - ) - tracked_code_phases[i] = get_code_phase(track_result) - tracked_carrier_phases[i] = get_carrier_phase(track_result) - tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz - tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz - tracked_prompts[i] = get_prompt(track_result) - code_phases[i] = comp_code_phase - carrier_phases[i] = comp_carrier_phase - end - @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4 - @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-2 - -# using PyPlot -# pygui(true) -# figure("carrier_phases") -# plot(tracked_carrier_phases) -# plot(carrier_phases) -# grid(true) -# figure("Code Phases") -# plot(300 * (tracked_code_phases .- code_phases)) -# figure("Carrier Doppler") -# plot(tracked_carrier_dopplers) -# figure("Code Doppler") -# plot(tracked_code_dopplers) -# figure("Prompt") -# plot(real.(tracked_prompts)) -# plot(imag.(tracked_prompts)) -end - -@testset "Track multiple signals with signal of type $type" for type in (Int16, Int32, Int64, Float32, Float64) - gpsl1 = GPSL1() - carrier_doppler = 200Hz - start_code_phase = 100 - code_frequency = carrier_doppler / 1540 + 1023e3Hz - sampling_frequency = 4e6Hz - prn = 1 - range = 0:3999 - start_carrier_phase = π / 2 - state = TrackingState(prn, gpsl1, carrier_doppler - 20Hz, start_code_phase, num_ants = NumAnts(3)) - - signal = cis.( - 2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase - ) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ start_code_phase, - prn - ) - signal_mat_temp = repeat(signal, outer = (1,3)) - scaling = 512 - signal_mat = type <: Integer ? complex.(floor.(type, real.(signal_mat_temp) * scaling), floor.(type, imag.(signal_mat_temp) * scaling)) : Complex{type}.(signal_mat_temp) - - @test_throws ArgumentError track(signal_mat', state, sampling_frequency) - - track_result = @inferred track(signal_mat, state, sampling_frequency) - - iterations = 2000 - code_phases = zeros(iterations) - carrier_phases = zeros(iterations) - tracked_code_phases = zeros(iterations) - tracked_carrier_phases = zeros(iterations) - tracked_code_dopplers = zeros(iterations) - tracked_carrier_dopplers = zeros(iterations) - for i = 1:iterations - carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency + - start_carrier_phase + π) - π - code_phase = mod( - code_frequency * 4000 * i / sampling_frequency + start_code_phase, - 1023 - ) - signal = cis.( - 2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase - ) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ code_phase, - prn - ) - signal_mat = repeat(signal, outer = (1,3)) - track_result = @inferred track( - signal_mat, - get_state(track_result), - sampling_frequency - ) - comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) / - sampling_frequency + start_carrier_phase + π) - π - comp_code_phase = mod( - code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, - 1023 - ) - tracked_code_phases[i] = get_code_phase(track_result) - tracked_carrier_phases[i] = get_carrier_phase(track_result) - tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz - tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz - code_phases[i] = comp_code_phase - carrier_phases[i] = comp_carrier_phase - end - @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4 - @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-2 - -# using PyPlot -# pygui(true) -# figure("carrier_phases") -# plot(tracked_carrier_phases) -# plot(carrier_phases) -# grid(true) -# figure("Code Phases") -# plot(300 * (tracked_code_phases .- code_phases)) -# figure("Carrier Doppler") -# plot(tracked_carrier_dopplers) -# figure("Code Doppler") -# plot(tracked_code_dopplers) -end diff --git a/test/tracking_results.jl b/test/tracking_results.jl deleted file mode 100644 index 21bc401..0000000 --- a/test/tracking_results.jl +++ /dev/null @@ -1,52 +0,0 @@ -@testset "Tracking results" begin - gpsl1 = GPSL1() - results = Tracking.TrackingResults( - TrackingState(1, gpsl1, 100Hz, 100), - EarlyPromptLateCorrelator(NumAnts(2)), - complex(1.2,1.3), - SVector(-1, 0, 1), - 1, - 1.0Hz, - 1.0, - true, - Tracking.BitBuffer(), - 45dBHz - ) - - @test @inferred(get_carrier_doppler(results)) == 100Hz - @test @inferred(get_carrier_phase(results)) == 0.0 - @test @inferred(get_code_doppler(results)) == 100Hz / 1540 - @test @inferred(get_code_phase(results)) == 100 - @test @inferred(get_correlator_carrier_phase(results)) ≈ 1.0 * 2π - @test @inferred(get_correlator_carrier_frequency(results)) == 1.0Hz - correlator = @inferred get_correlator(results) - @test correlator == EarlyPromptLateCorrelator(NumAnts(2)) - @test @inferred(get_correlator_sample_shifts(results)) == SVector(-1, 0, 1) - @test @inferred(get_early_late_index_shift(results)) == 1 - @test @inferred(get_early(results)) == [0.0, 0.0] - @test @inferred(get_prompt(results)) == [0.0, 0.0] - @test @inferred(get_late(results)) == [0.0, 0.0] - @test @inferred(get_filtered_prompt(results)) == complex(1.2,1.3) - @test @inferred(get_bits(results)) == 0 - @test @inferred(get_num_bits(results)) == 0 - @test @inferred(get_secondary_code_or_bit_found(results)) == false - @test @inferred(get_cn0(results)) == 45dBHz - - results = Tracking.TrackingResults( - TrackingState(1, gpsl1, 100Hz, 100), - EarlyPromptLateCorrelator(NumAnts(2)), - complex(1.2,1.3), - SVector(-1, 0, 1), - 1, - 1.0Hz, - 1.0, - false, - Tracking.BitBuffer(), - 45dBHz - ) - correlator = @inferred get_correlator(results) - @test @inferred(get_filtered_prompt(results)) == complex(1.2,1.3) - @test @inferred(get_early(results)) == [0, 0] - @test @inferred(get_prompt(results)) == [0, 0] - @test @inferred(get_late(results)) == [0, 0] -end diff --git a/test/tracking_state.jl b/test/tracking_state.jl index e6124a7..b7ee10a 100644 --- a/test/tracking_state.jl +++ b/test/tracking_state.jl @@ -1,73 +1,230 @@ @testset "Tracking state" begin + sampling_frequency = 5e6Hz + num_samples = 5000 + gpsl1 = GPSL1() + sat_states = [ + SatState(gpsl1, 1, sampling_frequency, 10.5, 10.0Hz), + SatState(gpsl1, 2, sampling_frequency, 11.5, 20.0Hz) + ] + + track_state = @inferred TrackState( + gpsl1, + sat_states; + num_samples, + ) + + @test get_system(track_state) isa GPSL1 + @test get_sat_state(track_state, 1).prn == 1 + @test get_sat_state(track_state, 2).prn == 2 + + @test track_state.doppler_estimator.plls_and_dlls[1][1].init_carrier_doppler == 10.0Hz + @test track_state.doppler_estimator.plls_and_dlls[1][2].init_carrier_doppler == 20.0Hz + + @test track_state.downconvert_and_correlator isa CPUDownconvertAndCorrelator + + system_sats_states = (SystemSatsState( + gpsl1, + [ + SatState(gpsl1, 1, sampling_frequency, 10.5, 10.0Hz), + SatState(gpsl1, 2, sampling_frequency, 11.5, 20.0Hz) + ] + ),) - carrier_doppler = 100Hz - code_phase = 100 + track_state = @inferred TrackState( + system_sats_states; + num_samples, + ) + + @test @inferred(get_system(track_state)) isa GPSL1 + @test @inferred(get_sat_state(track_state, 1)).prn == 1 + @test @inferred(get_sat_state(track_state, 2)).prn == 2 + + @test track_state.doppler_estimator.plls_and_dlls[1][1].init_carrier_doppler == 10.0Hz + @test track_state.doppler_estimator.plls_and_dlls[1][2].init_carrier_doppler == 20.0Hz + + @test track_state.downconvert_and_correlator isa CPUDownconvertAndCorrelator + + sat_states_num_ants2 = [ + SatState(gpsl1, 1, sampling_frequency, 10.5, 10.0Hz, num_ants = NumAnts(2)), + SatState(gpsl1, 2, sampling_frequency, 11.5, 20.0Hz, num_ants = NumAnts(2)) + ] + + track_state = @inferred TrackState( + gpsl1, + sat_states_num_ants2; + num_samples + ) + + @test track_state.downconvert_and_correlator isa CPUDownconvertAndCorrelator + @test track_state.downconvert_and_correlator.buffers[1][1].downconvert_signal_buffer isa AbstractMatrix + @test size(track_state.downconvert_and_correlator.buffers[1][1].downconvert_signal_buffer, 2) == 2 + @test size(track_state.downconvert_and_correlator.buffers[1][1].downconvert_signal_buffer, 1) == 5000 + @test track_state.downconvert_and_correlator.buffers[1][2].downconvert_signal_buffer isa AbstractMatrix + @test size(track_state.downconvert_and_correlator.buffers[1][2].downconvert_signal_buffer, 2) == 2 + @test size(track_state.downconvert_and_correlator.buffers[1][2].downconvert_signal_buffer, 1) == 5000 +end + +@testset "Add and remove satellite state to and from track state" begin + sampling_frequency = 5e6Hz + num_samples = 5000 gpsl1 = GPSL1() - state = TrackingState(1, gpsl1, carrier_doppler, code_phase) - - @test @inferred(Tracking.get_prn(state)) == 1 - @test @inferred(Tracking.get_code_phase(state)) == 100 - @test @inferred(Tracking.get_carrier_phase(state)) == 0.0 - @test @inferred(Tracking.get_init_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_init_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_correlator(state)) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_sc_bit_detector(state)) == SecondaryCodeOrBitDetector() - @test @inferred(Tracking.get_carrier_loop_filter(state)) == ThirdOrderBilinearLF() - @test @inferred(Tracking.get_code_loop_filter(state)) == SecondOrderBilinearLF() - @test @inferred(Tracking.get_prompt_accumulator(state)) == 0.0 - @test @inferred(Tracking.get_integrated_samples(state)) == 0 - @test @inferred(get_post_corr_filter(state)) == Tracking.DefaultPostCorrFilter() - - state = TrackingState( + sat_states = [ + SatState(gpsl1, 1, sampling_frequency, 10.5, 10.0Hz), + SatState(gpsl1, 2, sampling_frequency, 11.5, 20.0Hz) + ] + + track_state = @inferred TrackState( + gpsl1, + sat_states; + num_samples, + ) + + @inferred add_sats!( + track_state, 1, gpsl1, - carrier_doppler, - code_phase; - code_doppler = carrier_doppler * get_code_center_frequency_ratio(gpsl1), - carrier_phase = 0.0, - carrier_loop_filter = ThirdOrderBilinearLF(), - code_loop_filter = SecondOrderBilinearLF(), - sc_bit_detector = SecondaryCodeOrBitDetector(), - correlator = Tracking.get_default_correlator(gpsl1, NumAnts(1)), - integrated_samples = 0, - prompt_accumulator = zero(ComplexF64) + SatState(gpsl1, 3, sampling_frequency, 5.5, 80.0Hz), + ) + + @test length(@inferred(get_sat_states(track_state))) == 3 + @test @inferred(get_sat_state(track_state, 3)).prn == 3 + @test @inferred(get_sat_state(track_state, 1, 3)).prn == 3 + + @inferred remove_sats!( + track_state, + 1, + 2, + ) + @test length(@inferred(get_sat_states(track_state))) == 2 + @test @inferred(get_sat_state(track_state, 1)).prn == 1 + @test @inferred(get_sat_state(track_state, 2)).prn == 3 + + @inferred add_sats!( + track_state, + gpsl1, + [ + SatState(gpsl1, 6, sampling_frequency, 12.5, 60.0Hz), + SatState(gpsl1, 8, sampling_frequency, 67.5, 120.0Hz), + ] + ) + + @test length(@inferred(get_sat_states(track_state))) == 4 + @test @inferred(get_sat_state(track_state, 1)).prn == 1 + @test @inferred(get_sat_state(track_state, 2)).prn == 3 + @test @inferred(get_sat_state(track_state, 3)).prn == 6 + @test @inferred(get_sat_state(track_state, 4)).prn == 8 + + @inferred remove_sats!( + track_state, + [1, 6], ) - @test @inferred(Tracking.get_prn(state)) == 1 - @test @inferred(Tracking.get_code_phase(state)) == 100 - @test @inferred(Tracking.get_carrier_phase(state)) == 0.0 - @test @inferred(Tracking.get_init_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_init_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_correlator(state)) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_sc_bit_detector(state)) == SecondaryCodeOrBitDetector() - @test @inferred(Tracking.get_carrier_loop_filter(state)) == ThirdOrderBilinearLF() - @test @inferred(Tracking.get_code_loop_filter(state)) == SecondOrderBilinearLF() - @test @inferred(Tracking.get_prompt_accumulator(state)) == 0.0 - @test @inferred(Tracking.get_integrated_samples(state)) == 0 - - state = TrackingState(1, gpsl1, carrier_doppler, code_phase, num_ants = NumAnts(2)) - @test @inferred(Tracking.get_correlator(state)) == EarlyPromptLateCorrelator(NumAnts(2)) - - acq = AcquisitionResults(GPSL1(), 5, 5e6Hz, 100.0Hz, 524.6, 45.0, 1.0, randn(100,100), -500:100.0:500) - state = TrackingState(acq) - @test @inferred(Tracking.get_prn(state)) == 5 - @test @inferred(Tracking.get_code_phase(state)) == 524.6 - @test @inferred(Tracking.get_carrier_phase(state)) == 0.0 - @test @inferred(Tracking.get_init_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_init_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_code_doppler(state)) == 100Hz / 1540 - @test @inferred(Tracking.get_carrier_doppler(state)) == 100Hz - @test @inferred(Tracking.get_correlator(state)) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_sc_bit_detector(state)) == SecondaryCodeOrBitDetector() - @test @inferred(Tracking.get_carrier_loop_filter(state)) == ThirdOrderBilinearLF() - @test @inferred(Tracking.get_code_loop_filter(state)) == SecondOrderBilinearLF() - @test @inferred(Tracking.get_prompt_accumulator(state)) == 0.0 - @test @inferred(Tracking.get_integrated_samples(state)) == 0 - - changed_state = @inferred(TrackingState(state, post_corr_filter = Tracking.DefaultPostCorrFilter())) - @test get_post_corr_filter(changed_state) == Tracking.DefaultPostCorrFilter() + @test length(@inferred(get_sat_states(track_state))) == 2 + @test @inferred(get_sat_state(track_state, 1)).prn == 3 + @test @inferred(get_sat_state(track_state, 2)).prn == 8 + end + +@testset "Add and remove satellite state to track state with multiple systems" begin + sampling_frequency = 5e6Hz + num_samples = 5000 + gpsl1 = GPSL1() + galileo_e1b = GalileoE1B() + + system_sats_states = (gps = SystemSatsState( + gpsl1, + [ + SatState(gpsl1, 1, sampling_frequency, 10.5, 10.0Hz), + SatState(gpsl1, 2, sampling_frequency, 11.5, 20.0Hz) + ] + ), + gal = SystemSatsState( + galileo_e1b, + [ + SatState(galileo_e1b, 2, sampling_frequency, 10.5, 10.0Hz), + SatState(galileo_e1b, 3, sampling_frequency, 11.5, 20.0Hz) + ] + ),) + + track_state = @inferred TrackState( + system_sats_states; + num_samples, + ) + + @inferred add_sats!( + track_state, + 1, + gpsl1, + SatState(gpsl1, 3, sampling_frequency, 5.5, 80.0Hz), + ) + + @test length(@inferred(get_sat_states(track_state, Val(:gps)))) == 3 + @test length(@inferred(get_sat_states(track_state, Val(:gal)))) == 2 + @test length(@inferred(get_sat_states(track_state, Val(1)))) == 3 + @test length(@inferred(get_sat_states(track_state, Val(2)))) == 2 + @test @inferred(get_sat_state(track_state, Val(:gps), 3)).prn == 3 + @test @inferred(get_sat_state(track_state, Val(1), 3)).prn == 3 + + @inferred remove_sats!( + track_state, + 1, + 2, + ) + @test length(get_sat_states(track_state, :gps)) == 2 + @test length(get_sat_states(track_state, :gal)) == 2 + @test length(get_sat_states(track_state, 1)) == 2 + @test length(get_sat_states(track_state, 2)) == 2 + @test get_sat_state(track_state, :gps, 1).prn == 1 + @test get_sat_state(track_state, :gps, 2).prn == 3 + + @inferred add_sats!( + track_state, + :gps, + gpsl1, + SatState(gpsl1, 2, sampling_frequency, 5.5, 80.0Hz), + ) + @test length(get_sat_states(track_state, :gps)) == 3 + @test length(get_sat_states(track_state, :gal)) == 2 + @test get_sat_state(track_state, :gps, 1).prn == 1 + @test get_sat_state(track_state, :gps, 2).prn == 3 + @test get_sat_state(track_state, :gps, 3).prn == 2 + + @inferred remove_sats!( + track_state, + :gps, + 2, + ) + @test length(get_sat_states(track_state, :gps)) == 2 + @test length(get_sat_states(track_state, :gal)) == 2 + @test get_sat_state(track_state, :gps, 1).prn == 1 + @test get_sat_state(track_state, :gps, 2).prn == 3 + + @inferred add_sats!( + track_state, + :gps, + gpsl1, + [ + SatState(gpsl1, 6, sampling_frequency, 12.5, 60.0Hz), + SatState(gpsl1, 8, sampling_frequency, 67.5, 120.0Hz), + ] + ) + + @test length(get_sat_states(track_state, :gps)) == 4 + @test length(get_sat_states(track_state, :gal)) == 2 + @test get_sat_state(track_state, :gps, 1).prn == 1 + @test get_sat_state(track_state, :gps, 2).prn == 3 + @test get_sat_state(track_state, :gps, 3).prn == 6 + @test get_sat_state(track_state, :gps, 4).prn == 8 + + @inferred remove_sats!( + track_state, + :gps, + [1, 6], + ) + + @test length(get_sat_states(track_state, :gps)) == 2 + @test length(get_sat_states(track_state, :gal)) == 2 + @test get_sat_state(track_state, :gps, 1).prn == 3 + @test get_sat_state(track_state, :gps, 2).prn == 8 + +end \ No newline at end of file