106 lines
3.3 KiB
C
106 lines
3.3 KiB
C
#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;
|
|
}
|
|
|