algorithm
Snabb Fourier Transform

Introduktion

Radix 2 FFT

Den enklaste och kanske mest kända metoden för beräkning av FFT är Radix-2 Decimation in Time-algoritmen. Radix-2 FFT fungerar genom att sönderdela en N-punkts tiddomänssignal till N tiddomänssignaler som var och en består av en enda punkt .

Signalnedbrytning, eller "decimering i tid" uppnås genom att bitomvända index för uppsättningen av tidsdomändata. Således, för en sextonpunktssignal, byts prov 1 (binär 0001) med prov 8 (1000), prov 2 (0010) byts med 4 (0100) och så vidare. Provbyte med bit-reverse-teknik kan uppnås helt enkelt i mjukvara, men begränsar användningen av Radix 2 FFT till signaler med längden N = 2 ^ M.

Värdet på en 1-punktssignal i tidsdomänen är lika med dess värde i frekvensdomänen, varför denna uppsättning av sönderdelade enstaka tidsdomänspunkter inte kräver någon transformation för att bli en matris med frekvensdomänpunkter. N-poäng; måste emellertid rekonstrueras till ett frekvensspektra för N-punkter. Optimal rekonstruktion av det kompletta frekvensspektrumet utförs med hjälp av fjärilsberäkningar. Varje rekonstruktionssteg i Radix-2 FFT utför ett antal tvåpunktsfjärilar med användning av en liknande uppsättning exponentiella viktningsfunktioner, Wn ^ R.

FFT tar bort redundanta beräkningar i den diskreta Fourier-transformen genom att utnyttja periodiciteten för Wn ^ R. Spektral rekonstruktion är klar i log2 (N) stadier av fjärilsberäkningar som ger X [K]; de verkliga och imaginära frekvensdomändata i rektangulär form. För att konvertera till storlek och fas (polära koordinater) krävs att hitta det absoluta värdet, √ (Re2 + Im2) och argument, tan-1 (Im / Re).

Det kompletta flödesschemat för fjäril för en åttapunkts Radix 2 FFT visas nedan. Observera att insignalerna tidigare har sorterats om enligt tidigare decimering i tidsförfarandet.

FFT arbetar vanligtvis på komplexa ingångar och producerar en komplex utgång. För verkliga signaler kan den imaginära delen ställas in på noll och verklig del inställd på insignalen, x [n], men många optimeringar är möjliga som involverar omvandlingen av data som bara är verkliga. Värden för Wn ^ R som används under hela rekonstruktionen kan bestämmas med hjälp av den exponentiella viktningsekvationen.

Värdet på R (den exponentiella viktningseffekten) bestäms det aktuella steget i spektralrekonstruktionen och den aktuella beräkningen inom en viss fjäril.

Kodexempel (C / C ++)

AC / C ++ kodprov för beräkning av Radix 2 FFT kan hittas nedan. Detta är en enkel implementering som fungerar för alla storlekar N där N har en effekt på 2. Den är cirka 3x långsammare än den snabbaste FFTw-implementeringen, men fortfarande en mycket bra bas för framtida optimering eller för att lära sig hur denna algoritm fungerar.

#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

På grund av Fourier Transforms starka dualitet kan justering av utgången från en framåttransformering producera det omvända FFT. Data i frekvensdomänen kan konverteras till tidsdomänen med följande metod:

1. Hitta det komplexa konjugatet av frekvensdomändata genom att invertera den imaginära komponenten för alla fall av K.
2. Utför framåt FFT på den konjugerade frekvensdomändata.
3. Dela varje utgång från resultatet av denna FFT med N för att ge det sanna tidsdomänvärdet.
4. Hitta det komplexa konjugatet av utgången genom att invertera den imaginära komponenten av tidsdomändata för alla instanser av n.

Obs : både frekvens- och tidsdomändata är komplexa variabler. Typiskt är den imaginära komponenten av tidsdomänssignalen efter en omvänd FFT antingen noll eller ignoreras som avrundningsfel. Att öka precisionen för variabler från 32-bitars flottör till 64-bitars dubbel eller 128-bitars lång dubbel minskar avsevärt avrundningsfel som produceras av flera på varandra följande FFT-operationer.

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