6 minutos de lectura
Transformada discreta de Fourier (DFT) en JavaScript
Índice:
Este post explica cómo implementar la "Transformada discreta de Fourier" (DFT).
La implementación está basada en los siguiente recursos:
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.
Puntos aparte, existen versiones más preparadas de esta implementación como:
Sin embargo, esta implementación está más enfocada en el código.
Empecemos...
Transformada
La fórmula para la DFT es la siguiente:
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
arraysse 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 unarrayde números.k: Es la posición que representa el índice en la transformada.Xk: Representa la transformada enk.k: Es la posición que representa el índice en X (la serie de datos).Xn: Es el valor de la serieXen la posiciónn.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
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.
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 👇