Browse Source

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.
kleines Filmröllchen 3 years ago
parent
commit
00dd8f8fbe

+ 1 - 1
Userland/Applications/SoundPlayer/BarsVisualizationWidget.cpp

@@ -26,7 +26,7 @@ void BarsVisualizationWidget::paint_event(GUI::PaintEvent& event)
     if (m_sample_buffer.is_empty())
     if (m_sample_buffer.is_empty())
         return;
         return;
 
 
-    LibDSP::fft(m_sample_buffer, false);
+    LibDSP::fft(m_sample_buffer.span(), false);
     double max = AK::sqrt(m_sample_count * 2.);
     double max = AK::sqrt(m_sample_count * 2.);
 
 
     double freq_bin = m_samplerate / (double)m_sample_count;
     double freq_bin = m_samplerate / (double)m_sample_count;

+ 0 - 1
Userland/Libraries/LibDSP/CMakeLists.txt

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

+ 0 - 62
Userland/Libraries/LibDSP/FFT.cpp

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

+ 35 - 2
Userland/Libraries/LibDSP/FFT.h

@@ -7,10 +7,43 @@
 #pragma once
 #pragma once
 
 
 #include <AK/Complex.h>
 #include <AK/Complex.h>
-#include <AK/Vector.h>
+#include <AK/Math.h>
+#include <AK/Span.h>
 
 
 namespace LibDSP {
 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;
+    }
+}
 
 
 }
 }