diff options
Diffstat (limited to 'tools/spectral_editor/dct.js')
| -rw-r--r-- | tools/spectral_editor/dct.js | 201 |
1 files changed, 201 insertions, 0 deletions
diff --git a/tools/spectral_editor/dct.js b/tools/spectral_editor/dct.js index e48ce2b..deff8a9 100644 --- a/tools/spectral_editor/dct.js +++ b/tools/spectral_editor/dct.js @@ -29,3 +29,204 @@ function hanningWindow(size) { const hanningWindowArray = hanningWindow(dctSize); // Pre-calculate window +// ============================================================================ +// FFT-based DCT/IDCT Implementation +// ============================================================================ +// Fast Fourier Transform using Radix-2 Cooley-Tukey algorithm +// This implementation MUST match the C++ version in src/audio/fft.cc exactly + +// Bit-reversal permutation (in-place) +// Reorders array elements by reversing their binary indices +function bitReversePermute(real, imag, N) { + // Calculate number of bits needed + let temp_bits = N; + let num_bits = 0; + while (temp_bits > 1) { + temp_bits >>= 1; + num_bits++; + } + + for (let i = 0; i < N; i++) { + // Compute bit-reversed index + let j = 0; + let temp = i; + for (let 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 tmp_real = real[i]; + const 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 +function fftRadix2(real, imag, N, direction) { + const PI = Math.PI; + + // Butterfly operations + for (let stage_size = 2; stage_size <= N; stage_size *= 2) { + const half_stage = stage_size / 2; + const angle = direction * 2.0 * PI / stage_size; + + // Precompute twiddle factors for this stage + let wr = 1.0; + let wi = 0.0; + const wr_delta = Math.cos(angle); + const wi_delta = Math.sin(angle); + + for (let k = 0; k < half_stage; k++) { + // Apply butterfly to all groups at this stage + for (let group_start = k; group_start < N; group_start += stage_size) { + const i = group_start; + const j = group_start + half_stage; + + // Complex multiplication: (real[j] + i*imag[j]) * (wr + i*wi) + const temp_real = real[j] * wr - imag[j] * wi; + const 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 wr_old = wr; + wr = wr_old * wr_delta - wi * wi_delta; + wi = wr_old * wi_delta + wi * wr_delta; + } + } +} + +// Forward FFT: Time domain → Frequency domain +function fftForward(real, imag, N) { + bitReversePermute(real, imag, N); + fftRadix2(real, imag, N, +1); +} + +// Inverse FFT: Frequency domain → Time domain +function fftInverse(real, imag, N) { + bitReversePermute(real, imag, N); + fftRadix2(real, imag, N, -1); + + // Scale by 1/N + const scale = 1.0 / N; + for (let i = 0; i < N; i++) { + real[i] *= scale; + imag[i] *= scale; + } +} + +// DCT-II via FFT using double-and-mirror method (matches C++ dct_fft) +// This is a more robust algorithm that avoids reordering issues +function javascript_dct_fft(input, N) { + const PI = Math.PI; + + // Allocate arrays for 2N-point FFT + const M = 2 * N; + const real = new Float32Array(M); + const imag = new Float32Array(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 (let i = 0; i < N; i++) { + real[i] = input[i]; + } + for (let i = 0; i < N; i++) { + real[N + i] = input[N - 1 - i]; + } + // imag is already zeros (Float32Array default) + + // Apply 2N-point FFT + fftForward(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 + const output = new Float32Array(N); + for (let k = 0; k < N; k++) { + const angle = -PI * k / (2.0 * N); + const wr = Math.cos(angle); + const wi = Math.sin(angle); + + // Complex multiplication: (real + j*imag) * (wr + j*wi) + // Real part: real*wr - imag*wi + const 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 * Math.sqrt(1.0 / N) / 2.0; + } else { + output[k] = dct_value * Math.sqrt(2.0 / N) / 2.0; + } + } + + return output; +} + +// IDCT (Inverse DCT-II) via FFT using double-and-mirror method (matches C++ idct_fft) +function javascript_idct_fft(input, N) { + const PI = Math.PI; + + // Allocate arrays for 2N-point FFT + const M = 2 * N; + const real = new Float32Array(M); + const imag = new Float32Array(M); + + // Prepare FFT input from DCT coefficients + // IDCT = Re{IFFT[DCT * exp(j*pi*k/(2*N))]} * 2 + for (let k = 0; k < N; k++) { + const angle = PI * k / (2.0 * N); // Positive for inverse + const wr = Math.cos(angle); + const wi = Math.sin(angle); + + // Apply inverse normalization + let scaled_input; + if (k === 0) { + scaled_input = input[k] * Math.sqrt(N) * 2.0; + } else { + scaled_input = input[k] * Math.sqrt(N / 2.0) * 2.0; + } + + // 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 (let k = 1; k < N; k++) { + real[M - k] = real[k]; + imag[M - k] = -imag[k]; + } + + // Apply inverse FFT + fftInverse(real, imag, M); + + // Extract first N samples (real part only, imag should be ~0) + const output = new Float32Array(N); + for (let i = 0; i < N; i++) { + output[i] = real[i]; + } + + return output; +} + +// Convenience wrappers for dctSize = 512 (backward compatible) +function javascript_dct_512_fft(input) { + return javascript_dct_fft(input, dctSize); +} + +function javascript_idct_512_fft(input) { + return javascript_idct_fft(input, dctSize); +} + |