From 94bcda93f0bb4ac354505d68003adb2353d65910 Mon Sep 17 00:00:00 2001 From: zsoerenm Date: Wed, 3 Feb 2021 10:51:24 +0100 Subject: [PATCH] From fixedpoint to float and generic correlator (#17) * Implemented correlator for custom number of taps * Added required functions for GenericBOC signals * Bug fixes for generic correlator * Fix functions for GenericBOC signals * Fix wrong type inheritance * [GenericCorrelator] Fix multi antenna correlation * Fix integration time calculation for short codes * Fix order of early and late * Switch to full declaration of correlator shifts * Fixed integer based replica generation * Adapt generic correlator to new spacing definition * Performance improvements * Fix typo * Use SVector for GenericCorrelator * Performance improvements * Re-Added defaults for GenericBOC signals * Replace SVector{SVector} by Vector{SVector} * Adjust to changes in GNSSSignals * Change carrier fixed point calculation to float * Renamed `generic_boc.jl` to `boc.jl` * Fix infinite recursion error * Added tests for BOCcos Co-authored-by: Michael Niestroj Co-authored-by: Soeren Schoenbrod --- Project.toml | 8 +- src/Tracking.jl | 18 +- src/agc.jl | 67 ----- src/bit_buffer.jl | 6 +- src/boc.jl | 19 ++ src/carrier_replica.jl | 33 +-- src/code_replica.jl | 52 ++-- src/correlator.jl | 78 +++--- src/correlators/generic_correlator.jl | 361 +++++++++++++++++++++++++ src/discriminators.jl | 10 +- src/downconvert.jl | 69 ++--- src/galileo_e1b.jl | 4 +- src/gpsl1.jl | 4 +- src/gpsl5.jl | 4 +- src/secondary_code_or_bit_detector.jl | 4 +- src/tracking_loop.jl | 257 ++++++++++-------- src/tracking_state.jl | 97 +++++-- test/agc.jl | 43 --- test/bit_buffer.jl | 12 +- test/boc.jl | 18 ++ test/carrier_replica.jl | 17 +- test/cn0_estimation.jl | 22 +- test/code_replica.jl | 11 +- test/correlator.jl | 41 ++- test/discriminators.jl | 22 +- test/galileo_e1b.jl | 15 +- test/generic_correlator.jl | 184 +++++++++++++ test/gps_l1.jl | 11 +- test/gps_l5.jl | 11 +- test/runtests.jl | 2 +- test/secondary_code_or_bit_detector.jl | 11 +- test/tracking_loop.jl | 56 ++-- test/tracking_results.jl | 5 +- test/tracking_state.jl | 12 +- 34 files changed, 1060 insertions(+), 524 deletions(-) delete mode 100644 src/agc.jl create mode 100644 src/boc.jl create mode 100644 src/correlators/generic_correlator.jl delete mode 100644 test/agc.jl create mode 100644 test/boc.jl create mode 100644 test/generic_correlator.jl diff --git a/Project.toml b/Project.toml index bf62a8d..46b54f9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Tracking" uuid = "10b2438b-ffd4-5096-aa58-44041d5c8f3b" authors = ["Soeren Zorn "] -version = "0.13.0" +version = "0.14.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -15,9 +15,9 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] DocStringExtensions = "0.6, 0.7, 0.8" -GNSSSignals = "0.13" -LoopVectorization = "0.8" -StaticArrays = "0.9, 0.10, 0.11, 0.12" +GNSSSignals = "0.14" +LoopVectorization = "0.8, 0.9, 0.10, 0.11" +StaticArrays = "0.9, 0.10, 0.11, 0.12, 1.0" StructArrays = "0.4" TrackingLoopFilters = "0.1" Unitful = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 1.0" diff --git a/src/Tracking.jl b/src/Tracking.jl index eec87da..e542d40 100644 --- a/src/Tracking.jl +++ b/src/Tracking.jl @@ -13,12 +13,19 @@ module Tracking get_early, get_prompt, get_late, + get_tap, + get_taps, + get_num_taps, + get_early_index, + get_prompt_index, + get_late_index, get_correlator, get_carrier_doppler, get_carrier_phase, get_code_doppler, get_code_phase, - get_early_late_sample_shift, + get_correlator_sample_shifts, + get_early_late_sample_spacing, get_secondary_code_or_bit_found, get_correlator_carrier_phase, get_correlator_carrier_frequency, @@ -27,12 +34,14 @@ module Tracking track, TrackingState, NumAnts, + NumTaps, MomentsCN0Estimator, AbstractCN0Estimator, get_bits, get_num_bits, EarlyPromptLateCorrelator, VeryEarlyPromptLateCorrelator, + GenericCorrelator, SecondaryCodeOrBitDetector, GainControlledSignal @@ -41,7 +50,11 @@ module Tracking NumAnts(x) = NumAnts{x}() - include("agc.jl") + struct NumTaps{x} + end + + NumTaps(x) = NumTaps{x}() + include("code_replica.jl") include("carrier_replica.jl") include("downconvert.jl") @@ -56,4 +69,5 @@ module Tracking include("gpsl1.jl") include("gpsl5.jl") include("galileo_e1b.jl") + include("boc.jl") end diff --git a/src/agc.jl b/src/agc.jl deleted file mode 100644 index b1fdd9b..0000000 --- a/src/agc.jl +++ /dev/null @@ -1,67 +0,0 @@ -struct GainControlledSignal{ - S <: StructArray{Complex{Int16}}, - A <: Union{Real, Vector{<:Real}} -} - signal::S - attenuation::A - amplitude_power::Int -end - -get_signal(agc::GainControlledSignal) = agc.signal -get_attenuation(agc::GainControlledSignal) = agc.attenuation -get_amplitude_power(agc::GainControlledSignal) = agc.amplitude_power - -@inline function GainControlledSignal!( - agc_signal::StructArray{Complex{Int16}}, - signal::AbstractVector, - bits::Integer = 7 -) - size(agc_signal) == size(signal) || - throw(DimensionMismatch("size of AGC signal not equal to size of signal")) - max_ampl = find_max(signal) - amplification = 1 << bits / max_ampl - @inbounds @fastmath for i = eachindex(signal, agc_signal) - agc_signal.re[i] = floor(Int16, real(signal[i]) * amplification) - agc_signal.im[i] = floor(Int16, imag(signal[i]) * amplification) - end - GainControlledSignal(agc_signal, max_ampl, bits) -end - -@inline function GainControlledSignal!( - agc_signal::StructArray{Complex{Int16}}, - signal::AbstractMatrix, - bits::Integer = 7 -) - size(agc_signal) == size(signal) || - throw(DimensionMismatch("size of AGC signal not equal to size of signal")) - max_ampl = map(find_max, eachcol(signal)) - amplification = (1 << bits) ./ max_ampl - @inbounds @fastmath for a = axes(signal, 2), i = axes(signal, 1) - agc_signal.re[i, a] = floor(Int16, real(signal[i, a]) * amplification[a]) - agc_signal.im[i, a] = floor(Int16, imag(signal[i, a]) * amplification[a]) - end - GainControlledSignal(agc_signal, max_ampl, bits) -end - -@inline function GainControlledSignal(signal, bits::Integer = 7) - GainControlledSignal!( - StructArray{Complex{Int16}}(undef, size(signal)), - signal, - bits - ) -end - -@inline function find_max(signal) - max_real_value = 0.0 - max_imag_value = 0.0 - for i = 1:length(signal) - if real(signal[i]) > max_real_value - max_real_value = real(signal[i]) - end - if imag(signal[i]) > max_imag_value - max_imag_value = imag(signal[i]) - end - end - #sqrt(float(max_real_value)^2 + float(max_imag_value)^2) - max(max_real_value, max_imag_value) -end diff --git a/src/bit_buffer.jl b/src/bit_buffer.jl index f821cdb..73ce538 100644 --- a/src/bit_buffer.jl +++ b/src/bit_buffer.jl @@ -21,7 +21,7 @@ $(SIGNATURES) Buffer data bits based on the prompt accumulation and the current prompt value. """ function buffer( - ::Type{S}, + system::AbstractGNSS, bit_buffer, prompt_accumulator, secondary_code_or_bit_found, @@ -29,12 +29,12 @@ function buffer( code_phase, integration_time, prompt_correlator -) where S <: AbstractGNSSSystem +) prompt_accumulator = prompt_accumulator + secondary_code_or_bit_found * prompt_correlator if secondary_code_or_bit_found && - (code_phase - prev_code_phase < 0 || integration_time == 1 / get_data_frequency(S)) + (code_phase - prev_code_phase < 0 || integration_time == 1 / get_data_frequency(system)) bit = real(prompt_accumulator) > 0 bit_buffer = BitBuffer( get_bits(bit_buffer) << 1 + UInt64(bit), diff --git a/src/boc.jl b/src/boc.jl new file mode 100644 index 0000000..bac7c35 --- /dev/null +++ b/src/boc.jl @@ -0,0 +1,19 @@ +""" +$(SIGNATURES) + +Checks if upcoming integration is a new bit for generic BOC signal +""" +function is_upcoming_integration_new_bit( + boc::BOCcos, + prns, + num_prns +) + is_upcoming_integration_new_bit(boc.system, prns, num_prns) +end + +function get_default_correlator( + ::BOCcos, + numAnts::NumAnts +) + EarlyPromptLateCorrelator(numAnts) +end diff --git a/src/carrier_replica.jl b/src/carrier_replica.jl index 3ff8a3d..1ccaac8 100644 --- a/src/carrier_replica.jl +++ b/src/carrier_replica.jl @@ -1,21 +1,16 @@ function gen_carrier_replica!( - carrier_replica::StructArray{Complex{Int16}}, + carrier_replica::StructArray{Complex{T}}, carrier_frequency, sampling_frequency, start_phase, - carrier_amplitude_power::Val{N}, start_sample, num_samples -) where N - fpcarrier!( - carrier_replica, - carrier_frequency, - sampling_frequency, - start_phase, - start_sample = start_sample, - num_samples = num_samples, - bits = carrier_amplitude_power - ) +) where T + c_re = carrier_replica.re; c_im = carrier_replica.im + @avx for i in 0:num_samples - 1 + c_im[i + start_sample], c_re[i + start_sample] = + sincos(T(2π) * (i * T(upreferred(carrier_frequency / Hz)) / T(upreferred(sampling_frequency / Hz)) + T(start_phase))) + end carrier_replica end @@ -29,17 +24,9 @@ function update_carrier_phase( carrier_frequency, sampling_frequency, start_phase, - carrier_amplitude_power::Val{N} -) where N - n = N + 2 - fixed_point = 32 - n - 2 - delta = floor(Int32, carrier_frequency * 1 << (fixed_point + n) / sampling_frequency) - fixed_point_start_phase = floor(Int32, start_phase * 1 << (fixed_point + n)) - phase_fixed_point = delta * num_samples + fixed_point_start_phase - mod( - phase_fixed_point / 1 << (fixed_point + n) + 0.5, - 1 - ) - 0.5 +) + phase = num_samples * carrier_frequency / sampling_frequency + start_phase + mod(phase + 0.5, 1) - 0.5 end """ diff --git a/src/code_replica.jl b/src/code_replica.jl index 4affee7..0726a64 100644 --- a/src/code_replica.jl +++ b/src/code_replica.jl @@ -1,29 +1,45 @@ +""" +$(SIGNATURES) + +Generate a code replica for a signal from satellite system `S`. The +replica contains `num_samples` prompt samples as well as an additional +number of early and late samples specified by `correlator_sample_shifts`. +The codefrequency is specified by `code_frequency`, while the sampling +rate is given by `sampling_frequency`. The phase of the first prompt +sample is given by `start_code_phase`. The generated signal is returned +in the array `code_replica` with the first generated sample written to +index start_sample. +""" function gen_code_replica!( code_replica, - ::Type{S}, + system::AbstractGNSS, code_frequency, sampling_frequency, start_code_phase::AbstractFloat, start_sample::Integer, num_samples::Integer, - early_late_sample_shift, + correlator_sample_shifts::SVector, prn::Integer -) where S <: AbstractGNSSSystem - fixed_point = sizeof(Int) * 8 - 1 - min_bits_for_code_length(S) +) + most_early_sample_shift = correlator_sample_shifts[end] + most_late_sample_shift = correlator_sample_shifts[1] + num_early_late_samples = most_early_sample_shift - most_late_sample_shift + fixed_point = sizeof(Int) * 8 - 1 - min_bits_for_code_length(system) delta = floor(Int, code_frequency * 1 << fixed_point / sampling_frequency) + ((most_late_sample_shift * delta) >> fixed_point < -get_code_length(system) || + (most_early_sample_shift * delta) >> fixed_point > get_code_length(system)) && + throw(ArgumentError("The number of taps or the tab spacing is too large.")) modded_start_code_phase = mod( start_code_phase, - get_code_length(S) * get_secondary_code_length(S) + get_code_length(system) * get_secondary_code_length(system) ) fixed_point_start_code_phase = floor(Int, modded_start_code_phase * 1 << fixed_point) - max_sample_shift = maximum(early_late_sample_shift) # Assumes, that the number of early shifts is identical to the number of late shifts - early_late_samples = 2 * max_sample_shift - @inbounds for i = start_sample:num_samples + early_late_samples + start_sample - 1 - fixed_point_code_phase = (i - max_sample_shift - start_sample) * delta + + @inbounds for i = start_sample:num_samples + num_early_late_samples + start_sample - 1 + fixed_point_code_phase = (i + most_late_sample_shift - start_sample) * delta + fixed_point_start_code_phase code_index = fixed_point_code_phase >> fixed_point - code_replica[i] = get_code_unsafe(S, code_index, prn) + code_replica[i] = get_code_unsafe(system, code_index, prn) end code_replica end @@ -34,20 +50,20 @@ $(SIGNATURES) Updates the code phase. """ function update_code_phase( - ::Type{S}, + system::AbstractGNSS, num_samples, code_frequency, sampling_frequency, start_code_phase, secondary_code_or_bit_found -) where S <: AbstractGNSSSystem - if get_data_frequency(S) == 0Hz - secondary_code_or_bit_length = get_secondary_code_length(S) +) + if get_data_frequency(system) == 0Hz + secondary_code_or_bit_length = get_secondary_code_length(system) else secondary_code_or_bit_length = - Int(get_code_frequency(S) / (get_data_frequency(S) * get_code_length(S))) + Int(get_code_frequency(system) / (get_data_frequency(system) * get_code_length(system))) end - code_length = get_code_length(S) * + code_length = get_code_length(system) * (secondary_code_or_bit_found ? secondary_code_or_bit_length : 1) mod(code_frequency * num_samples / sampling_frequency + start_code_phase, code_length) # fixed_point = sizeof(Int) * 8 - 1 - min_bits_for_code_length(S) @@ -62,6 +78,6 @@ $(SIGNATURES) Calculates the current code frequency. """ -function get_current_code_frequency(::Type{S}, code_doppler) where S <: AbstractGNSSSystem - code_doppler + get_code_frequency(S) +function get_current_code_frequency(system::AbstractGNSS, code_doppler) + code_doppler + get_code_frequency(system) end diff --git a/src/correlator.jl b/src/correlator.jl index e8cc53c..19983ba 100644 --- a/src/correlator.jl +++ b/src/correlator.jl @@ -97,15 +97,21 @@ end """ $(SIGNATURES) -Calculate the shift between the early and late in samples. +Calculate the replica phase offset required for the correlator taps with +respect to the prompt correlator, expressed in samples. The shifts are +ordered from latest to earliest replica. """ -function get_early_late_sample_shift( - ::Type{S}, +function get_correlator_sample_shifts( + system::AbstractGNSS, correlator::EarlyPromptLateCorrelator, sampling_frequency, preferred_code_shift -) where S <: AbstractGNSSSystem - round(Int, preferred_code_shift * sampling_frequency / get_code_frequency(S)) +) + @SVector [i*round(Int, preferred_code_shift * sampling_frequency / get_code_frequency(system)) for i in -1:1] +end + +function get_early_late_sample_spacing(correlator::EarlyPromptLateCorrelator, correlator_sample_shifts::SVector{3}) + correlator_sample_shifts[3] - correlator_sample_shifts[1] end """ @@ -128,31 +134,28 @@ Perform a correlation. """ function correlate( correlator::EarlyPromptLateCorrelator, - downconverted_signal, + downconverted_signal::AbstractVector, code, - early_late_sample_shift, + correlator_sample_shifts::SVector{3}, start_sample, - num_samples_left, - agc_attenuation, - agc_bits, - carrier_bits::Val{NC} -) where NC - late = zero(Complex{Int32}) - prompt = zero(Complex{Int32}) - early = zero(Complex{Int32}) - @inbounds for i = start_sample:num_samples_left + start_sample - 1 + num_samples +) + late = zero(eltype(downconverted_signal)) + prompt = zero(eltype(downconverted_signal)) + early = zero(eltype(downconverted_signal)) + @inbounds for i = start_sample:num_samples + start_sample - 1 late = late + downconverted_signal[i] * code[i] end - @inbounds for i = start_sample:num_samples_left + start_sample - 1 - prompt = prompt + downconverted_signal[i] * code[i + early_late_sample_shift] + @inbounds for i = start_sample:num_samples + start_sample - 1 + prompt = prompt + downconverted_signal[i] * code[i + correlator_sample_shifts[2]-correlator_sample_shifts[1]] end - @inbounds for i = start_sample:num_samples_left + start_sample - 1 - early = early + downconverted_signal[i] * code[i + 2 * early_late_sample_shift] + @inbounds for i = start_sample:num_samples + start_sample - 1 + early = early + downconverted_signal[i] * code[i + correlator_sample_shifts[3]-correlator_sample_shifts[1]] end EarlyPromptLateCorrelator( - get_early(correlator) + early * agc_attenuation / 1 << (agc_bits + NC), - get_prompt(correlator) + prompt * agc_attenuation / 1 << (agc_bits + NC), - get_late(correlator) + late * agc_attenuation / 1 << (agc_bits + NC) + get_early(correlator) + early, + get_prompt(correlator) + prompt, + get_late(correlator) + late ) end @@ -160,32 +163,30 @@ function correlate( correlator::EarlyPromptLateCorrelator{<: SVector{N}}, downconverted_signal::AbstractMatrix, code, - early_late_sample_shift, + correlator_sample_shifts::SVector{3}, start_sample, num_samples_left, - agc_attenuation, - agc_bits, - carrier_bits::Val{NC} -) where {N, NC} - late = zero(MVector{N, Complex{Int32}}) - prompt = zero(MVector{N, Complex{Int32}}) - early = zero(MVector{N, Complex{Int32}}) +) where N + late = zero(MVector{N, eltype(downconverted_signal)}) + prompt = zero(MVector{N, eltype(downconverted_signal)}) + early = zero(MVector{N, eltype(downconverted_signal)}) @inbounds for j = 1:length(late), i = start_sample:num_samples_left + start_sample - 1 late[j] = late[j] + downconverted_signal[i,j] * code[i] end @inbounds for j = 1:length(late), i = start_sample:num_samples_left + start_sample - 1 - prompt[j] = prompt[j] + downconverted_signal[i,j] * code[i + early_late_sample_shift] + prompt[j] = prompt[j] + downconverted_signal[i,j] * code[i + correlator_sample_shifts[2] - correlator_sample_shifts[1]] end @inbounds for j = 1:length(late), i = start_sample:num_samples_left + start_sample - 1 - early[j] = early[j] + downconverted_signal[i,j] * code[i + 2 * early_late_sample_shift] + early[j] = early[j] + downconverted_signal[i,j] * code[i + correlator_sample_shifts[3] - correlator_sample_shifts[1]] end EarlyPromptLateCorrelator( - get_early(correlator) + early .* agc_attenuation / 1 << (agc_bits + NC), - get_prompt(correlator) + prompt .* agc_attenuation / 1 << (agc_bits + NC), - get_late(correlator) + late .* agc_attenuation / 1 << (agc_bits + NC) + get_early(correlator) + early, + get_prompt(correlator) + prompt, + get_late(correlator) + late ) end +#= struct VeryEarlyPromptLateCorrelator{T} <: AbstractCorrelator{T} very_early::T early::T @@ -255,6 +256,7 @@ function normalize(correlator::VeryEarlyPromptLateCorrelator, integrated_samples get_very_late(correlator) / integrated_samples ) end +=# # TODO: correlate and dump for very early, very late #= @@ -267,7 +269,7 @@ Base.@propagate_inbounds @inline function correlate_iteration( prn, total_code_length, prompt_code_phase -) where S <: AbstractGNSSSystem +) where S <: AbstractGNSS early_code_phase = prompt_code_phase + code_phase_delta * early_late_sample_shift early_code_phase += (early_code_phase < 0) * total_code_length late_code_phase = prompt_code_phase - code_phase_delta * early_late_sample_shift @@ -281,3 +283,5 @@ Base.@propagate_inbounds @inline function correlate_iteration( VeryEarlyPromptLateCorrelator(early, prompt, late) end =# + +include("correlators/generic_correlator.jl") \ No newline at end of file diff --git a/src/correlators/generic_correlator.jl b/src/correlators/generic_correlator.jl new file mode 100644 index 0000000..af111d9 --- /dev/null +++ b/src/correlators/generic_correlator.jl @@ -0,0 +1,361 @@ +""" +$(SIGNATURES) + +Generic correlator holding a user defined number of correlation values. +""" +struct GenericCorrelator{T, M} <: AbstractCorrelator{T} + taps::Vector{T} + num_taps::NumTaps{M} + early_index::Int64 + prompt_index::Int64 + late_index::Int64 +end + +""" +$(SIGNATURES) + +GenericCorrelator constructor without parameters assumes a single antenna +and tree correlator elements. +""" +function GenericCorrelator() + GenericCorrelator(NumAnts(1)) +end + +""" +$(SIGNATURES) + +GenericCorrelator constructor that allows for the configuration of the +number of antenna elements using `num_ants::NumAnts{N}` The number of +correlator taps is fixed to three. +""" +function GenericCorrelator(num_ants::NumAnts{N}) where N + GenericCorrelator( + num_ants, + NumTaps(3), + 2 + ) +end + +""" +$(SIGNATURES) + +GenericCorrelator constructor that allows for the configuration of the +number of correlator taps using `num_taps::NumTaps{N}` The number of +antenne elements is fixed to one. +""" +function GenericCorrelator(num_taps::NumTaps{M}) where M + GenericCorrelator( + NumAnts(1), + num_taps, + 2 + ) +end + +""" +$(SIGNATURES) + +GenericCorrelator constructor that allows for the configuration of the +number of antenna elements using `num_ants::NumAnts{N}`, the number of +correlator taps using `num_taps::NumTaps{M}`. The early-late spacing +is fixed to 2. +""" +function GenericCorrelator(num_ants::NumAnts{N}, num_taps::NumTaps{M}) where {M, N} + prompt_idx = ceil(M/2) + GenericCorrelator( + num_ants, + num_taps, + 2 + ) +end + +""" +$(SIGNATURES) + +GenericCorrelator constructor for single antenna correlation, that +allows for the configuration of the number of correlator taps using +`num_taps::NumTaps{M}` and the spacing of the early and late correlator +using `el_spacing`. +""" +function GenericCorrelator(num_ants::NumAnts{1}, num_taps::NumTaps{M}, el_spacing) where M + prompt_index = ceil(Int64, M/2) + early_index = prompt_index + ceil(Int64, el_spacing/2) + late_index = prompt_index - ceil(Int64, el_spacing/2) + @assert late_index <= M + GenericCorrelator( + [zero(ComplexF64) for i = 1:M], + num_taps, + early_index, + prompt_index, + late_index + ) +end + +""" +$(SIGNATURES) + +GenericCorrelator constructor that allows for the configuration of the +number of antenna elements using `num_ants::NumAnts{N}`, the number of +correlator taps using `num_taps::NumTaps{M}` and the spacing of the +early and late correlator using `el_spacing`. +""" +function GenericCorrelator(num_ants::NumAnts{N}, num_taps::NumTaps{M}, el_spacing) where {M, N} + prompt_index = ceil(Int64, M/2) + early_index = prompt_index + ceil(Int64, el_spacing/2) + late_index = prompt_index - ceil(Int64, el_spacing/2) + @assert late_index <= M + GenericCorrelator( + [zero(SVector{N, ComplexF64}) for i = 1:M], + num_taps, + early_index, + prompt_index, + late_index + ) +end + + +""" +$(SIGNATURES) + +Get number of antennas from correlator +""" +get_num_ants(correlator::GenericCorrelator{Complex{T},M}) where {T,M} = 1 +get_num_ants(correlator::GenericCorrelator{SVector{N,T},M}) where {T,M,N} = N + +""" +$(SIGNATURES) + +Get number of correlator taps +""" +get_num_taps(correlator::GenericCorrelator{T,M}) where {T,M} = M + +""" +$(SIGNATURES) + +Get early correlator index +""" +get_early_index(correlator::GenericCorrelator) = correlator.early_index + +""" +$(SIGNATURES) + +Get prompt correlator index +""" +get_prompt_index(correlator::GenericCorrelator) = correlator.prompt_index + +""" +$(SIGNATURES) + +Get late correlator index +""" +get_late_index(correlator::GenericCorrelator) = correlator.late_index + +""" +$(SIGNATURES) + +Get all correlator taps +""" +get_taps(correlator::GenericCorrelator) = correlator.taps + +""" +$(SIGNATURES) + +Get a specific correlator tap with `index` counted negative for late and +positive for early correlators. +""" +function get_tap(correlator::GenericCorrelator, index::Integer) + correlator.taps[index+get_prompt_index(correlator)] +end + +""" +$(SIGNATURES) + +Get the early correlator +""" +function get_early(correlator::GenericCorrelator) + correlator.taps[get_early_index(correlator)] +end + +""" +$(SIGNATURES) + +Get the prompt correlator +""" +function get_prompt(correlator::GenericCorrelator) + correlator.taps[get_prompt_index(correlator)] +end + +""" +$(SIGNATURES) + +Get the late correlator +""" +function get_late(correlator::GenericCorrelator) + correlator.taps[get_late_index(correlator)] +end + +""" +$(SIGNATURES) + +Reset the Correlator +""" +function zero(correlator::GenericCorrelator) + GenericCorrelator( + zero(correlator.taps), + NumTaps(get_num_taps(correlator)), + get_early_index(correlator), + get_prompt_index(correlator), + get_late_index(correlator) + ) +end + +""" +$(SIGNATURES) + +Filter the correlator by the function `post_corr_filter` +""" +function filter(post_corr_filter, correlator::GenericCorrelator) + GenericCorrelator( + map(x->post_corr_filter(x), get_taps(correlator)), + NumTaps(get_num_taps(correlator)), + get_early_index(correlator), + get_prompt_index(correlator), + get_late_index(correlator) + ) +end + +""" +$(SIGNATURES) + +Calculate the replica phase offset required for the correlator taps with +respect to the prompt correlator, expressed in samples. The shifts are +ordered from latest to earliest replica. +""" +function get_correlator_sample_shifts( + system::AbstractGNSS, + correlator::GenericCorrelator{T,M}, + sampling_frequency, + preferred_code_shift +) where {T,M} + numEl = floor(Int, M/2) + SVector{M}(-numEl:numEl) .* round(Int, preferred_code_shift * sampling_frequency / get_code_frequency(system)) +end + +""" +$(SIGNATURES) + +Calculate the total spacing between early and late correlator in samples. +""" +function get_early_late_sample_spacing( + correlator::GenericCorrelator{T,M}, + correlator_sample_shifts::SVector{M} +) where {T,M} + correlator_sample_shifts[get_early_index(correlator)] - + correlator_sample_shifts[get_late_index(correlator)] +end + +""" +$(SIGNATURES) + +Normalize the correlator +""" +function normalize(correlator::GenericCorrelator, integrated_samples) + GenericCorrelator( + map(x->x/integrated_samples, get_taps(correlator)), + NumTaps(get_num_taps(correlator)), + get_early_index(correlator), + get_prompt_index(correlator), + get_late_index(correlator) + ) +end +""" +$(SIGNATURES) + +Perform a correlation for multi antenna systems +""" +function correlate( + correlator::GenericCorrelator, + downconverted_signal::AbstractVector, + code, + correlator_sample_shifts, + start_sample, + num_samples +) + taps = map(Vector(correlator_sample_shifts)) do correlator_sample_shift + correlate_single_tap( + correlator_sample_shift - correlator_sample_shifts[1], + start_sample, + num_samples, + downconverted_signal, + code + ) + end + + return GenericCorrelator( + map(+, get_taps(correlator), taps), + NumTaps(get_num_taps(correlator)), + get_early_index(correlator), + get_prompt_index(correlator), + get_late_index(correlator) + ) +end + +function correlate_single_tap( + offset, + start_sample, + num_samples, + downconverted_signal, + code +) + tap = zero(eltype(downconverted_signal)) + @inbounds for i = start_sample:num_samples + start_sample - 1 + tap += downconverted_signal[i] * code[i + offset] + end + tap +end + +""" +$(SIGNATURES) + +Perform a correlation for multi antenna systems +""" +function correlate( + correlator::GenericCorrelator{<: SVector{N}}, + downconverted_signal::AbstractMatrix, + code, + correlator_sample_shifts, + start_sample, + num_samples, +) where N + taps = map(Vector(correlator_sample_shifts)) do correlator_sample_shift + correlate_single_tap( + NumAnts(N), + correlator_sample_shift - correlator_sample_shifts[1], + start_sample, + num_samples, + downconverted_signal, + code + ) + end + GenericCorrelator( + map(+, get_taps(correlator), taps), + NumTaps(get_num_taps(correlator)), + get_early_index(correlator), + get_prompt_index(correlator), + get_late_index(correlator) + ) +end + +function correlate_single_tap( + ::NumAnts{N}, + offset, + start_sample, + num_samples, + downconverted_signal, + code +) where N + tap = zero(MVector{N, eltype(downconverted_signal)}) + @inbounds for j = 1:length(tap), i = start_sample:num_samples + start_sample - 1 + tap[j] += downconverted_signal[i,j] * code[i + offset] + end + SVector(tap) +end diff --git a/src/discriminators.jl b/src/discriminators.jl index 76cd79c..708c3d0 100644 --- a/src/discriminators.jl +++ b/src/discriminators.jl @@ -4,14 +4,14 @@ $(SIGNATURES) Calculates the code phase error in chips. """ function dll_disc( - ::Type{S}, + system::AbstractGNSS, correlator, - early_late_sample_shift, + correlator_sample_shifts, code_phase_delta -) where S <: AbstractGNSSSystem +) E = abs(get_early(correlator)) L = abs(get_late(correlator)) - distance_between_early_and_late = 2 * early_late_sample_shift * code_phase_delta + distance_between_early_and_late = get_early_late_sample_spacing(correlator, correlator_sample_shifts) * code_phase_delta (E - L) / (E + L) / (2 * (2 - distance_between_early_and_late)) end @@ -20,7 +20,7 @@ $(SIGNATURES) Calculates the carrier phase error in radians. """ -function pll_disc(::Type{S}, correlator) where S <: AbstractGNSSSystem +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 f39998a..ef11ada 100644 --- a/src/downconvert.jl +++ b/src/downconvert.jl @@ -1,38 +1,37 @@ +# TODO: This function might not be needed! Without @avx it can be merged with +# the two dimensional case function downconvert!( - downconverted_signal_re::AbstractVector, - downconverted_signal_im::AbstractVector, - carrier_re, - carrier_im, - signal_re::AbstractVector, - signal_im::AbstractVector, + downconverted_signal::StructArray{Complex{T}, 1}, + signal::StructArray{Complex{T}, 1}, + carrier_replica::StructArray{Complex{T}, 1}, start_sample::Integer, - num_samples_left::Integer -) - @avx unroll = 3 for i = start_sample:num_samples_left + start_sample - 1 - downconverted_signal_re[i] = signal_re[i] * carrier_re[i] + - signal_im[i] * carrier_im[i] - downconverted_signal_im[i] = signal_im[i] * carrier_re[i] - - signal_re[i] * carrier_im[i] + num_samples::Integer +) where T + ds_re = downconverted_signal.re; ds_im = downconverted_signal.im + s_re = signal.re; s_im = signal.im + c_re = carrier_replica.re; c_im = carrier_replica.im + @avx for i = start_sample:num_samples + start_sample - 1 + ds_re[i] = s_re[i] * c_re[i] + s_im[i] * c_im[i] + ds_im[i] = s_im[i] * c_re[i] - s_re[i] * c_im[i] end + downconverted_signal end function downconvert!( - downconverted_signal_re::AbstractMatrix, - downconverted_signal_im::AbstractMatrix, - carrier_re, - carrier_im, - signal_re::AbstractMatrix, - signal_im::AbstractMatrix, + downconverted_signal::StructArray{Complex{T}, 2}, + signal::StructArray{Complex{T}, 2}, + carrier_replica::StructArray{Complex{T}, 1}, start_sample::Integer, - num_samples_left::Integer -) - @avx unroll = 3 for i = start_sample:num_samples_left + start_sample - 1, j = 1:size(signal_re, 2) - # Calculate signal * carrier' - downconverted_signal_re[i, j] = signal_re[i, j] * carrier_re[i] + - signal_im[i, j] * carrier_im[i] - downconverted_signal_im[i, j] = signal_im[i, j] * carrier_re[i] - - signal_re[i, j] * carrier_im[i] + num_samples::Integer +) where T + ds_re = downconverted_signal.re; ds_im = downconverted_signal.im + s_re = signal.re; s_im = signal.im + c_re = carrier_replica.re; c_im = carrier_replica.im + @avx for i = start_sample:num_samples + start_sample - 1, j = 1:size(s_re, 2) + ds_re[i, j] = s_re[i, j] * c_re[i] + s_im[i, j] * c_im[i] + ds_im[i, j] = s_im[i, j] * c_re[i] - s_re[i, j] * c_im[i] end + downconverted_signal end function downconvert!( @@ -40,17 +39,9 @@ function downconvert!( signal, carrier_replica, start_sample::Integer, - num_samples_left::Integer -) - downconvert!( - downconverted_signal.re, - downconverted_signal.im, - carrier_replica.re, - carrier_replica.im, - signal.re, - signal.im, - start_sample, - num_samples_left - ) + num_samples::Integer +) + sample_range = start_sample:start_sample + num_samples - 1 + downconverted_signal[sample_range] .= @view(signal[sample_range]) .* conj.(@view(carrier_replica[sample_range])) downconverted_signal end diff --git a/src/galileo_e1b.jl b/src/galileo_e1b.jl index d557960..e9cc74f 100644 --- a/src/galileo_e1b.jl +++ b/src/galileo_e1b.jl @@ -3,7 +3,7 @@ $(SIGNATURES) Checks if upcoming integration is a new bit for GalileoE1B. """ -function is_upcoming_integration_new_bit(::Type{GalileoE1B}, prns, num_prns) +function is_upcoming_integration_new_bit(galileo_e1b::GalileoE1B, prns, num_prns) num_prns < 8 && return false masked_bit_synchronizer = prns & 0xff # First 8 bits # Upcoming integration will be a new bit if masked_bit_synchronizer contains @@ -12,6 +12,6 @@ function is_upcoming_integration_new_bit(::Type{GalileoE1B}, prns, num_prns) end # TODO: Very early very late correlator? -function get_default_correlator(::Type{GalileoE1B}, num_ants::NumAnts{N}) where N +function get_default_correlator(galileo_e1b::GalileoE1B, num_ants::NumAnts{N}) where N EarlyPromptLateCorrelator(num_ants) end diff --git a/src/gpsl1.jl b/src/gpsl1.jl index 0d727d1..065ca03 100644 --- a/src/gpsl1.jl +++ b/src/gpsl1.jl @@ -3,7 +3,7 @@ $(SIGNATURES) Checks if upcoming integration is a new bit for GPSL1. """ -function is_upcoming_integration_new_bit(::Type{GPSL1}, prns, num_prns) +function is_upcoming_integration_new_bit(gpsl1::GPSL1, prns, num_prns) num_prns < 40 && return false masked_bit_synchronizer = prns & 0xffffffffff # First 40 bits # Upcoming integration will be a new bit if masked_bit_synchronizer contains @@ -11,6 +11,6 @@ function is_upcoming_integration_new_bit(::Type{GPSL1}, prns, num_prns) masked_bit_synchronizer == 0xfffff || masked_bit_synchronizer == 0xfffff00000 end -function get_default_correlator(::Type{GPSL1}, num_ants::NumAnts{N}) where N +function get_default_correlator(gpsl1::GPSL1, num_ants::NumAnts{N}) where N EarlyPromptLateCorrelator(num_ants) end diff --git a/src/gpsl5.jl b/src/gpsl5.jl index 32d2a7f..d9d9918 100644 --- a/src/gpsl5.jl +++ b/src/gpsl5.jl @@ -3,7 +3,7 @@ $(SIGNATURES) Checks if upcoming integration is a new bit for GPSL5. """ -function is_upcoming_integration_new_bit(::Type{GPSL5}, prns, num_prns) +function is_upcoming_integration_new_bit(gpsl5::GPSL5, prns, num_prns) num_prns < 10 && return false masked_bit_synchronizer = prns & 0x3ff # First 10 xored_bit_synchronizer = masked_bit_synchronizer ⊻ 0x35 # 0x35 == 0000110101 @@ -11,6 +11,6 @@ function is_upcoming_integration_new_bit(::Type{GPSL5}, prns, num_prns) xored_bit_synchronizer == 0 || xored_bit_synchronizer == 0x3ff end -function get_default_correlator(::Type{GPSL5}, num_ants::NumAnts{N}) where N +function get_default_correlator(gpsl5::GPSL5, num_ants::NumAnts{N}) where N EarlyPromptLateCorrelator(num_ants) end diff --git a/src/secondary_code_or_bit_detector.jl b/src/secondary_code_or_bit_detector.jl index 0e2f40b..eeecc6c 100644 --- a/src/secondary_code_or_bit_detector.jl +++ b/src/secondary_code_or_bit_detector.jl @@ -22,10 +22,10 @@ $(SIGNATURES) Finds the secondary code or the bit shift. """ -function find(::Type{S}, detector, prompt_correlator) where S <: AbstractGNSSSystem +function find(system::AbstractGNSS, detector, prompt_correlator) found(detector) && return detector prns = get_buffer(detector) << 1 + UInt64(real(prompt_correlator) > 0) num_prns = length(detector) + 1 - bit_found = is_upcoming_integration_new_bit(S, prns, num_prns) + bit_found = is_upcoming_integration_new_bit(system, prns, num_prns) SecondaryCodeOrBitDetector(prns, num_prns, bit_found) end diff --git a/src/tracking_loop.jl b/src/tracking_loop.jl index d922182..9a73dc8 100644 --- a/src/tracking_loop.jl +++ b/src/tracking_loop.jl @@ -11,50 +11,49 @@ Track the signal `signal` based on the current tracking `state`, the sampling fr - 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 between early and late `early_late_sample_shift` defaults to - `get_early_late_sample_shift(...)` +- 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( - gain_controlled_signal::GainControlledSignal, - state::TrackingState{S, C, CALF, COLF, CN, DS}, + signal, + state::TrackingState{S, C, CALF, COLF, CN, DS, CAR, COR}, prn::Integer, sampling_frequency; post_corr_filter = get_default_post_corr_filter(get_correlator(state)), intermediate_frequency = 0.0Hz, max_integration_time::typeof(1ms) = 1ms, min_integration_time::typeof(1.0ms) = 0.75ms, - early_late_sample_shift = get_early_late_sample_shift(S, + correlator_sample_shifts = get_correlator_sample_shifts(get_system(state), get_correlator(state), sampling_frequency, 0.5), carrier_loop_filter_bandwidth = 18Hz, code_loop_filter_bandwidth = 1Hz, velocity_aiding = 0Hz, - carrier_amplitude_power::Val{N} = Val(5) ) where { - S <: AbstractGNSSSystem, + S <: AbstractGNSS, C <: AbstractCorrelator, CALF <: AbstractLoopFilter, COLF <: AbstractLoopFilter, CN <: AbstractCN0Estimator, - DS <: StructArray, - N + DS, + CAR, + COR } - if get_data_frequency(S) != 0Hz - @assert rem(1 / get_data_frequency(S), max_integration_time) == 0ms + system = get_system(state) + if get_data_frequency(system) != 0Hz + @assert rem(1 / get_data_frequency(system), max_integration_time) == 0ms end - N > 7 && throw(ArgumentError("The carrier amplitude power should be less than 8 to stay within 16 bits.")) - (get_amplitude_power(gain_controlled_signal) + N) > 16 && throw(ArgumentError("The AGC amplitude + carrier replica amplitude should not exceed 16 bits")) - signal = get_signal(gain_controlled_signal) 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.")) - agc_amplitude_power = get_amplitude_power(gain_controlled_signal) - agc_attenuation = get_attenuation(gain_controlled_signal) - downconverted_signal = resize!(get_downconverted_signal(state), size(signal, 1)) - carrier_replica = resize!(get_carrier(state), size(signal, 1)) - code_replica = resize!(get_code(state), size(signal, 1) + 2 * maximum(early_late_sample_shift)) + 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) @@ -75,7 +74,7 @@ function track( got_correlator = false while true num_samples_left_to_integrate = get_num_samples_left_to_integrate( - S, + system, max_integration_time, sampling_frequency, code_doppler, @@ -88,44 +87,23 @@ function track( intermediate_frequency, carrier_doppler ) - code_frequency = get_current_code_frequency(S, code_doppler) - code_replica = gen_code_replica!( + code_frequency = get_current_code_frequency(system, code_doppler) + correlator = downconvert_and_correlate!( + system, + signal, + correlator, code_replica, - S, - code_frequency, - sampling_frequency, code_phase, - signal_start_sample, - num_samples_left, - early_late_sample_shift, - prn - ) - carrier_replica = gen_carrier_replica!( carrier_replica, - carrier_frequency, - sampling_frequency, carrier_phase, - carrier_amplitude_power, - signal_start_sample, - num_samples_left - ) - downconverted_signal = downconvert!( downconverted_signal, - signal, - carrier_replica, - signal_start_sample, - num_samples_left - ) - correlator = correlate( - correlator, - downconverted_signal, - code_replica, - early_late_sample_shift, + code_frequency, + correlator_sample_shifts, + carrier_frequency, + sampling_frequency, signal_start_sample, num_samples_left, - agc_attenuation, - agc_amplitude_power, - carrier_amplitude_power + prn ) integrated_samples += num_samples_left carrier_phase = update_carrier_phase( @@ -133,11 +111,10 @@ function track( carrier_frequency, sampling_frequency, carrier_phase, - carrier_amplitude_power ) prev_code_phase = code_phase code_phase = update_code_phase( - S, + system, num_samples_left, code_frequency, sampling_frequency, @@ -154,11 +131,11 @@ function track( valid_correlator_carrier_phase = carrier_phase valid_correlator_carrier_frequency = carrier_frequency filtered_correlator = filter(post_corr_filter, correlator) - pll_discriminator = pll_disc(S, filtered_correlator) + pll_discriminator = pll_disc(system, filtered_correlator) dll_discriminator = dll_disc( - S, + system, filtered_correlator, - early_late_sample_shift, + correlator_sample_shifts, code_frequency / sampling_frequency ) carrier_freq_update, carrier_loop_filter = filter_loop( @@ -174,7 +151,7 @@ function track( code_loop_filter_bandwidth ) carrier_doppler, code_doppler = aid_dopplers( - S, + system, init_carrier_doppler, init_code_doppler, carrier_freq_update, @@ -183,7 +160,7 @@ function track( ) cn0_estimator = update(cn0_estimator, get_prompt(filtered_correlator)) bit_buffer, prompt_accumulator = buffer( - S, + system, bit_buffer, prompt_accumulator, found(sc_bit_detector), @@ -192,7 +169,7 @@ function track( max_integration_time, get_prompt(filtered_correlator) ) - sc_bit_detector = find(S, sc_bit_detector, get_prompt(filtered_correlator)) + sc_bit_detector = find(system, sc_bit_detector, get_prompt(filtered_correlator)) correlator = zero(correlator) integrated_samples = 0 end @@ -200,7 +177,8 @@ function track( num_samples_left == signal_samples_left && break signal_start_sample += num_samples_left end - next_state = TrackingState{S, C, CALF, COLF, CN, DS}( + next_state = TrackingState{S, C, CALF, COLF, CN, DS, CAR, COR}( + system, init_carrier_doppler, init_code_doppler, carrier_doppler, @@ -230,50 +208,72 @@ function track( ) end -@inline function track( - signal::AbstractArray, - state::TrackingState{S, C, CALF, COLF, CN, DS}, - prn::Integer, - sampling_frequency; - post_corr_filter = get_default_post_corr_filter(get_correlator(state)), - intermediate_frequency = 0.0Hz, - max_integration_time::typeof(1ms) = 1ms, - min_integration_time::typeof(1.0ms) = 0.75ms, - early_late_sample_shift = get_early_late_sample_shift(S, - get_correlator(state), sampling_frequency, 0.5), - carrier_loop_filter_bandwidth = 18Hz, - code_loop_filter_bandwidth = 1Hz, - velocity_aiding = 0Hz, - carrier_amplitude_power::Val{N} = Val(5) -) where { - S <: AbstractGNSSSystem, - C <: AbstractCorrelator, - CALF <: AbstractLoopFilter, - COLF <: AbstractLoopFilter, - CN <: AbstractCN0Estimator, - DS <: StructArray, - N -} - 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.")) - track( - GainControlledSignal(signal), - state, - prn, +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, - post_corr_filter = post_corr_filter, - intermediate_frequency = intermediate_frequency, - max_integration_time = max_integration_time, - min_integration_time = min_integration_time, - early_late_sample_shift = early_late_sample_shift, - carrier_loop_filter_bandwidth = carrier_loop_filter_bandwidth, - code_loop_filter_bandwidth = code_loop_filter_bandwidth, - velocity_aiding = velocity_aiding, - carrier_amplitude_power = carrier_amplitude_power + 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 +function choose(replica::CarrierReplicaCPU, signal::AbstractArray{Complex{Float64}}) + replica.carrier_f64 +end +function choose(replica::CarrierReplicaCPU, signal::AbstractArray{Complex{Float32}}) + replica.carrier_f32 +end +function choose(replica::DownconvertedSignalCPU, signal::AbstractArray{Complex{Float64}}) + replica.downconverted_signal_f64 +end +function choose(replica::DownconvertedSignalCPU, signal::AbstractArray{Complex{Float32}}) + replica.downconverted_signal_f32 +end + """ $(SIGNATURES) @@ -294,15 +294,15 @@ Returns the appropiate integration time. It will be the maximum integration time secondary code or the bit shift has been found. """ function get_integration_time( - ::Type{S}, + system::AbstractGNSS, max_integration_time, secondary_code_or_bit_found::Bool -) where S <: AbstractGNSSSystem +) ifelse( secondary_code_or_bit_found, max_integration_time, min( - convert(typeof(1ms), get_code_length(S) / get_code_frequency(S)), + ceil(typeof(1ms), get_code_length(system) / get_code_frequency(system)), max_integration_time ) ) @@ -314,13 +314,13 @@ $(SIGNATURES) Calculates the number of chips to integrate. """ function get_num_chips_to_integrate( - ::Type{S}, + system::AbstractGNSS, max_integration_time, current_code_phase, secondary_code_or_bit_found -) where S <: AbstractGNSSSystem - max_phase = Int(upreferred(get_code_frequency(S) * - get_integration_time(S, max_integration_time, 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) max_phase - current_phase_mod_max_phase end @@ -331,20 +331,20 @@ $(SIGNATURES) Calculates the number of samples to integrate. """ function get_num_samples_left_to_integrate( - ::Type{S}, + system::AbstractGNSS, max_integration_time, sampling_frequency, current_code_doppler, current_code_phase, secondary_code_or_bit_found -) where S <: AbstractGNSSSystem +) phase_to_integrate = get_num_chips_to_integrate( - S, + system, max_integration_time, current_code_phase, secondary_code_or_bit_found ) - code_frequency = get_code_frequency(S) + current_code_doppler + code_frequency = get_code_frequency(system) + current_code_doppler ceil(Int, phase_to_integrate * sampling_frequency / code_frequency) end @@ -355,15 +355,15 @@ Aid dopplers. That is velocity aiding for the carrier doppler and carrier aiding for the code doppler. """ function aid_dopplers( - ::Type{S}, + system::AbstractGNSS, init_carrier_doppler, init_code_doppler, carrier_freq_update, code_freq_update, velocity_aiding -) where S <: AbstractGNSSSystem +) carrier_doppler = carrier_freq_update + velocity_aiding - code_doppler = code_freq_update + carrier_doppler * get_code_center_frequency_ratio(S) + code_doppler = code_freq_update + carrier_doppler * get_code_center_frequency_ratio(system) init_carrier_doppler + carrier_doppler, init_code_doppler + code_doppler end @@ -380,10 +380,27 @@ end size(signal, 1) end -function resize!(A::StructArray{Complex{T}, 2}, b::Integer) where T - if size(A, 1) == b - return A - end - num_ants = size(A, 2) - StructArray{Complex{T}}((Matrix{T}(undef, b, num_ants), Matrix{T}(undef, b, num_ants))) +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{Float32}}) + 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 newline at end of file diff --git a/src/tracking_state.jl b/src/tracking_state.jl index be8e4d3..daaa604 100644 --- a/src/tracking_state.jl +++ b/src/tracking_state.jl @@ -1,16 +1,73 @@ """ $(SIGNATURES) +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 +end + +function CarrierReplicaCPU() + CarrierReplicaCPU( + StructArray{Complex{Float32}}(undef, 0), + StructArray{Complex{Float64}}(undef, 0) + ) +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 +end + +function DownconvertedSignalCPU(num_ants::NumAnts{1}) + DownconvertedSignalCPU( + StructArray{Complex{Float32}}(undef, 0), + StructArray{Complex{Float64}}(undef, 0) + ) +end + +function DownconvertedSignalCPU(num_ants::NumAnts{N}) where N + DownconvertedSignalCPU( + StructArray{Complex{Float32}}(undef, 0), + StructArray{Complex{Float64}}(undef, 0, N) + ) +end + + +""" +$(SIGNATURES) + TrackingState that holds the tracking state. """ struct TrackingState{ - S <: AbstractGNSSSystem, + S <: AbstractGNSS, C <: AbstractCorrelator, CALF <: AbstractLoopFilter, COLF <: AbstractLoopFilter, CN <: AbstractCN0Estimator, - DS <: StructArray, + DS <: DownconvertedSignalCPU, # Union{DownconvertedSignalCPU, CuArray{ComplexF32}} + CAR <: CarrierReplicaCPU, # Union{CarrierReplicaCPU, CuArray{ComplexF32}} + COR <: Vector{Int8}, # Union{Vector{Int8}, CuArray{Float32}} } + system::S init_carrier_doppler::typeof(1.0Hz) init_code_doppler::typeof(1.0Hz) carrier_doppler::typeof(1.0Hz) @@ -25,8 +82,8 @@ struct TrackingState{ prompt_accumulator::ComplexF64 cn0_estimator::CN downconverted_signal::DS - carrier::StructArray{Complex{Int16},1,NamedTuple{(:re, :im),Tuple{Array{Int16,1},Array{Int16,1}}},Int64} - code::Vector{Int16} + carrier::CAR + code::COR end """ @@ -46,37 +103,38 @@ carrier doppler `carrier_doppler` and the code phase `code_phase`. Optional para - CN0 estimator `cn0_estimator`, that defaults to `MomentsCN0Estimator(20)` """ function TrackingState( - ::Type{S}, + system::S, carrier_doppler, code_phase; - code_doppler = carrier_doppler * get_code_center_frequency_ratio(S), + 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(S, num_ants), + correlator::C = get_default_correlator(system, num_ants), integrated_samples = 0, prompt_accumulator = zero(ComplexF64), cn0_estimator::CN = MomentsCN0Estimator(20) ) where { - S <: AbstractGNSSSystem, + S <: AbstractGNSS, C <: AbstractCorrelator, CALF <: AbstractLoopFilter, COLF <: AbstractLoopFilter, CN <: AbstractCN0Estimator } if found(sc_bit_detector) - code_phase = mod(code_phase, get_code_length(S) * - get_secondary_code_length(S)) + code_phase = mod(code_phase, get_code_length(system) * + get_secondary_code_length(system)) else - code_phase = mod(code_phase, get_code_length(S)) + code_phase = mod(code_phase, get_code_length(system)) end - downconverted_signal = init_downconverted_signal(num_ants) - carrier = StructArray{Complex{Int16}}(undef, 0) - code = Vector{Int16}(undef, 0) + downconverted_signal = DownconvertedSignalCPU(num_ants) + carrier = CarrierReplicaCPU() + code = Vector{Int8}(undef, 0) - TrackingState{S, C, CALF, COLF, CN, typeof(downconverted_signal)}( + TrackingState{S, C, CALF, COLF, CN, typeof(downconverted_signal), typeof(carrier), typeof(code)}( + system, carrier_doppler, code_doppler, carrier_doppler, @@ -96,14 +154,7 @@ function TrackingState( ) end -function init_downconverted_signal(num_ants::NumAnts{1}) - StructArray{Complex{Int16}}(undef, 0) -end - -function init_downconverted_signal(num_ants::NumAnts{N}) where N - StructArray{Complex{Int16}}(undef, 0, N) -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 diff --git a/test/agc.jl b/test/agc.jl deleted file mode 100644 index 29e6f4a..0000000 --- a/test/agc.jl +++ /dev/null @@ -1,43 +0,0 @@ -@testset "AGC" begin - @testset "Single signal" begin - signal = rand(ComplexF64, 100) - agc_signal = GainControlledSignal(signal, 5) - real_max = maximum(real(x) for x in signal) - imag_max = maximum(imag(x) for x in signal) - -# attenuation = sqrt(real_max^2 + imag_max^2) - attenuation = max(real_max, imag_max) - amplification = 1 << 5 / attenuation - @test Tracking.get_attenuation(agc_signal) ≈ attenuation - @test Tracking.get_amplitude_power(agc_signal) == 5 - @test real.(agc_signal.signal) ≈ floor.(Int16, real.(signal) .* amplification) - @test imag.(agc_signal.signal) ≈ floor.(Int16, imag.(signal) .* amplification) - - signal = rand(Complex{Int16}, 100) - agc_signal = GainControlledSignal(signal, 5) - real_max = maximum(real(x) for x in signal) - imag_max = maximum(imag(x) for x in signal) -# attenuation = sqrt(float(real_max)^2 + float(imag_max)^2) - attenuation = max(real_max, imag_max) - amplification = 1 << 5 / attenuation - @test Tracking.get_attenuation(agc_signal) ≈ attenuation - @test Tracking.get_amplitude_power(agc_signal) == 5 - @test real.(agc_signal.signal) ≈ floor.(Int16, real.(signal) .* amplification) - @test imag.(agc_signal.signal) ≈ floor.(Int16, imag.(signal) .* amplification) - end - - @testset "Multiple signals" begin - signal = randn(ComplexF64, 100, 4) - agc_signal = GainControlledSignal(signal, 5) - real_max = map(X -> maximum(real(x) for x in X), eachcol(signal)) - imag_max = map(X -> maximum(imag(x) for x in X), eachcol(signal)) - -# attenuation = sqrt.(real_max.^2 .+ imag_max.^2) - attenuation = max.(real_max, imag_max) - amplification = (1 << 5) ./ attenuation - @test Tracking.get_attenuation(agc_signal) ≈ attenuation - @test Tracking.get_amplitude_power(agc_signal) == 5 - @test real.(agc_signal.signal) ≈ floor.(Int16, real.(signal) .* amplification') - @test imag.(agc_signal.signal) ≈ floor.(Int16, imag.(signal) .* amplification') - end -end diff --git a/test/bit_buffer.jl b/test/bit_buffer.jl index 0382988..4d3ffe0 100644 --- a/test/bit_buffer.jl +++ b/test/bit_buffer.jl @@ -10,9 +10,9 @@ code_phase = 1023 prompt_correlator = 1.0 + 0.0im integration_time = 5ms - + gpsl1 = GPSL1() next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( - GPSL1, + gpsl1, bit_buffer, prompt_accumulator, secondary_code_or_bit_found, @@ -23,7 +23,7 @@ ) @test @allocated(Tracking.buffer( - GPSL1, + gpsl1, Tracking.BitBuffer(), 0.0 + 0.0im, false, @@ -45,7 +45,7 @@ integration_time = 5ms next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( - GPSL1, + gpsl1, bit_buffer, prompt_accumulator, secondary_code_or_bit_found, @@ -66,7 +66,7 @@ integration_time = 20ms next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( - GPSL1, + gpsl1, bit_buffer, prompt_accumulator, secondary_code_or_bit_found, @@ -87,7 +87,7 @@ integration_time = 20ms next_bit_buffer, next_prompt_accumulator = @inferred Tracking.buffer( - GPSL1, + gpsl1, bit_buffer, prompt_accumulator, secondary_code_or_bit_found, diff --git a/test/boc.jl b/test/boc.jl new file mode 100644 index 0000000..2ddd4d8 --- /dev/null +++ b/test/boc.jl @@ -0,0 +1,18 @@ +@testset "BOCcos" begin + system = BOCcos(GPSL1(), 1, 1) + @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0xfffff00000, 50)) == true + @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0xfffff, 10)) == false + @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0xfffff, 40)) == true + + @test @inferred(Tracking.get_default_correlator(system, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1)) + @test @inferred(Tracking.get_default_correlator(system, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3)) + + system = BOCcos(GPSL5(), 1, 1) + @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0x35, 50)) == true + @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0x35, 5)) == false + @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0x3ca, 10)) == true # 0x3ca == 1111001010 + + @test @inferred(Tracking.get_default_correlator(system, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1)) + @test @inferred(Tracking.get_default_correlator(system, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3)) + +end diff --git a/test/carrier_replica.jl b/test/carrier_replica.jl index 9e2f202..c75520a 100644 --- a/test/carrier_replica.jl +++ b/test/carrier_replica.jl @@ -1,19 +1,20 @@ @testset "Carrier replica" begin - carrier = StructArray(zeros(Complex{Int16}, 2500)) + carrier = StructArray(zeros(Complex{Float32}, 2500)) Tracking.gen_carrier_replica!( carrier, 1500Hz, 2.5e6Hz, 0.25, # π / 2 - Val(7), 111, 2390 ) - @test sqrt(mean(abs2.(carrier.re[111:2500] ./ 1 << 7 .- - cos.(2π * (0:2389) * 1500Hz / 2.5e6Hz .+ π / 2)))) < 8e-3 + @test carrier[111:2500] ≈ cis.(2π * (0:2389) * 1500Hz / 2.5e6Hz .+ π / 2) + + # @test sqrt(mean(abs2.(carrier.re[111:2500] ./ 1 << 7 .- + # cos.(2π * (0:2389) * 1500Hz / 2.5e6Hz .+ π / 2)))) < 8e-3 end @testset "Update carrier phase" begin @@ -26,9 +27,8 @@ end carrier_frequency, sampling_frequency, carrier_phase, - Val(7) ) - @test phase ≈ mod(0.25 + 0.1 * 2000 + 0.5, 1) - 0.5 atol = 1 / 1 << 30 * 2500 + @test phase ≈ mod(0.25 + 0.1 * 2000 + 0.5, 1) - 0.5 carrier_phase = 0.25 carrier_frequency = 10Hz @@ -38,8 +38,7 @@ end num_samples, carrier_frequency, sampling_frequency, - carrier_phase, - Val(7) + carrier_phase ) - @test phase ≈ mod(0.25 + 10 * 2500 / 2.5e6 + 0.5, 1) - 0.5 atol = 1 / 1 << 30 * 2500 + @test phase ≈ mod(0.25 + 10 * 2500 / 2.5e6 + 0.5, 1) - 0.5 end diff --git a/test/cn0_estimation.jl b/test/cn0_estimation.jl index 6dffe45..287f9bc 100644 --- a/test/cn0_estimation.jl +++ b/test/cn0_estimation.jl @@ -34,39 +34,31 @@ end range = 0:3999 start_carrier_phase = π / 2 cn0_estimator = MomentsCN0Estimator(20) - early_late_sample_shift = 2 + correlator_sample_shifts = SVector(-2, 0, 2) start_sample = 1 num_samples = 4000 - agc_bits = 5 + gpsl1 = GPSL1() for i = 1:20 signal = get_code.( - GPSL1, + gpsl1, code_frequency .* range ./ sampling_frequency .+ start_code_phase, prn ) .* 10^(45 / 20) .+ randn(ComplexF64, length(range)) .* sqrt(4e6) - attenuation = Tracking.find_max(signal) - agc_signal = StructArray{Complex{Int16}}(( - floor.(Int16, real.(signal * 1 << agc_bits / attenuation)), - floor.(Int16, imag.(signal * 1 << agc_bits / attenuation)) - )) code = get_code.( - GPSL1, + gpsl1, code_frequency .* (-2:4001) ./ sampling_frequency .+ start_code_phase, prn ) correlator = EarlyPromptLateCorrelator() correlator = Tracking.correlate( correlator, - agc_signal, + signal, code, - early_late_sample_shift, + correlator_sample_shifts, start_sample, num_samples, - 1.0, - agc_bits, - Val(1) ) cn0_estimator = Tracking.update(cn0_estimator, get_prompt(correlator)) end @@ -74,6 +66,6 @@ end @test @inferred(Tracking.length(cn0_estimator)) == 20 cn0_estimate = @inferred Tracking.estimate_cn0(cn0_estimator, 1ms) - @test cn0_estimate ≈ 45dBHz atol = 0.30dBHz + @test cn0_estimate ≈ 45dBHz atol = 1.05dBHz end diff --git a/test/code_replica.jl b/test/code_replica.jl index 311ceff..df88624 100644 --- a/test/code_replica.jl +++ b/test/code_replica.jl @@ -1,30 +1,31 @@ @testset "Code replica" begin code = zeros(Int16, 2502) - + gpsl1 = GPSL1() Tracking.gen_code_replica!( code, - GPSL1, + gpsl1, 1023e3, 2.5e6, 2.0, 11, 2480, - 1, + SVector(-1, 0, 1), 1 ) - @test code[11:2492] == get_code.(GPSL1, (-1:2480) * 1023e3 / 2.5e6 .+ 2.0, 1) + @test code[11:2492] == get_code.(gpsl1, (-1:2480) * 1023e3 / 2.5e6 .+ 2.0, 1) end @testset "Update code phase" begin + gpsl1 = GPSL1() code_phase = 10 code_frequency = 10Hz sampling_frequency = 100Hz num_samples = 2000 bit_found = true phase = @inferred Tracking.update_code_phase( - GPSL1, + gpsl1, num_samples, code_frequency, sampling_frequency, diff --git a/test/correlator.jl b/test/correlator.jl index 461c190..5d0cf84 100644 --- a/test/correlator.jl +++ b/test/correlator.jl @@ -1,7 +1,7 @@ @testset "Correlator" begin @testset "Early prompt late correlator" begin - + gpsl1 = GPSL1() correlator = @inferred EarlyPromptLateCorrelator() @test @inferred(get_early(correlator)) == 0.0 @@ -68,14 +68,16 @@ 1.0 + 0.0im ) - early_late_sample_shift = @inferred get_early_late_sample_shift( - GPSL1, + correlator_sample_shifts = @inferred get_correlator_sample_shifts( + gpsl1, correlator, 4e6Hz, 0.5 ) - @test early_late_sample_shift ≈ round(0.5 * 4e6 / 1023e3) - @test typeof(early_late_sample_shift) <: Integer + @test correlator_sample_shifts ≈ round(0.5 * 4e6 / 1023e3) * SVector(-1, 0, 1) + @test typeof(correlator_sample_shifts) <: SVector{3,<:Integer} + + @test get_early_late_sample_spacing(correlator, SVector(-2,0,2)) == 4 normalized_correlator = @inferred Tracking.normalize(correlator, 10) @test normalized_correlator == EarlyPromptLateCorrelator( @@ -96,14 +98,14 @@ SVector(0.1 + 0.0im, 0.1 + 0.0im) ) - signal = StructArray{Complex{Int16}}( - (get_code.(GPSL1, (1:2500) * 1023e3 / 2.5e6, 1) * Int16(1) << (7 + 2), - zeros(Int16, 2500)) + signal = StructArray{Complex{Float32}}( + (Float32.(get_code.(gpsl1, (1:2500) * 1023e3 / 2.5e6, 1)), + zeros(Float32, 2500)) ) - early_late_sample_shift = 1 + correlator_sample_shifts = SVector(-1, 0, 1) code = get_code.( - GPSL1, - (1 - early_late_sample_shift:2500 + early_late_sample_shift) * 1023e3 / 2.5e6, + gpsl1, + (1 + correlator_sample_shifts[1]:2500 + correlator_sample_shifts[end]) * 1023e3 / 2.5e6, 1 ) correlator = EarlyPromptLateCorrelator(0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im) @@ -112,35 +114,30 @@ correlator, signal, code, - early_late_sample_shift, + correlator_sample_shifts, 1, - 2500, - 1.0, - 2, - Val(7) + 2500 ) @test get_early(correlator_result) == get_late(correlator_result) @test get_prompt(correlator_result) == 2500 signal_mat = repeat(signal, outer = (1,3)) - early_late_sample_shift = 1 + correlator_sample_shifts = SVector(-1,0,1) correlator = EarlyPromptLateCorrelator(NumAnts(3)) correlator_result = Tracking.correlate( correlator, - signal, + signal_mat, code, - early_late_sample_shift, + correlator_sample_shifts, 1, 2500, - [1.0, 1.0, 1.0], - 2, - Val(7) ) @test get_early(correlator_result) == get_late(correlator_result) @test all(get_early(correlator_result) .== 1476) @test all(get_prompt(correlator_result) .== 2500) end + include("generic_correlator.jl") end diff --git a/test/discriminators.jl b/test/discriminators.jl index 93fd688..8e68108 100644 --- a/test/discriminators.jl +++ b/test/discriminators.jl @@ -15,37 +15,37 @@ 1 + sqrt(3) * 1im, 0.5 + sqrt(3) / 2im ) - - @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° + gpsl1 = GPSL1() + @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 - + gpsl1 = GPSL1() very_early_correlator = EarlyPromptLateCorrelator(0.0 + 0.0im, 0.5 + 0.0im, 1.0 + 0.0im) early_correlator = EarlyPromptLateCorrelator(0.25 + 0.0im, 0.75 + 0.0im, 0.75 + 0.0im) prompt_correlator = EarlyPromptLateCorrelator(0.5 + 0.0im, 1 + 0.0im, 0.5 + 0.0im) late_correlator = EarlyPromptLateCorrelator(0.75 + 0.0im, 0.75 + 0.0im, 0.25 + 0.0im) very_late_correlator = EarlyPromptLateCorrelator(1.0 + 0.0im, 0.5 + 0.0im, 0.0 + 0.0im) - sample_shift = 2 + sample_shift = SVector(-2, 0, 2) delta = 0.25 @test @inferred( - Tracking.dll_disc(GPSL1, very_early_correlator, sample_shift, delta) + Tracking.dll_disc(gpsl1, very_early_correlator, sample_shift, delta) ) == -0.5 @test @inferred( - Tracking.dll_disc(GPSL1, early_correlator, sample_shift, delta) + Tracking.dll_disc(gpsl1, early_correlator, sample_shift, delta) ) == -0.25 @test @inferred( - Tracking.dll_disc(GPSL1, prompt_correlator, sample_shift, delta) + Tracking.dll_disc(gpsl1, prompt_correlator, sample_shift, delta) ) == 0 @test @inferred( - Tracking.dll_disc(GPSL1, late_correlator, sample_shift, delta) + Tracking.dll_disc(gpsl1, late_correlator, sample_shift, delta) ) == 0.25 @test @inferred( - Tracking.dll_disc(GPSL1, very_late_correlator, sample_shift, delta) + Tracking.dll_disc(gpsl1, very_late_correlator, sample_shift, delta) ) == 0.5 end diff --git a/test/galileo_e1b.jl b/test/galileo_e1b.jl index 81389d2..46193b7 100644 --- a/test/galileo_e1b.jl +++ b/test/galileo_e1b.jl @@ -1,14 +1,15 @@ @testset "Galileo E1B" begin - @test @inferred(Tracking.is_upcoming_integration_new_bit(GalileoE1B, 0xf, 29)) == true + galileo_e1b = GalileoE1B() + @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf, 29)) == true - @test @inferred(Tracking.is_upcoming_integration_new_bit(GalileoE1B, 0xf, 6)) == false + @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf, 6)) == false - @test @inferred(Tracking.is_upcoming_integration_new_bit(GalileoE1B, 0xc, 40)) == false + @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xc, 40)) == false - @test @inferred(Tracking.is_upcoming_integration_new_bit(GalileoE1B, 0xf, 8)) == true + @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf, 8)) == true - @test @inferred(Tracking.is_upcoming_integration_new_bit(GalileoE1B, 0xf0, 8)) == true + @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf0, 8)) == true - @test @inferred(Tracking.get_default_correlator(GalileoE1B, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1)) - @test @inferred(Tracking.get_default_correlator(GalileoE1B, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3)) + @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/generic_correlator.jl b/test/generic_correlator.jl new file mode 100644 index 0000000..0145ca2 --- /dev/null +++ b/test/generic_correlator.jl @@ -0,0 +1,184 @@ + + +@testset "Generic Correlator" begin + +@testset "Empty constructor" begin + correlator = @inferred GenericCorrelator() + @test @inferred(get_early(correlator)) == 0.0 + @test @inferred(get_prompt(correlator)) == 0.0 + @test @inferred(get_late(correlator)) == 0.0 + @test length(get_taps(correlator)) == 3 + @test @inferred(get_num_taps(correlator)) == 3 + @test @inferred(get_early_index(correlator)) == 3 + @test @inferred(get_prompt_index(correlator)) == 2 + @test @inferred(get_late_index(correlator)) == 1 +end + +@testset "Single antenna" begin + correlator = @inferred GenericCorrelator(NumAnts(1)) + @test @inferred(get_early(correlator)) == 0.0 + @test @inferred(get_prompt(correlator)) == 0.0 + @test @inferred(get_late(correlator)) == 0.0 + @test length(get_taps(correlator)) == 3 + @test @inferred(get_num_taps(correlator)) == 3 + @test @inferred(get_early_index(correlator)) == 3 + @test @inferred(get_prompt_index(correlator)) == 2 + @test @inferred(get_late_index(correlator)) == 1 +end + +@testset "Multiple Antennas" begin + correlator = @inferred GenericCorrelator(NumAnts(2)) + @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 length(get_taps(correlator)) == 3 + @test @inferred(get_num_taps(correlator)) == 3 + @test @inferred(get_early_index(correlator)) == 3 + @test @inferred(get_prompt_index(correlator)) == 2 + @test @inferred(get_late_index(correlator)) == 1 +end + +@testset "Multiple antennas, multiple taps" begin + correlator = @inferred GenericCorrelator(NumAnts(2), NumTaps(5)) + for k in -2:2 + @test @inferred(get_tap(correlator, k)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + end + @test length(get_taps(correlator)) == 5 + @test @inferred(get_num_taps(correlator)) == 5 + @test @inferred(get_early_index(correlator)) == 4 + @test @inferred(get_prompt_index(correlator)) == 3 + @test @inferred(get_late_index(correlator)) == 2 +end + +@testset "Multiple antennas, multiple taps, custom spacing" begin + correlator = @inferred GenericCorrelator(NumAnts(2), NumTaps(7), 4) + for k in -3:3 + @test @inferred(get_tap(correlator, k)) == SVector(0.0 + 0.0im, 0.0 + 0.0im) + end + @test length(get_taps(correlator)) == 7 + @test @inferred(get_num_taps(correlator)) == 7 + @test @inferred(get_early_index(correlator)) == 6 + @test @inferred(get_prompt_index(correlator)) == 4 + @test @inferred(get_late_index(correlator)) == 2 +end + +@testset "Default constructor" begin + correlator = @inferred GenericCorrelator( + [SVector(i + 1.0i*im, i + 2.0im) for i = 1:7], + NumTaps(7), + 1, 4, 7 + ) + @test @inferred(get_early(correlator)) == SVector(1 + 1.0im, 1 + 2.0im) + @test @inferred(get_prompt(correlator)) == SVector(4 + 4.0im, 4 + 2.0im) + @test @inferred(get_late(correlator)) == SVector(7 + 7.0im, 7 + 2.0im) +end + +@testset "Correlator zeroing" begin + correlator = @inferred GenericCorrelator( + [SVector(i + 1im, 2i + 3im) for i = 1:3], + NumTaps(3), + 3, 2, 1 + ) + zero_correlator = @inferred Tracking.zero(correlator) + @test all([all(t .== 0) for t in get_taps(zero_correlator)]) +end + +@testset "Correlator filtering" begin + correlator = @inferred GenericCorrelator( + [SVector(1.0 + 0im, 1.0 + 0im) for i = 1:3], + NumTaps(3), + 1, 2, 3 + ) + filtered_correlator = @inferred Tracking.filter(x->x[1], correlator) + @test get_taps(filtered_correlator) == [1.0 + 0im for i = 1:3] + filtered_correlator = @inferred Tracking.filter(x -> x, correlator) + @test get_taps(filtered_correlator) == [SVector(1.0 + 0im, 1.0 + 0im) for i = 1:3] +end + +@testset "Sampleshift calculation" begin + gpsl1 = GPSL1() + correlator = @inferred GenericCorrelator(NumAnts(2), NumTaps(7)) + correlator_sample_shifts = @inferred get_correlator_sample_shifts( + gpsl1, + correlator, + 4e6Hz, + 0.5 + ) + @test correlator_sample_shifts ≈ SVector(-3, -2, -1, 0, 1, 2, 3) .* round(0.5 * 4e6 / 1.023e6) + @test typeof(correlator_sample_shifts) <: SVector +end + +@testset "Correlator normalization" begin + correlator = @inferred GenericCorrelator( + [SVector(10i + 0im, 20i + 0im) for i = 1:3], + NumTaps(3), + 1, 2, 3 + ) + normalized_correlator = @inferred Tracking.normalize(correlator, 10) + @test get_taps(normalized_correlator) == [SVector(1i + 0im, 2i + 0im) for i = 1:3] +end + +@testset "Single antenna correlation" begin + gpsl1 = GPSL1() + signal = StructArray{Complex{Float32}}( + (Float32.(get_code.(gpsl1, (1:2500) * 1023e3 / 2.5e6, 1)), + zeros(Float32, 2500)) + ) + correlator = GenericCorrelator(NumAnts(1), NumTaps(5)) + correlator_sample_shifts = get_correlator_sample_shifts( + gpsl1, + correlator, + 2.5e6Hz, + 0.5 + ) + code = get_code.( + gpsl1, + (1 + correlator_sample_shifts[1]:2500 + correlator_sample_shifts[end]) * 1023e3 / 2.5e6, + 1 + ) + correlator_result = Tracking.correlate( + correlator, + signal, + code, + correlator_sample_shifts, + 1, + 2500 + ) + @test get_prompt(correlator_result) == 2500 + @test get_early(correlator_result) == get_late(correlator_result) + @test get_taps(correlator_result) == reverse(get_taps(correlator_result)) +end + +@testset "Multi antenna correlation" begin + gpsl1 = GPSL1() + signal = StructArray{Complex{Float32}}( + (Float32.(get_code.(gpsl1, (1:2500) * 1023e3 / 2.5e6, 1)), + zeros(Float32, 2500)) + ) + correlator = GenericCorrelator(NumAnts(3), NumTaps(5)) + correlator_sample_shifts = get_correlator_sample_shifts( + gpsl1, + correlator, + 2.5e6Hz, + 0.5 + ) + code = get_code.( + gpsl1, + (1 + correlator_sample_shifts[1]:2500 + correlator_sample_shifts[end]) * 1023e3 / 2.5e6, + 1 + ) + signal_mat = repeat(signal, outer = (1,3)) + correlator_result = Tracking.correlate( + correlator, + signal_mat, + code, + correlator_sample_shifts, + 1, + 2500, + ) + @test all(get_prompt(correlator_result) .== 2500) + @test all(get_early(correlator_result) .== 1476) + @test get_taps(correlator) == reverse(get_taps(correlator)) +end + +end diff --git a/test/gps_l1.jl b/test/gps_l1.jl index 2e0c15e..79eacf0 100644 --- a/test/gps_l1.jl +++ b/test/gps_l1.jl @@ -1,10 +1,11 @@ @testset "GPS L1" begin - @test @inferred(Tracking.is_upcoming_integration_new_bit(GPSL1, 0xfffff00000, 50)) == true + gpsl1 = GPSL1() + @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, 10)) == false - @test @inferred(Tracking.is_upcoming_integration_new_bit(GPSL1, 0xfffff, 40)) == true + @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)) + @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/gps_l5.jl b/test/gps_l5.jl index 5b9da3c..755511e 100644 --- a/test/gps_l5.jl +++ b/test/gps_l5.jl @@ -1,10 +1,11 @@ @testset "GPS L5" begin - @test @inferred(Tracking.is_upcoming_integration_new_bit(GPSL5, 0x35, 50)) == true + gpsl5 = GPSL5() + @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, 0x35, 5)) == false - @test @inferred(Tracking.is_upcoming_integration_new_bit(GPSL5, 0x3ca, 10)) == true # 0x3ca == 1111001010 + @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)) + @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/runtests.jl b/test/runtests.jl index be92109..868e538 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,7 +10,6 @@ using import Unitful: MHz, kHz, Hz, s, ms, dBHz -include("agc.jl") include("code_replica.jl") include("carrier_replica.jl") include("downconvert.jl") @@ -18,6 +17,7 @@ include("discriminators.jl") include("gps_l1.jl") include("gps_l5.jl") include("galileo_e1b.jl") +include("boc.jl") include("secondary_code_or_bit_detector.jl") include("bit_buffer.jl") include("correlator.jl") diff --git a/test/secondary_code_or_bit_detector.jl b/test/secondary_code_or_bit_detector.jl index e62cc01..aef3819 100644 --- a/test/secondary_code_or_bit_detector.jl +++ b/test/secondary_code_or_bit_detector.jl @@ -1,21 +1,22 @@ @testset "Secondary code or bit detector" begin detector = SecondaryCodeOrBitDetector() + gpsl1 = GPSL1() @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) + 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) + 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) + 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 @@ -26,13 +27,13 @@ @test Tracking.length(detector) == 1 @test Tracking.found(detector) == true - next_detector = Tracking.find(GPSL1, detector, 1.0 + 0.0im) + 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) + 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 diff --git a/test/tracking_loop.jl b/test/tracking_loop.jl index d01a3bc..1d5d40a 100644 --- a/test/tracking_loop.jl +++ b/test/tracking_loop.jl @@ -9,20 +9,22 @@ 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) + 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) + 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( - GalileoE1B, + galileo_e1b, max_integration_time, bit_found ) @@ -31,7 +33,7 @@ end max_integration_time = 10ms bit_found = false time = @inferred Tracking.get_integration_time( - GalileoE1B, + galileo_e1b, max_integration_time, bit_found ) @@ -39,11 +41,12 @@ end 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, + gpsl1, max_integration_time, code_phase, bit_found @@ -52,7 +55,7 @@ end code_phase = 1200 chips = @inferred Tracking.get_num_chips_to_integrate( - GPSL1, + gpsl1, max_integration_time, code_phase, bit_found @@ -61,13 +64,14 @@ end 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, + gpsl1, max_integration_time, sampling_frequency, current_code_doppler, @@ -88,13 +92,14 @@ end end @testset "Code phase delta" begin + gpsl1 = GPSL1() code_doppler = 1Hz - frequency = @inferred Tracking.get_current_code_frequency(GPSL1, code_doppler) + 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 @@ -102,7 +107,7 @@ end velocity_aiding = 3Hz carrier_freq, code_freq = @inferred Tracking.aid_dopplers( - GPSL1, + gpsl1, init_carrier_doppler, init_code_doppler, carrier_freq_update, @@ -115,7 +120,7 @@ end end @testset "Tracking" begin - + gpsl1 = GPSL1() carrier_doppler = 200Hz start_code_phase = 100 code_frequency = carrier_doppler / 1540 + 1023kHz @@ -123,32 +128,17 @@ end prn = 1 range = 0:3999 start_carrier_phase = π / 2 - state = TrackingState(GPSL1, carrier_doppler - 20Hz, start_code_phase) + state = TrackingState(gpsl1, carrier_doppler - 20Hz, start_code_phase) signal = cis.( 2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase ) .* get_code.( - GPSL1, + gpsl1, code_frequency .* range ./ sampling_frequency .+ start_code_phase, prn ) - agc_signal = GainControlledSignal(signal, 11) - @test_throws ArgumentError track( - agc_signal, - state, - prn, - sampling_frequency, - carrier_amplitude_power = Val(6) - ) - @test_throws ArgumentError track( - signal, - state, - prn, - sampling_frequency, - carrier_amplitude_power = Val(8) - ) track_result = @inferred track(signal, state, prn, sampling_frequency) iterations = 2000 @@ -170,7 +160,7 @@ end 2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase ) .* get_code.( - GPSL1, + gpsl1, code_frequency .* range ./ sampling_frequency .+ code_phase, prn ) @@ -215,7 +205,7 @@ end end @testset "Track multiple signals" begin - + gpsl1 = GPSL1() carrier_doppler = 200Hz start_code_phase = 100 code_frequency = carrier_doppler / 1540 + 1023kHz @@ -223,13 +213,13 @@ end prn = 1 range = 0:3999 start_carrier_phase = π / 2 - state = TrackingState(GPSL1, carrier_doppler - 20Hz, start_code_phase, num_ants = NumAnts(3)) + state = TrackingState(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, + gpsl1, code_frequency .* range ./ sampling_frequency .+ start_code_phase, prn ) @@ -257,7 +247,7 @@ end 2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase ) .* get_code.( - GPSL1, + gpsl1, code_frequency .* range ./ sampling_frequency .+ code_phase, prn ) diff --git a/test/tracking_results.jl b/test/tracking_results.jl index 6bcd02a..d93d600 100644 --- a/test/tracking_results.jl +++ b/test/tracking_results.jl @@ -1,6 +1,7 @@ @testset "Tracking results" begin + gpsl1 = GPSL1() results = Tracking.TrackingResults( - TrackingState(GPSL1, 100Hz, 100), + TrackingState(gpsl1, 100Hz, 100), EarlyPromptLateCorrelator(NumAnts(2)), 1.0Hz, 1.0, @@ -26,7 +27,7 @@ @test @inferred(get_cn0(results)) == 45dBHz results = Tracking.TrackingResults( - TrackingState(GPSL1, 100Hz, 100), + TrackingState(gpsl1, 100Hz, 100), EarlyPromptLateCorrelator(NumAnts(2)), 1.0Hz, 1.0, diff --git a/test/tracking_state.jl b/test/tracking_state.jl index cd8eece..7053d92 100644 --- a/test/tracking_state.jl +++ b/test/tracking_state.jl @@ -2,8 +2,8 @@ carrier_doppler = 100Hz code_phase = 100 - - state = TrackingState(GPSL1, carrier_doppler, code_phase) + gpsl1 = GPSL1() + state = TrackingState(gpsl1, carrier_doppler, code_phase) @test @inferred(Tracking.get_code_phase(state)) == 100 @test @inferred(Tracking.get_carrier_phase(state)) == 0.0 @@ -19,15 +19,15 @@ @test @inferred(Tracking.get_integrated_samples(state)) == 0 state = TrackingState( - GPSL1, + gpsl1, carrier_doppler, code_phase; - code_doppler = carrier_doppler * get_code_center_frequency_ratio(GPSL1), + 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)), + correlator = Tracking.get_default_correlator(gpsl1, NumAnts(1)), integrated_samples = 0, prompt_accumulator = zero(ComplexF64) ) @@ -45,7 +45,7 @@ @test @inferred(Tracking.get_prompt_accumulator(state)) == 0.0 @test @inferred(Tracking.get_integrated_samples(state)) == 0 - state = TrackingState(GPSL1, carrier_doppler, code_phase, num_ants = NumAnts(2)) + state = TrackingState(gpsl1, carrier_doppler, code_phase, num_ants = NumAnts(2)) @test @inferred(Tracking.get_correlator(state)) == EarlyPromptLateCorrelator(NumAnts(2)) end