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 <michael.niestroj@nav.rwth-aachen.de> Co-authored-by: Soeren Schoenbrod <soeren.schoenbrod@rwth-aachen.de>
JuliaGNSS · Feb 3, 2021 · 94bcda9 · 94bcda9
1 parent 71d7c62
commit 94bcda9
Show file tree

Hide file tree

Showing 34 changed files with 1,060 additions and 524 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Tracking"
 uuid = "10b2438b-ffd4-5096-aa58-44041d5c8f3b"
 authors = ["Soeren Zorn <soeren.zorn@nav.rwth-aachen.de>"]
-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
@@ -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
diff --git a/src/bit_buffer.jl b/src/bit_buffer.jl
@@ -21,20 +21,20 @@ $(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,
     prev_code_phase,
     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
@@ -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
@@ -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
@@ -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