spectrogram : remove old FFT algorithm

Reuse the one embedded within ggwave through a static function

examples/spectrogram/fft.h

@@ -1,92 +0,0 @@
-#pragma once
-
-#include <cmath>
-#include <stdlib.h>
-
-// FFT routines taken from https://stackoverflow.com/a/37729648/4039976
-
-constexpr auto kMaxSamplesPerFrame = 1024;
-
-int log2(int N) {
-    int k = N, i = 0;
-    while(k) {
-        k >>= 1;
-        i++;
-    }
-    return i - 1;
-}
-
-int reverse(int N, int n) {
-    int j, p = 0;
-    for(j = 1; j <= log2(N); j++) {
-        if(n & (1 << (log2(N) - j)))
-            p |= 1 << (j - 1);
-    }
-    return p;
-}
-
-void ordina(float * f1, int N) {
-    float f2[2*kMaxSamplesPerFrame];
-    for (int i = 0; i < N; i++) {
-        int ir = reverse(N, i);
-        f2[2*i + 0] = f1[2*ir + 0];
-        f2[2*i + 1] = f1[2*ir + 1];
-    }
-    for (int j = 0; j < N; j++) {
-        f1[2*j + 0] = f2[2*j + 0];
-        f1[2*j + 1] = f2[2*j + 1];
-    }
-}
-
-void transform(float * f, int N) {
-    ordina(f, N);    //first: reverse order
-    float * W;
-    W = (float *)malloc(N*sizeof(float));
-    W[2*1 + 0] = cos(-2.*M_PI/N);
-    W[2*1 + 1] = sin(-2.*M_PI/N);
-    W[2*0 + 0] = 1;
-    W[2*0 + 1] = 0;
-    for (int i = 2; i < N / 2; i++) {
-        W[2*i + 0] = cos(-2.*i*M_PI/N);
-        W[2*i + 1] = sin(-2.*i*M_PI/N);
-    }
-    int n = 1;
-    int a = N / 2;
-    for(int j = 0; j < log2(N); j++) {
-        for(int i = 0; i < N; i++) {
-            if(!(i & n)) {
-                int wi = (i * a) % (n * a);
-                int fi = i + n;
-                float a = W[2*wi + 0];
-                float b = W[2*wi + 1];
-                float c = f[2*fi + 0];
-                float d = f[2*fi + 1];
-                float temp[2] = { f[2*i + 0], f[2*i + 1] };
-                float Temp[2] = { a*c - b*d, b*c + a*d };
-                f[2*i + 0]  = temp[0] + Temp[0];
-                f[2*i + 1]  = temp[1] + Temp[1];
-                f[2*fi + 0] = temp[0] - Temp[0];
-                f[2*fi + 1] = temp[1] - Temp[1];
-            }
-        }
-        n *= 2;
-        a = a / 2;
-    }
-    free(W);
-}
-
-void FFT(float * f, int N, float d) {
-    transform(f, N);
-    for (int i = 0; i < N; i++) {
-        f[2*i + 0] *= d;
-        f[2*i + 1] *= d;
-    }
-}
-
-void FFT(float * src, float * dst, int N, float d) {
-    for (int i = 0; i < N; ++i) {
-        dst[2*i + 0] = src[i];
-        dst[2*i + 1] = 0.0f;
-    }
-    FFT(dst, N, d);
-}

examples/spectrogram/main.cpp

@@ -1,5 +1,3 @@
-#include "fft.h"
-
 #include "ggwave/ggwave.h"
 #include "ggwave-common.h"
 
@@ -178,14 +176,34 @@ bool GGWave_mainLoop() {
         SDL_ClearQueuedAudio(g_devIdInp);
     }
 
-    int n = 0;
+    static bool isInitialzed = false;
+
     static float data[g_nSamplesPerFrame];
     static float out[2*g_nSamplesPerFrame];
+
+    static int   workI[2*g_nSamplesPerFrame];
+    static float workF[g_nSamplesPerFrame/2];
+
+    if (!isInitialzed) {
+        memset(data, 0, sizeof(data));
+        memset(out,  0, sizeof(out));
+
+        memset(workI, 0, sizeof(workI));
+        memset(workF, 0, sizeof(workF));
+
+        isInitialzed = true;
+    }
+
+    int n = 0;
+
     do {
         n = SDL_DequeueAudio(g_devIdInp, data, sizeof(float)*g_nSamplesPerFrame);
         if (n <= 0) break;
 
-        FFT(data, out, g_nSamplesPerFrame, 1.0);
+        if (GGWave::computeFFTR(data, out, g_nSamplesPerFrame, workI, workF) == false) {
+            fprintf(stderr, "Failed to compute FFT!\n");
+            return false;
+        }
 
         for (int i = 0; i < g_nSamplesPerFrame; ++i) {
             out[i] = std::sqrt(out[2*i + 0]*out[2*i + 0] + out[2*i + 1]*out[2*i + 1]);

include/ggwave/ggwave.h

@@ -765,6 +765,19 @@ public:
     //
     bool computeFFTR(const float * src, float * dst, int N);
 
+    // Compute FFT of real values (static)
+    //
+    //   src - input real-valued data, size is N
+    //   dst - output complex-valued data, size is 2*N
+    //   ip  - work buffer, with size 2*N
+    //   w   - work buffer, with size 3 + sqrt(N/2)
+    //
+    //   First time calling thid function, make sure that ip[0] == 0
+    //   This will initialize some internal coefficients and store them in ip and w for
+    //   future usage.
+    //
+    static bool computeFFTR(const float * src, float * dst, int N, int * ip, float * w);
+
     // Resample audio waveforms from one sample rate to another using sinc interpolation
     class Resampler {
     public:

src/ggwave.cpp

@@ -1300,6 +1300,12 @@ bool GGWave::computeFFTR(const float * src, float * dst, int N) {
     return true;
 }
 
+bool GGWave::computeFFTR(const float * src, float * dst, int N, int * ip, float * w) {
+    FFT(src, dst, N, ip, w);
+
+    return true;
+}
+
 //
 // GGWave::Resampler
 //