algorithm
Szybka transformata Fouriera

Wprowadzenie

Radix 2 FFT

Najprostszą i prawdopodobnie najbardziej znaną metodą obliczania FFT jest algorytm Radix-2 Decimation in Time. Radix-2 FFT działa poprzez rozkład N sygnału w dziedzinie czasu w N na sygnały w dziedzinie czasu, każdy złożony z jednego punktu .

Dekompozycja sygnału lub „dziesiątkowanie w czasie” osiągane jest przez odwrócenie bitów wskaźników dla tablicy danych w dziedzinie czasu. Zatem dla szesnastopunktowego sygnału próbka 1 (Binarna 0001) jest zamieniana próbką 8 (1000), próbka 2 (0010) zamieniana jest na 4 (0100) i tak dalej. Zamiana próbek przy użyciu techniki odwrotności bitów może być osiągnięta po prostu w oprogramowaniu, ale ogranicza użycie Radix 2 FFT do sygnałów o długości N = 2 ^ M.

Wartość 1-punktowego sygnału w dziedzinie czasu jest równa jego wartości w dziedzinie częstotliwości, a zatem ta tablica zdekomponowanych pojedynczych punktów w dziedzinie czasu nie wymaga transformacji, aby stać się tablicą punktów w dziedzinie częstotliwości. N pojedynczych punktów; należy jednak zrekonstruować je w jednym widmie częstotliwości w punktach N. Optymalną rekonstrukcję pełnego spektrum częstotliwości przeprowadza się za pomocą obliczeń motylkowych. Każdy etap rekonstrukcji w Radix-2 FFT wykonuje szereg dwupunktowych motyli, wykorzystując podobny zestaw funkcji ważenia wykładniczego, Wn ^ R.

FFT usuwa zbędne obliczenia w dyskretnej transformacie Fouriera, wykorzystując okresowość Wn ^ R. Rekonstrukcja widmowa jest wykonywana w log2 (N) etapach obliczeń motyla dających X [K]; rzeczywiste i urojone dane w dziedzinie częstotliwości w formie prostokątnej. Aby przekonwertować na wielkość i fazę (współrzędne biegunowe), należy znaleźć wartość bezwzględną √ (Re2 + Im2) i argument tan-1 (Im / Re).

Pełny schemat przepływu motyla dla ośmiopunktowego Radix 2 FFT pokazano poniżej. Uwaga: sygnały wejściowe zostały uprzednio ponownie uporządkowane zgodnie z procedurą dziesiętną w czasie opisaną wcześniej.

FFT zazwyczaj działa na skomplikowanych danych wejściowych i wytwarza złożone dane wyjściowe. W przypadku sygnałów rzeczywistych część urojoną można ustawić na zero, a część rzeczywistą ustawić na sygnał wejściowy x [n], jednak możliwe są liczne optymalizacje obejmujące transformację danych tylko rzeczywistych. Wartości Wn ^ R stosowane podczas rekonstrukcji można określić za pomocą wykładniczego równania wagowego.

Wartość R (wykładniczej mocy ważącej) określa się bieżący etap w rekonstrukcji widmowej i bieżące obliczenia w obrębie konkretnego motyla.

Przykład kodu (C / C ++)

Przykład kodu AC / C ++ do obliczania Radix 2 FFT można znaleźć poniżej. Jest to prosta implementacja, która działa dla dowolnego rozmiaru N, gdzie N jest potęgą 2. Jest około 3 razy wolniejsza niż najszybsza implementacja FFTw, ale nadal stanowi bardzo dobrą podstawę do przyszłej optymalizacji lub do nauki o tym, jak działa ten algorytm.

#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

Ze względu na silną dualność transformacji Fouriera, dostosowanie mocy wyjściowej transformacji do przodu może powodować odwrotne FFT. Dane w dziedzinie częstotliwości można przekonwertować na dziedzinę czasu za pomocą następującej metody:

1. Znajdź złożony koniugat danych w dziedzinie częstotliwości, odwracając urojoną składową dla wszystkich instancji K.
2. Wykonaj FFT przesyłania dalej na sprzężonych danych w dziedzinie częstotliwości.
3. Podziel każdy wynik wyniku tej analizy FFT przez N, aby dać prawdziwą wartość w dziedzinie czasu.
4. Znajdź złożony koniugat wyniku, odwracając urojoną składową danych w dziedzinie czasu dla wszystkich instancji n.

Uwaga : dane w dziedzinie częstotliwości i czasu są zmiennymi złożonymi. Zazwyczaj urojona składowa sygnału w dziedzinie czasu po odwrotnym FFT wynosi albo zero, albo jest ignorowana jako błąd zaokrąglenia. Zwiększenie precyzji zmiennych z 32-bitowego zmiennoprzecinkowego na 64-bitowy podwójny lub 128-bitowy podwójny znacząco zmniejsza błędy zaokrąglania generowane przez kilka kolejnych operacji FFT.

Przykład kodu (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
    }    
}