More

Context/c_FFT.c

@@ -0,0 +1,105 @@
+#include <c_FFT.h>
+#include <math.h>
+#include <string.h>
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+/* Private Helper: Verifies if an integer is a strict power of 2 */
+C_STATIC_FORCE_INLINE
+bool c_FFT_IsPowerOfTwo(c_size_t n) {
+    return (n > 0) && ((n & (n - 1)) == 0);
+}
+
+/* Private Helper: In-place O(N) bit-reversal permutation shuffling */
+static void c_FFT_BitReversePermutation(c_Complex_t* signal, c_size_t n) {
+    c_size_t j = 0;
+    for (c_size_t i = 0; i < n; ++i) {
+        if (i < j) {
+            c_Complex_t temp = signal[i];
+            signal[i] = signal[j];
+            signal[j] = temp;
+        }
+        c_size_t bit = n >> 1;
+        while (j & bit) {
+            j ^= bit;
+            bit >>= 1;
+        }
+        j ^= bit;
+    }
+}
+
+/* ------------------------------------------------------------------------------------------------------------------ */
+/*  */
+
+
+c_err_t c_FFT_Init(c_FFT_t* self, c_size_t N, c_Allocator_t* allocator) {
+    if (!self || N == 0) return C_ERR_PARAM;
+    if (!c_FFT_IsPowerOfTwo(N)) return C_ERR_PARAM; /* Algorithm constraint check */
+
+    self->allocator = allocator ? *allocator : c_DefaultAllocator;
+    self->N = N;
+
+    return C_SUCCESS;
+}
+
+void c_FFT_Destroy(c_FFT_t* self) {
+    if (!self) return;
+    self->N = 0;
+}
+
+/* ================================================================================================================== */
+/* Non-Recursive In-Place FFT Engine paired with your formal c_Complex_t APIs */
+
+c_err_t c_FFT_Execute(const c_FFT_t* self, c_Complex_t* signal, c_bool_t inverse) {
+    if (!self || !signal) return C_ERR_PARAM;
+    if (self->N == 0) return C_SUCCESS;
+
+    c_size_t n = self->N;
+
+    /* 1. Execute the in-place bit-reversal shuffle up-front */
+    c_FFT_BitReversePermutation(signal, n);
+
+    /* 2. Bottom-up butterfly merging loops */
+    /* len represents the current sub-problem size (2, 4, 8, ..., N) */
+    for (c_size_t len = 2; len <= n; len <<= 1) {
+        double angle = 2.0 * M_PI / (double)len * (inverse ? 1.0 : -1.0);
+
+        /* Calculate principal twiddle factor step block using your standard constructor factory hook */
+        c_Complex_t wlen = c_Complex_Make(cos(angle), sin(angle));
+
+        for (c_size_t i = 0; i < n; i += len) {
+            c_Complex_t w = c_Complex_Make(1.0, 0.0);
+            c_size_t half_len = len >> 1;
+
+            for (c_size_t j = 0; j < half_len; ++j) {
+                c_size_t idx_u = i + j;
+                c_size_t idx_v = i + j + half_len;
+
+                c_Complex_t u = signal[idx_u];
+
+                /* High-performance complex multiplication butterfly product segment: t = v * w */
+                c_Complex_t t = c_Complex_Mul(signal[idx_v], w);
+
+                /* Core butterfly update mappings leveraging clean addition and subtraction interfaces */
+                signal[idx_u] = c_Complex_Add(u, t);
+                signal[idx_v] = c_Complex_Sub(u, t);
+
+                /* Advance intermediate rotation matrix step: w = w * wlen */
+                w = c_Complex_Mul(w, wlen);
+            }
+        }
+    }
+
+    /* 3. Scale Inverse IFFT vectors uniformly by 1/N using scalar operations helper */
+    if (inverse) {
+        double scale = 1.0 / (double)n;
+        for (c_size_t i = 0; i < n; ++i) {
+            signal[i] = c_Complex_MulScalar(signal[i], scale);
+        }
+    }
+
+    return C_SUCCESS;
+}
+