algorithm
Trasformata di Fourier veloce

introduzione

Radix 2 FFT

Il metodo più semplice e forse più conosciuto per il calcolo della FFT è l'algoritmo di decimazione in tempo Radix-2. La FFT Radix-2 funziona decomponendo un segnale del dominio del tempo N in segnali del dominio del tempo N ciascuno composto da un singolo punto .

La scomposizione del segnale, o "decimazione nel tempo", viene ottenuta invertendo il bit degli indici per l'array di dati del dominio del tempo. Quindi, per un segnale a sedici punti, il campione 1 (Binario 0001) viene scambiato con il campione 8 (1000), il campione 2 (0010) viene scambiato con 4 (0100) e così via. Lo scambio di campioni usando la tecnica del bit reverse può essere ottenuto semplicemente con il software, ma limita l'uso della Radix 2 FFT a segnali di lunghezza N = 2 ^ M.

Il valore di un segnale a 1 punto nel dominio del tempo è uguale al suo valore nel dominio della frequenza, quindi questa matrice di punti del dominio del tempo decomposti non richiede alcuna trasformazione per diventare una matrice di punti del dominio di frequenza. I singoli punti N; tuttavia, devono essere ricostruiti in uno spettro di frequenza a N punti. La ricostruzione ottimale dello spettro di frequenza completo viene eseguita utilizzando i calcoli a farfalla. Ogni fase di ricostruzione della FFT Radix-2 esegue una serie di farfalle a due punti, utilizzando un set simile di funzioni di ponderazione esponenziale, Wn ^ R.

La FFT rimuove i calcoli ridondanti nella Trasformata di Fourier discreta sfruttando la periodicità di Wn ^ R. La ricostruzione spettrale è completata negli stadi log2 (N) dei calcoli a farfalla dando X [K]; i dati del dominio di frequenza reale e immaginario in forma rettangolare. Per convertire in grandezza e fase (coordinate polari) è necessario trovare il valore assoluto, √ (Re2 + Im2) e argomento, tan-1 (Im / Re).

Di seguito è mostrato il diagramma di flusso completo della farfalla per una FFT Radix 2 a otto punti. Nota che i segnali di ingresso sono stati precedentemente riordinati in base alla procedura di decimazione nel tempo descritta in precedenza.

La FFT in genere opera su input complessi e produce un output complesso. Per i segnali reali, la parte immaginaria può essere impostata su zero e la parte reale impostata sul segnale di ingresso, x [n], tuttavia sono possibili molte ottimizzazioni che implicano la trasformazione di dati solo reali. I valori di Wn ^ R utilizzati durante la ricostruzione possono essere determinati utilizzando l'equazione di ponderazione esponenziale.

Il valore di R (la potenza di ponderazione esponenziale) è determinato dallo stadio corrente nella ricostruzione spettrale e dal calcolo corrente all'interno di una particolare farfalla.

Esempio di codice (C / C ++)

Di seguito è riportato un esempio di codice AC / C ++ per il calcolo della FFT di Radix 2. Questa è un'implementazione semplice che funziona per qualsiasi dimensione N dove N è una potenza di 2. È approssimativamente 3 volte più lento dell'implementazione FFTw più veloce, ma è comunque una buona base per l'ottimizzazione futura o per imparare come funziona questo 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

A causa della forte dualità della Trasformata di Fourier, la regolazione dell'output di una trasformazione in avanti può produrre la FFT inversa. I dati nel dominio della frequenza possono essere convertiti nel dominio del tempo con il seguente metodo:

1. Trova il complesso coniugato dei dati del dominio della frequenza invertendo il componente immaginario per tutte le istanze di K.
2. Eseguire la FFT diretta sui dati del dominio di frequenza coniugato.
3. Dividere ogni uscita del risultato di questa FFT di N per fornire il valore del dominio del tempo reale.
4. Trova il complesso coniugato dell'output invertendo la componente immaginaria dei dati del dominio del tempo per tutte le istanze di n.

Nota : i dati relativi alla frequenza e al dominio del tempo sono variabili complesse. In genere, il componente immaginario del segnale del dominio del tempo che segue una FFT inversa è zero o ignorato come errore di arrotondamento. Aumentando la precisione delle variabili da 32-bit float a 64-bit double, o double lungo a 128-bit, si riducono significativamente gli errori di arrotondamento prodotti da diverse operazioni FFT consecutive.

Esempio di codice (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
    }    
}