diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 78d4f3e7..b10e4835 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -2,27 +2,40 @@ env: SECRET_CODECOV_TOKEN: "Q3fuMdJjaQy9h/uk43rwSqz8M6ulvlCedU2Ir0S3QLP4t9F8cf7pzrTkX+nVhkGycZ/r5FRtTOwPr445R3wK5v9mEAsJN5GMOgI5w/L8m2XDwLmW3PN8RMno+fm2JVxZyPMNNmIQqbYEmmQcBS6Q3nywW3xi0Cl5umJuwDB+NdOFbpq3wc2wrnbOAbwlBJoCJmlH+F4ncuVY6EMmsgNKAf9RqUNWQxIthG616X1cNwuYEpL4dO/PWY2GMXWXTQ8ndO/713p4b5yIlzDP0mr2MrO+1A5fhgPc7Vr+f9mUlIAx+9AsWQYPrqPTkr2L5+mfaTodVE3u2Cop877WJZQD7w==;U2FsdGVkX1/wk2jzfWlRZ66IWgionQK/5Fu0pg3u0b26hhmmMjAjOklyi7QZKhJHjjt4KjK/dJzhd3eK28S0qQ==" steps: - - label: "Julia v1.10 - GPU Tests" + - label: "Julia v1.10 - KA + CUDA Tests" plugins: - JuliaCI/julia#v1: version: "1.10" commands: - - julia --project=ext/TrackingCUDAExt/test -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - - julia --project=ext/TrackingCUDAExt/test ext/TrackingCUDAExt/test/runtests.jl + - julia --project -e 'using Pkg; Pkg.instantiate(); Pkg.add("CUDA")' + - julia --project test/ka_gpu_tests.jl agents: queue: "juliagpu" cuda: "*" if: build.message !~ /\[skip tests\]/ timeout_in_minutes: 60 - - label: "Julia v1.11 - GPU Tests" + - label: "Julia latest - KA + CUDA Tests" plugins: - JuliaCI/julia#v1: - version: "1.11" + version: "1" commands: - - julia --project=ext/TrackingCUDAExt/test -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - - julia --project=ext/TrackingCUDAExt/test ext/TrackingCUDAExt/test/runtests.jl + - julia --project -e 'using Pkg; Pkg.instantiate(); Pkg.add("CUDA")' + - julia --project test/ka_gpu_tests.jl agents: queue: "juliagpu" cuda: "*" if: build.message !~ /\[skip tests\]/ - timeout_in_minutes: 60 \ No newline at end of file + timeout_in_minutes: 60 + - label: "Benchmark: KA-CUDA vs CPU" + plugins: + - JuliaCI/julia#v1: + version: "1" + commands: + - julia --project=benchmark -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + - JULIA_NUM_THREADS=auto julia --project=benchmark benchmark/run_multi_sat_benchmarks.jl 2>&1 | tee benchmark_results.txt + - BODY=$$(awk '/^## Multi-Satellite/,0' benchmark_results.txt); if [ -n "$$BODY" ]; then echo "$$BODY" | buildkite-agent annotate --style "info" --context "benchmarks"; fi + agents: + queue: "juliagpu" + cuda: "*" + if: build.message !~ /\[skip tests\]/ + timeout_in_minutes: 60 diff --git a/.gitignore b/.gitignore index 2f01efcb..6f602baa 100644 --- a/.gitignore +++ b/.gitignore @@ -6,8 +6,10 @@ test/Manifest.toml ext/**/test/Manifest.toml docs/build/ +docs/plans/ docs/pdf/build/ docs/site/ docs/Manifest.toml +docs/plans .worktrees/ diff --git a/Project.toml b/Project.toml index 37f5729a..98b0edf6 100644 --- a/Project.toml +++ b/Project.toml @@ -10,35 +10,39 @@ Bumper = "8ce10254-0962-460f-a3d8-1f77fea1446e" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +FastSinCos = "eb71ef10-bbd6-487f-97cc-2efbd12ad3c4" GNSSSignals = "52c80523-2a4e-5c38-8979-05588f836870" +KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" +SIMD = "fdea26ae-647d-5447-a871-4b548cad5224" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" TrackingLoopFilters = "0814aff9-93cb-554c-9fff-9bf1cfdb5efa" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" - -[weakdeps] -Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" - -[extensions] -TrackingCUDAExt = ["Adapt", "CUDA"] +VectorizationBase = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" [compat] Accessors = "0.1" Acquisition = "0.1.1, 0.2, 0.3, 1" -LinearAlgebra = "1" -Adapt = "4.3" Bumper = "0.7.1" -CUDA = "5.0" ConstructionBase = "1.6.0" Dictionaries = "0.4.1" DocStringExtensions = "0.6, 0.7, 0.8, 0.9" +FastSinCos = "1" GNSSSignals = "0.17.3, 1" -LoopVectorization = "0.8, 0.9, 0.10, 0.11, 0.12" +KernelAbstractions = "0.9" +LinearAlgebra = "1" +SIMD = "3" StaticArrays = "0.9, 0.10, 0.11, 0.12, 1.0" -StructArrays = "0.4, 0.6.5, 0.7" TrackingLoopFilters = "0.2, 1" Unitful = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 1.0" +VectorizationBase = "0.21.72" +AMDGPU = "1" +CUDA = "5" julia = "1.10" + +[extras] +AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + +[targets] +test = ["AMDGPU", "CUDA"] diff --git a/ext/TrackingCUDAExt/test/Project.toml b/benchmark/Project.toml similarity index 59% rename from ext/TrackingCUDAExt/test/Project.toml rename to benchmark/Project.toml index 7aa5f3a7..49c7ad77 100644 --- a/ext/TrackingCUDAExt/test/Project.toml +++ b/benchmark/Project.toml @@ -1,8 +1,7 @@ [deps] -Bumper = "8ce10254-0962-460f-a3d8-1f77fea1446e" +AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" GNSSSignals = "52c80523-2a4e-5c38-8979-05588f836870" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Tracking = "10b2438b-ffd4-5096-aa58-44041d5c8f3b" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" diff --git a/benchmark/bench_cpu.jl b/benchmark/bench_cpu.jl new file mode 100644 index 00000000..923c3feb --- /dev/null +++ b/benchmark/bench_cpu.jl @@ -0,0 +1,139 @@ +# CPU downconvert-and-correlate + track benchmarks + +using BenchmarkTools +using GNSSSignals +using Unitful: Hz +using Tracking +using Tracking: + CPUDownconvertAndCorrelator, + EarlyPromptLateCorrelator, + NumAnts, + SystemSatsState, + SatState, + TrackState, + downconvert_and_correlate, + get_correlator_sample_shifts, + get_code_type, + gen_code_replica! +using StaticArrays + +# ── Helper: set up common benchmark state ────────────────────────────────── + +function setup_benchmark(; + signal_type = Float32, + num_samples = 2000, + sampling_frequency = 5e6Hz, + system = GPSL1(), + num_ants = 1, +) + code_phase = 10.5 + carrier_doppler = 1000.0Hz + code_doppler = carrier_doppler * GNSSSignals.get_code_center_frequency_ratio(system) + code_frequency = code_doppler + get_code_frequency(system) + + correlator = EarlyPromptLateCorrelator(; num_ants = NumAnts(num_ants)) + static_shifts = get_correlator_sample_shifts(correlator, sampling_frequency, code_frequency) + dynamic_shifts = collect(static_shifts) + + signal = num_ants == 1 ? + rand(Complex{signal_type}, num_samples) : + rand(Complex{signal_type}, num_samples, num_ants) + + code_replica = Vector{get_code_type(system)}( + undef, num_samples + maximum(static_shifts) - minimum(static_shifts), + ) + gen_code_replica!( + code_replica, system, code_frequency, sampling_frequency, + code_phase, 1, num_samples, static_shifts, 1, Val(sampling_frequency), + ) + + return (; + correlator, signal, code_replica, + static_shifts, dynamic_shifts, + sampling_frequency, carrier_doppler, code_phase, + system, num_samples, + ) +end + +# ── High-level downconvert_and_correlate (full pipeline) ─────────────────── + +function bench_downconvert_and_correlate(; + signal_type = Float32, + num_samples = 2000, + sampling_frequency = 5e6Hz, + system = GPSL1(), + num_ants = 1, +) + downconvert_and_correlator = CPUDownconvertAndCorrelator(Val(sampling_frequency)) + system_sats_state = SystemSatsState( + system, + [SatState(system, 1, 10.5, 1000.0Hz; num_ants = NumAnts(num_ants))], + ) + track_state = TrackState((system_sats_state,)) + signal = num_ants == 1 ? + rand(Complex{signal_type}, num_samples) : + rand(Complex{signal_type}, num_samples, num_ants) + + @benchmarkable Tracking.downconvert_and_correlate( + $downconvert_and_correlator, $signal, $track_state, 1, + $sampling_frequency, $(0.0Hz), + ) +end + +# ── Fused kernel microbenchmarks ─────────────────────────────────────────── + +function bench_fused_kernel(; + signal_type = Float32, + num_samples = 2000, + num_ants = 1, + shifts = :static, +) + s = setup_benchmark(; signal_type, num_samples, num_ants) + sample_shifts = shifts == :static ? s.static_shifts : s.dynamic_shifts + # Warmup to trigger compilation + Tracking.downconvert_and_correlate_fused!( + s.correlator, s.signal, s.code_replica, sample_shifts, + s.carrier_doppler, s.sampling_frequency, 0.0, 1, s.num_samples, + ) + @benchmarkable Tracking.downconvert_and_correlate_fused!( + $(s.correlator), $(s.signal), $(s.code_replica), $sample_shifts, + $(s.carrier_doppler), $(s.sampling_frequency), 0.0, 1, $(s.num_samples), + ) +end + +function cpu_suite() + suite = BenchmarkGroup() + + # Full pipeline: CPU, various signal types + foreach((Int16, Int32, Float32, Float64)) do signal_type + suite["downconvert and correlate"]["CPU"][string(signal_type)] = + bench_downconvert_and_correlate(; signal_type) + end + + # Full pipeline: multi-antenna + suite["downconvert and correlate"]["CPU"]["Float32 4ant"] = + bench_downconvert_and_correlate(; num_ants = 4) + suite["downconvert and correlate"]["CPU"]["Int16 4ant"] = + bench_downconvert_and_correlate(; signal_type = Int16, num_ants = 4) + + # Full pipeline: track + system = GPSL1() + sampling_frequency = 5e6Hz + downconvert_and_correlator = CPUDownconvertAndCorrelator(Val(sampling_frequency)) + track_state = TrackState(system, [SatState(system, 1, 0.0, 1000Hz)]) + signal = rand(ComplexF32, 2000) + suite["track"]["Float32"] = @benchmarkable track( + $signal, $track_state, $sampling_frequency; + downconvert_and_correlator = $downconvert_and_correlator, + ) + + # Fused kernel microbenchmarks (only available on branches with the fused kernel) + if isdefined(Tracking, :downconvert_and_correlate_fused!) + suite["fused kernel"]["1-ant static taps"] = bench_fused_kernel(; shifts = :static) + suite["fused kernel"]["1-ant dynamic taps"] = bench_fused_kernel(; shifts = :dynamic) + suite["fused kernel"]["4-ant static taps"] = bench_fused_kernel(; num_ants = 4, shifts = :static) + suite["fused kernel"]["4-ant dynamic taps"] = bench_fused_kernel(; num_ants = 4, shifts = :dynamic) + end + + suite +end diff --git a/benchmark/bench_multi_sat.jl b/benchmark/bench_multi_sat.jl new file mode 100644 index 00000000..eb181c74 --- /dev/null +++ b/benchmark/bench_multi_sat.jl @@ -0,0 +1,135 @@ +# Multi-satellite / multi-system downconvert-and-correlate benchmarks +# across all available backends (CPU, CPU-Threaded, KA-CUDA, KA-AMDGPU). +# +# Returns a BenchmarkGroup with entries for each available backend. + +using BenchmarkTools +using GNSSSignals: GPSL1, GalileoE1B +using Unitful: Hz +using Tracking +using Tracking: + CPUDownconvertAndCorrelator, + SystemSatsState, + SatState, + TrackState, + downconvert_and_correlate + +if !@isdefined(HAS_CUDA) + const HAS_CUDA = try + @eval using CUDA: CUDA, cu, CuArray + CUDA.functional() + catch + false + end + + const HAS_AMDGPU = try + @eval using AMDGPU: AMDGPU, ROCArray + AMDGPU.functional() + catch + false + end + + const AMD_GPU_CUS = HAS_AMDGPU ? try + AMDGPU.HIP.properties(AMDGPU.device()).multiProcessorCount + catch + 0 + end : 0 +end + +function multi_sat_suite() + suite = BenchmarkGroup() + gpsl1 = GPSL1() + gal = GalileoE1B() + + configs = [ + ("L1 4sat/5K", (gpsl1,), [4], 5e6Hz, 5000, 32, 1000.0), + ("L1 8sat/5K", (gpsl1,), [8], 5e6Hz, 5000, 32, 1000.0), + ("E1B 4sat/25K", (gal,), [4], 25e6Hz, 25000, 50, 100.0), + ("E1B 8sat/50K", (gal,), [8], 25e6Hz, 50000, 50, 100.0), + ("8L1+8E1B/25K", (gpsl1, gal), [8, 8], 25e6Hz, 25000, 50, 100.0), + ("8L1+8E1B/50K", (gpsl1, gal), [8, 8], 25e6Hz, 50000, 50, 100.0), + ] + + # Limit GPU configs on low-CU GPUs to avoid freezing the display + max_gpu_configs = AMD_GPU_CUS <= 8 ? 2 : length(configs) + if AMD_GPU_CUS > 0 && max_gpu_configs < length(configs) + @info "AMD GPU has only $AMD_GPU_CUS CUs — limiting GPU benchmarks to first $max_gpu_configs configs" + end + + for (cfg_idx, (label, systems, nsats_list, sfreq, nsamp, prn_max, code_dop)) in enumerate(configs) + all_sss = [] + total_sats = 0 + for (si, sys) in enumerate(systems) + ns = nsats_list[min(si, length(nsats_list))] + pm = sys isa GPSL1 ? 32 : prn_max + cd = sys isa GPSL1 ? 1000.0 : code_dop + sats = [ + SatState(sys, mod1(i, pm), 10.5 + i * 0.1, (cd + i * 10) * Hz) for i = 1:ns + ] + push!(all_sss, SystemSatsState(sys, sats)) + total_sats += ns + end + ts = TrackState(Tuple(all_sss)) + signal_cpu = rand(ComplexF32, nsamp) + + # CPU baseline + cpu_dc = CPUDownconvertAndCorrelator(Val(sfreq)) + suite[label]["CPU"] = @benchmarkable downconvert_and_correlate( + $cpu_dc, + $signal_cpu, + $ts, + 1, + $sfreq, + $(0.0Hz), + ) + + # CPU threaded (only available on branches with the threaded correlator) + if isdefined(Tracking, :CPUThreadedDownconvertAndCorrelator) + cpu_threaded_dc = Tracking.CPUThreadedDownconvertAndCorrelator( + systems, + Val(sfreq); + max_sats = max(total_sats, 4), + max_num_samples = nsamp, + ) + suite[label]["CPU-Threaded"] = @benchmarkable downconvert_and_correlate( + $cpu_threaded_dc, + $signal_cpu, + $ts, + 1, + $sfreq, + $(0.0Hz), + ) + end + + # KA + CUDA + if HAS_CUDA && isdefined(Tracking, :KADownconvertAndCorrelator) + signal_cu = cu(signal_cpu) + ka = Tracking.KADownconvertAndCorrelator(systems, CuArray; max_sats = max(total_sats, 4)) + suite[label]["KA-CUDA"] = @benchmarkable downconvert_and_correlate( + $ka, + $signal_cu, + $ts, + 1, + $sfreq, + $(0.0Hz), + ) + end + + # KA + AMDGPU (skip heavy configs on low-CU GPUs) + if HAS_AMDGPU && isdefined(Tracking, :KADownconvertAndCorrelator) && cfg_idx <= max_gpu_configs + signal_roc = ROCArray(signal_cpu) + ka = + Tracking.KADownconvertAndCorrelator(systems, ROCArray; max_sats = max(total_sats, 4)) + suite[label]["KA-AMD"] = @benchmarkable downconvert_and_correlate( + $ka, + $signal_roc, + $ts, + 1, + $sfreq, + $(0.0Hz), + ) + end + end + + suite +end diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 9727aa8a..5fd7235b 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -1,210 +1,14 @@ using BenchmarkTools -using GNSSSignals -using GNSSSignals: GalileoE1B -using Unitful: Hz -using CUDA -using Tracking -using Tracking: EarlyPromptLateCorrelator, get_correlator_sample_shifts, get_code_type, - NumAnts, gen_code_replica!, SystemSatsState, SatState, TrackState, - downconvert_and_correlate -using StaticArrays - -function bench_downconvert_and_correlate( - type; - signal_type = Float32, - num_samples_signal = 2000, - sampling_frequency = 5e6Hz, - system = GPSL1(), -) - type == :GPU && !CUDA.functional() && return - - code_phase = 10.5 - intermediate_frequency = 0.0Hz - - maximum_expected_sampling_frequency = Val(sampling_frequency) - - system_sats_state = - PACKAGE_VERSION <= v"0.15.6" ? - SystemSatsState( - system, - [SatState(system, 1, sampling_frequency, code_phase, 1000.0Hz)], - ) : SystemSatsState(system, [SatState(system, 1, code_phase, 1000.0Hz)]) - - multiple_system_sats_state = (system_sats_state,) - - downconvert_and_correlator = - PACKAGE_VERSION <= v"0.15.4" ? - (type == :CPU ? CPUDownconvertAndCorrelator() : GPUDownconvertAndCorrelator()) : - ( - PACKAGE_VERSION <= v"0.15.5" ? - ( - type == :CPU ? - CPUDownconvertAndCorrelator( - maximum_expected_sampling_frequency, - multiple_system_sats_state, - num_samples_signal, - ) : - GPUDownconvertAndCorrelator(multiple_system_sats_state, num_samples_signal) - ) : - ( - type == :CPU ? - CPUDownconvertAndCorrelator(maximum_expected_sampling_frequency) : - GPUDownconvertAndCorrelator() - ) - ) - - array_transform = type == :CPU ? Array : cu - - track_state = - PACKAGE_VERSION <= v"0.15.4" ? - TrackState( - multiple_system_sats_state; - num_samples = num_samples_signal, - downconvert_and_correlator, - ) : - ( - PACKAGE_VERSION <= v"0.15.5" ? - TrackState( - multiple_system_sats_state; - num_samples = num_samples_signal, - downconvert_and_correlator, - maximum_expected_sampling_frequency, - ) : TrackState(multiple_system_sats_state) - ) - - signal = array_transform(rand(Complex{signal_type}, num_samples_signal)) - - type == :GPU && !CUDA.functional() && CUDA.versioninfo() - if PACKAGE_VERSION <= v"0.15.3" - system_sats_sample_params = - Tracking.init_sample_params(multiple_system_sats_state, 1) - next_system_sats_sample_params = Tracking.calc_sample_params( - multiple_system_sats_state, - system_sats_sample_params, - num_samples_signal, - sampling_frequency, - 1, - ) - - return if PACKAGE_VERSION <= v"0.15.2" - @benchmarkable Tracking.downconvert_and_correlate( - $signal, - $track_state, - $next_system_sats_sample_params, - $sampling_frequency, - $intermediate_frequency, - $num_samples_signal, - ) - else - @benchmarkable Tracking.downconvert_and_correlate( - $signal, - $track_state, - $next_system_sats_sample_params, - $sampling_frequency, - $intermediate_frequency, - $num_samples_signal, - $(Val(sampling_frequency)), - ) - end - else - preferred_num_code_blocks_to_integrate = 1 - return if PACKAGE_VERSION <= v"0.15.4" - @benchmarkable Tracking.downconvert_and_correlate( - $signal, - $track_state, - $preferred_num_code_blocks_to_integrate, - $sampling_frequency, - $intermediate_frequency, - $num_samples_signal, - $(Val(sampling_frequency)), - ) - else - return if PACKAGE_VERSION <= v"0.15.5" - @benchmarkable Tracking.downconvert_and_correlate( - $signal, - $track_state, - $preferred_num_code_blocks_to_integrate, - $sampling_frequency, - $intermediate_frequency, - $num_samples_signal, - ) - else - @benchmarkable Tracking.downconvert_and_correlate( - $downconvert_and_correlator, - $signal, - $track_state, - $preferred_num_code_blocks_to_integrate, - $sampling_frequency, - $intermediate_frequency, - ) - end - end - end -end - -function bench_track(; - signal_type, - num_samples_signal = 2000, - sampling_frequency = 5e6Hz, - system = GPSL1(), -) - track_state = - PACKAGE_VERSION <= v"0.15.4" ? - TrackState( - system, - [SatState(system, 1, sampling_frequency, 0.0, 1000Hz)]; - num_samples = num_samples_signal, - ) : - ( - PACKAGE_VERSION <= v"0.15.5" ? - TrackState( - system, - [SatState(system, 1, sampling_frequency, 0.0, 1000Hz)]; - num_samples = num_samples_signal, - maximum_expected_sampling_frequency = Val(sampling_frequency), - ) : - PACKAGE_VERSION <= v"0.15.6" ? - TrackState(system, [SatState(system, 1, sampling_frequency, 0.0, 1000Hz)]) : - TrackState(system, [SatState(system, 1, 0.0, 1000Hz)]) - ) - signal = rand(Complex{signal_type}, num_samples_signal) - return if PACKAGE_VERSION <= v"0.15.2" - @benchmarkable track($signal, $track_state, $sampling_frequency) - else - return if PACKAGE_VERSION <= v"0.15.4" - @benchmarkable track( - $signal, - $track_state, - $sampling_frequency, - maximum_expected_sampling_frequency = $(Val(sampling_frequency)), - ) - else - return if PACKAGE_VERSION <= v"0.15.5" - @benchmarkable track($signal, $track_state, $sampling_frequency) - else - downconvert_and_correlator = - CPUDownconvertAndCorrelator(Val(sampling_frequency)) - @benchmarkable track( - $signal, - $track_state, - $sampling_frequency; - downconvert_and_correlator = $downconvert_and_correlator, - ) - end - end - end -end +include("bench_cpu.jl") const SUITE = BenchmarkGroup() -foreach((Int16, Int32, Float32, Float64)) do signal_type - SUITE["downconvert and correlate"]["CPU"][string(signal_type)] = - bench_downconvert_and_correlate(:CPU; signal_type) -end -if (CUDA.functional()) - SUITE["downconvert and correlate"]["GPU"] = bench_downconvert_and_correlate(:GPU) +for (k, v) in cpu_suite() + SUITE[k] = v end -SUITE["track"]["Float32"] = bench_track(; signal_type = Float32) + +# GPU vs CPU benchmarks are run separately via Buildkite (benchmark/run_gpu_benchmarks.jl), +# not through AirspeedVelocity, since they require a GPU runner. # ── Multi-satellite benchmarks (threaded if available, CPU fallback) ────── diff --git a/benchmark/run_multi_sat_benchmarks.jl b/benchmark/run_multi_sat_benchmarks.jl new file mode 100644 index 00000000..4ac30a82 --- /dev/null +++ b/benchmark/run_multi_sat_benchmarks.jl @@ -0,0 +1,62 @@ +# Run multi-satellite benchmarks across all backends and print results as a markdown table. +# Usage: JULIA_NUM_THREADS=auto julia --project=benchmark benchmark/run_multi_sat_benchmarks.jl + +using Pkg +Pkg.instantiate() + +include("bench_multi_sat.jl") + +suite = multi_sat_suite() +results = run(suite; verbose=true, seconds=5) + +configs = sort(collect(results), by=first) + +# Find all GPU backends (exclude CPU and CPU-Threaded) +gpu_backends = String[] +has_threaded = false +for (_, backends) in configs + for (name, _) in backends + if name == "CPU-Threaded" + global has_threaded = true + elseif name != "CPU" && name ∉ gpu_backends + push!(gpu_backends, name) + end + end +end +sort!(gpu_backends) + +function to_us(t_ns) + round(t_ns / 1e3, digits=1) +end + +println("\n## Multi-Satellite Benchmark Results\n") +println("Julia threads: $(Threads.nthreads())\n") + +for gpu in gpu_backends + if has_threaded + println("| Config | $gpu (μs) | CPU-Threaded (μs) | CPU (μs) | vs CPU | vs Threaded |") + println("|--------|--------:|--------:|--------:|--------:|--------:|") + else + println("| Config | $gpu (μs) | CPU (μs) | vs CPU |") + println("|--------|--------:|--------:|--------:|") + end + for (config, backends) in configs + haskey(backends, "CPU") || continue + haskey(backends, gpu) || continue + cpu_us = to_us(median(backends["CPU"]).time) + gpu_us = to_us(median(backends[gpu]).time) + vs_cpu = round(cpu_us / gpu_us, digits=2) + if has_threaded + if haskey(backends, "CPU-Threaded") + thr_us = to_us(median(backends["CPU-Threaded"]).time) + vs_thr = round(thr_us / gpu_us, digits=2) + println("| $config | $gpu_us | $thr_us | $cpu_us | $(vs_cpu)x | $(vs_thr)x |") + else + println("| $config | $gpu_us | — | $cpu_us | $(vs_cpu)x | — |") + end + else + println("| $config | $gpu_us | $cpu_us | $(vs_cpu)x |") + end + end + println() +end diff --git a/docs/src/index.md b/docs/src/index.md index a9753543..64a0dc75 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -234,20 +234,21 @@ suit better. For that, you can pass a custom `post_corr_filter` to the `SatState ## GPU Support -To use GPU acceleration features, you need to explicitly load CUDA: +GPU acceleration is available via [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl), +which supports CUDA, AMDGPU, and other backends: ```julia using Tracking -using CUDA # Activates GPU functionality - -# Access GPU types via the extension -ext = Base.get_extension(Tracking, :TrackingCUDAExt) -gpu_correlator = ext.GPUDownconvertAndCorrelator(...) +using Tracking: KADownconvertAndCorrelator +using CUDA: CuArray # or AMDGPU: ROCArray + +gpsl1 = GPSL1() +ka_dc = KADownconvertAndCorrelator((gpsl1,), CuArray) +track_state = track(signal_gpu, track_state, sampling_frequency; + downconvert_and_correlator = ka_dc, +) ``` -Note: The GPU implementation is available for reference but is not yet as optimized as the -CPU implementation which uses SIMD vectorization via LoopVectorization.jl. - ## Q/A - Why are the correlator values zero? diff --git a/ext/TrackingCUDAExt/TrackingCUDAExt.jl b/ext/TrackingCUDAExt/TrackingCUDAExt.jl deleted file mode 100644 index c55b1d27..00000000 --- a/ext/TrackingCUDAExt/TrackingCUDAExt.jl +++ /dev/null @@ -1,28 +0,0 @@ -module TrackingCUDAExt - -using Tracking -using CUDA - -# Import dependencies needed by GPU code -using GNSSSignals -using StaticArrays -using Dictionaries -using Accessors -using DocStringExtensions -using Unitful: Hz -import Adapt - -# Import functions from Tracking that we need to extend or use -import Tracking: downconvert_and_correlate, downconvert_and_correlate!, get_num_samples, update - -# Include the GPU implementation -include("downconvert_and_correlate_gpu.jl") - -# Export GPU types from the extension -# Users access these via: Base.get_extension(Tracking, :TrackingCUDAExt).GPUDownconvertAndCorrelator -export GPUSatDownconvertAndCorrelator, - GPUSystemDownconvertAndCorrelator, - GPUDownconvertAndCorrelator, - convert_code_to_texture_memory - -end diff --git a/ext/TrackingCUDAExt/downconvert_and_correlate_gpu.jl b/ext/TrackingCUDAExt/downconvert_and_correlate_gpu.jl deleted file mode 100644 index d6fa1d8c..00000000 --- a/ext/TrackingCUDAExt/downconvert_and_correlate_gpu.jl +++ /dev/null @@ -1,345 +0,0 @@ -""" -$(SIGNATURES) - -A buffer that holds GPU buffers for necessary replicas and downconverted -signal. -""" -struct GPUSatDownconvertAndCorrelator{T} - downconverted_and_decoded_signal::CuArray{Complex{T},3} -end - -""" -$(SIGNATURES) - -Convenient constructor to initialize buffers for the GPU with the correct lengths for a single -satellite. -""" -function GPUSatDownconvertAndCorrelator( - ::Type{T}, - correlator::AbstractCorrelator{M}, - num_samples, -) where {T,M} - GPUSatDownconvertAndCorrelator( - CuArray{ComplexF32}(undef, (num_samples, M, get_num_accumulators(correlator))), - ) -end - -""" -$(SIGNATURES) - -Convenient constructor to initialize buffers for the GPU with the correct lengths for a single -satellite. This constructor uses Float32 as the sample data type. -""" -function GPUSatDownconvertAndCorrelator( - correlator::AbstractCorrelator{M}, - num_samples, -) where {M} - GPUSatDownconvertAndCorrelator(Float32, correlator, num_samples) -end - -struct GPUSystemDownconvertAndCorrelator{I,T,S<:AbstractGNSS{<:CuTexture}} - buffers::Dictionary{I,GPUSatDownconvertAndCorrelator{T}} - textured_system::S -end - -struct GPUDownconvertAndCorrelator{N,I,T<:GPUSystemDownconvertAndCorrelator{I}} <: - AbstractDownconvertAndCorrelator - buffers::MultipleSystemType{N,T} -end - -function GPUDownconvertAndCorrelator( - multiple_system_sats_state::MultipleSystemSatsState, - num_samples::Int, -) - buffers = map(multiple_system_sats_state) do system_sats_state - buffers = map(system_sats_state.states) do sat_state - GPUSatDownconvertAndCorrelator(Float32, sat_state.correlator, num_samples) - end - GPUSystemDownconvertAndCorrelator( - buffers, - convert_code_to_texture_memory(system_sats_state.system), - ) - end - GPUDownconvertAndCorrelator(buffers) -end - -function merge_sats( - downconvert_and_correlator::GPUDownconvertAndCorrelator{N,I,G}, - system_idx, - sats_state::Dictionary{I,<:SatState}, - num_samples::Int, -) where {N,I,T,G<:GPUSystemDownconvertAndCorrelator{I,T}} - system_sats_state = get_system_sats_state(multiple_system_sats_state, system_idx) - new_buffers = map(sats_state) do sat_state - GPUSatDownconvertAndCorrelator{T}( - T, - system_sats_state.system, - sat_state.correlator, - num_samples, - ) - end - @set downconvert_and_correlator.buffers[system_idx].buffers = - merge(downconvert_and_correlator.buffers[system_idx].buffers, new_buffers) -end - -function filter_out_sats( - downconvert_and_correlator::GPUDownconvertAndCorrelator, - system_idx::Union{Symbol,Integer}, - identifiers, -) - filtered_buffers = map( - last, - filter( - ((id,),) -> !in(id, identifiers), - pairs(downconvert_and_correlator.buffers[system_idx].buffers), - ), - ) - @set downconvert_and_correlator.buffers[system_idx].buffers = filtered_buffers -end - -import Adapt - -Adapt.@adapt_structure GPSL1 -Adapt.@adapt_structure GPSL5 -Adapt.@adapt_structure GalileoE1B - -recreate_system_with_texture(system::GPSL1, texture) = GPSL1(texture) -recreate_system_with_texture(system::GPSL5, texture) = GPSL5(texture) -recreate_system_with_texture(system::GalileoE1B, texture) = GalileoE1B(texture) - -function convert_code_to_texture_memory(system::S) where {S<:AbstractGNSS} - # Get only base code without secondary code, since otherwise code might be too - # large for texture memory. Texture memory has a max size of 65536 in each - # 2D dimension. GPSL5 would have a size of 102300 with secondary code. - # Without secondary code GPSL5 has a code size of 10230. - # The secondary code is multiplied in the kernel instead. - # The same goes for any subcarrier code. - codes = get_codes(system)[1:get_code_length(system), :] - recreate_system_with_texture( - system, - CuTexture( - CuTextureArray(CuArray(Float32.(codes))); - address_mode = CUDA.ADDRESS_MODE_WRAP, - interpolation = CUDA.NearestNeighbour(), - ), - ) -end - -""" -$(SIGNATURES) - -Downconvert and correlate all available satellites on the GPU. -""" -function downconvert_and_correlate( - downconvert_and_correlator::GPUDownconvertAndCorrelator, - signal, - track_state::TrackState, - preferred_num_code_blocks_to_integrate::Int, - sampling_frequency, - intermediate_frequency, -) - num_samples_signal = get_num_samples(signal) - new_multiple_system_sats_state = map( - track_state.multiple_system_sats_state, - downconvert_and_correlator.buffers, - ) do system_sats_state, system_buffers - new_sat_states = - map(system_sats_state.states, system_buffers.buffers) do sat_state, buffer - signal_samples_to_integrate, is_integration_completed = - calc_signal_samples_to_integrate( - system_sats_state.system, - sat_state.signal_start_sample, - sampling_frequency, - sat_state.code_doppler, - sat_state.code_phase, - preferred_num_code_blocks_to_integrate, - has_bit_or_secondary_code_been_found(sat_state), - num_samples_signal, - ) - if signal_samples_to_integrate == 0 - return sat_state - end - carrier_frequency = sat_state.carrier_doppler + intermediate_frequency - code_frequency = - sat_state.code_doppler + get_code_frequency(system_sats_state.system) - sample_shifts = get_correlator_sample_shifts( - sat_state.correlator, - sampling_frequency, - code_frequency, - ) - new_correlator = downconvert_and_correlate!( - system_buffers.textured_system, - signal, - sat_state.correlator, - sample_shifts, - sat_state.code_phase, - sat_state.carrier_phase, - code_frequency, - carrier_frequency, - sampling_frequency, - sat_state.signal_start_sample, - signal_samples_to_integrate, - sat_state.prn, - buffer.downconverted_and_decoded_signal, - )::typeof(sat_state.correlator) - return update( - system_sats_state.system, - sat_state, - signal_samples_to_integrate, - intermediate_frequency, - sampling_frequency, - new_correlator, - is_integration_completed, - ) - end - return SystemSatsState(system_sats_state, new_sat_states) - end - return TrackState( - track_state; - multiple_system_sats_state = new_multiple_system_sats_state, - ) -end - -""" -$(SIGNATURES) - -Get code for the GPU -""" -function get_code(system::AbstractGNSS, phase, prn) - get_code(system, get_modulation(system), phase, prn) -end - -""" -$(SIGNATURES) - -Get code for the GPU -""" -function get_code(system::AbstractGNSS, modulation::GNSSSignals.LOC, phase, prn) - # Must add 0.5 because CUDA uses nearest neighbour instead of floor. - system.codes[phase+0.5f0, prn] * get_secondary_code(system, phase) -end - -""" -$(SIGNATURES) - -Get code for the GPU -""" -function get_code(system::AbstractGNSS, modulation::GNSSSignals.BOC, phase, prn) - # Must add 0.5 because CUDA uses nearest neighbour instead of floor. - system.codes[phase+0.5f0, prn] * - get_secondary_code(system, phase) * - GNSSSignals.get_subcarrier_code(modulation, phase) -end - -""" -$(SIGNATURES) - -Simple kernel to downconvert and decode the signal on the GPU. -The final step for correlation (summing the samples) is done within -the sum call outside this kernel. - -This kernel currently assumes that we have more threads than number of samples -to process -TODO: handle more samples than number of threads available -""" -function downconvert_and_decode_prn_kernel!( - downconverted_and_decoded_signal, - signal, - system::AbstractGNSS, - prn::Int32, - correlator_sample_shifts, - num_samples::Int32, - code_frequency, - carrier_frequency, - sampling_frequency, - start_code_phase::Float32, - start_carrier_phase::Float32, - start_sample::Int32, - num_ants::NumAnts{N}, -) where {N} - sample = ((blockIdx().x - 0x1) * blockDim().x + (threadIdx().x - 0x1)) - index = sample + 0x1 - if sample < num_samples - carrier_wipe_off = cis( - -Float32(2π) * - (sample * carrier_frequency / sampling_frequency + start_carrier_phase), - ) - for sample_shift_index in eachindex(correlator_sample_shifts) - sample_shift = correlator_sample_shifts[sample_shift_index] - code = get_code( - system, - (sample + sample_shift) * code_frequency / sampling_frequency + - start_code_phase, - prn, - ) - for antenna_index = 0x1:N - @inbounds downconverted_and_decoded_signal[ - index, - antenna_index, - sample_shift_index, - ] = signal[sample+start_sample, antenna_index] * carrier_wipe_off * code - end - end - end - return -end - -""" -$(SIGNATURES) - -Downconvert and correlator all available satellites on the GPU. -""" -function downconvert_and_correlate!( - code_buffer, - signal, - correlator::AbstractCorrelator{M}, - sample_shifts, - code_phase, - carrier_phase, - code_frequency, - carrier_frequency, - sampling_frequency, - signal_start_sample, - num_samples_left, - prn, - downconverted_and_decoded_signal, -) where {M} - # Assume 1024 to be the max number of threads - # TODO: Evaluate this at run time - threads = min(num_samples_left, 1024) - blocks = cld(num_samples_left, threads) - num_correlators = size(downconverted_and_decoded_signal, 3) - @cuda threads = threads blocks = blocks downconvert_and_decode_prn_kernel!( - downconverted_and_decoded_signal, - signal, - code_buffer, - Int32(prn), - sample_shifts, - Int32(num_samples_left), - Float32(code_frequency / Hz), - Float32(carrier_frequency / Hz), - Float32(sampling_frequency / Hz), - Float32(code_phase), - Float32(carrier_phase), - Int32(signal_start_sample), - NumAnts{M}(), - ) - correlated_signal = - sum(view(downconverted_and_decoded_signal, 1:num_samples_left, :, :); dims = 1) - result = reshape(Array(correlated_signal), M, num_correlators) - gpu_add_to_accumulators(correlator, result) -end - -function gpu_add_to_accumulators(correlator::AbstractCorrelator{1}, result) - update_accumulator( - correlator, - SVector(map((a, b) -> a + b[1], get_accumulators(correlator), eachcol(result))), - ) -end - -function gpu_add_to_accumulators(correlator::AbstractCorrelator{M}, result) where {M} - update_accumulator( - correlator, - SVector(map(+, get_accumulators(correlator), eachcol(result))), - ) -end diff --git a/ext/TrackingCUDAExt/test/runtests.jl b/ext/TrackingCUDAExt/test/runtests.jl deleted file mode 100644 index 4b36dadb..00000000 --- a/ext/TrackingCUDAExt/test/runtests.jl +++ /dev/null @@ -1,217 +0,0 @@ -module TrackingCUDAExtTests - -using Test: @test, @testset, @inferred -using Unitful: Hz -using Tracking -using CUDA: CUDA, cu -using GNSSSignals: - GPSL1, gen_code, get_code_frequency, get_code_center_frequency_ratio, get_code -using Pkg -using Bumper: SlabBuffer -using Tracking: - CPUDownconvertAndCorrelator, - SystemSatsState, - SatState, - TrackState, - downconvert_and_correlate, - get_correlator, - track, - get_code_phase, - get_carrier_phase, - get_carrier_doppler, - get_code_doppler, - get_last_fully_integrated_filtered_prompt, - NumAnts, - get_default_correlator - -# Access extension types via Base.get_extension -const TrackingCUDAExt = Base.get_extension(Tracking, :TrackingCUDAExt) -const GPUDownconvertAndCorrelator = TrackingCUDAExt.GPUDownconvertAndCorrelator - -# Helper function to check CUDA.jl version -function is_cuda_below_version(version_string::String) - cuda_version = Pkg.installed()["CUDA"] - return cuda_version < VersionNumber(version_string) -end - -@testset "GPU Downconvert and Correlate" begin - !CUDA.functional() && return - - gpsl1 = GPSL1() - sampling_frequency = 5e6Hz - code_phase = 10.5 - num_samples_signal = 5000 - intermediate_frequency = 0.0Hz - - system_sats_state = SystemSatsState( - gpsl1, - [SatState(gpsl1, 1, code_phase, 1000.0Hz), SatState(gpsl1, 2, 11.0, 500.0Hz)]; - ) - multiple_system_sats_state = (system_sats_state,) - - downconvert_and_correlator = - GPUDownconvertAndCorrelator(multiple_system_sats_state, num_samples_signal) - - track_state = TrackState(multiple_system_sats_state) - - preferred_num_code_blocks_to_integrate = 1 - - signal = cu( - gen_code( - num_samples_signal, - gpsl1, - 1, - sampling_frequency, - get_code_frequency(gpsl1) + 1000Hz * get_code_center_frequency_ratio(gpsl1), - code_phase, - ) .* cis.(2π * (0:(num_samples_signal-1)) * 1000.0Hz / sampling_frequency), - ) - - next_track_state = @inferred downconvert_and_correlate( - downconvert_and_correlator, - signal, - track_state, - preferred_num_code_blocks_to_integrate, - sampling_frequency, - intermediate_frequency, - ) - - # GPU uses floating point arithmetic and might differ a little with the fixed point arithmetic - if is_cuda_below_version("5.8.0") - @test real.(get_correlator(next_track_state, 1).accumulators) ≈ [2921, 4949, 2921] - else - # Unfortuntely, the texture index rounding is off with CUDA 5.9 TODO: report - accumulator = real.(get_correlator(next_track_state, 1).accumulators) - @test 2900 < accumulator[1] < 3000 - @test 4900 < accumulator[2] < 5000 - @test 2900 < accumulator[3] < 3000 - end - - signal = cu( - gen_code( - num_samples_signal, - gpsl1, - 2, - sampling_frequency, - get_code_frequency(gpsl1) + 500Hz * get_code_center_frequency_ratio(gpsl1), - 11.0, - ) .* cis.(2π * (0:(num_samples_signal-1)) * 500.0Hz / sampling_frequency), - ) - - next_track_state = @inferred downconvert_and_correlate( - downconvert_and_correlator, - signal, - track_state, - preferred_num_code_blocks_to_integrate, - sampling_frequency, - intermediate_frequency, - ) - - # GPU uses floating point arithmetic and might differ a little with the fixed point arithmetic - if is_cuda_below_version("5.8.0") - @test real.(get_correlator(next_track_state, 2).accumulators) ≈ [2919, 4947, 2919] - else - # Unfortuntely, the texture index rounding is off with CUDA 5.9 TODO: report - accumulator = real.(get_correlator(next_track_state, 2).accumulators) - @test 2900 < accumulator[1] < 3000 - @test 4900 < accumulator[2] < 5000 - @test 2900 < accumulator[3] < 3000 - end -end - -@testset "Track multiple signals with GPU" begin - !CUDA.functional() && return - gpsl1 = GPSL1() - carrier_doppler = 200Hz - start_code_phase = 100 - code_frequency = carrier_doppler / 1540 + get_code_frequency(gpsl1) - sampling_frequency = 4e6Hz - prn = 1 - range = 0:3999 - start_carrier_phase = π / 2 - - num_samples = 4000 - num_ants = NumAnts(3) - - correlator = get_default_correlator(gpsl1, num_ants) - sat_states = - [SatState(gpsl1, 1, start_code_phase, carrier_doppler - 20Hz; num_ants, correlator)] - - system_sats_state = SystemSatsState(gpsl1, sat_states) - - # TODO: Why doesn't @inferred work here? - downconvert_and_correlator = - GPUDownconvertAndCorrelator((system_sats_state,), num_samples) - - track_state = @inferred TrackState(system_sats_state) - - signal = - cis.(2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ start_code_phase, - prn, - ) - signal_mat = cu(repeat(signal; outer = (1, 3))) - - track_state = @inferred track( - signal_mat, - track_state, - sampling_frequency; - downconvert_and_correlator, - ) - - iterations = 2000 - code_phases = zeros(iterations) - carrier_phases = zeros(iterations) - tracked_code_phases = zeros(iterations) - tracked_carrier_phases = zeros(iterations) - tracked_code_dopplers = zeros(iterations) - tracked_carrier_dopplers = zeros(iterations) - tracked_prompts = zeros(ComplexF64, iterations) - for i = 1:iterations - carrier_phase = - mod2pi( - 2π * carrier_doppler * 4000 * i / sampling_frequency + - start_carrier_phase + - π, - ) - π - code_phase = - mod(code_frequency * 4000 * i / sampling_frequency + start_code_phase, 1023) - signal = - cis.(2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase) .* - get_code.( - gpsl1, - code_frequency .* range ./ sampling_frequency .+ code_phase, - prn, - ) - signal_mat = cu(repeat(signal; outer = (1, 3))) - track_state = @inferred track( - signal_mat, - track_state, - sampling_frequency; - downconvert_and_correlator, - ) - comp_carrier_phase = - mod2pi( - 2π * carrier_doppler * 4000 * (i + 1) / sampling_frequency + - start_carrier_phase + - π, - ) - π - comp_code_phase = mod( - code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase, - 1023, - ) - tracked_code_phases[i] = get_code_phase(track_state) - tracked_carrier_phases[i] = get_carrier_phase(track_state) - tracked_carrier_dopplers[i] = get_carrier_doppler(track_state) / Hz - tracked_code_dopplers[i] = get_code_doppler(track_state) / Hz - tracked_prompts[i] = get_last_fully_integrated_filtered_prompt(track_state) - code_phases[i] = comp_code_phase - carrier_phases[i] = comp_carrier_phase - end - @test tracked_code_phases[end] ≈ code_phases[end] atol = 1e-2 - @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-5 -end - -end diff --git a/src/Tracking.jl b/src/Tracking.jl index 5f58a481..aa7bdfd0 100644 --- a/src/Tracking.jl +++ b/src/Tracking.jl @@ -1,15 +1,16 @@ module Tracking using DocStringExtensions +using FastSinCos using GNSSSignals -using LoopVectorization +using SIMD using StaticArrays -using StructArrays using TrackingLoopFilters using Dictionaries using Accessors using ConstructionBase using Bumper +using KernelAbstractions using Acquisition: AcquisitionResults using Unitful: upreferred, Hz, dBHz, ms @@ -47,6 +48,7 @@ export get_early, SatState, SystemSatsState, CPUDownconvertAndCorrelator, + CPUThreadedDownconvertAndCorrelator, ConventionalPLLAndDLL, ConventionalAssistedPLLAndDLL, DefaultPostCorrFilter, @@ -71,12 +73,11 @@ export get_early, get_codes, get_modulation, get_secondary_code, - update_accumulator + update_accumulator, + KADownconvertAndCorrelator const Maybe{T} = Union{T,Nothing} -const StructVecOrMat{T} = Union{StructVector{T},StructArray{T,2}} - """ $(SIGNATURES) @@ -159,7 +160,9 @@ end include("sample_parameters.jl") include("downconvert_and_correlate.jl") +include("downconvert_and_correlate_fused.jl") include("downconvert_and_correlate_cpu.jl") +include("downconvert_and_correlate_ka.jl") include("conventional_pll_and_dll.jl") include("tracking_state.jl") include("track.jl") diff --git a/src/carrier_replica.jl b/src/carrier_replica.jl index 5e6b19b1..d010a3e0 100644 --- a/src/carrier_replica.jl +++ b/src/carrier_replica.jl @@ -1,30 +1,6 @@ """ $(SIGNATURES) -Fixed point CPU StructArray carrier replica generation -""" -function gen_carrier_replica!( - carrier_replica::StructArray{Complex{T}}, - carrier_frequency, - sampling_frequency, - start_phase, - start_sample, - num_samples, -) where {T} - c_re = carrier_replica.re - c_im = carrier_replica.im - carrier_freq = upreferred(carrier_frequency / Hz) - sampling_freq = upreferred(sampling_frequency / Hz) - @avx for i = 0:num_samples-1 - c_im[i+start_sample], c_re[i+start_sample] = - sincos(T(2π) * (i * T(carrier_freq) / T(sampling_freq) + T(start_phase))) - end - carrier_replica -end - -""" -$(SIGNATURES) - Updates the carrier phase. """ function update_carrier_phase( diff --git a/src/correlators/correlator.jl b/src/correlators/correlator.jl index cb593d91..d72abe30 100644 --- a/src/correlators/correlator.jl +++ b/src/correlators/correlator.jl @@ -142,88 +142,3 @@ function calc_preferred_code_shift_to_sample_shift( max(1, sample_shift) end -""" -$(SIGNATURES) - -Perform a correlation for single antenna systems -""" -function correlate( - correlator::AbstractCorrelator{1}, - downconverted_signal::AbstractVector, - sample_shifts, - code, - start_sample, - num_samples, -) - a_re = zero_accumulators(get_accumulators(correlator), downconverted_signal) - a_im = zero_accumulators(get_accumulators(correlator), downconverted_signal) - d_re = downconverted_signal.re - d_im = downconverted_signal.im - latest_sample_shift = minimum(sample_shifts) - @avx for i = start_sample:num_samples+start_sample-1 - for j = 1:length(a_re) - sample_shift = sample_shifts[j] - latest_sample_shift - a_re[j] += d_re[i] * code[i+sample_shift] - a_im[j] += d_im[i] * code[i+sample_shift] - end - end - accumulators_result = complex.(a_re, a_im) - update_accumulator( - correlator, - map(+, get_accumulators(correlator), accumulators_result), - ) -end - -function zero_accumulators(accumulators::SVector, signal) - zeros(MVector{length(accumulators),real(eltype(signal))}) -end -function zero_accumulators(accumulators::Vector, signal) - zeros(real(eltype(signal)), length(accumulators)) -end - -""" -$(SIGNATURES) - -Perform a correlation for multi antenna systems -""" -function correlate( - correlator::AbstractCorrelator{M}, - downconverted_signal::AbstractMatrix, - sample_shifts, - code, - start_sample, - num_samples, -) where {M} - a_re = zero_accumulators(get_accumulators(correlator), downconverted_signal) - a_im = zero_accumulators(get_accumulators(correlator), downconverted_signal) - d_re = downconverted_signal.re - d_im = downconverted_signal.im - latest_sample_shift = minimum(sample_shifts) - @avx for i = start_sample:num_samples+start_sample-1 - for k = 1:size(a_re, 2) - for j = 1:size(a_re, 1) - shift = sample_shifts[k] - latest_sample_shift - a_re[j, k] += d_re[i, j] * code[shift+i] - a_im[j, k] += d_im[i, j] * code[shift+i] - end - end - end - update_accumulator( - correlator, - add_to_previous(get_accumulators(correlator), a_re, a_im), - ) -end - -function add_to_previous(accumulators::SVector{NC,<:SVector}, a_re, a_im) where {NC} - SVector{NC,eltype(accumulators)}(map(+, accumulators, eachcol(complex.(a_re, a_im)))) -end -function add_to_previous(accumulators::Vector{<:SVector}, a_re, a_im) - map(+, accumulators, eachcol(complex.(a_re, a_im))) -end - -function zero_accumulators(accumulators::SVector{NC,<:SVector{NA}}, signal) where {NC,NA} - zeros(MMatrix{NA,NC,real(eltype(signal))}) -end -function zero_accumulators(accumulators::Vector{<:SVector{NA}}, signal) where {NA} - zeros(real(eltype(signal)), NA, length(accumulators)) -end diff --git a/src/downconvert.jl b/src/downconvert.jl index 1a12276d..81954293 100644 --- a/src/downconvert.jl +++ b/src/downconvert.jl @@ -1,166 +1,31 @@ -""" -$(SIGNATURES) +# SIMD helpers shared by the fused downconvert+correlate kernel. -Downconvert a signal that has only a single dimension, e.g. like -a single antenna channel. -This is for the case that the signal is of type StructArray. -""" -function downconvert!( - downconverted_signal::StructArray{Complex{T},1}, - signal::StructArray{Complex{ST},1}, - carrier_replica::StructArray{Complex{T},1}, - start_sample::Integer, - num_samples::Integer, -) where {T,ST} - 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 - -""" -$(SIGNATURES) - -Downconvert a signal that has multiple dimensions, e.g. like -multiple antenna channels. -This is for the case that the signal is of type StructArray. -""" -function downconvert!( - downconverted_signal::StructArray{Complex{T},2}, - signal::StructArray{Complex{ST},2}, - carrier_replica::StructArray{Complex{T},1}, - start_sample::Integer, - num_samples::Integer, -) where {T,ST} - 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 - -""" -$(SIGNATURES) +using Base.Cartesian: @nexprs +using VectorizationBase: pick_vector_width -Downconvert a signal that has only a single dimension, e.g. like -a single antenna channel. -This is for the case that the signal is a regular complex array. -""" -function downconvert!( - downconverted_signal::StructArray{Complex{T},1}, - signal::AbstractArray{Complex{ST},1}, - carrier_replica::StructArray{Complex{T},1}, - start_sample::Integer, - num_samples::Integer, -) where {T,ST} - signal_real = reinterpret(reshape, ST, signal) - ds_re = downconverted_signal.re - ds_im = downconverted_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] = signal_real[1, i] * c_re[i] + signal_real[2, i] * c_im[i] - ds_im[i] = signal_real[2, i] * c_re[i] - signal_real[1, i] * c_im[i] - end - downconverted_signal -end +# Optimal SIMD width for type T, determined at precompile time from CPU features. +_simd_width(::Type{T}) where {T} = Int(pick_vector_width(T)) -""" -$(SIGNATURES) - -Downconvert a signal that has multiple dimensions, e.g. like -multiple antenna channels. -This is for the case that the signal is a regular complex array. -""" -function downconvert!( - downconverted_signal::StructArray{Complex{T},2}, - signal::AbstractArray{Complex{ST},2}, - carrier_replica::StructArray{Complex{T},1}, - start_sample::Integer, - num_samples::Integer, -) where {T,ST} - signal_real = reinterpret(reshape, ST, signal) - ds_re = downconverted_signal.re - ds_im = downconverted_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(ds_re, 2) - ds_re[i, j] = signal_real[1, i, j] * c_re[i] + signal_real[2, i, j] * c_im[i] - ds_im[i, j] = signal_real[2, i, j] * c_re[i] - signal_real[1, i, j] * c_im[i] - end - downconverted_signal -end - -# This function creates the carrier and downconverts -# directly. It is a little bit faster than first creating the -# replica and then downconverting. This is not in use at the moment -# to be consistent (see next function). -function downconvert!( - downconverted_signal::StructArray{Complex{T}}, - signal::StructArray{Complex{TS}}, - carrier_frequency, - sampling_frequency, - start_phase, - start_sample, - num_samples, -) where {T,TS} - ds_re = downconverted_signal.re - ds_im = downconverted_signal.im - s_re = signal.re - s_im = signal.im - carrier_freq = upreferred(carrier_frequency / Hz) - sampling_freq = upreferred(sampling_frequency / Hz) - @avx for i = start_sample:start_sample+num_samples-1 - c_im, c_re = sincos( - T(2π) * - ((i - start_sample) * T(carrier_freq) / T(sampling_freq) + T(start_phase)), - ) - ds_re[i] = s_re[i] * c_re + s_im[i] * c_im - ds_im[i] = s_im[i] * c_re - s_re[i] * c_im - end - downconverted_signal +# Convert Vec{N,S} to Vec{N,T}. No-op when S == T. +# Uses @generated to unroll element-wise conversion, avoiding ntuple heap allocation +# that occurs for large N (e.g. N=16 on AVX-512). +@inline _to_vec(::Type{SIMD.Vec{N,T}}, v::SIMD.Vec{N,T}) where {N,T} = v +@inline @generated function _to_vec(::Type{SIMD.Vec{N,T}}, v::SIMD.Vec{N,S}) where {N,T,S} + elems = [:(T(v[$k])) for k in 1:N] + :(SIMD.Vec{$N,T}(($(elems...),))) end -# Same as above but for the multiple antenna case. It is faster -# for less than 4 antennas, but slower otherwise. This is not in -# use due to this circumstances. -function downconvert!( - downconverted_signal::StructArray{Complex{T},2}, - signal::StructArray{Complex{TS},2}, - carrier_frequency, - sampling_frequency, - start_phase, - start_sample, - num_samples, -) where {T,TS} - ds_re = downconverted_signal.re - ds_im = downconverted_signal.im - s_re = signal.re - s_im = signal.im - carrier_freq = upreferred(carrier_frequency / Hz) - sampling_freq = upreferred(sampling_frequency / Hz) - @avx for i = start_sample:start_sample+num_samples-1 - c_im, c_re = sincos( - T(2π) * - ((i - start_sample) * T(carrier_freq) / T(sampling_freq) + T(start_phase)), - ) - for j = 1:size(s_re, 2) - ds_re[i, j] = s_re[i, j] * c_re + s_im[i, j] * c_im - ds_im[i, j] = s_im[i, j] * c_re - s_re[i, j] * c_im - end +# SIMD deinterleave load: load N interleaved complex pairs [re1,im1,re2,im2,...] +# and separate into (re_vec, im_vec) using shufflevector. +@inline @generated function _deinterleave_load(::Type{SIMD.Vec{N,T}}, p::Ptr{ST}, byte_offset::Int) where {N,T,ST} + re_idx = ntuple(k -> 2(k - 1), N) + im_idx = ntuple(k -> 2(k - 1) + 1, N) + quote + base = Ptr{ST}(p + byte_offset) + lo = vload(SIMD.Vec{$N,ST}, base) + hi = vload(SIMD.Vec{$N,ST}, base + $N * sizeof(ST)) + re_raw = shufflevector(lo, hi, Val{$re_idx}()) + im_raw = shufflevector(lo, hi, Val{$im_idx}()) + _to_vec(SIMD.Vec{$N,T}, re_raw), _to_vec(SIMD.Vec{$N,T}, im_raw) end - downconverted_signal end diff --git a/src/downconvert_and_correlate_cpu.jl b/src/downconvert_and_correlate_cpu.jl index e2f48369..dddcdc1b 100644 --- a/src/downconvert_and_correlate_cpu.jl +++ b/src/downconvert_and_correlate_cpu.jl @@ -1,9 +1,8 @@ """ $(SIGNATURES) -CPU-based implementation of downconversion and correlation. Uses LoopVectorization -for SIMD-optimized processing. The MESF type parameter specifies the maximum -expected sampling frequency for buffer allocation. +CPU-based implementation of downconversion and correlation. The MESF type +parameter specifies the maximum expected sampling frequency. """ struct CPUDownconvertAndCorrelator{MESF,B} <: AbstractDownconvertAndCorrelator buffer::B @@ -16,19 +15,48 @@ function CPUDownconvertAndCorrelator( CPUDownconvertAndCorrelator{MESF,B}(buffer) end -function get_downconvert_signal_buffer( - ::Type{T}, - num_samples::Int, - correlator::AbstractCorrelator{1}, -) where {T} - StructVector{Complex{T}}(undef, num_samples) +""" +$(SIGNATURES) + +Multi-threaded CPU downconvert and correlate. Takes systems at construction to +pre-allocate correctly-typed code replica buffers per satellite slot. Spawns +one `Threads.@spawn` task per satellite. Does not use Bumper. + +The `max_sats` parameter controls how many satellite slots (and buffer sets) +are pre-allocated. Each slot gets its own code replica buffer. +""" +struct CPUThreadedDownconvertAndCorrelator{MESF,CT<:Tuple} <: AbstractDownconvertAndCorrelator + max_sats::Int + max_samples::Int + system_code_types::Dict{UInt64,Int} # objectid(system) => index into code_replica_buffers tuple + code_replica_buffers::CT # Tuple of Vector{Vector{CodeType}} per system end -function get_downconvert_signal_buffer( - ::Type{T}, - num_samples::Int, - correlator::AbstractCorrelator{M}, -) where {T,M} - StructArray{Complex{T}}(undef, num_samples, M) + +function CPUThreadedDownconvertAndCorrelator( + systems, + ::Val{MESF}; + max_sats::Int = 32, + max_sample_shift::Int = 20, + max_num_samples::Int = ceil(Int, upreferred(MESF / Hz) * 1e-3), +) where {MESF} + max_samples = max_num_samples + code_len = max_samples + 2 * max_sample_shift + + system_code_types = Dict{UInt64,Int}( + objectid(sys) => i for (i, sys) in enumerate(systems) + ) + + code_replica_buffers = Tuple( + [Vector{get_code_type(sys)}(undef, code_len) for _ in 1:max_sats] + for sys in systems + ) + + CPUThreadedDownconvertAndCorrelator{MESF,typeof(code_replica_buffers)}( + max_sats, + max_samples, + system_code_types, + code_replica_buffers, + ) end """ @@ -77,27 +105,13 @@ function downconvert_and_correlate( num_samples_signal + maximum(sample_shifts) - minimum(sample_shifts) ) - carrier_replica_buffer_re = @alloc(Float32, num_samples_signal) - carrier_replica_buffer_im = @alloc(Float32, num_samples_signal) - carrier_replica_buffer = StructArray{ComplexF32}(( - carrier_replica_buffer_re, - carrier_replica_buffer_im, - )) - downconvert_signal_buffer_re = @alloc(Float32, size(signal)...) - downconvert_signal_buffer_im = @alloc(Float32, size(signal)...) - downconvert_signal_buffer = StructArray{ComplexF32}(( - downconvert_signal_buffer_re, - downconvert_signal_buffer_im, - )) new_correlator = downconvert_and_correlate!( system_sats_state.system, signal, sat_state.correlator, code_replica_buffer, sat_state.code_phase, - carrier_replica_buffer, sat_state.carrier_phase, - downconvert_signal_buffer, code_frequency, carrier_frequency, sampling_frequency, @@ -128,17 +142,118 @@ end """ $(SIGNATURES) +Multi-threaded downconvert and correlate. Spawns one task per satellite across +all systems, using pre-allocated per-satellite code replica buffers and the +fused downconvert+correlate kernel. +""" +function downconvert_and_correlate( + dc::CPUThreadedDownconvertAndCorrelator{MESF}, + signal, + track_state::TrackState, + preferred_num_code_blocks_to_integrate::Int, + sampling_frequency, + intermediate_frequency, +) where {MESF} + num_samples_signal = get_num_samples(signal) + + buf_offset = 0 + new_multiple_system_sats_state = + map(track_state.multiple_system_sats_state) do system_sats_state + system = system_sats_state.system + sys_idx = dc.system_code_types[objectid(system)] + states = system_sats_state.states + n = length(states) + new_vals = Vector{valtype(states)}(undef, n) + + Threads.@threads for i in 1:n + buf_idx = buf_offset + i + sat_state = states.values[i] + + signal_samples_to_integrate, is_integration_completed = + calc_signal_samples_to_integrate( + system, + sat_state.signal_start_sample, + sampling_frequency, + sat_state.code_doppler, + sat_state.code_phase, + preferred_num_code_blocks_to_integrate, + has_bit_or_secondary_code_been_found(sat_state), + num_samples_signal, + ) + + if signal_samples_to_integrate == 0 + new_vals[i] = sat_state + continue + end + + carrier_frequency = sat_state.carrier_doppler + intermediate_frequency + code_frequency = sat_state.code_doppler + get_code_frequency(system) + + sample_shifts = get_correlator_sample_shifts( + sat_state.correlator, + sampling_frequency, + code_frequency, + ) + + code_replica = dc.code_replica_buffers[sys_idx][buf_idx] + gen_code_replica!( + code_replica, + system, + code_frequency, + sampling_frequency, + sat_state.code_phase, + sat_state.signal_start_sample, + signal_samples_to_integrate, + sample_shifts, + sat_state.prn, + Val{MESF}(), + ) + + new_correlator = downconvert_and_correlate_fused!( + sat_state.correlator, + signal, + code_replica, + sample_shifts, + carrier_frequency, + sampling_frequency, + sat_state.carrier_phase, + sat_state.signal_start_sample, + signal_samples_to_integrate, + )::typeof(sat_state.correlator) + + new_vals[i] = update( + system, + sat_state, + signal_samples_to_integrate, + intermediate_frequency, + sampling_frequency, + new_correlator, + is_integration_completed, + ) + end + + buf_offset += n + new_sat_states = Dictionary(keys(states), new_vals) + SystemSatsState(system_sats_state, new_sat_states) + end + TrackState( + track_state; + multiple_system_sats_state = new_multiple_system_sats_state, + ) +end + +""" +$(SIGNATURES) + Downconvert and correlate a single satellite on the CPU. """ function downconvert_and_correlate!( system, signal, - correlator, + correlator::AbstractCorrelator{M}, code_replica, code_phase, - carrier_replica, carrier_phase, - downconverted_signal, code_frequency, carrier_frequency, sampling_frequency, @@ -146,7 +261,7 @@ function downconvert_and_correlate!( num_samples_left, prn, maximum_expected_sampling_frequency, -) +) where {M} sample_shifts = get_correlator_sample_shifts(correlator, sampling_frequency, code_frequency) gen_code_replica!( @@ -161,60 +276,15 @@ function downconvert_and_correlate!( prn, maximum_expected_sampling_frequency, ) - gen_carrier_replica!( - carrier_replica, + downconvert_and_correlate_fused!( + correlator, + signal, + code_replica, + sample_shifts, carrier_frequency, sampling_frequency, carrier_phase, signal_start_sample, num_samples_left, ) - downconvert!( - downconverted_signal, - signal, - carrier_replica, - signal_start_sample, - num_samples_left, - ) - correlate( - correlator, - downconverted_signal, - sample_shifts, - code_replica, - signal_start_sample, - num_samples_left, - ) -end - -#= -# This is currently slower than splitting the loop. -# See https://github.com/JuliaSIMD/LoopVectorization.jl/issues/284 -function downconvert_and_correlate( - signal::StructArray{Complex{T}}, - correlator::C, - code, - correlator_sample_shifts, - carrier_frequency, - sampling_frequency, - start_phase, - start_sample, - num_samples -) where {T, C <: AbstractCorrelator} - s_re = signal.re; s_im = signal.im - accumulators = zero_accumulators(get_accumulators(correlator), signal) - a_re = real.(accumulators) - a_im = imag.(accumulators) - @avx for i = start_sample:start_sample + num_samples - 1 - c_im, c_re = sincos(T(2π) * ((i - start_sample) * T(upreferred(carrier_frequency / Hz)) / T(upreferred(sampling_frequency / Hz)) + T(start_phase))) - d_re = s_re[i] * c_re + s_im[i] * c_im - d_im = s_im[i] * c_re - s_re[i] * c_im - for j = 1:length(a_re) - sample_shift = correlator_sample_shifts[j] - correlator_sample_shifts[1] - a_re[j] += d_re * code[i + sample_shift] - a_im[j] += d_im * code[i + sample_shift] - end - end - accumulators_result = complex.(a_re, a_im) - C(map(+, get_accumulators(correlator), accumulators_result)) end -=# diff --git a/src/downconvert_and_correlate_fused.jl b/src/downconvert_and_correlate_fused.jl new file mode 100644 index 00000000..8da4ca17 --- /dev/null +++ b/src/downconvert_and_correlate_fused.jl @@ -0,0 +1,347 @@ +# Fused downconvert + correlate SIMD kernel. +# Generates the carrier on-the-fly, downconverts, and accumulates against the +# code replica in a single pass — keeping downconverted samples in registers +# instead of writing them to an intermediate buffer. + +using Base.Cartesian: @nexprs + +# Horizontal sum: reduce SIMD vector lanes to a scalar. +@inline _hsum(v::SIMD.Vec) = sum(Tuple(v)) + +# Build SIMD offset vector for unroll index u (1-based). +# Uses @generated to produce a literal tuple, avoiding ntuple heap allocation on AVX-512. +@inline @generated function _make_offset(::Type{T}, ::Val{W}, ::Val{U}) where {T,W,U} + elems = [:(T($(k - 1 + (U - 1) * W))) for k in 1:W] + :(SIMD.Vec{$W,T}(($(elems...),))) +end + +# Stack-friendly mutable copy / immutable conversion for accumulators. +# SVector → MVector (stack-allocated), Vector → copy (heap). +@inline _mutable_copy(v::SVector) = MVector(v) +@inline _mutable_copy(v::Vector) = copy(v) +@inline _to_immutable(v::MVector) = SVector(v) +@inline _to_immutable(v::Vector) = v + +# Add a single antenna's correlation result to accumulator element. +# M==1: scalar complex; M>1: SVector with one antenna updated. +@inline _add_antenna(acc::Complex, prev::Complex, j::Int, val::ComplexF64) = prev + val +@inline function _add_antenna(acc::SVector{M}, prev::SVector{M}, j::Int, val::ComplexF64) where {M} + Base.setindex(acc, acc[j] + val, j) +end + +""" + downconvert_and_correlate_fused!( + correlator, signal, code_replica, sample_shifts, + carrier_frequency, sampling_frequency, carrier_phase, + start_sample, num_samples, + ) + +Fused downconvert-and-correlate: generates the carrier replica on-the-fly, +downconverts the signal, and immediately accumulates the result against +shifted code replicas. Returns an updated correlator of the same type. + +Both antennas (M) and taps (NC) are fully unrolled at compile time via +`@generated`, so all accumulators live in named local variables (registers). +""" +@generated function downconvert_and_correlate_fused!( + correlator::AbstractCorrelator{M}, + signal::AbstractArray{Complex{ST}}, + code_replica, + sample_shifts::SVector{NC}, + carrier_frequency, + sampling_frequency, + carrier_phase, + start_sample::Integer, + num_samples::Integer, +) where {M,ST,NC} + W = _simd_width(Float32) # Compute at generation time, embed as literal + # ── Accumulator init: acc_re_j_k, acc_im_j_k (SIMD), s_re_j_k, s_im_j_k (scalar) ── + acc_init = Expr(:block) + for j in 1:M, k in 1:NC + push!(acc_init.args, :($(Symbol("acc_re_$(j)_$(k)")) = z)) + push!(acc_init.args, :($(Symbol("acc_im_$(j)_$(k)")) = z)) + push!(acc_init.args, :($(Symbol("s_re_$(j)_$(k)")) = 0.0)) + push!(acc_init.args, :($(Symbol("s_im_$(j)_$(k)")) = 0.0)) + end + + # ── Precompute per-tap shifted code pointers ── + shift_init = Expr(:block) + for k in 1:NC + push!(shift_init.args, :($(Symbol("p_code_$(k)")) = p_code + (sample_shifts[$k] - min_shift) * sizeof_CT)) + end + + # ── Main 4x-unrolled accumulate block ── + main_accum = Expr(:block) + for j in 1:M + push!(main_accum.args, :(col_byte_off = $(j - 1) * sig_col_bytes + row_byte_off)) + push!(main_accum.args, :(@nexprs 4 u -> (sr_u, si_u) = _deinterleave_load( + SIMD.Vec{W,T}, p_sig, + col_byte_off + (u - 1) * W * 2 * sizeof_ST, + ))) + push!(main_accum.args, :(@nexprs 4 u -> dre_u = sr_u * cr_u + si_u * ci_u)) + push!(main_accum.args, :(@nexprs 4 u -> dim_u = si_u * cr_u - sr_u * ci_u)) + for k in 1:NC + are = Symbol("acc_re_$(j)_$(k)") + aim = Symbol("acc_im_$(j)_$(k)") + pck = Symbol("p_code_$(k)") + push!(main_accum.args, quote + @nexprs 4 u -> begin + code_u = vload(SIMD.Vec{W,CT}, $pck + (i - 1 + (u - 1) * W) * sizeof_CT) + code_f_u = _to_vec(SIMD.Vec{W,T}, code_u) + $are = muladd(dre_u, code_f_u, $are) + $aim = muladd(dim_u, code_f_u, $aim) + end + end) + end + end + + # ── 1x cleanup accumulate block ── + cleanup_accum = Expr(:block) + for j in 1:M + sr_j = Symbol("sr_c_$(j)") + si_j = Symbol("si_c_$(j)") + dre_j = Symbol("dre_c_$(j)") + dim_j = Symbol("dim_c_$(j)") + push!(cleanup_accum.args, :(($sr_j, $si_j) = _deinterleave_load( + SIMD.Vec{W,T}, p_sig, + $(j - 1) * sig_col_bytes + row_byte_off, + ))) + push!(cleanup_accum.args, :($dre_j = $sr_j * c_re + $si_j * c_im)) + push!(cleanup_accum.args, :($dim_j = $si_j * c_re - $sr_j * c_im)) + for k in 1:NC + are = Symbol("acc_re_$(j)_$(k)") + aim = Symbol("acc_im_$(j)_$(k)") + cf_jk = Symbol("code_f_c_$(j)_$(k)") + pck = Symbol("p_code_$(k)") + push!(cleanup_accum.args, :($cf_jk = _to_vec(SIMD.Vec{W,T}, vload(SIMD.Vec{W,CT}, $pck + (i - 1) * sizeof_CT)))) + push!(cleanup_accum.args, :($are = muladd($dre_j, $cf_jk, $are))) + push!(cleanup_accum.args, :($aim = muladd($dim_j, $cf_jk, $aim))) + end + end + + # ── Scalar remainder accumulate block ── + scalar_accum = Expr(:block) + for j in 1:M + sig_j = Symbol("sig_s_$(j)") + sr_j = Symbol("sr_s_$(j)") + si_j = Symbol("si_s_$(j)") + dre_j = Symbol("dre_s_$(j)") + dim_j = Symbol("dim_s_$(j)") + push!(scalar_accum.args, :($sig_j = signal[i, $j])) + push!(scalar_accum.args, :($sr_j = T(real($sig_j)))) + push!(scalar_accum.args, :($si_j = T(imag($sig_j)))) + push!(scalar_accum.args, :($dre_j = $sr_j * c_re_s + $si_j * c_im_s)) + push!(scalar_accum.args, :($dim_j = $si_j * c_re_s - $sr_j * c_im_s)) + for k in 1:NC + sre = Symbol("s_re_$(j)_$(k)") + sim = Symbol("s_im_$(j)_$(k)") + cv_jk = Symbol("cv_s_$(j)_$(k)") + pck = Symbol("p_code_$(k)") + push!(scalar_accum.args, :($cv_jk = T(unsafe_load(Ptr{CT}($pck + (i - 1) * sizeof_CT))))) + push!(scalar_accum.args, :($sre += Float64($dre_j * $cv_jk))) + push!(scalar_accum.args, :($sim += Float64($dim_j * $cv_jk))) + end + end + + # ── Final reduction: build result tuple directly ── + if M == 1 + tap_exprs = [quote + complex( + Float64(_hsum($(Symbol("acc_re_1_$(k)")))) + $(Symbol("s_re_1_$(k)")), + Float64(_hsum($(Symbol("acc_im_1_$(k)")))) + $(Symbol("s_im_1_$(k)")), + ) + end for k in 1:NC] + else + tap_exprs = [begin + ant_exprs = [quote + complex( + Float64(_hsum($(Symbol("acc_re_$(j)_$(k)")))) + $(Symbol("s_re_$(j)_$(k)")), + Float64(_hsum($(Symbol("acc_im_$(j)_$(k)")))) + $(Symbol("s_im_$(j)_$(k)")), + ) + end for j in 1:M] + :(SVector(tuple($(ant_exprs...)))) + end for k in 1:NC] + end + + result_expr = quote + accumulators_result = SVector(tuple($(tap_exprs...))) + prev = get_accumulators(correlator) + new_acc = map(+, prev, accumulators_result) + update_accumulator(correlator, new_acc) + end + + quote + T = Float32 + CT = eltype(code_replica) + sizeof_CT = sizeof(CT) + sizeof_ST = sizeof(ST) + W = $W + min_shift = minimum(sample_shifts) + + carrier_freq = T(upreferred(carrier_frequency / Hz)) + sampling_freq = T(upreferred(sampling_frequency / Hz)) + phase0 = T(carrier_phase) + two_pi = T(2π) + freq_ratio = carrier_freq / sampling_freq + + num_samples_signal = size(signal, 1) + p_sig = Ptr{ST}(pointer(signal)) + sig_col_bytes = num_samples_signal * 2 * sizeof_ST + p_code = Ptr{CT}(pointer(code_replica)) + + z = zero(SIMD.Vec{W,T}) + $acc_init + $shift_init + + off_1 = _make_offset(T, Val{W}(), Val{1}()) + off_2 = _make_offset(T, Val{W}(), Val{2}()) + off_3 = _make_offset(T, Val{W}(), Val{3}()) + off_4 = _make_offset(T, Val{W}(), Val{4}()) + + # Phase via muladd: p_u = base_phase + init_u where + # base_phase = 2π·freq_ratio·(i - start_sample) (scalar, broadcast) + # init_u = 2π·(off_u·freq_ratio + phase0) (precomputed per-vector) + # One muladd per vector, no drift, no per-iteration multiply chain. + two_pi_fr = SIMD.Vec{W,T}(two_pi * freq_ratio) + init_1 = two_pi * (off_1 * freq_ratio + phase0) + init_2 = two_pi * (off_2 * freq_ratio + phase0) + init_3 = two_pi * (off_3 * freq_ratio + phase0) + init_4 = two_pi * (off_4 * freq_ratio + phase0) + + i = start_sample + last = start_sample + num_samples - 1 + + @inbounds while i + 4W - 1 <= last + base_phase = SIMD.Vec{W,T}(T(i - start_sample)) + p_1 = muladd(base_phase, two_pi_fr, init_1) + p_2 = muladd(base_phase, two_pi_fr, init_2) + p_3 = muladd(base_phase, two_pi_fr, init_3) + p_4 = muladd(base_phase, two_pi_fr, init_4) + ci_1, cr_1 = fast_sincos_u100k(p_1) + ci_2, cr_2 = fast_sincos_u100k(p_2) + ci_3, cr_3 = fast_sincos_u100k(p_3) + ci_4, cr_4 = fast_sincos_u100k(p_4) + row_byte_off = (i - 1) * 2 * sizeof_ST + $main_accum + i += 4W + end + + @inbounds while i + W - 1 <= last + base_phase = SIMD.Vec{W,T}(T(i - start_sample)) + phase = muladd(base_phase, two_pi_fr, init_1) + c_im, c_re = fast_sincos_u100k(phase) + row_byte_off = (i - 1) * 2 * sizeof_ST + $cleanup_accum + i += W + end + + @inbounds while i <= last + ph = two_pi * (T(i - start_sample) * freq_ratio + phase0) + c_im_s, c_re_s = sincos(ph) + $scalar_accum + i += 1 + end + + $result_expr + end +end + + +# ── Fallback for dynamic-length sample_shifts (AbstractVector) ──────── +# Fused downconvert + correlate using stack-allocated SoA tile buffers +# (via @no_escape) and @simd for the tap accumulation loop. +# No dependency on @avx / LoopVectorization. +function downconvert_and_correlate_fused!( + correlator::AbstractCorrelator{M}, + signal::AbstractArray{Complex{ST}}, + code_replica, + sample_shifts::AbstractVector, + carrier_frequency, + sampling_frequency, + carrier_phase, + start_sample::Integer, + num_samples::Integer, +) where {M,ST} + T = Float32 + sizeof_ST = sizeof(ST) + W = _simd_width(Float32) + num_taps = length(sample_shifts) + min_shift = minimum(sample_shifts) + + carrier_freq = T(upreferred(carrier_frequency / Hz)) + phase0 = T(carrier_phase) + two_pi = T(2π) + freq_ratio = carrier_freq / T(upreferred(sampling_frequency / Hz)) + + num_samples_signal = size(signal, 1) + p_sig = Ptr{ST}(pointer(signal)) + sig_col_bytes = num_samples_signal * 2 * sizeof_ST + + off_1 = _make_offset(T, Val{W}(), Val{1}()) + two_pi_fr = SIMD.Vec{W,T}(two_pi * freq_ratio) + init_1 = two_pi * (off_1 * freq_ratio + phase0) + + last = start_sample + num_samples - 1 + + @no_escape begin + # Flat buffer: M antennas × num_samples, SoA layout + tile_re = @alloc(T, num_samples * M) + tile_im = @alloc(T, num_samples * M) + + # Downconvert each antenna into its tile slice + i = start_sample + idx = 1 + @inbounds while i + W - 1 <= last + base_phase = SIMD.Vec{W,T}(T(i - start_sample)) + phase = muladd(base_phase, two_pi_fr, init_1) + ci, cr = fast_sincos_u100k(phase) + row_byte_off = (i - 1) * 2 * sizeof_ST + for j in 1:M + sr, si = _deinterleave_load( + SIMD.Vec{W,T}, p_sig, + (j - 1) * sig_col_bytes + row_byte_off, + ) + dre = sr * cr + si * ci + dim = si * cr - sr * ci + off = ((j - 1) * num_samples + idx - 1) * sizeof(T) + vstore(dre, pointer(tile_re) + off) + vstore(dim, pointer(tile_im) + off) + end + i += W + idx += W + end + @inbounds while i <= last + ph = two_pi * (T(i - start_sample) * freq_ratio + phase0) + c_im_s, c_re_s = sincos(ph) + for j in 1:M + sig = signal[i, j] + tile_re[(j-1)*num_samples + idx] = T(real(sig)) * c_re_s + T(imag(sig)) * c_im_s + tile_im[(j-1)*num_samples + idx] = T(imag(sig)) * c_re_s - T(real(sig)) * c_im_s + end + i += 1 + idx += 1 + end + + # Correlate: tap-outer, antenna-inner with @simd + prev = get_accumulators(correlator) + new_acc = _mutable_copy(prev) + @inbounds for k in 1:num_taps + shift_offset = sample_shifts[k] - min_shift + for j in 1:M + acc_r = zero(T) + acc_i = zero(T) + ant_off = (j - 1) * num_samples + @simd for n in 1:num_samples + c = code_replica[n + shift_offset] + acc_r += tile_re[ant_off + n] * c + acc_i += tile_im[ant_off + n] * c + end + corr_val = complex(Float64(acc_r), Float64(acc_i)) + new_acc[k] = _add_antenna(new_acc[k], prev[k], j, corr_val) + end + end + + update_accumulator(correlator, _to_immutable(new_acc)) + end +end + diff --git a/src/downconvert_and_correlate_ka.jl b/src/downconvert_and_correlate_ka.jl new file mode 100644 index 00000000..43ec9fd4 --- /dev/null +++ b/src/downconvert_and_correlate_ka.jl @@ -0,0 +1,523 @@ +using KernelAbstractions: @kernel, @index, get_backend, @Const, @uniform, @localmem, @synchronize, @groupsize, @private, synchronize +using GNSSSignals: LOC, BOC, CBOC, BOCsin, BOCcos, get_subcarrier_code, get_modulation + +""" + KADownconvertAndCorrelator + +GPU-accelerated downconvert-and-correlate using KernelAbstractions. + +Expanded code tables (with subcarrier baked in) live permanently on GPU. +Only a tiny param buffer is copied per call. + +Uses Int64 fixed-point arithmetic with 32 fractional bits for essentially +exact subcarrier indexing (zero quantization errors across 100K samples). +""" +struct KADownconvertAndCorrelator{CT,AR<:AbstractArray{ComplexF64,3},AP<:AbstractMatrix{Float32}} <: + AbstractDownconvertAndCorrelator + expanded_codes_gpu::CT # Tuple of GPU Float32 matrices: (code_length*sub_per_chip) × num_prns + sub_per_chip::Vector{Int32} # sub-phases per chip for each system (1 for LOC, 12 for CBOC, etc.) + system_indices::Dict{UInt64,Int32} # objectid(system) → system_idx in tuple + results_gpu::AR + results_cpu::Array{ComplexF64,3} + sat_params_gpu::AP # header_rows × max_sats Float32 (tiny param buffer) + sat_params_cpu::Matrix{Float32} + num_threads::Int + max_sats::Int + max_taps::Int + header_rows::Int +end + +""" + _get_sub_per_chip(modulation) + +Number of subcarrier sub-phases per code chip. +LOC = 1, BOCsin(m,n) = 2m, BOCcos(m,n) = 2m, CBOC = lcm(2m1, 2m2). +""" +_get_sub_per_chip(::LOC) = 1 +_get_sub_per_chip(mod::BOCsin) = 2 * mod.m +_get_sub_per_chip(mod::BOCcos) = 2 * mod.m +function _get_sub_per_chip(mod::CBOC) + lcm(2 * mod.boc1.m, 2 * mod.boc2.m) +end + +""" + _build_expanded_code_table(codes, modulation, code_length) + +Build expanded code table with subcarrier baked in. +Entry [(i-1)*sub_per_chip + j, prn] = codes[i, prn] * subcarrier((j-1)/sub_per_chip). +For LOC (sub_per_chip=1), this is just Float32.(codes). +""" +function _build_expanded_code_table(codes::AbstractMatrix, modulation::LOC, code_length::Int) + Float32.(codes) +end + +function _build_expanded_code_table(codes::AbstractMatrix, modulation, code_length::Int) + spc = _get_sub_per_chip(modulation) + num_prns = size(codes, 2) + expanded = Matrix{Float32}(undef, code_length * spc, num_prns) + for prn in 1:num_prns + for i in 1:code_length + chip_val = Float32(codes[i, prn]) + for j in 1:spc + sc_phase = (j - 1) / spc + sc_val = Float32(get_subcarrier_code(modulation, sc_phase)) + expanded[(i - 1) * spc + j, prn] = chip_val * sc_val + end + end + end + expanded +end + +""" + KADownconvertAndCorrelator(systems, ArrayType; kwargs...) + +Construct a GPU downconvert-and-correlator using Int64 fixed-point +arithmetic with 32 fractional bits for exact subcarrier indexing. +""" +function KADownconvertAndCorrelator( + systems, + ArrayType::Type{<:AbstractArray}; + num_threads::Int = 256, + max_ants::Int = 1, + max_taps::Int = 5, + max_sats::Int = 32, +) + expanded_codes_gpu = Tuple( + let mod = get_modulation(typeof(sys)), + codes = get_codes(sys), + code_length = get_code_length(sys) + ArrayType{Float32}(_build_expanded_code_table(codes, mod, code_length)) + end + for sys in systems + ) + sub_per_chip_vec = Int32[_get_sub_per_chip(get_modulation(typeof(sys))) for sys in systems] + system_indices = Dict{UInt64,Int32}( + objectid(sys) => Int32(i) for (i, sys) in enumerate(systems) + ) + + # Header layout per column (all Float32; integer values bit-packed via reinterpret): + # + # Row 1: carrier_freq_ratio Float32 + # Row 2: carrier_phase Float32 + # Row 3: num_samples Int32 reinterpreted + # Row 4: start_sample Int32 reinterpreted + # Row 5: prn Int32 reinterpreted + # Rows 6-7: code_phase_fixed Int64 as 2×Float32 + # Rows 8-9: delta_code_phase_fixed Int64 as 2×Float32 + # Rows 10..9+2*max_taps: tap_shift_fixed Int64 as 2×Float32 each + header_rows = 9 + 2 * max_taps + + results_gpu = ArrayType{ComplexF64}(undef, max_sats, max_ants, max_taps) + results_cpu = Array{ComplexF64}(undef, max_sats, max_ants, max_taps) + sat_params_gpu = ArrayType{Float32}(undef, header_rows, max_sats) + sat_params_cpu = Matrix{Float32}(undef, header_rows, max_sats) + KADownconvertAndCorrelator{typeof(expanded_codes_gpu),typeof(results_gpu),typeof(sat_params_gpu)}( + expanded_codes_gpu, + sub_per_chip_vec, + system_indices, + results_gpu, + results_cpu, + sat_params_gpu, + sat_params_cpu, + num_threads, + max_sats, + max_taps, + header_rows, + ) +end + +# Helper: store Int32 as Float32 bits (bit-exact round-trip) +@inline _f32_from_i32(x::Int32) = reinterpret(Float32, x) +# Helper: recover Int32 from Float32 bits +@inline _i32_from_f32(x::Float32) = reinterpret(Int32, x) + +# Helper: store Int64 as two Float32 values (lo, hi halves via UInt32) +@inline function _f32_pair_from_i64(x::Int64) + bits = reinterpret(UInt64, x) + lo = reinterpret(Float32, UInt32(bits & 0xFFFFFFFF)) + hi = reinterpret(Float32, UInt32((bits >> 32) & 0xFFFFFFFF)) + lo, hi +end +# Helper: recover Int64 from two Float32 values +@inline function _i64_from_f32_pair(lo::Float32, hi::Float32) + lo_bits = UInt64(reinterpret(UInt32, lo)) + hi_bits = UInt64(reinterpret(UInt32, hi)) + reinterpret(Int64, lo_bits | (hi_bits << 32)) +end + +# Fixed-point constants: 32 fractional bits +const FIXED_POINT_BITS = 32 +const FIXED_POINT_SCALE = Int64(1) << FIXED_POINT_BITS # 4294967296 + +# ============================================================ +# Unified kernel: handles single and multi-antenna via num_ants +# ============================================================ +@kernel function ka_dc_kernel!( + results, + @Const(signal), + @Const(expanded_codes), + ::Val{num_taps}, + ::Val{combined}, + expanded_code_length::Int32, + num_ants::Int32, + @Const(sat_params), + result_offset::Int32, +) where {num_taps, combined} + N = @uniform @groupsize()[1] + local_tid = @index(Local, Linear) + group_idx = @index(Group, Linear) + + shmem_re = @localmem Float64 (combined ? num_taps : 1, 256) + shmem_im = @localmem Float64 (combined ? num_taps : 1, 256) + sat_fp = @localmem Float32 (2,) + sat_ip = @localmem Int32 (3,) + sat_phase = @localmem Int64 (2,) + sat_sh = @localmem Int64 (num_taps,) + + if local_tid == 1 + col = result_offset + Int32(group_idx) + @inbounds sat_fp[1] = sat_params[1, col] + @inbounds sat_fp[2] = sat_params[2, col] + @inbounds sat_ip[1] = _i32_from_f32(sat_params[3, col]) + @inbounds sat_ip[2] = _i32_from_f32(sat_params[4, col]) + @inbounds sat_ip[3] = _i32_from_f32(sat_params[5, col]) + @inbounds sat_phase[1] = _i64_from_f32_pair(sat_params[6, col], sat_params[7, col]) + @inbounds sat_phase[2] = _i64_from_f32_pair(sat_params[8, col], sat_params[9, col]) + @inbounds for t in 1:num_taps + sat_sh[t] = _i64_from_f32_pair(sat_params[9 + 2*(t-1) + 1, col], sat_params[9 + 2*(t-1) + 2, col]) + end + end + @synchronize + + acc_re = @private Float32 (num_taps,) + acc_im = @private Float32 (num_taps,) + + for ant in Int32(1):num_ants + for tap in 1:num_taps + @inbounds acc_re[tap] = zero(Float32) + @inbounds acc_im[tap] = zero(Float32) + end + + @inbounds begin + initial_code_fixed = sat_phase[1] + delta_code_fixed = sat_phase[2] + prn = sat_ip[3] + + ecl_shifted = Int64(expanded_code_length) << Int64(FIXED_POINT_BITS) + delta_stride = delta_code_fixed * Int64(N) + + code_phase = mod((Int64(local_tid) - Int64(1)) * delta_code_fixed + initial_code_fixed, ecl_shifted) + + @fastmath begin + step_phase = Float32(6.283185307179586) * sat_fp[1] * Float32(N) + s_step, c_step = sincos(step_phase) + init_phase = Float32(6.283185307179586) * (Float32(Int32(local_tid) - Int32(1)) * sat_fp[1] + sat_fp[2]) + s, c = sincos(init_phase) + end + + sample = Int32(local_tid) + while sample <= sat_ip[1] + sig = signal[sample + sat_ip[2] - Int32(1), ant] + sig_re = real(sig) + sig_im = imag(sig) + dc_re = sig_re * c + sig_im * s + dc_im = sig_im * c - sig_re * s + + for tap in 1:num_taps + tap_phase = code_phase + sat_sh[tap] + tap_phase -= (tap_phase >= ecl_shifted) * ecl_shifted + tap_phase += (tap_phase < Int64(0)) * ecl_shifted + idx = Int32(tap_phase >> Int64(FIXED_POINT_BITS)) + Int32(1) + code_val = expanded_codes[idx, prn] + acc_re[tap] += dc_re * code_val + acc_im[tap] += dc_im * code_val + end + + code_phase += delta_stride + code_phase -= (code_phase >= ecl_shifted) * ecl_shifted + + @fastmath begin + c_new = c * c_step - s * s_step + s_new = s * c_step + c * s_step + c = c_new + s = s_new + end + sample += Int32(N) + end + end + + # Reduction + if combined + @inbounds for tap in 1:num_taps + shmem_re[tap, local_tid] = Float64(acc_re[tap]) + shmem_im[tap, local_tid] = Float64(acc_im[tap]) + end + for half in (128, 64, 32, 16, 8, 4, 2) + @synchronize + if local_tid <= half + @inbounds for tap in 1:num_taps + shmem_re[tap, local_tid] += shmem_re[tap, local_tid + half] + shmem_im[tap, local_tid] += shmem_im[tap, local_tid + half] + end + end + end + @synchronize + if local_tid == 1 + col = result_offset + Int32(group_idx) + @inbounds for tap in 1:num_taps + results[col, ant, tap] = complex( + shmem_re[tap, 1] + shmem_re[tap, 2], + shmem_im[tap, 1] + shmem_im[tap, 2], + ) + end + end + else + for tap in 1:num_taps + @inbounds shmem_re[1, local_tid] = Float64(acc_re[tap]) + @inbounds shmem_im[1, local_tid] = Float64(acc_im[tap]) + for half in (128, 64, 32, 16, 8, 4, 2) + @synchronize + if local_tid <= half + @inbounds shmem_re[1, local_tid] += shmem_re[1, local_tid + half] + @inbounds shmem_im[1, local_tid] += shmem_im[1, local_tid + half] + end + end + @synchronize + if local_tid == 1 + col = result_offset + Int32(group_idx) + @inbounds results[col, ant, tap] = complex( + shmem_re[1, 1] + shmem_re[1, 2], + shmem_im[1, 1] + shmem_im[1, 2], + ) + end + @synchronize + end + end + end +end + +# ============================================================ +# CPU-side result update helpers +# ============================================================ +function ka_update_from_results( + correlator::AbstractCorrelator{1}, + results, + batch_idx, +) + accumulators = get_accumulators(correlator) + new_accums = map(enumerate(accumulators)) do (tap, prev) + @inbounds prev + results[batch_idx, 1, tap] + end + update_accumulator(correlator, new_accums) +end + +function ka_update_from_results( + correlator::AbstractCorrelator{M}, + results, + batch_idx, +) where {M} + accumulators = get_accumulators(correlator) + new_accums = map(enumerate(accumulators)) do (tap, prev) + ant_vals = ntuple(M) do ant + @inbounds results[batch_idx, ant, tap] + end + prev + SVector{M}(ant_vals) + end + update_accumulator(correlator, new_accums) +end + +# ============================================================ +# Param packing (Int64 expanded units) +# ============================================================ +function _pack_sat_params!( + sat_params_cpu, batch_idx, + carrier_freq_ratio, carrier_phase, + signal_samples_to_integrate, signal_start_sample, prn, + code_freq_ratio, norm_phase, code_length, spc, + sample_shifts, num_taps, +) + expanded_freq_ratio = code_freq_ratio * spc + expanded_norm_phase = norm_phase * spc + + delta_code_fixed = round(Int64, expanded_freq_ratio * FIXED_POINT_SCALE) + code_phase_fixed = round(Int64, expanded_norm_phase * FIXED_POINT_SCALE) + + @inbounds sat_params_cpu[1, batch_idx] = carrier_freq_ratio + @inbounds sat_params_cpu[2, batch_idx] = carrier_phase + @inbounds sat_params_cpu[3, batch_idx] = _f32_from_i32(Int32(signal_samples_to_integrate)) + @inbounds sat_params_cpu[4, batch_idx] = _f32_from_i32(Int32(signal_start_sample)) + @inbounds sat_params_cpu[5, batch_idx] = _f32_from_i32(Int32(prn)) + + lo, hi = _f32_pair_from_i64(code_phase_fixed) + @inbounds sat_params_cpu[6, batch_idx] = lo + @inbounds sat_params_cpu[7, batch_idx] = hi + lo, hi = _f32_pair_from_i64(delta_code_fixed) + @inbounds sat_params_cpu[8, batch_idx] = lo + @inbounds sat_params_cpu[9, batch_idx] = hi + + for t in 1:num_taps + tap_shift_expanded = Float64(sample_shifts[t]) * expanded_freq_ratio + tap_shift_fixed = round(Int64, tap_shift_expanded * FIXED_POINT_SCALE) + lo, hi = _f32_pair_from_i64(tap_shift_fixed) + @inbounds sat_params_cpu[9 + 2*(t-1) + 1, batch_idx] = lo + @inbounds sat_params_cpu[9 + 2*(t-1) + 2, batch_idx] = hi + end +end + +# ============================================================ +# Main dispatch: downconvert_and_correlate +# ============================================================ +function downconvert_and_correlate( + dc::KADownconvertAndCorrelator, + signal, + track_state::TrackState, + preferred_num_code_blocks_to_integrate::Int, + sampling_frequency, + intermediate_frequency, +) + num_samples_signal = get_num_samples(signal) + backend = get_backend(first(dc.expanded_codes_gpu)) + + batch_idx = 0 + M = get_num_ants(track_state) + header_rows = dc.header_rows + + system_data = map(track_state.multiple_system_sats_state) do system_sats_state + system = system_sats_state.system + sys_idx = dc.system_indices[objectid(system)] + states = system_sats_state.states + code_length = get_code_length(system) + spc = dc.sub_per_chip[sys_idx] + + batch_offset = batch_idx + + sat_infos = map(states) do sat_state + signal_samples_to_integrate, is_integration_completed = + calc_signal_samples_to_integrate( + system, + sat_state.signal_start_sample, + sampling_frequency, + sat_state.code_doppler, + sat_state.code_phase, + preferred_num_code_blocks_to_integrate, + has_bit_or_secondary_code_been_found(sat_state), + num_samples_signal, + ) + carrier_frequency = sat_state.carrier_doppler + intermediate_frequency + code_frequency = sat_state.code_doppler + get_code_frequency(system) + sample_shifts = get_correlator_sample_shifts( + sat_state.correlator, + sampling_frequency, + code_frequency, + ) + (; + signal_samples_to_integrate, + is_integration_completed, + carrier_frequency, + code_frequency, + sample_shifts, + ) + end + + key_to_batch = Dict{eltype(keys(states)),Int}() + system_num_taps = 0 + + for key in keys(states) + info = sat_infos[key] + info.signal_samples_to_integrate == 0 && continue + + batch_idx += 1 + key_to_batch[key] = batch_idx + sat_state = states[key] + num_taps = length(info.sample_shifts) + system_num_taps = num_taps + + code_freq_ratio = Float64(upreferred(info.code_frequency / Hz) / upreferred(sampling_frequency / Hz)) + norm_phase = mod(Float64(sat_state.code_phase), code_length) + + carrier_freq_ratio = Float32( + upreferred(info.carrier_frequency / Hz) / + upreferred(sampling_frequency / Hz), + ) + + _pack_sat_params!( + dc.sat_params_cpu, batch_idx, + carrier_freq_ratio, Float32(sat_state.carrier_phase), + info.signal_samples_to_integrate, sat_state.signal_start_sample, sat_state.prn, + code_freq_ratio, norm_phase, code_length, spc, + info.sample_shifts, num_taps, + ) + end + + num_active = batch_idx - batch_offset + + (; system_sats_state, system, sys_idx, sat_infos, key_to_batch, batch_offset, num_active, system_num_taps, + code_length, spc) + end + + num_total_active = batch_idx + + if num_total_active == 0 + return track_state + end + + # Tiny copy to GPU + n = num_total_active + elems = header_rows * n + copyto!(vec(dc.sat_params_gpu), 1, vec(dc.sat_params_cpu), 1, elems) + + num_threads = dc.num_threads + + # Per-system kernel launches + for sd in system_data + sd.num_active == 0 && continue + + codes = dc.expanded_codes_gpu[sd.sys_idx] + combined = Val(sd.system_num_taps <= 8) + expanded_code_length = Int32(sd.code_length * sd.spc) + + ka_dc_kernel!(backend, (num_threads,))( + dc.results_gpu, signal, codes, + Val(sd.system_num_taps), combined, + expanded_code_length, Int32(M), + dc.sat_params_gpu, Int32(sd.batch_offset); + ndrange = num_threads * sd.num_active, + ) + end + + # Synchronize + copy results + synchronize(backend) + copyto!(dc.results_cpu, dc.results_gpu) + + # Distribute results back + new_multiple_system_sats_state = map(system_data) do sd + states = sd.system_sats_state.states + new_sat_states = map(pairs(states)) do (key, sat_state) + batch_pos = get(sd.key_to_batch, key, nothing) + if batch_pos === nothing + return sat_state + end + info = sd.sat_infos[key] + new_correlator = ka_update_from_results( + sat_state.correlator, + dc.results_cpu, + batch_pos, + )::typeof(sat_state.correlator) + return update( + sd.system, + sat_state, + info.signal_samples_to_integrate, + intermediate_frequency, + sampling_frequency, + new_correlator, + info.is_integration_completed, + ) + end + return SystemSatsState(sd.system_sats_state, new_sat_states) + end + + return TrackState( + track_state; + multiple_system_sats_state = new_multiple_system_sats_state, + ) +end diff --git a/src/sat_state.jl b/src/sat_state.jl index 1f155a3f..465218af 100644 --- a/src/sat_state.jl +++ b/src/sat_state.jl @@ -28,6 +28,7 @@ Get the PRN (Pseudo-Random Noise) number of the satellite. """ get_prn(s::SatState) = s.prn get_num_ants(s::SatState{<:AbstractCorrelator{M}}) where {M} = M +get_num_ants(::Type{<:SatState{<:AbstractCorrelator{M}}}) where {M} = M """ $(SIGNATURES) @@ -213,6 +214,8 @@ struct SystemSatsState{S<:AbstractGNSS,SS<:SatState,I} states::Dictionary{I,SS} end +get_num_ants(::SystemSatsState{<:Any,SS}) where {SS} = get_num_ants(SS) + """ Type alias for a tuple or named tuple of `SystemSatsState` objects, representing tracking state across multiple GNSS systems. diff --git a/src/tracking_state.jl b/src/tracking_state.jl index a1a2e8a5..f41871f3 100644 --- a/src/tracking_state.jl +++ b/src/tracking_state.jl @@ -210,6 +210,7 @@ end # Convenient methods get_prn(s::TrackState, id...) = get_prn(get_sat_state(s, id...)) +get_num_ants(s::TrackState) = get_num_ants(first(s.multiple_system_sats_state)) get_num_ants(s::TrackState, id...) = get_num_ants(get_sat_state(s, id...)) get_code_phase(s::TrackState, id...) = get_code_phase(get_sat_state(s, id...)) get_code_doppler(s::TrackState, id...) = get_code_doppler(get_sat_state(s, id...)) diff --git a/test/carrier_replica.jl b/test/carrier_replica.jl index 7efc555b..18f1b923 100644 --- a/test/carrier_replica.jl +++ b/test/carrier_replica.jl @@ -2,26 +2,7 @@ module CarrierReplicaTest using Test: @test, @testset, @inferred using Unitful: Hz -using StructArrays: StructArray -using Tracking: gen_carrier_replica!, update_carrier_phase - -@testset "Carrier replica" begin - carrier = StructArray(zeros(Complex{Float32}, 2500)) - - gen_carrier_replica!( - carrier, - 1500Hz, - 2.5e6Hz, - 0.25, # π / 2 - 111, - 2390, - ) - - @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 +using Tracking: update_carrier_phase @testset "Update carrier phase" begin carrier_phase = 0.25 diff --git a/test/cn0_estimation.jl b/test/cn0_estimation.jl index d070c627..60241702 100644 --- a/test/cn0_estimation.jl +++ b/test/cn0_estimation.jl @@ -4,15 +4,14 @@ using Test: @test, @testset, @inferred using Random: Random using Unitful: kHz, MHz, Hz, ms, dBHz using StaticArrays: SVector -using StructArrays: StructArray using GNSSSignals: GPSL1, get_code +import Tracking using Tracking: MomentsCN0Estimator, get_prompt_buffer, get_current_index, update, EarlyPromptLateCorrelator, - correlate, get_prompt, get_correlator_sample_shifts, estimate_cn0, @@ -72,12 +71,14 @@ end correlator = EarlyPromptLateCorrelator() sample_shifts = get_correlator_sample_shifts(correlator, sampling_frequency, code_frequency) - signal_struct = StructArray(signal) - correlator = correlate( + correlator = Tracking.downconvert_and_correlate_fused!( correlator, - signal_struct, - sample_shifts, + signal, code, + sample_shifts, + 0.0Hz, + sampling_frequency, + 0.0, start_sample, num_samples, ) diff --git a/test/correlator.jl b/test/correlator.jl index bdc76dcd..aacf56a8 100644 --- a/test/correlator.jl +++ b/test/correlator.jl @@ -3,14 +3,12 @@ module CorrelatorTest using Test: @test, @testset, @inferred using Unitful: Hz using GNSSSignals: GPSL1, GalileoE1B, get_code, get_code_frequency -using StaticArrays: SVector, MVector, MMatrix -using StructArrays: StructArray +using StaticArrays: SVector using Tracking: EarlyPromptLateCorrelator, VeryEarlyPromptLateCorrelator, NumAnts, NumAccumulators, - add_to_previous, get_very_early, get_early, get_prompt, @@ -22,10 +20,8 @@ using Tracking: DefaultPostCorrFilter, apply, normalize, - correlate, get_correlator_sample_shifts, - get_initial_accumulator, - zero_accumulators + get_initial_accumulator @testset "Correlator" begin @testset "Get initial accumulators" begin @@ -38,66 +34,6 @@ using Tracking: Vector{SVector{4,ComplexF64}} end - @testset "Zero accumulators" begin - accumulator = - @inferred(zero_accumulators(ones(SVector{3,ComplexF64}), zeros(ComplexF32, 1))) - @test accumulator isa MVector{3,Float32} - @test accumulator == [0, 0, 0] - - accumulator = @inferred( - zero_accumulators( - SVector(ones(SVector{3,ComplexF64}), ones(SVector{3,ComplexF64})), - zeros(ComplexF32, 1), - ) - ) - @test accumulator isa MMatrix{3,2,Float32} - @test accumulator == [0.0 0.0; 0.0 0.0; 0.0 0.0] - - accumulator = - @inferred(zero_accumulators(ones(ComplexF64, 3), zeros(ComplexF32, 1))) - @test accumulator isa Vector{Float32} - @test accumulator == [0, 0, 0] - - accumulator = @inferred( - zero_accumulators( - [ones(SVector{3,ComplexF64}), ones(SVector{3,ComplexF64})], - zeros(ComplexF32, 1), - ) - ) - @test accumulator isa Matrix{Float32} - @test accumulator == [0.0 0.0; 0.0 0.0; 0.0 0.0] - end - - @testset "Add to previous accumulators" begin - accumulator = @inferred( - add_to_previous( - get_initial_accumulator(NumAnts(2), NumAccumulators(3)), - ones(MMatrix{2,3,Float32}), - ones(MMatrix{2,3,Float32}), - ) - ) - @test accumulator isa SVector{3,SVector{2,ComplexF64}} - @test accumulator == [ - [1.0 + 1.0im, 1.0 + 1.0im], - [1.0 + 1.0im, 1.0 + 1.0im], - [1.0 + 1.0im, 1.0 + 1.0im], - ] - - accumulator = @inferred( - add_to_previous( - get_initial_accumulator(NumAnts(2), 3), - ones(Float32, 2, 3), - ones(Float32, 2, 3), - ) - ) - @test accumulator isa Vector{SVector{2,ComplexF64}} - @test accumulator == [ - [1.0 + 1.0im, 1.0 + 1.0im], - [1.0 + 1.0im, 1.0 + 1.0im], - [1.0 + 1.0im, 1.0 + 1.0im], - ] - end - @testset "Correlator constructor" begin gpsl1 = GPSL1() sampling_frequency = 5e6Hz @@ -280,36 +216,6 @@ using Tracking: ) end - @testset "Correlate" begin - gpsl1 = GPSL1() - sampling_frequency = 2.5e6Hz - signal = StructArray{Complex{Float32}}(( - Float32.(get_code.(gpsl1, (1:2500) * 1023e3 / (sampling_frequency / Hz), 1)), - zeros(Float32, 2500), - )) - code = get_code.(gpsl1, ((1+-1):(2500+1)) * 1023e3 / 2.5e6, 1) - - @testset "↳ $na antennas" for na ∈ [1, 4] - signal_mat = na == 1 ? signal : repeat(signal; outer = (1, na)) - - correlator = EarlyPromptLateCorrelator(; num_ants = NumAnts(na)) - - sample_shifts = get_correlator_sample_shifts( - correlator, - sampling_frequency, - get_code_frequency(gpsl1), - ) - - correlator_result = - @inferred correlate(correlator, signal_mat, sample_shifts, code, 1, 2500) - early = get_early(correlator_result) - prompt = get_prompt(correlator_result) - late = get_late(correlator_result) - @test early == late - @test all(early .== 1476) - @test all(prompt .== 2500) - end - end end end diff --git a/test/downconvert_and_correlate.jl b/test/downconvert_and_correlate.jl index f08a239d..5f63f845 100644 --- a/test/downconvert_and_correlate.jl +++ b/test/downconvert_and_correlate.jl @@ -2,15 +2,25 @@ module DownconvertAndCorrelateTest using Test: @test, @testset, @inferred using Unitful: Hz -using GNSSSignals: GPSL1, gen_code, get_code_frequency, get_code_center_frequency_ratio +using GNSSSignals: + GPSL1, GalileoE1B, gen_code, get_code_frequency, get_code_center_frequency_ratio, get_code_type using Bumper: SlabBuffer +import Tracking using Tracking: + AbstractCorrelator, CPUDownconvertAndCorrelator, + CPUThreadedDownconvertAndCorrelator, + EarlyPromptLateCorrelator, + KADownconvertAndCorrelator, + NumAnts, SystemSatsState, SatState, TrackState, downconvert_and_correlate, - get_correlator + get_accumulators, + get_correlator, + get_correlator_sample_shifts, + update_accumulator @testset "Downconvert and Correlator" begin downconvert_and_correlator = CPUDownconvertAndCorrelator(Val(5e6Hz)) @@ -58,7 +68,7 @@ end intermediate_frequency, ) - @test real.(get_correlator(next_track_state, 1).accumulators) ≈ [2921, 4949, 2917] + @test real.(get_correlator(next_track_state, 1).accumulators) ≈ [2921, 4949, 2917] rtol=1e-3 signal = gen_code( @@ -79,7 +89,641 @@ end intermediate_frequency, ) - @test real.(get_correlator(next_track_state, 2).accumulators) ≈ [2919, 4947, 2915] + @test real.(get_correlator(next_track_state, 2).accumulators) ≈ [2919, 4947, 2915] rtol=1e-3 +end + +@testset "Fused downconvert and correlate" begin + gpsl1 = GPSL1() + sampling_frequency = 5e6Hz + code_phase = 10.5 + num_samples_signal = 5000 + carrier_doppler = 1000.0Hz + prn = 1 + + carrier_frequency = carrier_doppler + 0.0Hz + code_doppler = carrier_doppler * get_code_center_frequency_ratio(gpsl1) + code_frequency = code_doppler + get_code_frequency(gpsl1) + + correlator = EarlyPromptLateCorrelator() + sample_shifts = get_correlator_sample_shifts(correlator, sampling_frequency, code_frequency) + + signal = + gen_code( + num_samples_signal, + gpsl1, + prn, + sampling_frequency, + code_frequency, + code_phase, + ) .* cis.(2π * (0:(num_samples_signal - 1)) * carrier_doppler / sampling_frequency) + + code_replica_length = num_samples_signal + maximum(sample_shifts) - minimum(sample_shifts) + code_replica = Vector{get_code_type(gpsl1)}(undef, code_replica_length) + Tracking.gen_code_replica!( + code_replica, gpsl1, code_frequency, sampling_frequency, + code_phase, 1, num_samples_signal, sample_shifts, prn, + Val(sampling_frequency), + ) + + # Scalar reference: plain Julia loop, no SIMD, no tricks + NC = length(sample_shifts) + min_shift = minimum(sample_shifts) + ref_accumulators = zeros(ComplexF64, NC) + carrier_step = 2π * Float64(carrier_frequency / sampling_frequency) + for i in 1:num_samples_signal + carrier = cis(-carrier_step * (i - 1)) + dc_sample = signal[i] * carrier + for tap in 1:NC + code_idx = i + sample_shifts[tap] - min_shift + ref_accumulators[tap] += dc_sample * code_replica[code_idx] + end + end + + # Static (@generated) path + fused_static = Tracking.downconvert_and_correlate_fused!( + correlator, signal, code_replica, sample_shifts, + carrier_frequency, sampling_frequency, 0.0, 1, num_samples_signal, + ) + + @test get_accumulators(fused_static) ≈ ref_accumulators rtol = 1e-4 + + # Dynamic (AbstractVector) path + dynamic_shifts = collect(sample_shifts) + fused_dynamic = Tracking.downconvert_and_correlate_fused!( + correlator, signal, code_replica, dynamic_shifts, + carrier_frequency, sampling_frequency, 0.0, 1, num_samples_signal, + ) + + @test get_accumulators(fused_dynamic) ≈ ref_accumulators rtol = 1e-4 + + # Static and dynamic paths should produce nearly identical results + @test get_accumulators(fused_static) ≈ get_accumulators(fused_dynamic) rtol = 1e-4 +end + +@testset "Fused downconvert and correlate multi-antenna" begin + using StaticArrays: SVector + + gpsl1 = GPSL1() + sampling_frequency = 5e6Hz + code_phase = 10.5 + num_samples_signal = 5000 + carrier_doppler = 1000.0Hz + prn = 1 + num_ants = 2 + + carrier_frequency = carrier_doppler + 0.0Hz + code_doppler = carrier_doppler * get_code_center_frequency_ratio(gpsl1) + code_frequency = code_doppler + get_code_frequency(gpsl1) + + correlator = EarlyPromptLateCorrelator(; num_ants = NumAnts(num_ants)) + sample_shifts = get_correlator_sample_shifts(correlator, sampling_frequency, code_frequency) + + # Multi-antenna signal: each antenna gets slightly different data + signal_ant1 = + gen_code( + num_samples_signal, gpsl1, prn, sampling_frequency, + code_frequency, code_phase, + ) .* cis.(2π * (0:(num_samples_signal - 1)) * carrier_doppler / sampling_frequency) + signal_ant2 = signal_ant1 .* cis(0.3) # phase-shifted copy + signal = hcat(signal_ant1, signal_ant2) + + code_replica_length = num_samples_signal + maximum(sample_shifts) - minimum(sample_shifts) + code_replica = Vector{get_code_type(gpsl1)}(undef, code_replica_length) + Tracking.gen_code_replica!( + code_replica, gpsl1, code_frequency, sampling_frequency, + code_phase, 1, num_samples_signal, sample_shifts, prn, + Val(sampling_frequency), + ) + + # Scalar reference for multi-antenna + NC = length(sample_shifts) + min_shift = minimum(sample_shifts) + ref_accumulators = [zeros(ComplexF64, num_ants) for _ in 1:NC] + carrier_step = 2π * Float64(carrier_frequency / sampling_frequency) + for i in 1:num_samples_signal + carrier = cis(-carrier_step * (i - 1)) + for j in 1:num_ants + dc_sample = signal[i, j] * carrier + for tap in 1:NC + code_idx = i + sample_shifts[tap] - min_shift + ref_accumulators[tap][j] += dc_sample * code_replica[code_idx] + end + end + end + + # Static (@generated) path + fused_static = Tracking.downconvert_and_correlate_fused!( + correlator, signal, code_replica, sample_shifts, + carrier_frequency, sampling_frequency, 0.0, 1, num_samples_signal, + ) + for tap in 1:NC + @test collect(get_accumulators(fused_static)[tap]) ≈ ref_accumulators[tap] rtol = 1e-4 + end + + # Dynamic (AbstractVector) path + dynamic_shifts = collect(sample_shifts) + fused_dynamic = Tracking.downconvert_and_correlate_fused!( + correlator, signal, code_replica, dynamic_shifts, + carrier_frequency, sampling_frequency, 0.0, 1, num_samples_signal, + ) + for tap in 1:NC + @test collect(get_accumulators(fused_dynamic)[tap]) ≈ ref_accumulators[tap] rtol = 1e-4 + end + + # Static and dynamic should match + for k in 1:NC + @test collect(get_accumulators(fused_static)[k]) ≈ collect(get_accumulators(fused_dynamic)[k]) rtol = 1e-4 + end + + # Dynamic path with non-SIMD-aligned signal length (covers scalar remainder loop) + num_samples_odd = 5003 # not a multiple of SIMD width + signal_odd_ant1 = + gen_code( + num_samples_odd, gpsl1, prn, sampling_frequency, + code_frequency, code_phase, + ) .* cis.(2π * (0:(num_samples_odd - 1)) * carrier_doppler / sampling_frequency) + signal_odd = hcat(signal_odd_ant1, signal_odd_ant1 .* cis(0.3)) + code_replica_odd = Vector{get_code_type(gpsl1)}(undef, num_samples_odd + maximum(sample_shifts) - minimum(sample_shifts)) + Tracking.gen_code_replica!( + code_replica_odd, gpsl1, code_frequency, sampling_frequency, + code_phase, 1, num_samples_odd, sample_shifts, prn, + Val(sampling_frequency), + ) + ref_odd = [zeros(ComplexF64, num_ants) for _ in 1:NC] + for i in 1:num_samples_odd + carrier = cis(-carrier_step * (i - 1)) + for j in 1:num_ants + dc_sample = signal_odd[i, j] * carrier + for tap in 1:NC + ref_odd[tap][j] += dc_sample * code_replica_odd[i + sample_shifts[tap] - min_shift] + end + end + end + fused_odd = Tracking.downconvert_and_correlate_fused!( + correlator, signal_odd, code_replica_odd, dynamic_shifts, + carrier_frequency, sampling_frequency, 0.0, 1, num_samples_odd, + ) + for tap in 1:NC + @test collect(get_accumulators(fused_odd)[tap]) ≈ ref_odd[tap] rtol = 1e-4 + end +end + +@testset "Fused downconvert and correlate with Vector accumulators" begin + # Minimal correlator with Vector (not SVector) accumulators + struct DynamicCorrelator <: AbstractCorrelator{1} + accumulators::Vector{ComplexF64} + shifts::Vector{Int} + end + Tracking.get_accumulators(c::DynamicCorrelator) = c.accumulators + Tracking.update_accumulator(c::DynamicCorrelator, acc) = + DynamicCorrelator(collect(acc), c.shifts) + + gpsl1 = GPSL1() + sampling_frequency = 5e6Hz + code_phase = 10.5 + num_samples_signal = 5000 + carrier_doppler = 1000.0Hz + prn = 1 + + carrier_frequency = carrier_doppler + 0.0Hz + code_doppler = carrier_doppler * get_code_center_frequency_ratio(gpsl1) + code_frequency = code_doppler + get_code_frequency(gpsl1) + + # Use the same shifts as EarlyPromptLateCorrelator but as a Vector + epl = EarlyPromptLateCorrelator() + static_shifts = get_correlator_sample_shifts(epl, sampling_frequency, code_frequency) + dynamic_shifts = collect(static_shifts) + + dynamic_correlator = DynamicCorrelator(zeros(ComplexF64, 3), dynamic_shifts) + + signal = + gen_code( + num_samples_signal, + gpsl1, + prn, + sampling_frequency, + code_frequency, + code_phase, + ) .* cis.(2π * (0:(num_samples_signal - 1)) * carrier_doppler / sampling_frequency) + + code_replica_length = num_samples_signal + maximum(dynamic_shifts) - minimum(dynamic_shifts) + code_replica = Vector{get_code_type(gpsl1)}(undef, code_replica_length) + Tracking.gen_code_replica!( + code_replica, gpsl1, code_frequency, sampling_frequency, + code_phase, 1, num_samples_signal, static_shifts, prn, + Val(sampling_frequency), + ) + + # Scalar reference + NC = length(dynamic_shifts) + min_shift = minimum(dynamic_shifts) + ref_accumulators = zeros(ComplexF64, NC) + carrier_step = 2π * Float64(carrier_frequency / sampling_frequency) + for i in 1:num_samples_signal + carrier = cis(-carrier_step * (i - 1)) + dc_sample = signal[i] * carrier + for tap in 1:NC + code_idx = i + dynamic_shifts[tap] - min_shift + ref_accumulators[tap] += dc_sample * code_replica[code_idx] + end + end + + # Test: use the dynamic Vector accumulator path + fused_dynamic = Tracking.downconvert_and_correlate_fused!( + dynamic_correlator, signal, code_replica, dynamic_shifts, + carrier_frequency, sampling_frequency, 0.0, 1, num_samples_signal, + ) + + @test get_accumulators(fused_dynamic) ≈ ref_accumulators rtol = 1e-4 +end + +@testset "Downconvert and correlate with CPUThreaded" begin + gpsl1 = GPSL1() + sampling_frequency = 5e6Hz + code_phase = 10.5 + num_samples_signal = 5000 + intermediate_frequency = 0.0Hz + + system_sats_state = SystemSatsState( + gpsl1, + [SatState(gpsl1, 1, code_phase, 1000.0Hz), SatState(gpsl1, 2, 11.0, 500.0Hz)]; + ) + multiple_system_sats_state = (system_sats_state,) + + downconvert_and_correlator = + CPUThreadedDownconvertAndCorrelator((gpsl1,), Val(sampling_frequency)) + + track_state = TrackState(multiple_system_sats_state) + + preferred_num_code_blocks_to_integrate = 1 + + signal = + gen_code( + num_samples_signal, + gpsl1, + 1, + sampling_frequency, + get_code_frequency(gpsl1) + 1000Hz * get_code_center_frequency_ratio(gpsl1), + code_phase, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 1000.0Hz / sampling_frequency) + + next_track_state = downconvert_and_correlate( + downconvert_and_correlator, + signal, + track_state, + preferred_num_code_blocks_to_integrate, + sampling_frequency, + intermediate_frequency, + ) + + @test real.(get_correlator(next_track_state, 1).accumulators) ≈ [2921, 4949, 2917] rtol=1e-3 + + signal = + gen_code( + num_samples_signal, + gpsl1, + 2, + sampling_frequency, + get_code_frequency(gpsl1) + 500Hz * get_code_center_frequency_ratio(gpsl1), + 11.0, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 500.0Hz / sampling_frequency) + + next_track_state = downconvert_and_correlate( + downconvert_and_correlator, + signal, + track_state, + preferred_num_code_blocks_to_integrate, + sampling_frequency, + intermediate_frequency, + ) + + @test real.(get_correlator(next_track_state, 2).accumulators) ≈ [2919, 4947, 2915] rtol=1e-3 + + # Test early return when signal_start_sample is past the signal end + sat_past_end = SatState( + SatState(gpsl1, 1, code_phase, 1000.0Hz); + signal_start_sample = num_samples_signal + 1, + ) + sss_skip = SystemSatsState(gpsl1, [sat_past_end]) + ts_skip = TrackState((sss_skip,)) + result_skip = downconvert_and_correlate( + downconvert_and_correlator, signal, ts_skip, 1, sampling_frequency, intermediate_frequency, + ) + @test get_correlator(result_skip, 1).accumulators == sat_past_end.correlator.accumulators +end + +@testset "Downconvert and correlate with KA (CPU backend)" begin + gpsl1 = GPSL1() + sampling_frequency = 5e6Hz + code_phase = 10.5 + num_samples_signal = 5000 + intermediate_frequency = 0.0Hz + + system_sats_state = SystemSatsState( + gpsl1, + [SatState(gpsl1, 1, code_phase, 1000.0Hz), SatState(gpsl1, 2, 11.0, 500.0Hz)]; + ) + multiple_system_sats_state = (system_sats_state,) + + downconvert_and_correlator = KADownconvertAndCorrelator((gpsl1,), Array) + + track_state = TrackState(multiple_system_sats_state) + + preferred_num_code_blocks_to_integrate = 1 + + signal = + gen_code( + num_samples_signal, + gpsl1, + 1, + sampling_frequency, + get_code_frequency(gpsl1) + 1000Hz * get_code_center_frequency_ratio(gpsl1), + code_phase, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 1000.0Hz / sampling_frequency) + + next_track_state = downconvert_and_correlate( + downconvert_and_correlator, + signal, + track_state, + preferred_num_code_blocks_to_integrate, + sampling_frequency, + intermediate_frequency, + ) + + @test real.(get_correlator(next_track_state, 1).accumulators) ≈ [2921, 4949, 2917] atol = 25 + + signal = + gen_code( + num_samples_signal, + gpsl1, + 2, + sampling_frequency, + get_code_frequency(gpsl1) + 500Hz * get_code_center_frequency_ratio(gpsl1), + 11.0, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 500.0Hz / sampling_frequency) + + next_track_state = downconvert_and_correlate( + downconvert_and_correlator, + signal, + track_state, + preferred_num_code_blocks_to_integrate, + sampling_frequency, + intermediate_frequency, + ) + + @test real.(get_correlator(next_track_state, 2).accumulators) ≈ [2919, 4947, 2915] atol = 25 +end + +@testset "Downconvert and correlate GalileoE1B: KA vs CPU" begin + gal = GalileoE1B() + sampling_frequency = 25e6Hz + code_phase = 10.5 + num_samples_signal = 100000 + intermediate_frequency = 0.0Hz + + system_sats_state = SystemSatsState( + gal, + [SatState(gal, 1, code_phase, 100.0Hz), SatState(gal, 2, 11.0, 200.0Hz)]; + ) + track_state = TrackState((system_sats_state,)) + + preferred_num_code_blocks_to_integrate = 1 + + signal = + gen_code( + num_samples_signal, + gal, + 1, + sampling_frequency, + get_code_frequency(gal) + 100Hz * get_code_center_frequency_ratio(gal), + code_phase, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 100.0Hz / sampling_frequency) + + cpu_dc = CPUDownconvertAndCorrelator(Val(sampling_frequency)) + cpu_result = downconvert_and_correlate( + cpu_dc, signal, track_state, preferred_num_code_blocks_to_integrate, + sampling_frequency, intermediate_frequency, + ) + + ka_dc = KADownconvertAndCorrelator((gal,), Array) + ka_result = downconvert_and_correlate( + ka_dc, signal, track_state, preferred_num_code_blocks_to_integrate, + sampling_frequency, intermediate_frequency, + ) + + cpu_accum = real.(get_correlator(cpu_result, 1).accumulators) + ka_accum = real.(get_correlator(ka_result, 1).accumulators) + @test ka_accum ≈ cpu_accum rtol = 0.01 + + # Prompt correlator (index 3 for VeryEarlyPromptLate) should be peak + prompt_idx = length(cpu_accum) ÷ 2 + 1 + @test abs(cpu_accum[prompt_idx]) == maximum(abs.(cpu_accum)) +end + +@testset "Downconvert and correlate multi-system GPSL1+GalileoE1B (KA CPU)" begin + gpsl1 = GPSL1() + gal = GalileoE1B() + sampling_frequency = 25e6Hz + num_samples_signal = 25000 + intermediate_frequency = 0.0Hz + + sss_l1 = SystemSatsState( + gpsl1, + [SatState(gpsl1, 1, 10.5, 1000.0Hz)]; + ) + sss_gal = SystemSatsState( + gal, + [SatState(gal, 1, 10.5, 100.0Hz)]; + ) + track_state = TrackState((sss_l1, sss_gal)) + + # Generate GPSL1 signal for PRN 1 + signal = + gen_code( + num_samples_signal, + gpsl1, + 1, + sampling_frequency, + get_code_frequency(gpsl1) + 1000Hz * get_code_center_frequency_ratio(gpsl1), + 10.5, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 1000.0Hz / sampling_frequency) + + ka_dc = KADownconvertAndCorrelator((gpsl1, gal), Array) + result = downconvert_and_correlate( + ka_dc, signal, track_state, 1, sampling_frequency, intermediate_frequency, + ) + + # GPSL1 sat should correlate well (prompt peak at index 2 for EarlyPromptLate) + l1_accum = real.(get_correlator(result, 1, 1).accumulators) + l1_prompt_idx = length(l1_accum) ÷ 2 + 1 + @test l1_accum[l1_prompt_idx] > l1_accum[1] + @test l1_accum[l1_prompt_idx] > l1_accum[end] + @test l1_accum[l1_prompt_idx] > 4000 # strong correlation + + # GalileoE1B sat should NOT correlate with GPSL1 signal (noise-level) + gal_accum = real.(get_correlator(result, 2, 1).accumulators) + gal_prompt_idx = length(gal_accum) ÷ 2 + 1 + @test abs(gal_accum[gal_prompt_idx]) < abs(l1_accum[l1_prompt_idx]) / 10 +end + +const CUDA_AVAILABLE = try + @eval using CUDA: CUDA, CuArray, cu + CUDA.functional() +catch + false +end + +if CUDA_AVAILABLE + @testset "Downconvert and correlate with KA (CUDA backend)" begin + gpsl1 = GPSL1() + sampling_frequency = 5e6Hz + code_phase = 10.5 + num_samples_signal = 5000 + intermediate_frequency = 0.0Hz + + system_sats_state = SystemSatsState( + gpsl1, + [ + SatState(gpsl1, 1, code_phase, 1000.0Hz), + SatState(gpsl1, 2, 11.0, 500.0Hz), + ]; + ) + multiple_system_sats_state = (system_sats_state,) + + downconvert_and_correlator = KADownconvertAndCorrelator((gpsl1,), CuArray) + + track_state = TrackState(multiple_system_sats_state) + + preferred_num_code_blocks_to_integrate = 1 + + signal_cpu = + gen_code( + num_samples_signal, + gpsl1, + 1, + sampling_frequency, + get_code_frequency(gpsl1) + + 1000Hz * get_code_center_frequency_ratio(gpsl1), + code_phase, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 1000.0Hz / sampling_frequency) + signal = CuArray(ComplexF32.(signal_cpu)) + + next_track_state = downconvert_and_correlate( + downconvert_and_correlator, + signal, + track_state, + preferred_num_code_blocks_to_integrate, + sampling_frequency, + intermediate_frequency, + ) + + @test real.(get_correlator(next_track_state, 1).accumulators) ≈ + [2921, 4949, 2917] atol = 25 + + signal_cpu = + gen_code( + num_samples_signal, + gpsl1, + 2, + sampling_frequency, + get_code_frequency(gpsl1) + + 500Hz * get_code_center_frequency_ratio(gpsl1), + 11.0, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 500.0Hz / sampling_frequency) + signal = CuArray(ComplexF32.(signal_cpu)) + + next_track_state = downconvert_and_correlate( + downconvert_and_correlator, + signal, + track_state, + preferred_num_code_blocks_to_integrate, + sampling_frequency, + intermediate_frequency, + ) + + @test real.(get_correlator(next_track_state, 2).accumulators) ≈ + [2919, 4947, 2915] atol = 25 + end +end + +const AMDGPU_AVAILABLE = try + @eval using AMDGPU + AMDGPU.functional() +catch + false +end + +if AMDGPU_AVAILABLE + @testset "Downconvert and correlate with KA (AMDGPU backend)" begin + gpsl1 = GPSL1() + sampling_frequency = 5e6Hz + code_phase = 10.5 + num_samples_signal = 5000 + intermediate_frequency = 0.0Hz + + system_sats_state = SystemSatsState( + gpsl1, + [ + SatState(gpsl1, 1, code_phase, 1000.0Hz), + SatState(gpsl1, 2, 11.0, 500.0Hz), + ]; + ) + multiple_system_sats_state = (system_sats_state,) + + downconvert_and_correlator = KADownconvertAndCorrelator((gpsl1,), ROCArray) + + track_state = TrackState(multiple_system_sats_state) + + preferred_num_code_blocks_to_integrate = 1 + + signal_cpu = + gen_code( + num_samples_signal, + gpsl1, + 1, + sampling_frequency, + get_code_frequency(gpsl1) + + 1000Hz * get_code_center_frequency_ratio(gpsl1), + code_phase, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 1000.0Hz / sampling_frequency) + signal = ROCArray(ComplexF32.(signal_cpu)) + + next_track_state = downconvert_and_correlate( + downconvert_and_correlator, + signal, + track_state, + preferred_num_code_blocks_to_integrate, + sampling_frequency, + intermediate_frequency, + ) + + @test real.(get_correlator(next_track_state, 1).accumulators) ≈ + [2921, 4949, 2917] atol = 25 + + signal_cpu = + gen_code( + num_samples_signal, + gpsl1, + 2, + sampling_frequency, + get_code_frequency(gpsl1) + + 500Hz * get_code_center_frequency_ratio(gpsl1), + 11.0, + ) .* cis.(2π * (0:(num_samples_signal-1)) * 500.0Hz / sampling_frequency) + signal = ROCArray(ComplexF32.(signal_cpu)) + + next_track_state = downconvert_and_correlate( + downconvert_and_correlator, + signal, + track_state, + preferred_num_code_blocks_to_integrate, + sampling_frequency, + intermediate_frequency, + ) + + @test real.(get_correlator(next_track_state, 2).accumulators) ≈ + [2919, 4947, 2915] atol = 25 + end end end diff --git a/test/ka_gpu_tests.jl b/test/ka_gpu_tests.jl new file mode 100644 index 00000000..e5b66e97 --- /dev/null +++ b/test/ka_gpu_tests.jl @@ -0,0 +1,10 @@ +# KernelAbstractions GPU tests only — run on Buildkite with CUDA/AMDGPU agents. +# Usage: julia --project -e 'using Pkg; Pkg.test()' +# Or directly: julia --project test/ka_gpu_tests.jl + +using Test + +@testset "KA GPU Tests" begin + include("downconvert_and_correlate.jl") + include("track.jl") +end diff --git a/test/runtests.jl b/test/runtests.jl index 58d831f4..7df64986 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,6 @@ include("downconvert_and_correlate.jl") include("post_corr_filter.jl") include("code_replica.jl") include("carrier_replica.jl") -include("downconvert.jl") include("discriminators.jl") include("gps_l1.jl") include("gps_l5.jl") diff --git a/test/track.jl b/test/track.jl index 6d4653b0..27f8e766 100644 --- a/test/track.jl +++ b/test/track.jl @@ -30,7 +30,8 @@ using Tracking: NumAnts, get_num_ants, ConventionalPLLAndDLL, - ConventionalAssistedPLLAndDLL + ConventionalAssistedPLLAndDLL, + KADownconvertAndCorrelator @testset "Tracking with signal of type $type" for type in (Int16, Int32, Int64, Float32, Float64) @@ -552,4 +553,239 @@ end # plot(imag.(tracked_prompts)) end +@testset "Tracking GPSL1 with KADownconvertAndCorrelator (CPU backend)" begin + gpsl1 = GPSL1() + carrier_doppler = 200Hz + start_code_phase = 100 + code_frequency = + carrier_doppler * get_code_center_frequency_ratio(gpsl1) + get_code_frequency(gpsl1) + sampling_frequency = 4e6Hz + prn = 1 + range = 0:3999 + start_carrier_phase = π / 2 + + ka_dc = KADownconvertAndCorrelator((gpsl1,), Array) + + track_state = TrackState( + gpsl1, + [SatState(gpsl1, prn, start_code_phase, carrier_doppler - 20Hz)], + ) + + signal = Complex{Float32}.( + cis.(2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase) .* + gen_code(4000, gpsl1, prn, sampling_frequency, code_frequency, start_code_phase) + ) + track_state = track(signal, track_state, sampling_frequency; + downconvert_and_correlator = ka_dc, + ) + + iterations = 2000 + for i = 1:iterations + carrier_phase = + mod2pi( + 2π * carrier_doppler * 4000 * i / sampling_frequency + + start_carrier_phase + π, + ) - π + code_phase = + mod(code_frequency * 4000 * i / sampling_frequency + start_code_phase, 1023) + signal = Complex{Float32}.( + cis.(2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase) .* + gen_code(4000, gpsl1, prn, sampling_frequency, code_frequency, code_phase) + ) + track_state = track(signal, track_state, sampling_frequency; + downconvert_and_correlator = ka_dc, + ) + end + comp_code_phase = mod( + code_frequency * 4000 * (iterations + 1) / sampling_frequency + start_code_phase, + 1023, + ) + comp_carrier_phase = + mod2pi( + 2π * carrier_doppler * 4000 * (iterations + 1) / sampling_frequency + + start_carrier_phase + π, + ) - π + @test get_code_phase(track_state) ≈ comp_code_phase atol = 5e-5 + @test get_carrier_phase(track_state) + π ≈ comp_carrier_phase atol = 1e-3 +end + +@testset "Track GPSL1+GalileoE1B with KADownconvertAndCorrelator (CPU backend)" begin + gpsl1 = GPSL1() + galileo_e1b = GalileoE1B() + carrier_doppler_gps = 200Hz + carrier_doppler_gal = 1200Hz + start_code_phase = 100 + code_frequency_gps = + carrier_doppler_gps * get_code_center_frequency_ratio(gpsl1) + + get_code_frequency(gpsl1) + code_frequency_gal = + carrier_doppler_gal * get_code_center_frequency_ratio(galileo_e1b) + + get_code_frequency(galileo_e1b) + sampling_frequency = 15e6Hz + prn = 1 + range = 0:3999 + start_carrier_phase = π / 2 + + ka_dc = KADownconvertAndCorrelator((gpsl1, galileo_e1b), Array) + + track_state = TrackState(( + gps = SystemSatsState( + gpsl1, + [SatState(gpsl1, prn, start_code_phase, carrier_doppler_gps)], + ), + gal = SystemSatsState( + galileo_e1b, + [SatState(galileo_e1b, prn, start_code_phase, carrier_doppler_gal)], + ), + );) + + signal_temp = + cis.( + 2π .* carrier_doppler_gps .* range ./ sampling_frequency .+ start_carrier_phase, + ) .* gen_code( + 4000, gpsl1, prn, sampling_frequency, code_frequency_gps, start_code_phase, + ) .+ + cis.( + 2π .* carrier_doppler_gal .* range ./ sampling_frequency .+ start_carrier_phase, + ) .* gen_code( + 4000, galileo_e1b, prn, sampling_frequency, code_frequency_gal, start_code_phase, + ) + signal = ComplexF32.(signal_temp) + + track_state = track(signal, track_state, sampling_frequency; + downconvert_and_correlator = ka_dc, + ) + + iterations = 2000 + for i = 1:iterations + carrier_phase_gps = + mod2pi( + 2π * carrier_doppler_gps * 4000 * i / sampling_frequency + + start_carrier_phase + π, + ) - π + code_phase_gps = mod( + code_frequency_gps * 4000 * i / sampling_frequency + start_code_phase, + get_code_length(gpsl1), + ) + carrier_phase_gal = + mod2pi( + 2π * carrier_doppler_gal * 4000 * i / sampling_frequency + + start_carrier_phase + π, + ) - π + code_phase_gal = mod( + code_frequency_gal * 4000 * i / sampling_frequency + start_code_phase, + get_code_length(galileo_e1b), + ) + signal_temp = + cis.( + 2π .* carrier_doppler_gps .* range ./ sampling_frequency .+ + carrier_phase_gps, + ) .* gen_code( + 4000, gpsl1, prn, sampling_frequency, code_frequency_gps, code_phase_gps, + ) .+ + cis.( + 2π .* carrier_doppler_gal .* range ./ sampling_frequency .+ + carrier_phase_gal, + ) .* gen_code( + 4000, galileo_e1b, prn, sampling_frequency, code_frequency_gal, code_phase_gal, + ) + signal = ComplexF32.(signal_temp) + track_state = track(signal, track_state, sampling_frequency; + downconvert_and_correlator = ka_dc, + ) + end + comp_code_phase_gps = mod( + code_frequency_gps * 4000 * (iterations + 1) / sampling_frequency + + start_code_phase, + get_code_length(gpsl1), + ) + comp_carrier_phase_gps = + mod2pi( + 2π * carrier_doppler_gps * 4000 * (iterations + 1) / sampling_frequency + + start_carrier_phase + π, + ) - π + comp_code_phase_gal = mod( + code_frequency_gal * 4000 * (iterations + 1) / sampling_frequency + + start_code_phase, + get_code_length(galileo_e1b), + ) + comp_carrier_phase_gal = + mod2pi( + 2π * carrier_doppler_gal * 4000 * (iterations + 1) / sampling_frequency + + start_carrier_phase + π, + ) - π + @test get_code_phase(track_state, :gps, 1) ≈ comp_code_phase_gps atol = 5e-3 + @test mod(get_carrier_phase(track_state, :gps, 1), π) ≈ mod(comp_carrier_phase_gps, π) atol = + 2e-2 + @test get_code_phase(track_state, :gal, 1) ≈ comp_code_phase_gal atol = 5e-3 + @test mod(get_carrier_phase(track_state, :gal, 1), π) ≈ mod(comp_carrier_phase_gal, π) atol = + 5e-3 +end + +const CUDA_AVAILABLE = try + @eval using CUDA: CUDA, CuArray + CUDA.functional() +catch + false +end + +if CUDA_AVAILABLE + @testset "Tracking GPSL1 with KADownconvertAndCorrelator (CUDA backend)" begin + gpsl1 = GPSL1() + carrier_doppler = 200Hz + start_code_phase = 100 + code_frequency = + carrier_doppler * get_code_center_frequency_ratio(gpsl1) + + get_code_frequency(gpsl1) + sampling_frequency = 4e6Hz + prn = 1 + range = 0:3999 + start_carrier_phase = π / 2 + + ka_dc = KADownconvertAndCorrelator((gpsl1,), CuArray) + + track_state = TrackState( + gpsl1, + [SatState(gpsl1, prn, start_code_phase, carrier_doppler - 20Hz)], + ) + + signal = CuArray(Complex{Float32}.( + cis.(2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase) .* + gen_code(4000, gpsl1, prn, sampling_frequency, code_frequency, start_code_phase) + )) + track_state = track(signal, track_state, sampling_frequency; + downconvert_and_correlator = ka_dc, + ) + + iterations = 2000 + for i = 1:iterations + carrier_phase = + mod2pi( + 2π * carrier_doppler * 4000 * i / sampling_frequency + + start_carrier_phase + π, + ) - π + code_phase = + mod(code_frequency * 4000 * i / sampling_frequency + start_code_phase, 1023) + signal = CuArray(Complex{Float32}.( + cis.(2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase) .* + gen_code(4000, gpsl1, prn, sampling_frequency, code_frequency, code_phase) + )) + track_state = track(signal, track_state, sampling_frequency; + downconvert_and_correlator = ka_dc, + ) + end + comp_code_phase = mod( + code_frequency * 4000 * (iterations + 1) / sampling_frequency + start_code_phase, + 1023, + ) + comp_carrier_phase = + mod2pi( + 2π * carrier_doppler * 4000 * (iterations + 1) / sampling_frequency + + start_carrier_phase + π, + ) - π + @test get_code_phase(track_state) ≈ comp_code_phase atol = 5e-5 + @test get_carrier_phase(track_state) + π ≈ comp_carrier_phase atol = 1e-3 + end +end + end