algorithm
Transformada rápida de Fourier

Introducción

Radix 2 FFT

El método más simple y quizás el más conocido para calcular la FFT es el algoritmo de diezmado en el tiempo Radix-2. El Radix-2 FFT funciona al descomponer una señal de dominio de tiempo de punto N en N señales de dominio de tiempo cada una compuesta de un solo punto .

La descomposición de la señal, o 'decimación en el tiempo' se logra invirtiendo los índices para la matriz de datos del dominio del tiempo. Así, para una señal de dieciséis puntos, la muestra 1 (0001 binario) se intercambia con la muestra 8 (1000), la muestra 2 (0010) se intercambia con 4 (0100) y así sucesivamente. El intercambio de muestras con la técnica de inversión de bits se puede lograr simplemente en el software, pero limita el uso de Radix 2 FFT a señales de longitud N = 2 ^ M.

El valor de una señal de 1 punto en el dominio del tiempo es igual a su valor en el dominio de la frecuencia, por lo que esta matriz de puntos del dominio del tiempo único descompuestos no requiere transformación para convertirse en una matriz de puntos del dominio de la frecuencia. Los N puntos individuales; sin embargo, deben ser reconstruidos en un espectro de frecuencia de N puntos. La reconstrucción óptima del espectro de frecuencia completo se realiza mediante cálculos de mariposa. Cada etapa de reconstrucción en el Radix-2 FFT realiza una serie de mariposas de dos puntos, utilizando un conjunto similar de funciones de ponderación exponencial, Wn ^ R.

La FFT elimina los cálculos redundantes en la Transformada Discreta de Fourier explotando la periodicidad de Wn ^ R. La reconstrucción espectral se completa en las etapas log2 (N) de los cálculos de mariposa que dan X [K]; Los datos del dominio de frecuencia real e imaginario en forma rectangular. Para convertir a magnitud y fase (coordenadas polares) se requiere encontrar el valor absoluto, √ (Re2 + Im2), y el argumento, tan-1 (Im / Re).

El diagrama de flujo completo de la mariposa para un Radix 2 FFT de ocho puntos se muestra a continuación. Tenga en cuenta que las señales de entrada se han reordenado previamente de acuerdo con el procedimiento de reducción de tiempo descrito anteriormente.

La FFT normalmente funciona con entradas complejas y produce una salida compleja. Para señales reales, la parte imaginaria se puede establecer en cero y la parte real se establece en la señal de entrada, x [n], sin embargo, son posibles muchas optimizaciones que involucran la transformación de datos solo reales. Los valores de Wn ^ R utilizados a lo largo de la reconstrucción pueden determinarse utilizando la ecuación de ponderación exponencial.

El valor de R (la potencia de ponderación exponencial) se determina la etapa actual en la reconstrucción espectral y el cálculo actual dentro de una mariposa particular.

Ejemplo de código (C / C ++)

El ejemplo de código AC / C ++ para calcular la FFT de Radix 2 se puede encontrar a continuación. Esta es una implementación simple que funciona para cualquier tamaño N donde N es una potencia de 2. Es aproximadamente 3 veces más lenta que la implementación más rápida de FFTw, pero aún así es una muy buena base para futuras optimizaciones o para aprender cómo funciona este algoritmo.

#include <math.h>

#define PI       3.1415926535897932384626433832795    // PI for sine/cos calculations
#define TWOPI    6.283185307179586476925286766559     // 2*PI for sine/cos calculations
#define Deg2Rad  0.017453292519943295769236907684886  // Degrees to Radians factor
#define Rad2Deg  57.295779513082320876798154814105    // Radians to Degrees factor
#define log10_2  0.30102999566398119521373889472449   // Log10 of 2 
#define log10_2_INV 3.3219280948873623478703194294948 // 1/Log10(2)

// complex variable structure (double precision)
struct complex
{
public:
    double  Re, Im;        // Not so complicated after all
};

// Returns true if N is a power of 2
bool isPwrTwo(int N, int *M)
{
    *M = (int)ceil(log10((double)N) * log10_2_INV);// M is number of stages to perform. 2^M = N
    int NN = (int)pow(2.0, *M);
    if ((NN != N) || (NN == 0)) // Check N is a power of 2. 
        return false;

    return true;
}

void rad2FFT(int N, complex *x, complex *DFT)
{
    int M = 0;

    // Check if power of two. If not, exit        
    if (!isPwrTwo(N, &M))
        throw "Rad2FFT(): N must be a power of 2 for Radix FFT";

    // Integer Variables

    int BSep;                  // BSep is memory spacing between butterflies
    int BWidth;                // BWidth is memory spacing of opposite ends of the butterfly
    int P;                     // P is number of similar Wn's to be used in that stage
    int j;                     // j is used in a loop to perform all calculations in each stage
    int stage = 1;             // stage is the stage number of the FFT. There are M stages in total (1 to M).
    int HiIndex;               // HiIndex is the index of the DFT array for the top value of each butterfly calc 
    unsigned int iaddr;        // bitmask for bit reversal 
    int ii;                    // Integer bitfield for bit reversal (Decimation in Time)
    int MM1 = M - 1;

    unsigned int i;
    int l;
    unsigned int nMax = (unsigned int)N;

    // Double Precision Variables
    double TwoPi_N = TWOPI / (double)N;    // constant to save computational time.  = 2*PI / N
    double TwoPi_NP;

    // complex Variables (See 'struct complex')
    complex WN;               // Wn is the exponential weighting function in the form a + jb
    complex TEMP;             // TEMP is used to save computation in the butterfly calc
    complex *pDFT = DFT;      // Pointer to first elements in DFT array
    complex *pLo;             // Pointer for lo / hi value of butterfly calcs
    complex *pHi;
    complex *pX;              // Pointer to x[n]


    // Decimation In Time - x[n] sample sorting
    for (i = 0; i < nMax; i++, DFT++)
    {
        pX = x + i;             // Calculate current x[n] from base address *x and index i.
        ii = 0;                 // Reset new address for DFT[n]
        iaddr = i;              // Copy i for manipulations
        for (l = 0; l < M; l++) // Bit reverse i and store in ii...
        {
            if (iaddr & 0x01)     // Detemine least significant bit
                ii += (1 << (MM1 - l));    // Increment ii by 2^(M-1-l) if lsb was 1
            iaddr >>= 1;                // right shift iaddr to test next bit. Use logical operations for speed increase
            if (!iaddr)
                break;
        }
        DFT = pDFT + ii;        // Calculate current DFT[n] from base address *pDFT and bit reversed index ii    
        DFT->Re = pX->Re;       // Update the complex array with address sorted time domain signal x[n]
        DFT->Im = pX->Im;       // NB: Imaginary is always zero
    }

    // FFT Computation by butterfly calculation
    for (stage = 1; stage <= M; stage++) // Loop for M stages, where 2^M = N
    {
        BSep = (int)(pow(2, stage)); // Separation between butterflies = 2^stage
        P = N / BSep;             // Similar Wn's in this stage = N/Bsep
        BWidth = BSep / 2;     // Butterfly width (spacing between opposite points) = Separation / 2.

        TwoPi_NP = TwoPi_N*P;

        for (j = 0; j < BWidth; j++) // Loop for j calculations per butterfly
        {
            if (j != 0)              // Save on calculation if R = 0, as WN^0 = (1 + j0)
            {
                //WN.Re = cos(TwoPi_NP*j)
                WN.Re = cos(TwoPi_N*P*j);     // Calculate Wn (Real and Imaginary)
                WN.Im = -sin(TwoPi_N*P*j);
            }

            for (HiIndex = j; HiIndex < N; HiIndex += BSep) // Loop for HiIndex Step BSep butterflies per stage
            {
                pHi = pDFT + HiIndex;                  // Point to higher value
                pLo = pHi + BWidth;                    // Point to lower value (Note VC++ adjusts for spacing between elements)

                if (j != 0)                            // If exponential power is not zero...
                {
                    //CMult(pLo, &WN, &TEMP);          // Perform complex multiplication of Lovalue with Wn
                    TEMP.Re = (pLo->Re * WN.Re) - (pLo->Im * WN.Im);
                    TEMP.Im = (pLo->Re * WN.Im) + (pLo->Im * WN.Re);

                    //CSub (pHi, &TEMP, pLo);
                    pLo->Re = pHi->Re - TEMP.Re;       // Find new Lovalue (complex subtraction)
                    pLo->Im = pHi->Im - TEMP.Im;

                    //CAdd (pHi, &TEMP, pHi);          // Find new Hivalue (complex addition)
                    pHi->Re = (pHi->Re + TEMP.Re);
                    pHi->Im = (pHi->Im + TEMP.Im);
                }
                else
                {
                    TEMP.Re = pLo->Re;
                    TEMP.Im = pLo->Im;

                    //CSub (pHi, &TEMP, pLo);
                    pLo->Re = pHi->Re - TEMP.Re;       // Find new Lovalue (complex subtraction)
                    pLo->Im = pHi->Im - TEMP.Im;

                    //CAdd (pHi, &TEMP, pHi);          // Find new Hivalue (complex addition)
                    pHi->Re = (pHi->Re + TEMP.Re);
                    pHi->Im = (pHi->Im + TEMP.Im);
                }
            }
        }
    }


    pLo = 0;    // Null all pointers
    pHi = 0;
    pDFT = 0;
    DFT = 0;
    pX = 0;
}

Radix 2 Inverse FFT

Debido a la fuerte dualidad de la Transformada de Fourier, el ajuste de la salida de una transformación directa puede producir la FFT inversa. Los datos en el dominio de frecuencia se pueden convertir al dominio de tiempo mediante el siguiente método:

1. Encuentre el complejo conjugado de los datos del dominio de la frecuencia invirtiendo el componente imaginario para todas las instancias de K.
2. Realice la FFT directa en los datos del dominio de frecuencia conjugada.
3. Divida cada salida del resultado de esta FFT por N para obtener el valor real del dominio de tiempo.
4. Encuentre el complejo conjugado de la salida invirtiendo el componente imaginario de los datos del dominio de tiempo para todas las instancias de n.

Nota : los datos de dominio de frecuencia y tiempo son variables complejas. Normalmente, el componente imaginario de la señal de dominio de tiempo después de una FFT inversa es cero o se ignora como error de redondeo. El aumento de la precisión de las variables de flotación de 32 bits a doble de 64 bits, o doble de 128 bits de longitud reduce significativamente los errores de redondeo producidos por varias operaciones FFT consecutivas.

Ejemplo de código (C / C ++)

#include <math.h>

#define PI       3.1415926535897932384626433832795    // PI for sine/cos calculations
#define TWOPI    6.283185307179586476925286766559     // 2*PI for sine/cos calculations
#define Deg2Rad  0.017453292519943295769236907684886  // Degrees to Radians factor
#define Rad2Deg  57.295779513082320876798154814105    // Radians to Degrees factor
#define log10_2  0.30102999566398119521373889472449   // Log10 of 2 
#define log10_2_INV 3.3219280948873623478703194294948 // 1/Log10(2)

// complex variable structure (double precision)
struct complex
{
public:
    double  Re, Im;        // Not so complicated after all
}; 

void rad2InverseFFT(int N, complex *x, complex *DFT)
{
    // M is number of stages to perform. 2^M = N
    double Mx = (log10((double)N) / log10((double)2));
    int a = (int)(ceil(pow(2.0, Mx)));
    int status = 0;
    if (a != N) // Check N is a power of 2
    {
        x = 0;
        DFT = 0;
        throw "rad2InverseFFT(): N must be a power of 2 for Radix 2 Inverse FFT";
    }

    complex *pDFT = DFT;        // Reset vector for DFT pointers
    complex *pX = x;            // Reset vector for x[n] pointer
    double NN = 1 / (double)N;  // Scaling factor for the inverse FFT        

    for (int i = 0; i < N; i++, DFT++)
        DFT->Im *= -1;          // Find the complex conjugate of the Frequency Spectrum

    DFT = pDFT;                 // Reset Freq Domain Pointer
    rad2FFT(N, DFT, x); // Calculate the forward FFT with variables switched (time & freq)

    int i;
    complex* x;
    for ( i = 0, x = pX; i < N; i++, x++){
        x->Re *= NN;    // Divide time domain by N for correct amplitude scaling
        x->Im *= -1;    // Change the sign of ImX
    }    
}