Increase performance (#5)

* Increase performance

* Add coarse fine acquire to Readme.md

* Change random seed

* Further increase performance

* Fix failing test

* Add JuliaFormatter

.JuliaFormatter.toml

@@ -0,0 +1,8 @@
+# Configuration file for JuliaFormatter.jl
+# For more information, see: https://domluna.github.io/JuliaFormatter.jl/stable/config/
+
+margin = 92
+remove_extra_newlines = true
+separate_kwargs_with_semicolon = true
+format_docstrings = true
+format_markdown = true

Project.toml

@@ -8,6 +8,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 GNSSSignals = "52c80523-2a4e-5c38-8979-05588f836870"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
@@ -16,8 +17,9 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 DocStringExtensions = "0.6, 0.7, 0.8, 0.9"
 FFTW = "1.0"
 GNSSSignals = "0.15, 0.16"
-Unitful = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 1.0"
+LoopVectorization = "0.12"
 RecipesBase = "1.0"
+Unitful = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 1.0"
 julia = "1.6"
 
 [extras]

README.md

@@ -17,11 +17,12 @@ pkg> add Acquisition
 ```julia
 using Acquisition, Plots
 import Acquisition: GPSL1, Hz
-stream1 = open("signal.dat")
+stream = open("signal.dat")
 signal = Vector{Complex{Int16}}(undef, 10000)
-read!(stream1, signal)
+read!(stream, signal)
 gpsl1 = GPSL1()
-acq_res = [acquire(gpsl1, signal, 5e6Hz, prn) for prn = 1:32]
+acq_res = acquire(gpsl1, signal, 5e6Hz, 1:32)
+# or acq_res = coarse_fine_acquire(gpsl1, signal, 5e6Hz, 1:32)
 plot(acq_res[1])
 ```
 

src/Acquisition.jl

@@ -1,86 +1,38 @@
 module Acquisition
 
-    using DocStringExtensions, GNSSSignals, RecipesBase, FFTW, Statistics, LinearAlgebra
-    import Unitful: s, Hz
-
-    export acquire, plot_acquisition_results, coarse_fine_acquire
-
-    struct AcquisitionResults{S<:AbstractGNSS}
-        system::S
-        prn::Int
-        sampling_frequency::typeof(1.0Hz)
-        carrier_doppler::typeof(1.0Hz)
-        code_phase::Float64
-        CN0::Float64
-        power_bins::Array{Float64, 2}
-        dopplers::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
-    end
-
-    @recipe function f(acq_res::AcquisitionResults, log_scale::Bool = false;)
-        seriestype := :surface
-        seriescolor --> :viridis
-        yguide --> "Code phase"
-        xguide --> "Dopplers [Hz]"
-        zguide --> (log_scale ? "Magnitude [dB]" : "Magnitude")
-        y = (1:size(acq_res.power_bins, 1)) ./ acq_res.sampling_frequency .* get_code_frequency(acq_res.system)
-        acq_res.dopplers, y, log_scale ? 10 * log10.(acq_res.power_bins) : acq_res.power_bins
-    end
-
-    """
-    $(SIGNATURES)
-    Perform the aquisition of the satellite `sat_prn` in System `S` in signal `signal`
-    sampled at rate `sampling_freq`. The aquisition is performed as parallel code phase
-    search using the doppler frequencies `dopplers`.
-    """
-    function acquire(S::AbstractGNSS, signal, sampling_freq, sat_prn; interm_freq = 0.0Hz, max_doppler = 7000Hz, dopplers = -max_doppler:1 / 3 / (length(signal) / sampling_freq):max_doppler)
-        code_period = get_code_length(S) / get_code_frequency(S)
-        powers = power_over_doppler_and_code(S, signal, sat_prn, dopplers, sampling_freq, interm_freq)
-        signal_power, noise_power, code_index, doppler_index = est_signal_noise_power(powers, sampling_freq, get_code_frequency(S))
-        CN0 = 10 * log10(signal_power / noise_power / code_period / 1.0Hz)
-        doppler = (doppler_index - 1) * step(dopplers) + first(dopplers)
-        code_phase = (code_index - 1) / (sampling_freq / get_code_frequency(S))
-        AcquisitionResults(S, sat_prn, sampling_freq, doppler, code_phase, CN0, powers, dopplers / 1.0Hz)
-    end
-
-    """
-    $(SIGNATURES)
-    Performs a coarse aquisition and fine acquisition of the satellite `sat_prn` in System `S` in signal `signal`
-    sampled at rate `sampling_freq`. The aquisition is performed as parallel code phase
-    search using the doppler frequencies `dopplers`.
-    """
-    function coarse_fine_acquire(S::AbstractGNSS, signal, sampling_freq, sat_prn; interm_freq = 0.0Hz, max_doppler = 7000Hz)
-        coarse_step = 1 / 3 / (length(signal) / sampling_freq)
-        acq_res = acquire(S, signal, sampling_freq, sat_prn; interm_freq, dopplers = -max_doppler:coarse_step:max_doppler)
-        fine_step = 1 / 12 / (length(signal) / sampling_freq)
-        dopplers = acq_res.carrier_doppler - 2 * coarse_step:fine_step:acq_res.carrier_doppler + 2 * coarse_step
-        acquire(S, signal, sampling_freq, sat_prn; interm_freq, dopplers = dopplers)
-    end
-
-    function power_over_doppler_and_code(S::AbstractGNSS, signal, sat_prn, doppler_steps, sampling_freq, interm_freq)
-        code = get_code.(S, (1:length(signal)) .* get_code_frequency(S) ./ sampling_freq, sat_prn)
-        fft_plan = plan_fft(code)
-        code_freq_domain = fft_plan * code
-        mapreduce(doppler -> power_over_code(S, signal, fft_plan, code_freq_domain, doppler, sampling_freq, interm_freq), hcat, doppler_steps)
-    end
-
-    function power_over_code(S::AbstractGNSS, signal, fft_plan, code_freq_domain, doppler, sampling_freq, interm_freq)
-        Δt = length(signal) / sampling_freq
-        code_interval = get_code_length(S) / get_code_frequency(S)
-        signal_baseband = signal .* cis.(-2π .* (1:length(signal)) .* (interm_freq + doppler) ./ sampling_freq)
-        signal_baseband_freq_domain = fft_plan * signal_baseband
-        code_freq_baseband_freq_domain = code_freq_domain .* conj(signal_baseband_freq_domain)
-        powers = abs2.(fft_plan \ code_freq_baseband_freq_domain)
-        powers[1:convert(Int, sampling_freq * min(Δt, code_interval))]
-    end
+using DocStringExtensions,
+    GNSSSignals, RecipesBase, FFTW, Statistics, LinearAlgebra, LoopVectorization, Unitful
+
+import Unitful: s, Hz
+
+export acquire,
+    plot_acquisition_results,
+    coarse_fine_acquire,
+    coarse_fine_acquire!,
+    acquire!,
+    AcquisitionPlan,
+    CoarseFineAcquisitionPlan
+
+struct AcquisitionResults{S<:AbstractGNSS,T}
+    system::S
+    prn::Int
+    sampling_frequency::typeof(1.0Hz)
+    carrier_doppler::typeof(1.0Hz)
+    code_phase::Float64
+    CN0::Float64
+    noise_power::T
+    power_bins::Matrix{T}
+    dopplers::StepRangeLen{
+        Float64,
+        Base.TwicePrecision{Float64},
+        Base.TwicePrecision{Float64},
+    }
+end
 
-    function est_signal_noise_power(power_bins, sampling_freq, code_freq)
-        samples_per_chip = floor(Int, sampling_freq / code_freq)
-        signal_noise_power, index = findmax(power_bins)
-        lower_code_phases = 1:index[1] - samples_per_chip
-        upper_code_phases = index[1] + samples_per_chip:size(power_bins, 1)
-        samples = (length(lower_code_phases) + length(upper_code_phases)) * size(power_bins, 2)
-        noise_power = (sum(power_bins[lower_code_phases,:]) + sum(power_bins[upper_code_phases,:])) / samples
-        signal_power = signal_noise_power - noise_power
-        signal_power, noise_power, index[1], index[2]
-    end
+include("plan_acquire.jl")
+include("downconvert.jl")
+include("plot.jl")
+include("calc_powers.jl")
+include("est_signal_noise_power.jl")
+include("acquire.jl")
 end

src/acquire.jl

@@ -0,0 +1,187 @@
+"""
+$(SIGNATURES)
+Perform the aquisition of multiple satellites with `prns` in system `system` with signal `signal`
+sampled at rate `sampling_freq`. Optional arguments are the intermediate frequency `interm_freq`
+(default 0Hz), the maximum expected Doppler `max_doppler` (default 7000Hz). If the maximum Doppler
+is too unspecific you can instead pass a Doppler range with with your individual step size using
+the argument `dopplers`.
+"""
+function acquire(
+    system::AbstractGNSS,
+    signal,
+    sampling_freq,
+    prns::AbstractVector{<:Integer};
+    interm_freq = 0.0Hz,
+    max_doppler = 7000Hz,
+    dopplers = -max_doppler:1/3/(length(signal)/sampling_freq):max_doppler,
+)
+    acq_plan = AcquisitionPlan(
+        system,
+        length(signal),
+        sampling_freq;
+        dopplers,
+        prns,
+        fft_flag = FFTW.ESTIMATE,
+    )
+    acquire!(acq_plan, signal, prns; interm_freq)
+end
+
+"""
+$(SIGNATURES)
+This acquisition function uses a predefined acquisition plan to accelerate the computing time.
+This will be useful, if you have to calculate the acquisition multiple time in a row.
+"""
+function acquire!(
+    acq_plan::AcquisitionPlan,
+    signal,
+    prns::AbstractVector{<:Integer};
+    interm_freq = 0.0Hz,
+    doppler_offset = 0.0Hz,
+    noise_power = nothing,
+)
+    all(map(prn -> prn in acq_plan.avail_prn_channels, prns)) ||
+        throw(ArgumentError("You'll need to plan every PRN"))
+    code_period = get_code_length(acq_plan.system) / get_code_frequency(acq_plan.system)
+    powers_per_sats =
+        power_over_doppler_and_codes!(acq_plan, signal, prns, interm_freq, doppler_offset)
+    map(powers_per_sats, prns) do powers, prn
+        signal_power, noise_power, code_index, doppler_index = est_signal_noise_power(
+            powers,
+            acq_plan.sampling_freq,
+            get_code_frequency(acq_plan.system),
+            noise_power,
+        )
+        CN0 = 10 * log10(signal_power / noise_power / code_period / 1.0Hz)
+        doppler =
+            (doppler_index - 1) * step(acq_plan.dopplers) +
+            first(acq_plan.dopplers) +
+            doppler_offset
+        code_phase =
+            (code_index - 1) /
+            (acq_plan.sampling_freq / get_code_frequency(acq_plan.system))
+        AcquisitionResults(
+            acq_plan.system,
+            prn,
+            acq_plan.sampling_freq,
+            doppler,
+            code_phase,
+            CN0,
+            noise_power,
+            powers,
+            (acq_plan.dopplers .+ doppler_offset) / 1.0Hz,
+        )
+    end
+end
+
+function acquire!(
+    acq_plan::CoarseFineAcquisitionPlan,
+    signal,
+    prns::AbstractVector{<:Integer};
+    interm_freq = 0.0Hz,
+)
+    acq_res = acquire!(acq_plan.coarse_plan, signal, prns; interm_freq)
+    map(acq_res, prns) do res, prn
+        acquire!(
+            acq_plan.fine_plan,
+            signal,
+            prn;
+            interm_freq,
+            doppler_offset = res.carrier_doppler,
+            noise_power = res.noise_power,
+        )
+    end
+end
+
+"""
+$(SIGNATURES)
+Perform the aquisition of a single satellite `prn` in system `system` with signal `signal`
+sampled at rate `sampling_freq`. Optional arguments are the intermediate frequency `interm_freq`
+(default 0Hz), the maximum expected Doppler `max_doppler` (default 7000Hz). If the maximum Doppler
+is too unspecific you can instead pass a Doppler range with with your individual step size using
+the argument `dopplers`.
+"""
+function acquire(
+    system::AbstractGNSS,
+    signal,
+    sampling_freq,
+    prn::Integer;
+    interm_freq = 0.0Hz,
+    max_doppler = 7000Hz,
+    dopplers = -max_doppler:1/3/(length(signal)/sampling_freq):max_doppler,
+)
+    only(acquire(system, signal, sampling_freq, [prn]; interm_freq, dopplers))
+end
+
+function acquire!(
+    acq_plan::AcquisitionPlan,
+    signal,
+    prn;
+    interm_freq = 0.0Hz,
+    doppler_offset = 0.0Hz,
+    noise_power = nothing,
+)
+    only(acquire!(acq_plan, signal, [prn]; interm_freq, doppler_offset, noise_power))
+end
+
+function acquire!(acq_plan::CoarseFineAcquisitionPlan, signal, prn; interm_freq = 0.0Hz)
+    only(acquire!(acq_plan, signal, [prn]; interm_freq))
+end
+
+"""
+$(SIGNATURES)
+Performs a coarse aquisition and fine acquisition of multiple satellites `prns` in system `system` with signal `signal`
+sampled at rate `sampling_freq`. The aquisition is performed as parallel code phase
+search using the Doppler frequencies with coarse step size `coarse_step` and fine step size `fine_step`.
+"""
+function coarse_fine_acquire(
+    system::AbstractGNSS,
+    signal,
+    sampling_freq,
+    prns::AbstractVector{<:Integer};
+    interm_freq = 0.0Hz,
+    max_doppler = 7000Hz,
+    coarse_step = 1 / 3 / (length(signal) / sampling_freq),
+    fine_step = 1 / 12 / (length(signal) / sampling_freq),
+)
+    acq_plan = CoarseFineAcquisitionPlan(
+        system,
+        length(signal),
+        sampling_freq;
+        max_doppler,
+        coarse_step,
+        fine_step,
+        prns,
+        fft_flag = FFTW.ESTIMATE,
+    )
+    acquire!(acq_plan, signal, prns; interm_freq)
+end
+
+"""
+$(SIGNATURES)
+Performs a coarse aquisition and fine acquisition of a single satellite with PRN `prn` in system `system` with signal `signal`
+sampled at rate `sampling_freq`. The aquisition is performed as parallel code phase
+search using the Doppler frequencies with coarse step size `coarse_step` and fine step size `fine_step`.
+"""
+function coarse_fine_acquire(
+    system::AbstractGNSS,
+    signal,
+    sampling_freq,
+    prn::Integer;
+    interm_freq = 0.0Hz,
+    max_doppler = 7000Hz,
+    coarse_step = 1 / 3 / (length(signal) / sampling_freq),
+    fine_step = 1 / 12 / (length(signal) / sampling_freq),
+)
+    only(
+        coarse_fine_acquire(
+            system,
+            signal,
+            sampling_freq,
+            [prn];
+            interm_freq,
+            max_doppler,
+            coarse_step,
+            fine_step,
+        ),
+    )
+end