La FFT (Transformada Rápida de Fourier) es un algoritmo que calcula la DFT (Transformada Discreta de Fourier) con un número mucho menor de operaciones aritméticas que las requeridas por la aplicación directa de la fórmula de la DFT. La FFT es una operación esencial en DSP (Procesamiento Digital de Señales), ya que a menudo resulta ser la forma más eficiente de implementar filtros digitales, realizar análisis espectral de señales periódicas y aleatorias, correlacionar secuencias, etc. A continuación se entrega el código fuente de una función escrita en C que calcula la FFT. La función fue traducida desde el lenguaje BASIC y validada en un microcontrolador ARM7. int N = 16; /* Número de muestras del bloque */ float REX[16]; /* Parte real de datos y resultados */ float IMX[16]; /* Parte imaginaria de datos y resultados */ float Mag[16]; /* Magnitud de resultados */ float PI = 3.14159265; int I, IP, JM1, K, L, LE, LE2, NM1, ND2, M, J; float TR, TI, UR, UI, SR, SI; int main (void) { NM1 = N-1; ND2 = N/2; M = ceil(log(N)/log(2)); J = ND2; /* Generación de la secuencia de prueba */ for(L=0; L<=(N-1); L++) { REX[L] = 2+2*cos(2*PI*L/N)+4*sin(2*PI*3*L/N); } /* Fin de la generación de la secuencia de prueba */ for(I=1; I<=N-2; I++) { if(I>=J) { goto L1; } TR = REX[J]; TI = IMX[J]; REX[J] = REX[I]; IMX[J] = IMX[I]; REX[I] = TR; IMX[I] = TI; L1: K = ND2; L2: if(K>J) { goto L3; } J = J-K; K = K/2; goto L2; L3: J = J+K; } for(L=1; L<=M; L++) { LE = ceil(pow(2,L)); LE2 = LE/2; UR = 1; UI = 0; SR = cos(PI/LE2); SI = -sin(PI/LE2); for(J=1; J<=LE2; J++) { JM1 = J-1; for(I=JM1; I<=NM1; I=I+LE) { IP = I+LE2; TR = REX[IP]*UR-IMX[IP]*UI; TI = REX[IP]*UI+IMX[IP]*UR; REX[IP] = REX[I]-TR; IMX[IP] = IMX[I]-TI; REX[I] = REX[I]+TR; IMX[I] = IMX[I]+TI; } TR = UR; UR = TR*SR-UI*SI; UI = TR*SI+UI*SR; } } REX[0]=REX[0]/2; IMX[0]=IMX[0]/2; for(L=0; L<=N-1; L++) { REX[L] = REX[L]/ND2; IMX[L] = IMX[L]/ND2; } for(L=0; L<=((N/2)-1); L++) { Mag[L] = sqrt(REX[L]*REX[L] + IMX[L]*IMX[L]); } } Julio de 2011 |
