algorithm
Snelle Fourier-transformatie

Invoering

Radix 2 FFT

De eenvoudigste en misschien wel bekendste methode voor het berekenen van de FFT is het Radix-2 Decimation in Time-algoritme. De Radix-2 FFT werkt door een N-punttijddomeinsignaal op te splitsen in N-tijddomeinsignalen die elk bestaan uit een enkel punt .

Signaalontleding, of 'decimering in de tijd' wordt bereikt door bit de indices voor de reeks tijddomeingegevens om te keren. Voor een zestienpuntsignaal wordt monster 1 (binair 0001) dus verwisseld met monster 8 (1000), monster 2 (0010) wordt verwisseld met 4 (0100) enzovoort. Monsteruitwisseling met behulp van de bit reverse techniek kan eenvoudig worden bereikt in software, maar beperkt het gebruik van de Radix 2 FFT tot signalen met een lengte N = 2 ^ M.

De waarde van een 1-punts signaal in het tijddomein is gelijk aan zijn waarde in het frequentiedomein, dus vereist deze reeks ontlede afzonderlijke tijddomeinpunten geen transformatie om een reeks frequentiedomeinpunten te worden. De N enkele punten; moeten echter worden gereconstrueerd in één N-punt frequentiespectra. Optimale reconstructie van het volledige frequentiespectrum wordt uitgevoerd met behulp van vlinderberekeningen. Elke reconstructiefase in de Radix-2 FFT voert een aantal tweepuntsvlinders uit, met behulp van een vergelijkbare set exponentiële wegingsfuncties, Wn ^ R.

De FFT verwijdert overbodige berekeningen in de Discrete Fourier Transform door gebruik te maken van de periodiciteit van Wn ^ R. Spectrale reconstructie wordt voltooid in log2 (N) stadia van vlinderberekeningen die X [K] geven; de echte en denkbeeldige frequentiedomeingegevens in rechthoekige vorm. Om te converteren naar grootte en fase (poolcoördinaten) moet de absolute waarde, √ (Re2 + Im2), en argument, tan-1 (Im / Re) worden gevonden.

Het volledige vlinderstroomdiagram voor een achtpuntige Radix 2 FFT wordt hieronder weergegeven. Merk op dat de ingangssignalen eerder opnieuw zijn gerangschikt volgens de eerder beschreven decimatie in tijdprocedure.

De FFT werkt meestal op complexe ingangen en produceert een complexe uitgang. Voor echte signalen kan het denkbeeldige deel worden ingesteld op nul en het reële deel worden ingesteld op het ingangssignaal, x [n], echter veel optimalisaties zijn mogelijk met betrekking tot de transformatie van alleen-echte gegevens. De waarden van Wn ^ R die tijdens de reconstructie worden gebruikt, kunnen worden bepaald met behulp van de exponentiële wegingsvergelijking.

De waarde van R (het exponentiële weegvermogen) wordt bepaald de huidige fase in de spectrale reconstructie en de huidige berekening binnen een bepaalde vlinder.

Codevoorbeeld (C / C ++)

AC / C ++ codevoorbeeld voor het berekenen van de Radix 2 FFT vindt u hieronder. Dit is een eenvoudige implementatie die werkt voor elke grootte N waarbij N een macht van 2 is. Het is ongeveer 3x langzamer dan de snelste FFTw-implementatie, maar nog steeds een zeer goede basis voor toekomstige optimalisatie of om te leren hoe dit algoritme werkt.

#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

Vanwege de sterke dualiteit van de Fourier-transformatie kan het aanpassen van de uitvoer van een voorwaartse transformatie de omgekeerde FFT produceren. Gegevens in het frequentiedomein kunnen op de volgende manier worden geconverteerd naar het tijdsdomein:

1. Vind de complexe conjugaat van de gegevens van het frequentiedomein door de denkbeeldige component voor alle instanties van K om te keren.
2. Voer de voorwaartse FFT uit op de geconjugeerde frequentiedomeingegevens.
3. Deel elke uitgang van het resultaat van deze FFT door N om de werkelijke tijdsdomeinwaarde te geven.
4. Vind de complexe conjugaat van de uitvoer door de denkbeeldige component van de tijddomeingegevens voor alle instanties van n om te keren.

Opmerking : zowel frequentie- als tijddomeingegevens zijn complexe variabelen. Typisch is de denkbeeldige component van het tijdsdomeinsignaal na een inverse FFT nul of wordt deze genegeerd als afrondingsfout. Het vergroten van de precisie van variabelen van 32-bit float naar 64-bit double, of 128-bit long double vermindert aanzienlijk de afrondingsfouten die worden geproduceerd door verschillende opeenvolgende FFT-bewerkingen.

Codevoorbeeld (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
    }    
}