moloxe.io / blog / transformada-discreta-de-fourier-dft-en-javascript

Marzo 22, 2024

6 minutos de lectura

Transformada discreta de Fourier (DFT) en JavaScript


Este post explica cómo implementar la "Transformada discreta de Fourier" (DFT).

La implementación está basada en los siguiente recursos:

https://www.youtube.com/watch?v=MY4luNgGfms

En el video la transformada discreta adelanta los cálculos de la inversa y se ordenan las transformadas por amplitudes. En esta implementación se hará por separado.

Sin embargo, esta implementación está más enfocada en el código.

Empecemos...


Transformada

La fórmula para la DFT es la siguiente:

discrete fourier transform

TRANQUI, la intención es transformar eso a código. Es complicada, no entiendo a la perfección el funcionamiento (honestidad). Pero, si se divide en partes pequeñas, cobra sentido.

No se puede explicar de izquierda a derecha, pero sí de adentro hacia afuera.

Iré rápido. Si la intención es dibujar caritas, no veo tan grave saltarse algunas cosas.

e**( -(2 * PI * i) / N * k * n)

Gracias a Euler, esta parte equivale a:

cos( (2 * PI * k * n) / N ) - i * sin( (2 * PI * k * n) / N )

Ahora, esta es una fórmula que trabaja con ángulos. Precisamente es la parte que se repite, por lo que sirve más darle un nombre:

const angle = (2 * PI * k * n) / N

Finalmente se reduce a:

cos(angle) - i * sin(angle)

IMPORTANTE

Por qué son necesarios todos estos cálculos? El problema es i, es necesario mantenerlo para poder salir del plano complejo. Si no entiendes esto, no te preocupes! es como teletransportarse, llegar de un lugar a otro de una forma inexplicable: sqrt(-1).

Con estos cambios, la fórmula sería equivalente a:

Xk = Σ Xn * (cos(angle) - i * sin(angle))

2 cosas en cuenta:

  • Aún no me acerco 100% a código.
  • La Σ sería de 0 a N-1, ya no desde 1, ya que los arrays se indexan desde 0.

Si llegaste hasta este punto, puedes sentirte más relajado. Ahora es más sencillo explicar las:

Variables

  • X: Es la serie de datos, en otros términos es un array de números.
  • k: Es la posición que representa el índice en la transformada.
  • Xk: Representa la transformada en k.
  • k: Es la posición que representa el índice en X (la serie de datos).
  • Xn: Es el valor de la serie X en la posición n.
  • angle: Esta parte podría complicar la explicación, de momento es eso, un ángulo.
  • i: Es la raíz de -1 y de mis problemas.

Entrada

La entrada será la variable X. En el código X[n] representará a Xn.

Salida

Aquí es donde puede haber confusiones, no se puede llamar X a la entrada y a la salida. Se representará los valores de Xk con transform[k].

transform[k] = Σ X[n] * (cos(angle) - i * sin(angle))


Algoritmo

Empezemos con angle:

const angle = (2 * PI * k * n) / N

Algo a notar es que k y n son índices, el código sería el siguiente:

function dft(X) {
  const N = X.length
  for (let k = 0; k < N; k++) {
    for (let n = 0; n < N; n++) {
      const angle = (2 * PI * k * n) / N
    }
  }
}

Ahora la parte complicada (compleja 👉 👉). Un número complejo se puede representar de la siguiente forma: a + bi, donde a es la parte real y bi es la parte imaginaria.

Entonces aquí está la magia: Se almacena el producto (la parte imaginaria) en una variable llamada bi. Se va suponer que es el producto de b * i pero en una sola variable.

a = Σ X[n] * cos(angle)
bi = Σ -(X[n] * sin(angle))

Este sería el código hasta calcular a y bi:

function dft(X) {
  const N = X.length
  for (let k = 0; k < N; k++) {
    // a + b*i
    let a = 0
    let bi = 0
    for (let n = 0; n < N; n++) {
      const angle = (2 * PI * k * n) / N
      a += X[n] * cos(angle)
      bi += -(X[n] * sin(angle))
    }
  }
}

Ahora faltaría almacenar los valores en transform[k].

Ya no es necesario el término de la sumatoria Σ, ahora está representado por las sumas en segundo bucle:

function dft(X) {
  const transform = []
  const N = X.length
  for (let k = 0; k < N; k++) {
    // transform[k] = Σ X[n] * cos(angle) - i * sin(angle)
    // transform[k] = a + b*i
    let a = 0
    let bi = 0
    for (let n = 0; n < N; n++) {
      const angle = (2 * PI * k * n) / N
      a += X[n] * cos(angle)
      bi += -(X[n] * sin(angle))
    }
    transform.push({
      a,
      bi,
    })
  }
  return transform
}

Listo! Eso es casi todo.

Spoiler: Se puede hacer 3d.


Transformada inversa

inverse discrete fourier transform

No se explicará la implementación de esta formula, es prácticamente el mismo procedimiento. Se los dejo de tarea, pero este sería el código:

function idft(transform) {
  const iTransform = []
  const N = transform.length
  for (let n = 0; n < N; n++) {
    // X[n] = Σ transform[k] * cos(angle) + i * sin(angle)
    let x = 0
    let y = 0
    for (let k = 0; k < N; k++) {
      const { a, bi } = transform[k]
      const f = k
      const A = sqrt(a ** 2 + bi ** 2)
      const hPhase = atan2(bi, a)
      const angle = -f * ((2 * PI * n) / N) + hPhase
      x += (A * cos(angle)) / N
      y += (A * sin(angle)) / N
    }
    iTransform.push({ x, y })
  }
  return iTransform
}

Nota

Si se desea trabajar con varios ejes (y sí se desea) es necesario rotarlo en el eje correspondiente. Lo único necesario sería pasar un parámetro rotation a la función idft y sumarlo al ángulo:

const angle = -f * ((2 * PI * n) / N) + hPhase + rotation

Simple! 🥲


Epicycle

Así se llaman los círculos que dibujan el rostro. Para lograr esto se tiene que modificar la función inversa (idft) para no calcular todo el rostro, solo una parte evaluada e ir trazando cada epiciclo.

function epicycle(x, y, transform, rotation) {
  const N = transform.length;
  for (let k = 0; k < N; k++) {
    const prevX = x;
    const prevY = y;
    const { a, bi } = transform[k];
    const f = k;
    const A = sqrt(a ** 2 + bi ** 2);
    const hPhase = atan2(bi, a);
    x += (A * cos(f * time + hPhase + rotation)) / N;
    y += (A * sin(f * time + hPhase + rotation)) / N;
    // Esto dibuja cada epicycle
    const radius = A / N;
    stroke("#ffffff50");
    noFill();
    ellipse(prevX, prevY, radius * 2);
    line(prevX, prevY, x, y);
  }
  return { x, y };
}

Es posible tener 4 transformadas, 2 para el plano XY y uno para el plano ZY. Podría ser cualquiera, pero no voy a mentir que esto explota la mente.

Puedes mover la perspectiva arrastrando con el mouse/touch, así se pueden percibir mejor las 3 dimensiones.

gif-block

Peor que probar cosas cuestionables

También es posible utilizar sólo 3 transformadas, 2 transformadas se proyectarían en el plano XY y 1 en el plano XZ.


Usando el plano complejo

Es posible ajustar la trayectoria de la transformada para realizar trazos 2D sin la necesidad de realizar 2 transformadas. Para ello la serie tiene que ser de números complejos, en el que la parte real es un eje y la parte imaginaria el otro eje.

El único cambio es en el cálculo de la transformada, ya que en la inversa ya se calcula un eje y.

function dft(X) {
  const fourier = [];
  const N = X.length;
  for (let k = 0; k < N; k++) {
    const Xk = { a: 0, bi: 0, f: k };
    for (let n = 0; n < N; n++) {
      const angle = (2 * PI * k * n) / N;
      const a = X[n].x;
      const bi = X[n].y;
      const c = cos(angle);
      const di = -sin(angle);
      // FOIL: (a + bi) * (c + di) = (ac - bd) + (ad + bc)i
      Xk.a += a * c - bi * di;
      Xk.bi += a * di + bi * c;
    }
    Xk.A = sqrt(Xk.a ** 2 + Xk.bi ** 2);
    fourier.push(Xk);
  }

  // Ordenar las amplitudes ayuda a la visualización de los epiciclos
  fourier.sort((f1, f2) => f2.A - f1.A);

  return fourier;
}

Ejemplo del rostro utilizando una sola transformada:


FFT

Fast Fourier Transform permite calcular la transformada con menos operaciones (reduciendo la complejidad de O(N^2) a O(N log(N)).

Un caso específico en el que es útil (que hay miles) es el interactuar con la formula en tiempo real.

En este ejemplo se pasa una imagen a un dominio de frecuencias y se realiza un shift para poner las frecuencias más "relevantes" al centro. Si solo se consideran las del centro, se obtiene una imagen comprimida (una representación con menos datos).

El algoritmo de la FFT tiene una implementación diferente, en este ejemplo se usó la de:


Fin

Al principio buscaba implementaciones para entenderla, pero no lo logré hasta tratar directamente con la fórmula. Lo cierto es que la expresión resume bastante todos los cálculos, así son las mates.

Quisiera probarla con otras formas, dimensiones o cualquier idea. Si tienes alguna, compártela en los comentarios 👇