fft arm

Appli/Core/Src/audio_analysis.c

@@ -1,4 +1,5 @@
 #include "audio_analysis.h"
+#include "arm_math.h"
 #include <math.h>
 #include <string.h>
 
@@ -15,14 +16,14 @@ static volatile uint8_t g_capture_ready_buf = 0xFFU; /* 0xFF = aucun buffer comp
 
 static bool g_analysis_enabled = false;
 
-/* --- Tables precalculees (a l'init uniquement, jamais recalculees en boucle) --- */
+/* --- Fenetre precalculee (a l'init uniquement, jamais recalculee en boucle) --- */
 static float g_hann_window[FFT_SIZE];
-static float g_twiddle_cos[FFT_SIZE / 2U];
-static float g_twiddle_sin[FFT_SIZE / 2U];
 
-/* --- Buffers de travail pour la FFT (reels/imaginaires en place) et resultat --- */
-static float g_fft_real[FFT_SIZE];
-static float g_fft_imag[FFT_SIZE];
+/* --- FFT reelle CMSIS-DSP : instance initialisee une fois, buffers de travail --- */
+static arm_rfft_fast_instance_f32 g_rfft_instance;
+static float g_fft_work[FFT_SIZE];   /* echantillons fenetres ; ecrase (scratch) par arm_rfft_fast_f32 */
+static float g_fft_out[FFT_SIZE];    /* spectre complexe empaquete, format CMSIS (voir plus bas) */
+
 static float g_magnitudes_db[NUM_BINS]; /* seule la moitie 0..Fe/2 est conservee */
 static float g_peak_freq_hz = 0.0f;
 
@@ -33,67 +34,6 @@ static void AudioAnalysis_BuildTables(void)
         /* Fenetre de Hann : reduit les fuites spectrales dues a la troncature du bloc. */
         g_hann_window[i] = 0.5f - 0.5f * cosf(2.0f * (float)M_PI * (float)i / (float)(FFT_SIZE - 1U));
     }
-    for (uint32_t i = 0; i < FFT_SIZE / 2U; i++)
-    {
-        float angle = -2.0f * (float)M_PI * (float)i / (float)FFT_SIZE;
-        g_twiddle_cos[i] = cosf(angle);
-        g_twiddle_sin[i] = sinf(angle);
-    }
-}
-
-/* Permutation en bit-reverse, etape preliminaire classique d'une FFT radix-2
- * "decimation in time" calculee en place. */
-static void AudioAnalysis_BitReverse(float *real, float *imag)
-{
-    uint32_t j = 0;
-    for (uint32_t i = 0; i < FFT_SIZE - 1U; i++)
-    {
-        if (i < j)
-        {
-            float tmp = real[i]; real[i] = real[j]; real[j] = tmp;
-            tmp = imag[i]; imag[i] = imag[j]; imag[j] = tmp;
-        }
-        uint32_t k = FFT_SIZE >> 1U;
-        while (k <= j)
-        {
-            j -= k;
-            k >>= 1U;
-        }
-        j += k;
-    }
-}
-
-/* FFT radix-2 iterative (Cooley-Tukey), en place, a partir des tables de twiddle
- * precalculees. Complexite O(N log N), largement dans les capacites d'un Cortex-M7
- * a FPU materielle pour N=1024, execute hors du chemin audio temps reel. */
-static void AudioAnalysis_ComputeFFT(float *real, float *imag)
-{
-    AudioAnalysis_BitReverse(real, imag);
-
-    for (uint32_t size = 2U; size <= FFT_SIZE; size <<= 1U)
-    {
-        uint32_t half = size >> 1U;
-        uint32_t table_step = FFT_SIZE / size;
-        for (uint32_t i = 0; i < FFT_SIZE; i += size)
-        {
-            for (uint32_t k = 0; k < half; k++)
-            {
-                uint32_t twiddle_idx = k * table_step;
-                float c = g_twiddle_cos[twiddle_idx];
-                float s = g_twiddle_sin[twiddle_idx];
-                uint32_t idx1 = i + k;
-                uint32_t idx2 = i + k + half;
-
-                float tr = real[idx2] * c - imag[idx2] * s;
-                float ti = real[idx2] * s + imag[idx2] * c;
-
-                real[idx2] = real[idx1] - tr;
-                imag[idx2] = imag[idx1] - ti;
-                real[idx1] += tr;
-                imag[idx1] += ti;
-            }
-        }
-    }
 }
 
 void AudioAnalysis_Init(void)
@@ -107,6 +47,7 @@ void AudioAnalysis_Init(void)
     memset(g_magnitudes_db, 0, sizeof(g_magnitudes_db));
 
     AudioAnalysis_BuildTables();
+    arm_rfft_fast_init_f32(&g_rfft_instance, FFT_SIZE);
 }
 
 void AudioAnalysis_PushSample(float sample)
@@ -143,11 +84,14 @@ bool AudioAnalysis_Update(void)
 
     for (uint32_t i = 0; i < FFT_SIZE; i++)
     {
-        g_fft_real[i] = g_capture_buf[ready][i] * g_hann_window[i];
-        g_fft_imag[i] = 0.0f;
+        g_fft_work[i] = g_capture_buf[ready][i] * g_hann_window[i];
     }
 
-    AudioAnalysis_ComputeFFT(g_fft_real, g_fft_imag);
+    /* arm_rfft_fast_f32 (forward, ifftFlag=0) empaquete le spectre d'un signal reel
+     * en NUM_BINS paires {re,im}, sauf la premiere qui contient les deux seules
+     * valeurs reelles (DC et Nyquist) : g_fft_out[0]=X[0] (DC), g_fft_out[1]=X[N/2]
+     * (Nyquist, non utilise ici). Pour k=1..NUM_BINS-1 : re=g_fft_out[2k], im=g_fft_out[2k+1]. */
+    arm_rfft_fast_f32(&g_rfft_instance, g_fft_work, g_fft_out, 0);
 
     /* Magnitude normalisee (une sinusoide plein-echelle donne ~0dB) puis convertie
      * en dB, avec un plancher pour eviter log(0). Destine a un affichage visuel,
@@ -161,8 +105,8 @@ bool AudioAnalysis_Update(void)
 
     for (uint32_t i = 0; i < NUM_BINS; i++)
     {
-        float re = g_fft_real[i];
-        float im = g_fft_imag[i];
+        float re = (i == 0U) ? g_fft_out[0] : g_fft_out[2U * i];
+        float im = (i == 0U) ? 0.0f : g_fft_out[2U * i + 1U];
         float mag = sqrtf(re * re + im * im) / ((float)FFT_SIZE / 4.0f);
 
         if (i >= 1U && mag > peak_mag)