Implementé un FFT radix-4 de 4 puntos y descubrí que necesito hacer alguna manipulación de los términos de salida para que coincida con un dft.
Mi código es una implementación bastante directa de la formulación de matriz, así que no tengo claro cuál es el problema.
// |
// radix-4 butterfly matrix form | complex multiplication
// |
// +- -+ +- -+ | a+ib
// X[0] = | 1 1 1 1 | |x[0]| | * c+id
// X[1] = | 1 -i -1 i | |x[1]| | -------
// X[2] = | 1 -1 1 -1 | |x[2]| | ac + ibc
// X[3] = | 1 i -1 -i | |x[3]| | iad - bd
// +- -+ +- -+ | ------------------
// | (ac-bd) + i(bc+ad)
// |
¿Alguien puede ver dónde me equivoqué?
Gracias,
-David
typedef double fp; // base floating-point type
// naiive N-point DFT implementation as reference to check fft implementation against
//
void dft(int inv, struct cfp *x, struct cfp *y, int N) {
long int i, j;
struct cfp w;
fp ang;
for(i=0; i<N; i++) { // do N-point FFT/IFFT
y[i].r = y[i].i = 0;
if (inv) ang = 2*PI*(fp)i/(fp)N;
else ang = -2*PI*(fp)i/(fp)N;
for (j=0; j<N; j++) {
w.r = cos(j*ang);
w.i = sin(j*ang);
y[i].r += (x[j].r * w.r - x[j].i * w.i);
y[i].i += (x[j].r * w.i + x[j].i * w.r);
}
}
// scale output in the case of an IFFT
if (inv) {
for (i=0; i<N; i++) {
y[i].r = y[i].r/(fp)N;
y[i].i = y[i].i/(fp)N;
}
}
} // dft()
void r4fft4(int inv, int reorder, struct cfp *x, struct cfp *y) {
struct cfp x1[4], w[4];
fp ang, temp;
int i;
// |
// radix-4 butterfly matrix form | complex multiplication
// |
// +- -+ +- -+ | a+ib
// y[0] = | 1 1 1 1 | |x[0]| | * c+id
// y[1] = | 1 -i -1 i | |x[1]| | -------
// y[2] = | 1 -1 1 -1 | |x[2]| | ac + ibc
// y[3] = | 1 i -1 -i | |x[3]| | iad - bd
// +- -+ +- -+ | ------------------
// | (ac-bd) + i(bc+ad)
// |
if (inv) ang = 2*PI/(fp)4; // invert sign for IFFT
else ang = -2*PI/(fp)4;
//
w[1].r = cos(ang*1); w[1].i = sin(ang*1); // twiddle1 = exp(-2*pi/4 * 1);
w[2].r = cos(ang*2); w[2].i = sin(ang*2); // twiddle2 = exp(-2*pi/4 * 2);
w[3].r = cos(ang*3); w[3].i = sin(ang*3); // twiddle3 = exp(-2*pi/4 * 3);
// *1 *1 *1 *1
y[0].r = x[0].r + x[1].r + x[2].r + x[3].r;
y[0].i = x[0].i + x[1].i + x[2].i + x[3].i;
// *1 *-i *-1 *i
x1[1].r = x[0].r + x[1].i - x[2].r - x[3].i;
x1[1].i = x[0].i - x[1].r - x[2].i + x[3].r;
// *1 *-1 *1 *-1
x1[2].r = x[0].r - x[1].r + x[2].r - x[3].r;
x1[2].i = x[0].i - x[1].i + x[2].i - x[3].i;
// *1 *i *-1 *-i
x1[3].r = x[0].r - x[1].i - x[2].r + x[3].i;
x1[3].i = x[0].i + x[1].r - x[2].i - x[3].r;
//
y[1].r = x1[1].r*w[1].r - x1[1].i*w[1].i; // scale radix-4 output
y[1].i = x1[1].i*w[1].r + x1[1].r*w[1].i;
//
y[2].r = x1[2].r*w[2].r - x1[2].i*w[2].i; // scale radix-4 output
y[2].i = x1[2].i*w[2].r + x1[2].r*w[2].i;
//
y[3].r = x1[3].r*w[3].r - x1[3].i*w[3].i; // scale radix-4 output
y[3].i = x1[3].i*w[3].r + x1[3].r*w[3].i;
// reorder output stage ... mystery as to why I need this
if (reorder) {
temp = y[1].r;
y[1].r = -1*y[1].i;
y[1].i = temp;
//
y[2].r = -1*y[2].r;
//
temp = y[3].r;
y[3].r = y[3].i;
y[3].i = -1*temp;
}
// scale output for inverse FFT
if (inv) {
for (i=0; i<4; i++) { // scale output by 1/N for IFFT
y[i].r = y[i].r/(fp)4;
y[i].i = y[i].i/(fp)4;
}
}
} // r4fft4()
Respuestas:
Acabo de portar un radix-4 DIF fft del código de S. Burrus Fortran a Java. En realidad, carece de varias optimizaciones, en primer lugar del factor twiddle impulsado por la tabla (los factores sen y cos deben calcularse previamente). Esto debería acelerar el fft un poco más (tal vez 50%). Tengo que hackear un poco para eso, pero si alguien tiene la respuesta correcta, estaré muy feliz y agradecido. Publicaré el código optimizado lo antes posible, espero que tal vez con algunas pruebas de velocidad vs algoritmo radix-2.
Más, las multiplicaciones por 1 y sqrt (-1) no se eliminan. Eliminarlos acelerará un poco más. Pero, en general, la IMHO radix-4 parece no ser más de un 25% más rápida que una radix-2, por lo que no sé si la relación velocidad / complejidad realmente vale la pena. Tenga en cuenta que las bibliotecas muy optimizadas como FFTW están ampliamente disponibles y utilizadas, por lo que este esfuerzo podría ser solo un "desvío" personal.
Aquí está el código de Java. Portarlo a C, C ++ o C # debería ser muy fácil.
Aquí está el código original Radix-4 DIF FORTRAN de Sidney Burrus:
Radix-4, DIF, una mariposa FFT
fuente
Primero, su supuesta 'mariposa radix-4' es un DFT de 4 puntos, no un FFT. Tiene 16 operaciones complejas (es decir, N al cuadrado). Una FFT típica de 4 puntos solo tendría Nlog (base 2) N (= 8 para N = 4). En segundo lugar, tiene algunos supuestos factores de 'escala' w [] .r y w [] .i que no pertenecen. Tal vez los obtuvo de una mariposa radix-4 que se muestra en un gráfico más grande. Tal mariposa tendría algunos giros entre etapas, pero en realidad no son parte de la mariposa. Un FFT de 4 puntos solo tiene una mariposa interna de -j cuando está diseñado para un exponente negativo FFT.
En lugar de intentar arreglar su código, es igual de fácil escribir el mío, como se muestra a continuación (compilador DevC ++; salidas agregadas al final del código):
Primero, se imprimen los datos de entrada (4 reales, 4 imaginarios). Luego se toma un DFT de 4 puntos. Los resultados (año [] e yi [] más amplificador / fase) se imprimen. Como los datos originales r [] e i [] no se sobrescribieron al hacer el DFT, esas entradas se reutilizan como entradas al FFT de 4 puntos. Tenga en cuenta que este último tiene menos operaciones +/- que el DFT.
El código para el FFT no es particularmente elegante ni eficiente: hay muchas formas de hacer mariposas. El código anterior corresponde a las cuatro mariposas radix-2 que se muestran en el libro de Rabiner y Gold "Teoría y aplicación del procesamiento de señales digitales" (p. 580, Fig. 10.9), con twiddles modificados para reflejar un exponente negativo (los utilizados para el figura en el libro fueron positivos). Tenga en cuenta que solo hay un twiddle de -j en el código, y esto no requiere una multiplicación (es un cambio de intercambio / signo).
Después de la FFT, se imprimen los resultados. Son lo mismo que el DFT
Y finalmente, los resultados escalados de la FFT se usan como entradas para una FFT inversa. Esto se logra mediante el método de 'intercambio' o 'invertir la lista' (es decir: si FFT (r, i) es una FFT directa, entonces FFT (i, r) es una inversa, siempre que, por supuesto, la FFT sea capaz de manejar entradas / salidas complejas, en otras palabras, no hay rutinas 'solo reales', que generalmente suponen que las entradas imaginarias son cero). Este método fue descrito hace casi 25 años en:
P. Duhamel, B. Piron, JM Etcheto, "Sobre el cómputo de la DFT inversa", IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 36, febrero de 1988, págs. 285-286.
Luego se imprime el resultado del inverso. Es lo mismo que los datos de entrada originales.
fuente