3 files changed, 473 insertions, 0 deletions
diff --git a/src/audio/fft.cc b/src/audio/fft.cc
new file mode 100644
index 0000000..25477b9
--- /dev/null
+++ b/src/audio/fft.cc
@@ -0,0 +1,194 @@
+// Fast Fourier Transform (FFT) implementation
+// Radix-2 Cooley-Tukey algorithm
+// Reference: Numerical Recipes, Press et al.
+
+#include "audio/fft.h"
+
+#include <cmath>
+#include <cstring>
+
+// Bit-reversal permutation (in-place)
+// Reorders array elements by reversing their binary indices
+static void bit_reverse_permute(float* real, float* imag, size_t N) {
+  const size_t bits = 0;
+  size_t temp_bits = N;
+  size_t num_bits = 0;
+  while (temp_bits > 1) {
+    temp_bits >>= 1;
+    num_bits++;
+  }
+
+  for (size_t i = 0; i < N; i++) {
+    // Compute bit-reversed index
+    size_t j = 0;
+    size_t temp = i;
+    for (size_t b = 0; b < num_bits; b++) {
+      j = (j << 1) | (temp & 1);
+      temp >>= 1;
+    }
+
+    // Swap if j > i (to avoid swapping twice)
+    if (j > i) {
+      const float tmp_real = real[i];
+      const float tmp_imag = imag[i];
+      real[i] = real[j];
+      imag[i] = imag[j];
+      real[j] = tmp_real;
+      imag[j] = tmp_imag;
+    }
+  }
+}
+
+// In-place radix-2 FFT (after bit-reversal)
+// direction: +1 for forward FFT, -1 for inverse FFT
+static void fft_radix2(float* real, float* imag, size_t N, int direction) {
+  const float PI = 3.14159265358979323846f;
+
+  // Butterfly operations
+  for (size_t stage_size = 2; stage_size <= N; stage_size *= 2) {
+    const size_t half_stage = stage_size / 2;
+    const float angle = direction * 2.0f * PI / stage_size;
+
+    // Precompute twiddle factors for this stage
+    float wr = 1.0f;
+    float wi = 0.0f;
+    const float wr_delta = cosf(angle);
+    const float wi_delta = sinf(angle);
+
+    for (size_t k = 0; k < half_stage; k++) {
+      // Apply butterfly to all groups at this stage
+      for (size_t group_start = k; group_start < N; group_start += stage_size) {
+        const size_t i = group_start;
+        const size_t j = group_start + half_stage;
+
+        // Complex multiplication: (real[j] + i*imag[j]) * (wr + i*wi)
+        const float temp_real = real[j] * wr - imag[j] * wi;
+        const float temp_imag = real[j] * wi + imag[j] * wr;
+
+        // Butterfly operation
+        real[j] = real[i] - temp_real;
+        imag[j] = imag[i] - temp_imag;
+        real[i] = real[i] + temp_real;
+        imag[i] = imag[i] + temp_imag;
+      }
+
+      // Update twiddle factor for next k (rotation)
+      const float wr_old = wr;
+      wr = wr_old * wr_delta - wi * wi_delta;
+      wi = wr_old * wi_delta + wi * wr_delta;
+    }
+  }
+}
+
+void fft_forward(float* real, float* imag, size_t N) {
+  bit_reverse_permute(real, imag, N);
+  fft_radix2(real, imag, N, +1);
+}
+
+void fft_inverse(float* real, float* imag, size_t N) {
+  bit_reverse_permute(real, imag, N);
+  fft_radix2(real, imag, N, -1);
+
+  // Scale by 1/N
+  const float scale = 1.0f / N;
+  for (size_t i = 0; i < N; i++) {
+    real[i] *= scale;
+    imag[i] *= scale;
+  }
+}
+
+// DCT-II via FFT using double-and-mirror method
+// This is a more robust algorithm that avoids reordering issues
+// Reference: Numerical Recipes, Press et al.
+void dct_fft(const float* input, float* output, size_t N) {
+  const float PI = 3.14159265358979323846f;
+
+  // Allocate temporary arrays for 2N-point FFT
+  const size_t M = 2 * N;
+  float* real = new float[M];
+  float* imag = new float[M];
+
+  // Pack input: [x[0], x[1], ..., x[N-1], x[N-1], x[N-2], ..., x[1]]
+  // This creates even symmetry for real-valued DCT
+  for (size_t i = 0; i < N; i++) {
+    real[i] = input[i];
+  }
+  for (size_t i = 0; i < N; i++) {
+    real[N + i] = input[N - 1 - i];
+  }
+  memset(imag, 0, M * sizeof(float));
+
+  // Apply 2N-point FFT
+  fft_forward(real, imag, M);
+
+  // Extract DCT coefficients
+  // DCT[k] = Re{FFT[k] * exp(-j*pi*k/(2*N))} * normalization
+  // Note: Need to divide by 2 because we doubled the signal length
+  for (size_t k = 0; k < N; k++) {
+    const float angle = -PI * k / (2.0f * N);
+    const float wr = cosf(angle);
+    const float wi = sinf(angle);
+
+    // Complex multiplication: (real + j*imag) * (wr + j*wi)
+    // Real part: real*wr - imag*wi
+    const float dct_value = real[k] * wr - imag[k] * wi;
+
+    // Apply DCT-II normalization (divide by 2 for double-length FFT)
+    if (k == 0) {
+      output[k] = dct_value * sqrtf(1.0f / N) / 2.0f;
+    } else {
+      output[k] = dct_value * sqrtf(2.0f / N) / 2.0f;
+    }
+  }
+
+  delete[] real;
+  delete[] imag;
+}
+
+// IDCT (Inverse DCT-II) via FFT using double-and-mirror method
+// This is the inverse of the DCT-II (used for synthesis)
+void idct_fft(const float* input, float* output, size_t N) {
+  const float PI = 3.14159265358979323846f;
+
+  // Allocate temporary arrays for 2N-point FFT
+  const size_t M = 2 * N;
+  float* real = new float[M];
+  float* imag = new float[M];
+
+  // Prepare FFT input from DCT coefficients
+  // IDCT = Re{IFFT[DCT * exp(j*pi*k/(2*N))]} * 2
+  for (size_t k = 0; k < N; k++) {
+    const float angle = PI * k / (2.0f * N);  // Positive for inverse
+    const float wr = cosf(angle);
+    const float wi = sinf(angle);
+
+    // Apply inverse normalization
+    float scaled_input;
+    if (k == 0) {
+      scaled_input = input[k] * sqrtf(N) * 2.0f;
+    } else {
+      scaled_input = input[k] * sqrtf(N / 2.0f) * 2.0f;
+    }
+
+    // Complex multiplication: DCT[k] * exp(j*theta)
+    real[k] = scaled_input * wr;
+    imag[k] = scaled_input * wi;
+  }
+
+  // Fill second half with conjugate symmetry (for real output)
+  for (size_t k = 1; k < N; k++) {
+    real[M - k] = real[k];
+    imag[M - k] = -imag[k];
+  }
+
+  // Apply inverse FFT
+  fft_inverse(real, imag, M);
+
+  // Extract first N samples (real part only, imag should be ~0)
+  for (size_t i = 0; i < N; i++) {
+    output[i] = real[i];
+  }
+
+  delete[] real;
+  delete[] imag;
+}
diff --git a/src/audio/fft.h b/src/audio/fft.h
new file mode 100644
index 0000000..81a12d4
--- /dev/null
+++ b/src/audio/fft.h
@@ -0,0 +1,34 @@
+// Fast Fourier Transform (FFT) implementation
+// Radix-2 Cooley-Tukey algorithm for power-of-2 sizes
+// This implementation matches the JavaScript version in tools/spectral_editor/dct.js
+
+#ifndef AUDIO_FFT_H_
+#define AUDIO_FFT_H_
+
+#include <cstddef>
+
+// Forward FFT: Time domain → Frequency domain
+// Input: real[] (length N), imag[] (length N, can be zeros)
+// Output: real[] and imag[] contain complex frequency bins
+// N must be a power of 2
+void fft_forward(float* real, float* imag, size_t N);
+
+// Inverse FFT: Frequency domain → Time domain
+// Input: real[] (length N), imag[] (length N)
+// Output: real[] and imag[] contain complex time samples (scaled by 1/N)
+// N must be a power of 2
+void fft_inverse(float* real, float* imag, size_t N);
+
+// DCT Type-II via FFT (matches existing dct_512 signature)
+// Input: input[] (length N)
+// Output: output[] (length N) containing DCT-II coefficients
+// N must be a power of 2
+void dct_fft(const float* input, float* output, size_t N);
+
+// IDCT (Inverse DCT Type-II) via FFT (matches existing idct_512 signature)
+// Input: input[] (length N)
+// Output: output[] (length N) containing time-domain samples
+// N must be a power of 2
+void idct_fft(const float* input, float* output, size_t N);
+
+#endif  /* AUDIO_FFT_H_ */
diff --git a/src/tests/test_fft.cc b/src/tests/test_fft.cc
new file mode 100644
index 0000000..948090a
--- /dev/null
+++ b/src/tests/test_fft.cc
@@ -0,0 +1,245 @@
+// Tests for FFT-based DCT/IDCT implementation
+// Verifies correctness against reference O(N²) implementation
+
+#include "audio/fft.h"
+
+#include <cassert>
+#include <cmath>
+#include <cstdio>
+#include <cstring>
+
+// Reference O(N²) DCT-II implementation (from original code)
+static void dct_reference(const float* input, float* output, size_t N) {
+  const float PI = 3.14159265358979323846f;
+
+  for (size_t k = 0; k < N; k++) {
+    float sum = 0.0f;
+    for (size_t n = 0; n < N; n++) {
+      sum += input[n] * cosf((PI / N) * k * (n + 0.5f));
+    }
+
+    // Apply DCT-II normalization
+    if (k == 0) {
+      output[k] = sum * sqrtf(1.0f / N);
+    } else {
+      output[k] = sum * sqrtf(2.0f / N);
+    }
+  }
+}
+
+// Reference O(N²) IDCT implementation (from original code)
+static void idct_reference(const float* input, float* output, size_t N) {
+  const float PI = 3.14159265358979323846f;
+
+  for (size_t n = 0; n < N; ++n) {
+    float sum = input[0] / 2.0f;
+    for (size_t k = 1; k < N; ++k) {
+      sum += input[k] * cosf((PI / N) * k * (n + 0.5f));
+    }
+    output[n] = sum * (2.0f / N);
+  }
+}
+
+// Compare two arrays with tolerance
+// Note: FFT-based DCT accumulates slightly more rounding error than O(N²) direct method
+// A tolerance of 1e-3 is still excellent for audio applications (< -60 dB error)
+static bool arrays_match(const float* a,
+                         const float* b,
+                         size_t N,
+                         float tolerance = 1e-3f) {
+  for (size_t i = 0; i < N; i++) {
+    const float diff = fabsf(a[i] - b[i]);
+    if (diff > tolerance) {
+      fprintf(stderr,
+              "Mismatch at index %zu: %.6f vs %.6f (diff=%.6e)\n",
+              i,
+              a[i],
+              b[i],
+              diff);
+      return false;
+    }
+  }
+  return true;
+}
+
+// Test 1: DCT correctness (FFT-based vs reference)
+static void test_dct_correctness() {
+  printf("Test 1: DCT correctness (FFT vs reference O(N²))...\n");
+
+  const size_t N = 512;
+  float input[N];
+  float output_ref[N];
+  float output_fft[N];
+
+  // Test case 1: Impulse at index 0
+  memset(input, 0, N * sizeof(float));
+  input[0] = 1.0f;
+
+  dct_reference(input, output_ref, N);
+  dct_fft(input, output_fft, N);
+
+  assert(arrays_match(output_ref, output_fft, N));
+  printf("  ✓ Impulse test passed\n");
+
+  // Test case 2: Impulse at middle
+  memset(input, 0, N * sizeof(float));
+  input[N / 2] = 1.0f;
+
+  dct_reference(input, output_ref, N);
+  dct_fft(input, output_fft, N);
+
+  assert(arrays_match(output_ref, output_fft, N));
+  printf("  ✓ Middle impulse test passed\n");
+
+  // Test case 3: Sinusoidal input
+  for (size_t i = 0; i < N; i++) {
+    input[i] = sinf(2.0f * 3.14159265358979323846f * 5.0f * i / N);
+  }
+
+  dct_reference(input, output_ref, N);
+  dct_fft(input, output_fft, N);
+
+  assert(arrays_match(output_ref, output_fft, N));
+  printf("  ✓ Sinusoidal input test passed\n");
+
+  // Test case 4: Random-ish input (deterministic)
+  for (size_t i = 0; i < N; i++) {
+    input[i] = sinf(i * 0.1f) * cosf(i * 0.05f);
+  }
+
+  dct_reference(input, output_ref, N);
+  dct_fft(input, output_fft, N);
+
+  assert(arrays_match(output_ref, output_fft, N));
+  printf("  ✓ Complex input test passed\n");
+
+  printf("Test 1: PASSED ✓\n\n");
+}
+
+// Test 2: IDCT correctness (FFT-based vs reference)
+static void test_idct_correctness() {
+  printf("Test 2: IDCT correctness (FFT vs reference O(N²))...\n");
+
+  const size_t N = 512;
+  float input[N];
+  float output_ref[N];
+  float output_fft[N];
+
+  // Test case 1: DC component only
+  memset(input, 0, N * sizeof(float));
+  input[0] = 1.0f;
+
+  idct_reference(input, output_ref, N);
+  idct_fft(input, output_fft, N);
+
+  assert(arrays_match(output_ref, output_fft, N));
+  printf("  ✓ DC component test passed\n");
+
+  // Test case 2: Single frequency bin
+  memset(input, 0, N * sizeof(float));
+  input[10] = 1.0f;
+
+  idct_reference(input, output_ref, N);
+  idct_fft(input, output_fft, N);
+
+  assert(arrays_match(output_ref, output_fft, N));
+  printf("  ✓ Single bin test passed\n");
+
+  // Test case 3: Mixed frequencies
+  for (size_t i = 0; i < N; i++) {
+    input[i] = (i % 10 == 0) ? 1.0f : 0.0f;
+  }
+
+  idct_reference(input, output_ref, N);
+  idct_fft(input, output_fft, N);
+
+  assert(arrays_match(output_ref, output_fft, N));
+  printf("  ✓ Mixed frequencies test passed\n");
+
+  printf("Test 2: PASSED ✓\n\n");
+}
+
+// Test 3: Round-trip (DCT → IDCT should recover original)
+static void test_roundtrip() {
+  printf("Test 3: Round-trip (DCT → IDCT = identity)...\n");
+
+  const size_t N = 512;
+  float input[N];
+  float dct_output[N];
+  float reconstructed[N];
+
+  // Test case 1: Sinusoidal input
+  for (size_t i = 0; i < N; i++) {
+    input[i] = sinf(2.0f * 3.14159265358979323846f * 3.0f * i / N);
+  }
+
+  dct_fft(input, dct_output, N);
+  idct_fft(dct_output, reconstructed, N);
+
+  assert(arrays_match(input, reconstructed, N));
+  printf("  ✓ Sinusoidal round-trip passed\n");
+
+  // Test case 2: Complex signal
+  for (size_t i = 0; i < N; i++) {
+    input[i] = sinf(i * 0.1f) * cosf(i * 0.05f) + cosf(i * 0.03f);
+  }
+
+  dct_fft(input, dct_output, N);
+  idct_fft(dct_output, reconstructed, N);
+
+  assert(arrays_match(input, reconstructed, N));
+  printf("  ✓ Complex signal round-trip passed\n");
+
+  printf("Test 3: PASSED ✓\n\n");
+}
+
+// Test 4: Output known values for JavaScript comparison
+static void test_known_values() {
+  printf("Test 4: Known values (for JavaScript verification)...\n");
+
+  const size_t N = 512;
+  float input[N];
+  float output[N];
+
+  // Simple test case: impulse at index 0
+  memset(input, 0, N * sizeof(float));
+  input[0] = 1.0f;
+
+  dct_fft(input, output, N);
+
+  printf("  DCT of impulse at 0:\n");
+  printf("    output[0] = %.8f (expected ~0.04419417)\n", output[0]);
+  printf("    output[1] = %.8f (expected ~0.04419417)\n", output[1]);
+  printf("    output[10] = %.8f (expected ~0.04419417)\n", output[10]);
+
+  // IDCT test
+  memset(input, 0, N * sizeof(float));
+  input[0] = 1.0f;
+
+  idct_fft(input, output, N);
+
+  printf("  IDCT of DC component:\n");
+  printf("    output[0] = %.8f (expected ~0.04419417)\n", output[0]);
+  printf("    output[100] = %.8f (expected ~0.04419417)\n", output[100]);
+  printf("    output[511] = %.8f (expected ~0.04419417)\n", output[511]);
+
+  printf("Test 4: PASSED ✓\n");
+  printf("(Copy these values to JavaScript test for verification)\n\n");
+}
+
+int main() {
+  printf("===========================================\n");
+  printf("FFT-based DCT/IDCT Test Suite\n");
+  printf("===========================================\n\n");
+
+  test_dct_correctness();
+  test_idct_correctness();
+  test_roundtrip();
+  test_known_values();
+
+  printf("===========================================\n");
+  printf("All tests PASSED ✓\n");
+  printf("===========================================\n");
+
+  return 0;
+}