¿Cuáles son algunas de las diferencias entre DFT y FFT que hacen que FFT sea tan rápido?

Estoy tratando de entender los FFT, esto es lo que tengo hasta ahora:

Para encontrar la magnitud de las frecuencias en una forma de onda, uno debe sondearlas multiplicando la onda por la frecuencia que están buscando, en dos fases diferentes (sin y cos) y promediando cada una. La fase se encuentra por su relación con los dos, y el código para eso es algo como esto:

//simple pseudocode
var wave = [...];                //an array of floats representing amplitude of wave
var numSamples = wave.length;
var spectrum = [1,2,3,4,5,6...]  //all frequencies being tested for.  

function getMagnitudesOfSpectrum() {
   var magnitudesOut = [];
   var phasesOut = [];

   for(freq in spectrum) {
       var magnitudeSin = 0;
       var magnitudeCos = 0;

       for(sample in numSamples) {
          magnitudeSin += amplitudeSinAt(sample, freq) * wave[sample];
          magnitudeCos += amplitudeCosAt(sample, freq) * wave[sample];
       }

       magnitudesOut[freq] = (magnitudeSin + magnitudeCos)/numSamples;
       phasesOut[freq] = //based off magnitudeSin and magnitudeCos
   }

   return magnitudesOut and phasesOut;
}

Para hacer esto para muchas frecuencias muy rápidamente, las FFT usan muchos trucos.

¿Cuáles son algunos de los trucos utilizados para hacer FFTs mucho más rápido que DFT?

PD: He intentado ver algoritmos FFT completos en la web, pero todos los trucos tienden a condensarse en un hermoso código sin mucha explicación. Lo que necesito primero, antes de que pueda entender todo, es una introducción a cada uno de estos cambios eficientes como conceptos.

Gracias.

La implementación ingenua de un DFT de punto es básicamente una multiplicación por una matrizEsto resulta en una complejidad de .N × N O ( N 2$N$$N \times N$$\mathcal{O}(N^2)$

Uno de los algoritmos más comunes de la Transformada rápida de Fourier (FFT) es el algoritmo FFT de decimación en el tiempo Radix-2 de Cooley-Tukey. Este es un enfoque básico de divide y vencerás.

Primero defina el "factor twiddle" como: donde es la unidad imaginaria, luego la DFT de viene dado por Si es par (y es un entero), la suma se puede dividir en dos sumas de la siguiente manera donde la primera suma se ocupa de las muestras pares de y el segundo con las muestras impares de . Definiendo y j≜

$$W_N \triangleq e^{-j\frac{2\pi}{N}}$$ X[k]x[n]X[k]= N - 1 ∑ n = 0 x[n$j \triangleq \sqrt{-1}$$X[k]$$x[n]$N $$X[k] = \sum_{n=0}^{N-1} x[n] \, W_N^{kn} \, .$$ $N$ X[k]= N / 2 - 1 ∑ n=0x[2n]W 2 k n N + N / 2 - 1 ∑ n=0x[2n+1]W k ( 2 n + 1 ) N x[n]x[n]x$\tfrac{N}{2}$ $$X[k] = \sum_{n=0}^{N/2-1} x[2n]W_N^{2kn} + \sum_{n=0}^{N/2-1} x[2n+1]W_N^{k(2n+1)}$$ $x[n]$$x[n]$x o [ n ] ≜ x [ 2 n + 1 $x_e[n] \triangleq x[2n]$$x_o[n] \triangleq x[2n+1]$ y usando el hecho de que

1. $W_N^{k(2n+1)} = W_N^{2kn}W_N^k$ , y
2. $W_N^{2kn} = W_{N/2}^{kn}$ ,

esto puede reescribirse como donde y son las transformaciones DFT de puntos de las muestras pares e impares de respectivamente. Así que simplemente transformamos un único DFT de puntos en dos DFT más pequeños de puntos. Esto reduce el costo computacional porque cuando .

$$\begin{align} X[k] &= \sum_{n=0}^{N/2-1} x_e[n] W_{N/2}^{kn} + W_N^k\sum_{n=0}^{N/2-1} x_o[n]W_{N/2}^{kn} \\ & = X_e[k] + W_N^k X_o[k] \end{align}$$ $X_e[k]$$X_o[k]$$\tfrac{N}{2}$$x[n]$$N$$\tfrac{N}{2}$ $$2 \left( \frac{N}{2} \right)^2 + N < N^2$$ $N > 2$

Luego podemos repetir el mismo proceso en esos dos DFT más pequeños. Este enfoque de divide y vencerás permite alcanzar la complejidad de , que es mucho mejor que el que teníamos con la ingenua implementación de DFT (como es ilustrado en gran medida por la respuesta de leftaroundabout ).$\mathcal{O}(N\log N)$$\mathcal{O}(N^2)$

http://nbviewer.jupyter.org/gist/leftaroundabout/83df89a7d3bdc24373ea470fb50be629

DFT, talla 16

FFT, talla 16

La diferencia en complejidad es bastante evidente a partir de eso, ¿no?

---

Así es como entiendo FFT.

En primer lugar, siempre pensaría en las transformaciones de Fourier principalmente como transformaciones de funciones continuas , es decir, un mapeo biyectivo . A la luz de esto, está claro que no puede ser realmente necesario ir al "nivel más profundo" y recorrer los elementos individuales , porque los "elementos individuales" son puntos únicos en la línea real, de los cuales hay infinitamente infinitos .$\operatorname{FT} : \mathcal{L}^2(\mathbb{R}) \to \mathcal{L}^2(\mathbb{R})$

Entonces, ¿cómo es que esta transformación todavía está bien definida? Bueno, es crucial que no funcione en el espacio de funciones generales sino solo en el espacio de funciones integrables (Lebesgue-, square-) . Ahora, esta integrabilidad no es una propiedad muy fuerte (mucho más débil que la diferenciabilidad, etc.), pero sí exige que la función se convierta en "localmente identificable con información contable". Tal descripción está dada por los coeficientes de una Transformada de Fourier de corto tiempo . ^†$\mathbb{R}\to\mathbb{C}$El caso más simple es que su función es continua y la divide en regiones tan pequeñas que es básicamente constante en cada una de ellas. Entonces cada uno de los STFT tiene más fuertemente un término cero. Si ignora los otros coeficientes (de todos modos en descomposición), cada dominio es solo un punto de datos. De todos estos coeficientes de límite de LF a corto plazo, podría tomar una transformada de Fourier discreta. De hecho, ¡eso es exactamente lo que haces al realizar cualquier FT en datos medidos del mundo real!

Sin embargo, los datos medidos no necesariamente corresponden a una cantidad física fundamental. Por ejemplo, cuando mides algo de intensidad de luz , en realidad solo estás midiendo la amplitud de una onda electromagnética cuya frecuencia es demasiado alta para ser muestreada con un ADC. Pero claramente también puede calcular el DFT de una señal de intensidad de luz muestreada, y de manera económica, a pesar de la frecuencia loca de la onda de luz.

Esto podría entenderse como la razón más importante por la que FFT es barata:

No se moleste en tratar de ver los ciclos de oscilación individuales desde el nivel más alto. En cambio, transforme solo información de alto nivel que ya haya sido preprocesada localmente.

Sin embargo, eso no es todo. Lo mejor de FFT es que aún le brinda toda la información que le daría un DFT completo . Es decir, toda la información que también obtendría al muestrear la onda electromagnética exacta de un haz de luz. ¿Se puede lograr esto transformando una señal de fotodiodo? - ¿Puedes medir la frecuencia de luz exacta a partir de eso?

Bueno, la respuesta es no, no puedes. Es decir, a menos que apliques trucos adicionales.
En primer lugar, debe medir al menos aproximadamente la frecuencia en los bloques de tiempo cortos. Bueno, eso es posible con un espectrógrafo. Pero solo es posible hasta una precisión de , una relación de incertidumbre típica ^‡ .$\Delta \nu = 1/{\Delta t}$

Al tener un período de tiempo más largo en general, también deberíamos poder reducir la incertidumbre de frecuencia. Y esto es realmente posible si se mide localmente no solo la frecuencia aproximada sino también la fase de la onda. Sabes que una señal de 1000 Hz tendrá exactamente la misma fase si la miras un segundo después. Mientras que una señal de 1000.5 Hz, aunque no se puede distinguir en la escala corta, habrá invertido la fase un segundo después.

Afortunadamente, esa información de fase se puede almacenar muy bien en un único número complejo. ¡Y así es como funciona FFT! Comienza con muchas pequeñas transformaciones locales. Estos son baratos, por un lado, obviamente, porque solo usan una pequeña cantidad de datos, pero en segundo lugar porque saben que, debido al corto lapso de tiempo, de todos modos no pueden resolver la frecuencia con mucha precisión, por lo que aún es asequible a pesar de que usted hacer muchas de esas transformaciones.

Sin embargo, estos también registran la fase y, a partir de eso, puede hacer que la resolución de frecuencia sea más exacta en el nivel superior. La transformación requerida es nuevamente barata, porque no se molesta con las oscilaciones de alta frecuencia sino solo con los datos de baja frecuencia preprocesados.

---

^† _{Sí, mi argumentación es un poco circular en este punto. Vamos a llamarlo recursivo y estamos bien ...}

^‡ _{Esta relación no es mecánica cuántica, pero la incertidumbre de Heisenberg tiene en realidad la misma razón fundamental.}

$W_N^{nk}$$e^{j2\pi \frac{nk}{N}}$

Tenga en cuenta la ruta que se muestra y la ecuación debajo muestra el resultado para el bin de frecuencia X (1), según lo dado por la ecuación de Robert.

Las líneas discontinuas no son diferentes a las líneas continuas solo para dejar en claro dónde están las uniones de suma.

esencialmente, al calcular el DFT ingenuo directamente de la sumatoria:

$$X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{j 2 \pi \frac{nk}{N}}$$

$N$$e^{j 2 \pi \frac{nk}{N}}$$N$$N-1$$X[k]$$k$$X[k+1]$

1. entonces la FFT conserva algunos datos intermedios.
2. la FFT también hará uso de factorizar un poco el factor twiddle para que el mismo factor pueda usarse para una combinación intermedia de datos.

Soy una persona visual Prefiero imaginar el FFT como un truco matricial en lugar de un truco de suma.

Para explicar a un alto nivel:

Un ingenuo DFT calcula cada muestra de salida de forma independiente y utiliza cada muestra de entrada en cada cálculo (algoritmo clásico N²).

Una FFT común usa simetrías y patrones en la definición de DFT para hacer el cálculo en "capas" (log N capas), cada capa con un requisito de tiempo constante por muestra creando un algoritmo N log N.

Más detalles:

Una forma de visualizar estas simetrías es mirar el DFT como una entrada de matriz 1 × N multiplicada por una matriz NxN de todos sus exponenciales complejos. Comencemos con el caso "radix 2". Vamos a dividir las filas pares e impares de la matriz (correspondientes a las muestras de entrada pares e impares) y considerarlas como dos multiplicaciones de matriz separadas que se suman para obtener el mismo resultado final.

Ahora observe estas matrices: en la primera, la mitad izquierda es idéntica a la mitad derecha. En el otro, la mitad derecha es la mitad izquierda x −1. Esto significa que solo tenemos que usar la mitad izquierda de estas matrices para la multiplicación y crear la mitad derecha a bajo costo multiplicando por 1 o −1. Luego, observe que la segunda matriz difiere de la primera matriz por factores que son iguales en cada columna, por lo que podemos factorizar eso y multiplicarlo en la entrada, por lo que ahora las muestras pares e impares usan la misma matriz, pero requieren un multiplicador primero. Y el paso final es observar que esta matriz N / 2 × N / 2 resultante es idéntica a una matriz N / 2 DFT y podemos hacer esto una y otra vez hasta llegar a una matriz 1 × 1 donde la DFT es una función de identidad.

Para generalizar más allá de la raíz 2, puede ver dividir cada tercera fila y ver tres trozos de columnas, o cada 4ta, etc.

En el caso de entradas de tamaño principal, existe un método para poner a cero correctamente, FFT y truncar, pero eso está más allá del alcance de esta respuesta.

Ver: http://whoiskylefinn.com/MatrixFFT.html

El DFT multiplica la matriz de fuerza bruta N ^ 2.

Los FFT hacen trucos ingeniosos, explotando las propiedades de la matriz (desgeneralizando la matriz de multiplicación) para reducir el costo computacional.

Veamos primero un pequeño DFT:

W = fft (ojo (4));

x = rand (4,1) + 1j * rand (4,1);

X_ref = fft (x);

X = W * x;

afirmar (max (abs (X-X_ref)) <1e-7)

Genial, por lo que podemos sustituir la llamada de MATLAB a la biblioteca FFTW por una pequeña multiplicación de matriz 4x4 (compleja) al llenar una matriz desde la función FFT. Entonces, ¿cómo se ve esta matriz?

N = 4

Wn = exp (-1j * 2 * pi / N),

f = ((0: N-1) '* (0: N-1))

f =

 0     0     0     0
 0     1     2     3
 0     2     4     6
 0     3     6     9

W = Wn. ^ F

W =

1 1 1 1

1 -i -1 i

1 -1 1 -1

1 i -1 -i

Cada elemento es +1, -1, + 1j o -1j. Obviamente, esto significa que podemos evitar multiplicaciones complejas completas. Además, la primera columna es idéntica, lo que significa que estamos multiplicando el primer elemento de x una y otra vez por el mismo factor.

Resulta que los productos de tensor de Kronecker, los "factores twiddle" y una matriz de permutación donde el índice cambia según la representación binaria invertida son compactos y ofrecen una perspectiva alternativa sobre cómo se calculan las FFT como un conjunto de operaciones de matriz dispersas.

Las líneas a continuación es una simple Decimación en Frecuencia (DIF) radix 2 hacia adelante FFT. Si bien los pasos pueden parecer engorrosos, es conveniente reutilizarlos para FFT directa / inversa, radix4 / split-radix o decimación en el tiempo, a la vez que es una representación justa de cómo las FFT en el lugar tienden a implementarse en el mundo real, Yo creo.

N = 4;

x = randn (N, 1) + 1j * randn (N, 1);

T1 = exp (-1j * 2 * pi * ([ceros (1, N / 2), 0: (N / 2-1)]). '/ N),

M0 = corona (ojo (2), fft (ojo (2))),

M1 = kron (fft (ojo (2)), ojo (2)),

X = bitrevorder (x. '* M1 * diag (T1) * M0),

X_ref = fft (x)

afirmar (max (abs (X (:) - X_ref (:))) <1e-6)

CF Van Loan tiene un gran libro sobre este tema.

Si quieres beber de la manguera de bomberos de la sabiduría, te sugiero:

"Transformaciones rápidas - Algoritmos, análisis, aplicaciones" por Douglas F. Elliott, K. Ramamohan Rao

Cubre FFT, Hartley, Winograd y aplicaciones.

Un punto fuerte es que muestra cómo la FFT es un conjunto de factorizaciones de matriz dispersas con un orden de inversión de bits.