algorithm
Transformée de Fourier Rapide

Introduction

Radix 2 FFT

La méthode la plus simple et peut-être la mieux connue pour calculer la FFT est l'algorithme Radix-2 Decimation in Time. Le Radix-2 FFT fonctionne en décomposant un signal de domaine temporel à N points en N signaux de domaine temporel composés chacun d’un seul point .

La décomposition du signal, ou «décimation dans le temps», est obtenue en inversant les indices pour le tableau des données du domaine temporel. Ainsi, pour un signal à 16 points, l'échantillon 1 (binaire 0001) est échangé avec l'échantillon 8 (1000), l'échantillon 2 (0010) est échangé avec 4 (0100), etc. L'échange d'échantillons en utilisant la technique d'inversion de bits peut être réalisé simplement dans un logiciel, mais limite l'utilisation de la FFT Radix 2 aux signaux de longueur N = 2 ^ M.

La valeur d'un signal à 1 point dans le domaine temporel est égale à sa valeur dans le domaine fréquentiel. Par conséquent, ce tableau de points temporels décomposés ne nécessite aucune transformation pour devenir un tableau de points de domaine fréquentiel. Les N points uniques; cependant, doivent être reconstruits dans un spectre de fréquence de N-point. La reconstruction optimale du spectre de fréquences complet est réalisée à l'aide de calculs papillons. Chaque étape de reconstruction de la FFT Radix-2 exécute un certain nombre de papillons à deux points, en utilisant un ensemble similaire de fonctions de pondération exponentielle, Wn ^ R.

La FFT supprime les calculs redondants dans la transformée de Fourier discrète en exploitant la périodicité de Wn ^ R. La reconstruction spectrale est complétée par les étapes log2 (N) des calculs papillons donnant X [K]; les données de domaine de fréquence réelles et imaginaires sous forme rectangulaire. Pour convertir en amplitude et en phase (coordonnées polaires), il faut trouver la valeur absolue √ (Re2 + Im2) et l'argument tan-1 (Im / Re).

L'organigramme complet d'un papillon Radix 2 FFT à huit points est présenté ci-dessous. Notez que les signaux d'entrée ont déjà été réorganisés en fonction de la procédure de décimation dans le temps décrite précédemment.

La FFT opère généralement sur des entrées complexes et produit une sortie complexe. Pour les signaux réels, la partie imaginaire peut être mise à zéro et la partie réelle définie sur le signal d'entrée, x [n], mais de nombreuses optimisations sont possibles, impliquant la transformation de données uniquement réelles. Les valeurs de Wn ^ R utilisées tout au long de la reconstruction peuvent être déterminées en utilisant l'équation de pondération exponentielle.

La valeur de R (la puissance de pondération exponentielle) est déterminée au stade actuel de la reconstruction spectrale et du calcul en cours dans un papillon particulier.

Exemple de code (C / C ++)

Vous trouverez ci-dessous un exemple de code AC / C ++ pour le calcul de la FFT Radix 2. Ceci est une implémentation simple qui fonctionne pour toute taille N où N est une puissance de 2. Il est environ 3 fois plus lent que l'implémentation FFTw la plus rapide, mais reste une très bonne base pour une optimisation future ou pour apprendre comment fonctionne cet algorithme.

#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;
}

FFT Radix 2 Inverse

En raison de la forte dualité de la transformée de Fourier, l'ajustement de la sortie d'une transformée directe peut produire l'inverse de la FFT. Les données du domaine fréquentiel peuvent être converties dans le domaine temporel par la méthode suivante:

1. Trouver le conjugué complexe des données du domaine de fréquence en inversant la composante imaginaire pour toutes les instances de K.
2. Effectuez la FFT avant sur les données du domaine de fréquence conjuguées.
3. Divisez chaque sortie du résultat de cette FFT par N pour obtenir la valeur réelle du domaine temporel.
4. Trouvez le conjugué complexe de la sortie en inversant la composante imaginaire des données du domaine temporel pour toutes les instances de n.

Remarque : les données de fréquence et de domaine temporel sont des variables complexes. Typiquement, la composante imaginaire du signal de domaine temporel suivant une FFT inverse est soit nulle, soit ignorée comme erreur d'arrondi. L'augmentation de la précision des variables de 32 bits à 64 bits ou de 128 bits double réduit de manière significative les erreurs d'arrondi générées par plusieurs opérations FFT consécutives.

Exemple de code (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
    }    
}