Esencia de la Transformada Discreta de Fourier (DFT)

La Transformada Discreta de Fourier (DFT) representa una transformación lineal que mapea una secuencia de números complejos \(\left<a_{0..n-1}\right>\) a otra secuencia \(\left<b_{0..n-1}\right>\). Matemáticamente, el proceso se define mediante la siguiente sumatoria, donde \(i\) denota la unidad imaginaria:

\[b_j = \sum_{k=0}^{n-1} a_k e^{-\frac{i2\pi jk}{n}}\]

El término \(e^{-\frac{i2\pi jk}{n}}\) corresponde a la \(jk\)-ésima potencia de la raíz \(n\)-ésima de la unidad, denotada frecuentemente como \(\omega_n^{jk}\). Esta operación puede expresarse como un producto matricial utilizando una matriz de Vandermonde:

\[\begin{bmatrix} \omega_n^{0\times 0} & \omega_n^{0\times 1} & \cdots & \omega_n^{0\times (n-1)} \\ \omega_n^{1\times 0} & \omega_n^{1\times 1} & \cdots & \omega_n^{1\times (n-1)} \\ \vdots & \vdots & \ddots & \vdots \\ \omega_n^{(n-1)\times 0} & \omega_n^{(n-1)\times 1} & \cdots & \omega_n^{(n-1)\times (n-1)} \end{bmatrix}\]

Una de las propiedades fundamentales de la DFT es su relación con la convolución. Si realizamos el producto punto de dos secuencias transformadas \(\hat{f}\) y \(\hat{g}\), el resultado obtenido en el dominio del tiempo (tras aplicar la transformada inversa) equivale a la convolución cíclica de las seceuncias originales:

\[h_k = \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} f_i g_j [(i+j) \bmod n = k]\]

Para recuperar la secuencia original a partir de su espectro de frecuencias, empleamos la Transformada Inversa (IDFT). Esta se puede calcular invirtiendo el signo del exponente en las raíces de la unidad o reordenando la secuencia de entrada (del índice 1 al \(n-1\)) y aplicando una DFT estándar dividida por \(n\).

Algoritmo de Transformada Rápida (FFT)

La implemetnación más eficiente de la DFT es la FFT, la cual reduce la complejidad de \(\mathcal{O}(n^2)\) a \(\mathcal{O}(n \log n)\) mediante una estrategia de "divide y vencerás". El requisito principal es que \(n\) sea una potencia de 2. Al separar los términos pares e impares de una secuencia, podemos derivar la estructura recursiva:

\[\hat{f}_i = \hat{f}_{par, i} + \omega_n^i \hat{f}_{impar, i}\]

Dado que los valores se repiten periódicamente, obtenemos las identidades de mariposa:

\[\hat{f}_i = \hat{f}_{par, i} + \omega_n^i \hat{f}_{impar, i}\] \[\hat{f}_{i+n/2} = \hat{f}_{par, i} - \omega_n^i \hat{f}_{impar, i}\]

Para evitar la sobrecarga de la recursividad, se utiliza una implementación iterativa basada en la permutaicón de inversión de bits (bit-reversal). A continuación, se presenta una estructura refinada del algoritmo:

void perform_fft(complex<double>* data, int n, bool invert) {
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >> 1) j ^= bit;
        j ^= bit;
        if (i < j) swap(data[i], data[j]);
    }

    for (int len = 2; len <= n; len <<= 1) {
        double ang = 2 * M_PI / len * (invert ? -1 : 1);
        complex<double> wlen(cos(ang), sin(ang));
        for (int i = 0; i < n; i += len) {
            complex<double> w(1);
            for (int j = 0; j < len / 2; j++) {
                complex<double> u = data[i + j], v = data[i + j + len / 2] * w;
                data[i + j] = u + v;
                data[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
    if (invert) {
        for (int i = 0; i < n; i++) data[i] /= n;
    }
}

Para casos donde el tamaño de la secuencia no es una potencia de 2, se puede recurrir al algoritmo de Bluestein, que transforma la DFT en una convolución lineal mediante una sustitución algebraica: \(jk = \binom{i+j}{2} - \binom{i}{2} - \binom{j}{2}\).

Optimización mediante el Principio de Transposición

El principio de transposición establece que si un algoritmo realiza una transformación lineal \(\mathbf{y} = \mathbf{M}\mathbf{x}\), existe un algoritmo para calcular \(\mathbf{z} = \mathbf{M}^T\mathbf{x}\) con la misma complejidad computacional. Dado que la matriz de la DFT es simétrica, su transpuesta es idéntica a la original.

Sin embargo, al descomponer la matriz FFT en factores elementales \(\mathbf{A} = \mathbf{A}_n \mathbf{A}_{n-1} \dots \mathbf{A}_1\), podemos obtener la transpuesta invirtiendo el orden de las operaciones y transponiendo cada factor: \(\mathbf{A}^T = \mathbf{A}_1^T \dots \mathbf{A}_n^T\).

Esta técnica permite derivar implementaciones que eliminan la necesidad de la etapa de inversión de bits, la cual suele ser costosa debido a los patrones de acceso a memoria no contiguos. En una Transformada de Teoría de Números (NTT) transpuesta, la estructura de mariposa se modifica de la siguiente manera:

// Implementación optimizada de NTT (DIF - Decimation in Frequency)
void forward_ntt(int* poly, int n, const int* roots) {
    for (int len = n, half = n >> 1; len >= 2; len >>= 1, half >>= 1) {
        for (int i = 0; i < n; i += len) {
            for (int j = 0; j < half; j++) {
                int left = poly[i + j];
                int right = poly[i + j + half];
                poly[i + j] = (left + right) % MOD;
                poly[i + j + half] = 1LL * (left - right + MOD) * roots[half + j] % MOD;
            }
        }
    }
}

// Implementación transpuesta (DIT - Decimation in Time)
void inverse_ntt(int* poly, int n, const int* inv_roots) {
    for (int len = 2, half = 1; len <= n; len <<= 1, half <<= 1) {
        for (int i = 0; i < n; i += len) {
            for (int j = 0; j < half; j++) {
                int left = poly[i + j];
                int right = 1LL * poly[i + j + half] * inv_roots[half + j] % MOD;
                poly[i + j] = (left + right) % MOD;
                poly[i + j + half] = (left - right + MOD) % MOD;
            }
        }
    }
    // Normalización final omitida por brevedad
}

Esta simetría entre los algoritmos de Diezmo en Frecuencia (DIF) y Diezmo en Tiempo (DIT) permite construir pipelines de procesamiento de señales mucho más eficientes, donde la salida de una etapa coincide con el orden esperado de la siguiente, minimizando los movimientos de datos innecesarios en la memoria caché.