c++ - make - image fft matlab
FFTW vs Matlab FFT (4)
Lo publiqué en el centro de Matlab pero no obtuve ninguna respuesta, así que pensé que volvería a publicar aquí.
Recientemente escribí una rutina simple en Matlab que usa una FFT en for-loop; la FFT domina los cálculos. Escribí la misma rutina en mex solo con fines de experimentación y llama a la biblioteca FFTW 3.3. Resulta que la rutina matlab se ejecuta más rápido que la rutina mex para arreglos muy grandes (aproximadamente el doble de rápido). La rutina mex usa sabiduría y realiza los mismos cálculos FFT. También sé que matlab utiliza FFTW, pero ¿es posible que su versión esté un poco más optimizada? Incluso utilicé el indicador FFTW_EXHAUSTIVE y todavía es aproximadamente el doble de lenta para arreglos grandes que MATLAB. Además, me aseguré de que el matlab que utilicé tuviera una sola hebra con el distintivo "-singleCompThread" y el archivo mex que utilicé no estaba en modo de depuración. Solo curiosidad por saber si este fue el caso, o si hay algunas optimizaciones que está usando matlab bajo el capó y que yo desconozco. Gracias.
Aquí está la porción mex:
void class_cg_toeplitz::analysis() {
// This method computes CG iterations using FFTs
// Check for wisdom
if(fftw_import_wisdom_from_filename("cd.wis") == 0) {
mexPrintf("wisdom not loaded./n");
} else {
mexPrintf("wisdom loaded./n");
}
// Set FFTW Plan - use interleaved FFTW
fftw_plan plan_forward_d_buffer;
fftw_plan plan_forward_A_vec;
fftw_plan plan_backward_Ad_buffer;
fftw_complex *A_vec_fft;
fftw_complex *d_buffer_fft;
A_vec_fft = fftw_alloc_complex(n);
d_buffer_fft = fftw_alloc_complex(n);
// CREATE MASTER PLAN - Do this on an empty vector as creating a plane
// with FFTW_MEASURE will erase the contents;
// Use d_buffer
// This is somewhat dangerous because Ad_buffer is a vector; but it does not
// get resized so &Ad_buffer[0] should work
plan_forward_d_buffer = fftw_plan_dft_r2c_1d(d_buffer.size(),&d_buffer[0],d_buffer_fft,FFTW_EXHAUSTIVE);
plan_forward_A_vec = fftw_plan_dft_r2c_1d(A_vec.height,A_vec.value,A_vec_fft,FFTW_WISDOM_ONLY);
// A_vec_fft.*d_buffer_fft will overwrite d_buffer_fft
plan_backward_Ad_buffer = fftw_plan_dft_c2r_1d(Ad_buffer.size(),d_buffer_fft,&Ad_buffer[0],FFTW_EXHAUSTIVE);
// Get A_vec_fft
fftw_execute(plan_forward_A_vec);
// Find initial direction - this is the initial residual
for (int i=0;i<n;i++) {
d_buffer[i] = b.value[i];
r_buffer[i] = b.value[i];
}
// Start CG iterations
norm_ro = norm(r_buffer);
double fft_reduction = (double)Ad_buffer.size(); // Must divide by size of vector because inverse FFT does not do this
while (norm(r_buffer)/norm_ro > relativeresidual_cutoff) {
// Find Ad - use fft
fftw_execute(plan_forward_d_buffer);
// Get A_vec_fft.*fft(d) - A_vec_fft is only real, but d_buffer_fft
// has complex elements; Overwrite d_buffer_fft
for (int i=0;i<n;i++) {
d_buffer_fft[i][0] = d_buffer_fft[i][0]*A_vec_fft[i][0]/fft_reduction;
d_buffer_fft[i][1] = d_buffer_fft[i][1]*A_vec_fft[i][0]/fft_reduction;
}
fftw_execute(plan_backward_Ad_buffer);
// Calculate r''*r
rtr_buffer = 0;
for (int i=0;i<n;i++) {
rtr_buffer = rtr_buffer + r_buffer[i]*r_buffer[i];
}
// Calculate alpha
alpha = 0;
for (int i=0;i<n;i++) {
alpha = alpha + d_buffer[i]*Ad_buffer[i];
}
alpha = rtr_buffer/alpha;
// Calculate new x
for (int i=0;i<n;i++) {
x[i] = x[i] + alpha*d_buffer[i];
}
// Calculate new residual
for (int i=0;i<n;i++) {
r_buffer[i] = r_buffer[i] - alpha*Ad_buffer[i];
}
// Calculate beta
beta = 0;
for (int i=0;i<n;i++) {
beta = beta + r_buffer[i]*r_buffer[i];
}
beta = beta/rtr_buffer;
// Calculate new direction vector
for (int i=0;i<n;i++) {
d_buffer[i] = r_buffer[i] + beta*d_buffer[i];
}
*total_counter = *total_counter+1;
if(*total_counter >= iteration_cutoff) {
// Set total_counter to -1, this indicates failure
*total_counter = -1;
break;
}
}
// Store Wisdom
fftw_export_wisdom_to_filename("cd.wis");
// Free fft alloc''d memory and plans
fftw_destroy_plan(plan_forward_d_buffer);
fftw_destroy_plan(plan_forward_A_vec);
fftw_destroy_plan(plan_backward_Ad_buffer);
fftw_free(A_vec_fft);
fftw_free(d_buffer_fft);
};
Aquí está la porción de matlab:
% Take FFT of A_vec.
A_vec_fft = fft(A_vec); % Take fft once
% Find initial direction - this is the initial residual
x = zeros(n,1); % search direction
r = zeros(n,1); % residual
d = zeros(n+(n-2),1); % search direction; pad to allow FFT
for i = 1:n
d(i) = b(i);
r(i) = b(i);
end
% Enter CG iterations
total_counter = 0;
rtr_buffer = 0;
alpha = 0;
beta = 0;
Ad_buffer = zeros(n+(n-2),1); % This holds the product of A*d - calculate this once per iteration and using FFT; only 1:n is used
norm_ro = norm(r);
while(norm(r)/norm_ro > 10^-6)
% Find Ad - use fft
Ad_buffer = ifft(A_vec_fft.*fft(d));
% Calculate rtr_buffer
rtr_buffer = r''*r;
% Calculate alpha
alpha = rtr_buffer/(d(1:n)''*Ad_buffer(1:n));
% Calculate new x
x = x + alpha*d(1:n);
% Calculate new residual
r = r - alpha*Ad_buffer(1:n);
% Calculate beta
beta = r''*r/(rtr_buffer);
% Calculate new direction vector
d(1:n) = r + beta*d(1:n);
% Update counter
total_counter = total_counter+1;
end
En términos de tiempo, para N = 50000 yb = 1: n toma aproximadamente 10.5 segundos con mex y 4.4 segundos con matlab. Estoy usando R2011b. Gracias
Algunas observaciones en lugar de una respuesta definitiva, ya que no conozco ninguno de los detalles de la implementación de MATLAB FFT:
- Basado en el código que tienes, puedo ver dos explicaciones para la diferencia de velocidad:
- la diferencia de velocidad se explica por las diferencias en los niveles de optimización de la FFT
- el ciclo while en MATLAB se ejecuta un número significativamente menor de veces
Asumiré que ya investigó el segundo problema y que el número de iteraciones es comparable. (Si no es así, es muy probable que se deba a algunos problemas de precisión y que valga la pena investigar más).
Ahora, con respecto a la comparación de velocidad de FFT:
- Sí, la teoría es que FFTW es más rápido que otras implementaciones de FFT de alto nivel, pero solo es relevante siempre que compares manzanas con manzanas: aquí estás comparando implementaciones en un nivel más abajo, en el nivel de ensamblaje, donde no solo el selección del algoritmo pero su optimización real para un procesador específico y por desarrolladores de software con diferentes habilidades entra en juego
- He optimizado o revisado FFT optimizados en ensamblado en muchos procesadores a lo largo del año (estuve en la industria de evaluación comparativa) y los grandes algoritmos son solo una parte de la historia. Hay consideraciones que son muy específicas para la arquitectura que está codificando (contabilidad de latencias, programación de instrucciones, optimización del uso de registros, disposición de datos en la memoria, contabilidad de latencias tomadas / no tomadas, etc.) y que hacen diferencias tan importante como la selección del algoritmo.
- Con N = 500000, también estamos hablando de grandes búferes de memoria: otra puerta para más optimizaciones que pueden volverse bastante específicas para la plataforma en la que ejecuta su código: qué tan bien se las arregla para evitar fallas de caché no será dictado por el Algoritmo tanto como por el flujo de datos y las optimizaciones que un desarrollador de software puede haber usado para traer datos dentro y fuera de la memoria de manera eficiente.
- Aunque no conozco los detalles de la implementación de MATLAB FFT, estoy bastante seguro de que un ejército de ingenieros de DSP ha estado (y sigue) perfeccionando su optimización, ya que es clave para tantos diseños. Esto bien podría significar que MATLAB tenía la combinación correcta de desarrolladores para producir una FFT mucho más rápida.
Esta es una ganancia de rendimiento clásica gracias a la optimización de bajo nivel y específica de la arquitectura.
Matlab utiliza FFT del binario Intel MKL (Math Kernel Library) (mkl.dll). Estas son rutinas optimizadas (a nivel de ensamblaje) por Intel para procesadores Intel. Incluso en AMD parece dar buenos aumentos de rendimiento.
FFTW parece una biblioteca c normal que no está tan optimizada. De ahí la ganancia de rendimiento para usar el MKL.
He encontrado el siguiente comentario en el sitio web de MathWorks [1]:
Nota sobre potencias grandes de 2: para dimensiones de FFT que son potencias de 2, entre 2 ^ 14 y 2 ^ 22, el software MATLAB utiliza información precargada especial en su base de datos interna para optimizar el cálculo de FFT. No se realiza ninguna afinación cuando la dimensión del FTT es una potencia de 2, a menos que borre la base de datos usando el comando fftw (''sabiduría'', []).
Aunque se relaciona con los poderes de 2, puede insinuar que MATLAB emplea su propia "sabiduría especial" cuando usa FFTW para ciertos tamaños de matriz (grandes). Considere: 2 ^ 16 = 65536.
[1] R2013b Documentación disponible en http://www.mathworks.de/de/help/matlab/ref/fftw.html (consultado el 29 de octubre de 2013)
EDITAR: la respuesta de @wakjah a esta respuesta es precisa: FFTW admite almacenamiento de memoria real e imaginaria dividida a través de su interfaz Guru. Por lo tanto, mi afirmación sobre la piratería no es precisa, pero puede aplicarse si la interfaz Guru de FFTW no se utiliza, que es el caso por defecto, así que ¡cuidado!
En primer lugar, perdón por haber llegado un año tarde. No estoy convencido de que el aumento de velocidad que ve proviene de MKL u otras optimizaciones. Hay algo completamente diferente entre FFTW y Matlab, y así es como los datos complejos se almacenan en la memoria.
En Matlab, las partes real e imaginaria de un vector complejo X son matrices separadas Xre [i] y Xim [i] (lineal en memoria, eficiente cuando se opera en cualquiera de ellas por separado).
En FFTW, las partes real e imaginaria están entrelazadas como dobles [2] por defecto, es decir, X [i] [0] es la parte real, y X [i] [1] es la parte imaginaria.
Por lo tanto, para usar la biblioteca FFTW en archivos mex uno no puede usar la matriz Matlab directamente, sino que primero debe asignar nueva memoria, luego empacar la entrada de Matlab en formato FFTW y luego descomprimir la salida de FFTW en formato Matlab. es decir
X = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
Y = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
entonces
for (size_t i=0; i<N; ++i) {
X[i][0] = Xre[i];
X[i][1] = Xim[i];
}
entonces
for (size_t i=0; i<N; ++i) {
Yre[i] = Y[i][0];
Yim[i] = Y[i][1];
}
Por lo tanto, esto requiere 2x asignaciones de memoria + 4x lecturas + 4x escrituras - todas de tamaño N. Esto tiene un costo en cuanto a la velocidad en grandes problemas.
Tengo la corazonada de que Mathworks ha pirateado el código FFTW3 para permitirle leer vectores de entrada directamente en el formato Matlab, lo que evita todo lo anterior.
En este escenario, uno solo puede asignar X y usar X para Y ejecutar FFTW en el lugar (como fftw_plan_*(N, X, X, ...) lugar de fftw_plan_*(N, X, Y, ...) ), ya que se copiará al vector Yre y Yim Matlab, a menos que la aplicación requiera / se beneficie de mantener separadas X e Y.
EDITAR : Viendo el consumo de memoria en tiempo real al ejecutar fft2 () de Matlab y mi código basado en la biblioteca fftw3, muestra que Matlab solo asigna una matriz compleja adicional (la salida), mientras que mi código necesita dos de tales matrices ( el *fftw_complex buffer más la salida de Matlab). No es posible realizar una conversión in situ entre los formatos Matlab y fftw porque las matrices real e imaginaria de Matlab no son consecutivas en la memoria. Esto sugiere que Mathworks hackeó la biblioteca fftw3 para leer / escribir los datos usando el formato Matlab.
Otra optimización para llamadas múltiples es asignar persistentemente (usando mexMakeMemoryPersistent() ). No estoy seguro de si la implementación de Matlab también lo hace.
Aclamaciones.
ps Como nota al margen, el formato de almacenamiento de datos complejo de Matlab es más eficiente para operar en los vectores reales o imaginarios por separado. En el formato de FFTW tendrías que hacer ++ 2 lecturas de memoria.