From 6de633dbef574202b14796cdb614cfe8062e75c0 Mon Sep 17 00:00:00 2001
From: zymon <s@zymon.org>
Date: Sat, 18 Nov 2023 23:43:21 +0100
Subject: [PATCH] sinc vectorization test

---
 random/sinc.jl | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)
 create mode 100644 random/sinc.jl

diff --git a/random/sinc.jl b/random/sinc.jl
new file mode 100644
index 0000000..08d6f4f
--- /dev/null
+++ b/random/sinc.jl
@@ -0,0 +1,28 @@
+# @INFO: there is no SIMD version for sinc (sinpi etc.)
+# We need a sinc kernel that can be SIMD vectorized
+# Below is a sinple comparision of accuracy between a few approaches
+using GLMakie
+
+pix(x) = π * x
+sinc2(x) = x |> pix |> (x -> iszero(x) ? 1.0 : sin(x)/x)
+sinc3(x) = x |> pix |> (x -> x + eps(x)) |> (x -> sin(x)/x)
+
+
+fs = 480e3
+Δt = 1 / fs
+t = -0.5:Δt:0.5
+y_sinc = t .|> sinc;
+y_sinc2 = t .|> sinc2;
+y_sinc3 = t .|> sinc3;
+
+
+@show maximum(abs, y_sinc .- y_sinc2)
+@show maximum(abs, y_sinc .- y_sinc3)
+
+fig = Figure()
+ax = Axis(fig[1, 1])
+lines!(ax, t, abs.(y_sinc - y_sinc2) .|> log10; label="if sinc")
+lines!(ax, t, abs.(y_sinc - y_sinc3) .|> log10; label="eps sinc")
+ylims!(ax, -16, -12)
+xlims!(ax, (t |> extrema)...)
+axislegend(ax)