From 00dd8f8fbef2952aa4db166b8a2ebf965f094775 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?kleines=20Filmr=C3=B6llchen?= <filmroellchen@serenityos.org>
Date: Thu, 3 Mar 2022 21:16:30 +0100
Subject: [PATCH] LibDSP: Generalize & improve FFT

Several related improvements to our Fast Fourier Transform
implementation:
- FFT now operates on spans, allowing it to use many more container
  types other than Vector. It's intended anyways that FFT transmutes the
  input data.
- FFT is now constexpr, moving the implementation to the header and
  removing the cpp file. This means that if we have static collections
  of samples, we can transform them at compile time.
- sample_data.data() weirdness is now gone.
---
 .../SoundPlayer/BarsVisualizationWidget.cpp   |  2 +-
 Userland/Libraries/LibDSP/CMakeLists.txt      |  1 -
 Userland/Libraries/LibDSP/FFT.cpp             | 62 -------------------
 Userland/Libraries/LibDSP/FFT.h               | 37 ++++++++++-
 4 files changed, 36 insertions(+), 66 deletions(-)
 delete mode 100644 Userland/Libraries/LibDSP/FFT.cpp

diff --git a/Userland/Applications/SoundPlayer/BarsVisualizationWidget.cpp b/Userland/Applications/SoundPlayer/BarsVisualizationWidget.cpp
index 0c87ced853..3e6ca3d75b 100644
--- a/Userland/Applications/SoundPlayer/BarsVisualizationWidget.cpp
+++ b/Userland/Applications/SoundPlayer/BarsVisualizationWidget.cpp
@@ -26,7 +26,7 @@ void BarsVisualizationWidget::paint_event(GUI::PaintEvent& event)
     if (m_sample_buffer.is_empty())
         return;
 
-    LibDSP::fft(m_sample_buffer, false);
+    LibDSP::fft(m_sample_buffer.span(), false);
     double max = AK::sqrt(m_sample_count * 2.);
 
     double freq_bin = m_samplerate / (double)m_sample_count;
diff --git a/Userland/Libraries/LibDSP/CMakeLists.txt b/Userland/Libraries/LibDSP/CMakeLists.txt
index f7dc72fcd1..a5de6f120e 100644
--- a/Userland/Libraries/LibDSP/CMakeLists.txt
+++ b/Userland/Libraries/LibDSP/CMakeLists.txt
@@ -3,7 +3,6 @@ set(SOURCES
     Track.cpp
     Effects.cpp
     Synthesizers.cpp
-    FFT.cpp
 )
 
 serenity_lib(LibDSP dsp)
diff --git a/Userland/Libraries/LibDSP/FFT.cpp b/Userland/Libraries/LibDSP/FFT.cpp
deleted file mode 100644
index 8083333f96..0000000000
--- a/Userland/Libraries/LibDSP/FFT.cpp
+++ /dev/null
@@ -1,62 +0,0 @@
-/*
- * Copyright (c) 2021, Cesar Torres <shortanemoia@protonmail.com>
- *
- * SPDX-License-Identifier: BSD-2-Clause
- */
-
-#include "FFT.h"
-#include <AK/Complex.h>
-#include <AK/Math.h>
-
-namespace LibDSP {
-
-// This function uses the input vector as output too. therefore, if you wish to
-// leave it intact, pass a copy to this function
-//
-// The sampling frequency must be more than twice the frequency to resolve.
-// The sample window must be at least large enough to reflect the periodicity
-// of the smallest frequency to be resolved.
-//
-// For example, to resolve a 10 KHz and a 2 Hz sine waves we need at least
-// a samplerate of 20 KHz and a window of 0.5 seconds
-//
-// If invert is true, this function computes the inverse discrete fourier transform.
-//
-// The data vector must be a power of 2
-// Adapted from https://cp-algorithms.com/algebra/fft.html
-void fft(Vector<Complex<double>>& sample_data, bool invert)
-{
-    int n = sample_data.size();
-    auto data = sample_data.data();
-
-    for (int i = 1, j = 0; i < n; i++) {
-        int bit = n >> 1;
-        for (; j & bit; bit >>= 1)
-            j ^= bit;
-        j ^= bit;
-
-        if (i < j)
-            swap(data[i], data[j]);
-    }
-
-    for (int len = 2; len <= n; len <<= 1) {
-        double ang = 2 * AK::Pi<double> / len * (invert ? -1 : 1);
-        Complex<double> wlen(AK::cos(ang), AK::sin(ang));
-        for (int i = 0; i < n; i += len) {
-            Complex<double> w = { 1., 0. };
-            for (int j = 0; j < len / 2; j++) {
-                Complex<double> u = data[i + j], v = data[i + j + len / 2] * w;
-                data[i + j] = u + v;
-                data[i + j + len / 2] = u - v;
-                w *= wlen;
-            }
-        }
-    }
-
-    if (invert) {
-        for (int i = 0; i < n; i++)
-            data[i] /= n;
-    }
-}
-
-}
diff --git a/Userland/Libraries/LibDSP/FFT.h b/Userland/Libraries/LibDSP/FFT.h
index d13d358af2..c9db991b1d 100644
--- a/Userland/Libraries/LibDSP/FFT.h
+++ b/Userland/Libraries/LibDSP/FFT.h
@@ -7,10 +7,43 @@
 #pragma once
 
 #include <AK/Complex.h>
-#include <AK/Vector.h>
+#include <AK/Math.h>
+#include <AK/Span.h>
 
 namespace LibDSP {
 
-void fft(Vector<Complex<double>>& sample_data, bool invert = false);
+constexpr void fft(Span<Complex<double>> sample_data, bool invert = false)
+{
+    int n = sample_data.size();
+
+    for (int i = 1, j = 0; i < n; i++) {
+        int bit = n >> 1;
+        for (; j & bit; bit >>= 1)
+            j ^= bit;
+        j ^= bit;
+
+        if (i < j)
+            swap(sample_data[i], sample_data[j]);
+    }
+
+    for (int len = 2; len <= n; len <<= 1) {
+        double ang = 2 * AK::Pi<double> / len * (invert ? -1 : 1);
+        Complex<double> wlen(AK::cos(ang), AK::sin(ang));
+        for (int i = 0; i < n; i += len) {
+            Complex<double> w = { 1., 0. };
+            for (int j = 0; j < len / 2; j++) {
+                Complex<double> u = sample_data[i + j], v = sample_data[i + j + len / 2] * w;
+                sample_data[i + j] = u + v;
+                sample_data[i + j + len / 2] = u - v;
+                w *= wlen;
+            }
+        }
+    }
+
+    if (invert) {
+        for (int i = 0; i < n; i++)
+            sample_data[i] /= n;
+    }
+}
 
 }