path: root/tools/spectral_editor/dct.js

const dctSize = 512; // Default DCT size, read from header

// --- Utility Functions for Audio Processing ---
// JavaScript equivalent of C++ idct_512
function javascript_idct_512(input) {
    const output = new Float32Array(dctSize);
    const PI = Math.PI;
    const N = dctSize;

    for (let n = 0; n < N; ++n) {
        let sum = input[0] / 2.0;
        for (let k = 1; k < N; ++k) {
            sum += input[k] * Math.cos((PI / N) * k * (n + 0.5));
        }
        output[n] = sum * (2.0 / N);
    }
    return output;
}

// Hanning window for smooth audio transitions (JavaScript equivalent)
function hanningWindow(size) {
    const window = new Float32Array(size);
    const PI = Math.PI;
    for (let i = 0; i < size; i++) {
        window[i] = 0.5 * (1 - Math.cos((2 * PI * i) / (size - 1)));
    }
    return window;
}

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);
}