FFT: Qué es realmente?

Cursando en la facultad, escuché varias veces hablar sobre la FFT e incluso la utilicé en una materia. Si bien tuvimos una preparación razonable en cuanto a las series y transformadas de Fourier, en ningún momento nos explicaron nada acerca de este algoritmo.

Esto me llevó naturalmente a suponer que debía ser un algoritmo complejo y, en consecuencia, a sorprenderme al ver que la idea básica tras él es muy simple y bien conocida por los que estudiamos Informática: “divide & conquer”.

Para empezar consideremos primero la definición de la transformada discreta de Fourier:

X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n} \quad \quad k = 0, \dots, N-1

Si lo implementáramos directamente obtendríamos un algoritmo como el siguiente:

def naive_dft(X):
    n = len(X)
    basic_root = math.e ** (-2.0 * math.pi * 1j / n)
    K = [sum(x * root ** i for i, x in enumerate(X))
         for root in (basic_root ** j for j in range(n))]
    return K

Como tenemos esencialmente una iteración completa sobre los N elementos de x_n para calcular cada uno de los N elementos X_k, la complejidad del algoritmo será \Theta(N^2) (puede recortarse el número de operaciones aprovechando ciertas identidades pero el número de operaciones seguirá siendo cuadrático).

Pero la expresión dada anteriormente para la DFT puede reescribirse del siguiente modo (asumiendo que N es par y tomando M = N/2):

X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n}

X_k = \sum_{m=0}^{M-1} (x_{2m} e^{-\frac{2 \pi i}{2M} k (2m)} + x_{2m+1} e^{-\frac{2 \pi i}{2M} k (2m + 1)})

X_k = \sum_{m=0}^{M-1} x_{2m} e^{-\frac{2 \pi i}{M} k m} + \sum_{m=0}^{M-1} x_{2m+1} e^{-\frac{2 \pi i}{M} k m} e^{-\frac{2 \pi i}{2M} k}

X_k = \sum_{m=0}^{M-1} x_{2m} e^{-\frac{2 \pi i}{M} k m} + e^{-\frac{2 \pi i}{N} k} \left( \sum_{m=0}^{M-1} x_{2m+1} e^{-\frac{2 \pi i}{M} k m} \right)

Puede observarse que obtuvimos dos expresiones muy similares a la expresión original para la DFT, pero tratando con secuencias con una longitud de un medio de la original. La única diferencia es que en la secuencia original podíamos garantizar que k < N, mientras que k puede ser mayor que M. Pero como

e^{-\frac{2 \pi i}{M} k m} = e^{-\frac{2 \pi i}{M} (k - a M) m} e^{-\frac{2 \pi i}{M} a M m}

e^{-\frac{2 \pi i}{M} k m} = e^{-\frac{2 \pi i}{M} (k - a M) m} e^{-2 \pi i a m}

e^{-\frac{2 \pi i}{M} k m} = e^{-\frac{2 \pi i}{M} (k - a M) m},

los valores de las DFT parciales simplemente se repetirán para K = M, \dots, N-1.

Expresando esta definición en Python tendremos:

def fft(X):
    n = len(X)
    assert n & (n - 1) == 0 # n must be a power of two
    if n == 1:
        return X
    even_fft = fft(X[0::2])
    odd_fft = fft(X[1::2])
    w = math.e ** (-2.0 * math.pi * 1j / n)
    first_half = [ke + (w ** i ) * ko
                  for i, ke, ko in zip(itertools.count(),
                                       even_fft, odd_fft)]
    second_half = [ke - (w ** i ) * ko
                   for i, ke, ko in zip(itertools.count(),
                                        even_fft, odd_fft)]
    return first_half + second_half

donde el signo “-” entre las partes par e impar en la segunda mitad es debido a que

e^{-\frac{2 \pi i}{N} k} = e^{-\frac{2 \pi i}{2 M} (k - M)} e^{-\frac{2 \pi i}{2 M} M} = e^{-\frac{2 \pi i}{2 M} (k - M)} e^{-\pi i} = -e^{-\frac{2 \pi i}{2 M} (k - M)}.

Obviamente esta implementación en Python no es competitiva con FFTW 😀 pero se observan claramente las diferencias que el comportamiento asintótico hace en la performance de los algoritmos. En este caso tenemos una complejidad total lineal en cada nivel de recursión, dando una complejidad general \Theta(n \log n).

Sería interesante comparar la performance entre una implementación en C del algoritmo DFT básico y esta implementación en Python, de modo de ver para que valor de n el comportamiento asintótico pasa a dominar sobre las constantes…

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s