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.

Author: kleinesfilmroellchen

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;

Userland/Libraries/LibDSP/CMakeLists.txt

@@ -3,7 +3,6 @@ set(SOURCES
     Track.cpp
     Effects.cpp
     Synthesizers.cpp
-    FFT.cpp
 )
 
 serenity_lib(LibDSP dsp)

Userland/Libraries/LibDSP/FFT.cpp

@@ -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;
+    }
+}
 
 }