From d3d096ec2db9147c1a339838db25316a9c944186 Mon Sep 17 00:00:00 2001 From: Georgi Gerganov Date: Mon, 11 Jul 2022 18:21:20 +0300 Subject: [PATCH] spectrogram : remove old FFT algorithm Reuse the one embedded within ggwave through a static function --- examples/spectrogram/fft.h | 92 ----------------------------------- examples/spectrogram/main.cpp | 26 ++++++++-- include/ggwave/ggwave.h | 13 +++++ src/ggwave.cpp | 6 +++ 4 files changed, 41 insertions(+), 96 deletions(-) delete mode 100644 examples/spectrogram/fft.h diff --git a/examples/spectrogram/fft.h b/examples/spectrogram/fft.h deleted file mode 100644 index 31b6da7..0000000 --- a/examples/spectrogram/fft.h +++ /dev/null @@ -1,92 +0,0 @@ -#pragma once - -#include -#include - -// 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); -} diff --git a/examples/spectrogram/main.cpp b/examples/spectrogram/main.cpp index be15c05..08f61bd 100644 --- a/examples/spectrogram/main.cpp +++ b/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]); diff --git a/include/ggwave/ggwave.h b/include/ggwave/ggwave.h index 2543bec..5cdd01c 100644 --- a/include/ggwave/ggwave.h +++ b/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: diff --git a/src/ggwave.cpp b/src/ggwave.cpp index f2ceaf7..3abed39 100644 --- a/src/ggwave.cpp +++ b/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 //