📲 🎐 🏺 Transformada de Fourier. El rápido y el furioso 🍖 💋 🍞

A menudo, al desarrollar algoritmos, nos topamos con el límite de la complejidad computacional, que, al parecer, es imposible de superar. La transformada de Fourier tiene complejidad

O (n^{2})

$O (n ^ 2)$ , y una versión rápida, propuesta alrededor de 1805 por la Casa ¹ (y reinventada en 1965 por James Cooley y John Tukey)

O (n l o g (n))

$O (nlog (n))$ . En este artículo quiero mostrarle que puede obtener los resultados de conversión en tiempo lineal

O (n)

$O (n)$ o incluso lograr dificultad constante

O (1)

$O (1)$ bajo ciertas condiciones que se encuentran en problemas reales.
Análisis de Fourier

Cuando me enfrenté a la tarea de escribir un programa para analizar las funciones de transferencia de sistemas de sonido en tiempo real, como todos los demás, primero recurrí a la conversión rápida. Todo estaba bien, pero con el gran tamaño de la ventana de tiempo, la carga de la CPU se volvió indecentemente grande y había que hacer algo. Se decidió pausar y estudiar la transformación nuevamente, y al mismo tiempo buscar formas de resolver el problema. Volvamos a la transformación original de Joseph Fourier ² :

f (x) = s u m l i m i t s_{- i n f t y}^{+ i n f t y} c_{k} e^{2 p i i k x / T} c_{k} = f r a c 1 T i n t l i m i t s_{0}^{T} f (x) e^{- 2 p i i k x / T} d x

$f (x) = \ sum \ limits _ {- \ infty} ^ {+ \ infty} c_ke ^ {2 \ pi ikx / T} \\ c_k = \ frac {1} {T} \ int \ limits_0 ^ Tf ( x) e ^ {- 2 \ pi ikx / T} dx$

Analizaremos cuidadosamente lo que está sucediendo aquí. Cada valor de salida en el dominio de frecuencia

c_{k}

$c_k$ es la suma de todos los valores de entrada de la señal

f (x)

$f (x)$ multiplicado por

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ . Para hacer cálculos, necesitamos revisar todos los datos de entrada para cada valor de salida, es decir. para cumplir esos

n^{2}

$n ^ 2$ operaciones

Deshacerse de n

Permítame recordarle que inicialmente la tarea era analizar datos de sonido en tiempo real. Para hacer esto, la ventana de tiempo seleccionada (esencialmente un búfer) de tamaño N se llena con datos con una frecuencia f _d correspondiente a la frecuencia de muestreo. Con el período T, los datos de entrada se convierten de la ventana de tiempo a la ventana de frecuencia. Si observa los números reales, entonces N varía de 2 ¹⁴ (16 384) a 2 ¹⁶ (65 536) muestras (los valores se heredan de la FFT, donde el tamaño de la ventana debe ser una potencia de dos). Tiempo T = 80ms (12.5 veces por segundo), lo que le permite ver los cambios de manera conveniente y no sobrecargar la CPU y la GPU. La frecuencia de muestreo f _{d es} estándar y es de 48 kHz. Calculemos cuántos datos en la ventana de tiempo cambian entre dimensiones. Durante el tiempo T, ingresa al búfer

48000 * 0.08 = $ 384

$48000 * 0.08 = $ 384$ muestras Por lo tanto, solo del 5% al 23% de los datos se actualizan en la ventana. En el peor de los casos, el 95% (¡y en el mejor de los casos el 73%, que también es mucho!) De las muestras procesadas caerán en la conversión una y otra vez, a pesar de que ya se han procesado en iteraciones anteriores.

El atento lector en ese momento levantará la mano y dirá: "espera, pero ¿qué pasa con el coeficiente?

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ ? Después de todo, con cada nueva transformación, los mismos datos se ubicarán en las nuevas posiciones de la serie y, como resultado, ¿tienen coeficientes diferentes? Por cada cinco para su cuidado, recordemos un detalle importante en la transformación que a menudo se olvida. En el estudio de los valores de las funciones.

f (t)

$f (t)$ durante el intervalo de 0 a t, la función se considera periódica, lo que le permite desplazar sin dolor la función hacia la izquierda o hacia la derecha en el tiempo. Como resultado, tenemos todo el derecho de no insertar un nuevo valor al final y eliminar el valor anterior desde el principio, sino de reemplazar cíclicamente los datos en el búfer.

Para mayor claridad, puede escribir en forma tabular cómo cambiará el búfer:

t = 0	f (0)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 1	f (10)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 2	f (10)	f (11)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 3	f (10)	f (11)	f (12)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 4	f (10)	f (11)	f (12)	f (13)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)

Puede escribir cómo la transformación en el tiempo cambia de t ₁ a t ₂ :

F_{t} = F_{t - 1} + D e l t a F D e l t a F : D e l t a c_{k} = f r a c 1 T i n t l i m i t s_{t_{1}}^{t_{2}} (f_{t} (x) - f_{t - 1} (x)) e^{- 2 p i i k x / T} d x

$F_t = F_ {t-1} + \ Delta F \\ \ Delta F: \ Delta c_k = \ frac {1} {T} \ int \ limits_ {t_1} ^ {t_2} (f_t (x) - f_ { t-1} (x)) e ^ {- 2 \ pi ikx / T} dx$

Valor

F_{t - 1} (x)

$F_ {t-1} (x)$ es el resultado de la conversión anterior y la complejidad del cálculo

D e l t a f (x)

$\ Delta f (x)$ no depende del tamaño de la ventana de tiempo y, por lo tanto, es constante. Como resultado, la complejidad de la conversión será

O (n)

$O (n)$ ^* porque todo lo que nos queda es pasar por la ventana de frecuencia una vez y aplicar los cambios para las muestras T que han cambiado a lo largo del tiempo. También me gustaría llamar su atención sobre el hecho de que las probabilidades

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ se puede calcular por adelantado, lo que proporciona una ganancia adicional en productividad, y solo quedan dos operaciones en el ciclo: restar números reales y multiplicar un número real por uno complejo, en la práctica ambas operaciones son simples y económicas.

Para completar la imagen, solo queda indicar el estado inicial, pero aquí todo es simple:

F_{0} (x) = 0

$F_0 (x) = 0$

* - por supuesto, la complejidad final de toda la transformación seguirá siendo así

O (n^{2})

$O (n ^ 2)$ , pero se ejecutará gradualmente, durante n iteraciones, mientras se actualiza el búfer.

O (n)

$O (n)$ - esta es la complejidad de actualizar los datos, pero esto es exactamente lo que necesitamos (cuando usamos el FFT, la complejidad de cada transformación

O (n l o g (n))

$O (nlog (n))$ )

Pero, ¿y si cavas más profundo? O deshacerse de la segunda n

Quiero hacer una reserva de inmediato para que los siguientes pasos sean aplicables solo si no planea realizar la transformación inversa para el resultado (para corregir la señal u obtener una respuesta de impulso). Para comenzar, quiero recordarles que como resultado de la conversión, obtenemos una matriz de datos simétrica, que nos permite reducir a la mitad el número de conversiones de inmediato.

Ahora analicemos el conjunto de datos resultante, dadas las condiciones del problema. Tenemos un conjunto de números complejos, cada uno de los cuales describe la amplitud y la fase de las oscilaciones a una frecuencia particular. La frecuencia se puede determinar mediante la fórmula:

f [j] = j f r a c f d N

$f [j] = j \ frac {fd} {N}$ para

j < f r a c N 2

$j <\ frac {N} {2}$ . Evaluemos el paso de la ventana de frecuencia en nuestros datos:

D e l t a f = f r a c f d N

$\ Delta f = \ frac {fd} {N}$ Para N = 2 ¹⁴ : 2.93 Hz (y para 2 ¹⁶ : 0.73 Hz). Por lo tanto, en el rango de 1 kHz a 2 kHz obtenemos 341 resultados. Intente evaluar de forma independiente cuántos datos estarán en el rango de 8 kHz a 16 kHz para N = 65536. Mucho, ¿verdad? Mucho! ¿Necesitamos tantos datos? Por supuesto, en los problemas de visualización de las características de frecuencia de los sistemas de sonido, la respuesta es no. Y, por otro lado, para el análisis en la región de baja frecuencia, un pequeño paso es muy útil. No olvide que todavía hay un cronograma por delante que debería estos volúmenes (

f r a c N 2

$\ frac {N} {2}$ ) se convierten en una forma legible para humanos (promediación, spline o suavizado) y se muestran en la pantalla. Y a altas frecuencias, incluso teniendo una pantalla 4K y mostrando el gráfico en modo de pantalla completa con el eje de frecuencia logarítmica, el tamaño del paso rápidamente resulta ser mucho menor que 1 píxel.

Por experiencia, puede descubrir que es suficiente tener solo 48 puntos por octava, y tener los datos un poco más suaves y promediados, sugiero detenerse en 96. En el rango de frecuencia de audio de 20 Hz a 20 kHz, es fácil contar solo 10 octavas: 20, 40, 80 , 160, 320, 640, 1280, 2560, 5120, 10240, 20480, cada uno de los cuales puede dividirse en un número determinado de subbandas (no olvide que la partición debe hacerse geométricamente y no aritméticamente), por lo tanto, es más que suficiente para realizar la conversión solo para 960 frecuencias para obtener Performan que en 16 ... 65 veces más pequeño que la versión original.

Por lo tanto, combinando ambos enfoques, obtenemos la complejidad constante del algoritmo de actualización de datos

O (1)

$O (1)$ .

Miel cuadrada y una cucharada de alquitrán

Ahora podemos decirlo con seguridad desde la complejidad

O (n^{2})

$O (n ^ 2)$ llegamos a la complejidad

O (1)

$O (1)$ usando dos simples trucos de la vida:

Después de analizar el problema, notamos que los datos se agregan gradualmente, y el período de actualización completa de la ventana de tiempo es mucho mayor que el período de transformaciones y luego calculamos la diferencia de la transformación de Fourier.
pasado del paso aritmético en la ventana de frecuencia a limitado solo por los valores especificados, lo que puede reducir drásticamente el número de conversiones.

Pero, por supuesto, la vida sería realmente un cuento de hadas, si no uno, pero. La aplicación de estos dos enfoques nos permitió descargar realmente la CPU para adivinar que calcula la transformación de Fourier y muestra los resultados en la pantalla incluso con

N = 2^{20}

$N = 2 ^ {20}$ Fue dificil. Pero el castigo no se hizo esperar cuando sus señales en realidad no son periódicas (y esto es necesario para obtener los resultados de conversión correctos) y no es posible seleccionar el tamaño de ventana apropiado, se hace necesario usar varias funciones de ventana, lo que ya no le permite usar el primer paso por completo. La práctica ha demostrado que el uso de funciones de ventana es crítico en el estudio de señales con una frecuencia inferior a

0, 1 f_{d}

$0,1f_d$ . A altas frecuencias, el número de períodos que caen en la ventana de tiempo atenúa significativamente las distorsiones que surgen de la presencia de una brecha de primer orden (entre f (0) yf (N-1)) en la función original.

Al final, también rechacé el segundo paso y volví a la FFT, porque la ganancia en esta tarea ya era pequeña.

En conclusión

El primer enfoque se puede aplicar si sus datos tienen un carácter periódico pronunciado y deben analizarse a lo largo del tiempo utilizando una ventana de tiempo grande, que, recuerdo, no tiene que ser de grado 2, es decir. Cualquier número natural.
El segundo enfoque es aplicable (incluso teniendo en cuenta las funciones de la ventana) si solo se analiza un cierto conjunto pequeño de frecuencias en los datos.

Por desgracia, para mí en este problema, esto siguió siendo solo un pequeño entretenimiento matemático, pero espero que lo inspire a estudiar otros algoritmos en vacaciones en términos de cambios en los datos de entrada a lo largo del tiempo :)

Literatura

Imagen tomada del manga por Michio Shibuya. “EXCITING MATHEMATICS. Análisis de Fourier

Transformada de Fourier. El rápido y el furioso

Deshacerse de n

Pero, ¿y si cavas más profundo? O deshacerse de la segunda n

Miel cuadrada y una cucharada de alquitrán

En conclusión

Literatura

More articles: