Mejorando la velocidad de implementación de FFT

Question

Soy un principiante en programación y actualmente estoy tratando de trabajar en un proyecto que requiera la implementación de la Transformada rápida de Fourier.

que he conseguido hasta ahora poner en práctica lo siguiente:

cualquier persona tiene alternativas y sugerencias para mejorar la velocidad del programa sin renunciar a la precisión.

short FFTMethod::FFTcalc(short int dir,long m,double *x,double *y) 
{ 
long n,i,i1,j,k,i2,l,l1,l2; 
double c1,c2,tx,ty,t1,t2,u1,u2,z; 

/* Calculate the number of points */ 
n = 1; 
for (i=0;i<m;i++) 
    n *= 2; 

/* Do the bit reversal */ 
i2 = n >> 1; 
j = 0; 
for (i=0;i<n-1;i++) { 
    if (i < j) { 
    tx = x[i]; 
    ty = y[i]; 
    x[i] = x[j]; 
    y[i] = y[j]; 
    x[j] = tx; 
    y[j] = ty; 
    } 
    k = i2; 
    while (k <= j) { 
    j -= k; 
    k >>= 1; 
    } 
    j += k; 
} 

/* Compute the FFT */ 
c1 = -1.0; 
c2 = 0.0; 
l2 = 1; 
for (l=0;l<m;l++) { 
    l1 = l2; 
    l2 <<= 1; 
    u1 = 1.0; 
    u2 = 0.0; 
    for (j=0;j<l1;j++) { 
    for (i=j;i<n;i+=l2) { 
     i1 = i + l1; 
     t1 = u1 * x[i1] - u2 * y[i1]; 
     t2 = u1 * y[i1] + u2 * x[i1]; 
     x[i1] = x[i] - t1; 
     y[i1] = y[i] - t2; 
     x[i] += t1; 
     y[i] += t2; 
    } 
    z = u1 * c1 - u2 * c2; 
    u2 = u1 * c2 + u2 * c1; 
    u1 = z; 
    } 
    c2 = sqrt((1.0 - c1)/2.0); 
    if (dir == 1) 
    c2 = -c2; 
    c1 = sqrt((1.0 + c1)/2.0); 
    } 

/* Scaling for forward transform */ 
if (dir == 1) { 
    for (i=0;i<n;i++) { 
     x[i] /= n; 
     y[i] /= n; 
    } 
} 


    return(1); 
} 

Answer 1

Recientemente encontré este excelente PDF en el Construction of a high performance FFTs de Eric Postpischil. Después de haber desarrollado varias FFT sé lo difícil que es competir con las bibliotecas comerciales. Créame que lo está haciendo bien si su FFT es solo 4 veces más lenta que Intel o FFTW, ¡no 40 veces! Sin embargo, puedes competir, y así es cómo.

Para resumir ese artículo, el autor afirma que las FFT de Radix2 son simples pero ineficientes, la construcción más eficiente es la FFT radix4. Un método aún más eficiente es el Radix8, sin embargo, esto a menudo no encaja en los registros de una CPU, por lo que se prefiere Radix4.

Las FFT se pueden construir en etapas, por lo que para calcular una FFT de 1024 puntos se pueden realizar 10 etapas de la FFT Radix2 (como 2^10 - 1024) o 5 etapas de la FFT Radix4 (4^5 = 1024) . Incluso puede calcular una FFT de 1024 puntos en etapas de 8 * 4 * 4 * 4 * 2 si así lo desea.Menos etapas significan menos lecturas y escrituras en la memoria (el cuello de botella para el rendimiento de FFT es el ancho de banda de la memoria) por lo tanto, elegir dinámicamente radix 4, 8 o superior es una necesidad. La etapa Radix4 es particularmente eficiente ya que todos los pesos son 1 + 0i, 0 + 1i, -1 + 0i, 0-1i y el código de mariposa Radix4 se puede escribir para que quepa completamente en la memoria caché.

En segundo lugar, cada etapa en la FFT no es la misma. La primera etapa, los pesos son todos iguales a 1 + 0i. no tiene sentido calcular este peso e incluso multiplicarlo ya que es un complejo multiplicado por 1, por lo que la primera etapa se puede realizar sin pesos. La etapa final también se puede tratar de manera diferente y se puede usar para realizar la decimación en el tiempo (inversión de bit). El documento de Eric Postpischil cubre todo esto.

Los pesos se pueden precalcular y almacenar en una tabla. Los cálculos Sin/Cos demoran alrededor de 100-150 ciclos cada uno en el hardware x86, por lo que precomputarlos puede ahorrar un 10-20% del tiempo total de cálculo, ya que el acceso a la memoria es, en este caso, más rápido que los cálculos de la CPU. El uso de algoritmos rápidos para calcular sincos de una sola vez es particularmente beneficioso (tenga en cuenta que cos es igual a sqrt (1.0 - seno senoidal), o usando tablas de búsqueda, cos es solo un cambio de fase de seno).

Finalmente, una vez que tenga su implementación aerodinámica FFT puede utilizar la vectorización SIMD para calcular 4x coma flotante o 2x operaciones de punto flotante doble por ciclo dentro de la rutina de mariposa para otra mejora de 100-300% de velocidad. ¡Tomando todo lo anterior, te harás una FFT bastante resbaladiza y rápida!

Para ir más lejos puede realizar la optimización sobre la marcha al proporcionar diferentes implementaciones de las etapas de FFT dirigidas a arquitecturas de procesador específicas. El tamaño de la caché, el recuento de registros, los conjuntos de instrucciones SSE/SSE2/3/4, etc., difieren según la máquina, por lo que la elección de un enfoque único para todos es a menudo superada por las rutinas específicas. En FFTW, por ejemplo, muchas FFT de tamaño más pequeño son implementaciones altamente optimizadas desenrolladas (sin bucles) dirigidas a una arquitectura específica. Al combinar estos constructos más pequeños (como las rutinas RadixN) puede elegir la rutina más rápida y mejor para la tarea en cuestión.

Answer 2

Mientras que no puedo dar una pista de rendimiento en este momento, me gustaría dar algunos consejos para su optimización que es demasiado largo para un comentario:

1. Si no lo ha hecho hecho así, escriba una serie de pruebas de corrección para su código en este momento. Pruebas simples como "hacer una FFT de esta matriz y ver si los resultados coinciden con los que he proporcionado" son suficientes, pero antes de optimizar el código, necesita una prueba de la unidad firme y automatizada que confirme que su código optimizado es correcto.
2. A continuación, perfile su código para ver dónde está el cuello de botella real. Aunque sospecho que el bucle más interno es for (i=j;i<n;i+=l2) {, ver es mejor que creer.

Answer 3

Esto parece una implementación básica de FFT radix-2 directamente de un libro de texto antiguo. Hay muchas docenas de documentos de hace décadas sobre la optimización de las FFT de varias maneras, dependiendo de muchos factores. Por ejemplo, ¿sus datos son más pequeños que el caché de la CPU? Agregado: Por ejemplo, si el vector de datos más una tabla de coeficientes encajarán en la CPU dcache y/o si multiplicaciones son mucho más lentas que los accesos a memoria en su CPU, precomputar una tabla de factores puede reducir el ciclo total contar para el uso repetido de la FFT. Pero si no, la precomputación podría ser más lenta. Punto de referencia. YMMV.

Answer 4

hay varias cosas que puedo recomendar tratar:

1. No intercambie los elementos de entrada, en lugar calcular el índice de bits invertida. Esto le ahorrará una cantidad de lecturas y escrituras de memoria.
2. Precalcula los coeficientes si estás haciendo muchas FFT del mismo tamaño. Esto ahorrará algunos cálculos.
3. Utilice radix-4 FFT en lugar de radix-2. Esto dará como resultado menos iteraciones en los bucles internos.

La respuesta definitiva se puede encontrar, por supuesto, perfilando el código.