algorithm
Schnelle Fourier-Transformation

Einführung

Radix 2 FFT

Das einfachste und vielleicht bekannteste Verfahren zur Berechnung der FFT ist der Radix-2-Algorithmus Decimation in Time. Die Radix-2-FFT zerlegt ein N-Punkt-Zeitbereichssignal in N-Zeitbereichssignale, die jeweils aus einem einzelnen Punkt bestehen .

Die Signalzerlegung oder zeitliche Dezimierung wird durch Bitumkehrung der Indizes für das Array der Zeitbereichsdaten erreicht. Bei einem 16-Punkt-Signal wird also Probe 1 (Binär 0001) mit Probe 8 (1000) und Probe 2 (0010) mit 4 (0100) usw. ausgetauscht. Das Austauschen von Mustern unter Verwendung der Bit-Reverse-Technik kann einfach in Software erreicht werden, beschränkt jedoch die Verwendung der Radix 2-FFT auf Signale der Länge N = 2 ^ M.

Der Wert eines 1-Punkt-Signals im Zeitbereich ist gleich seinem Wert im Frequenzbereich. Daher erfordert dieses Array zerlegter Einzelzeitpunktpunkte keine Transformation, um ein Array von Frequenzdomänenpunkten zu werden. Die N Einzelpunkte; Sie müssen jedoch in ein N-Punkt-Frequenzspektrum rekonstruiert werden. Die optimale Rekonstruktion des gesamten Frequenzspektrums wird mithilfe von Schmetterlingsberechnungen durchgeführt. Jede Rekonstruktionsstufe in der Radix-2-FFT führt eine Anzahl von Zweipunkt-Schmetterlingen aus, wobei ein ähnlicher Satz von exponentiellen Gewichtungsfunktionen Wn ^ R verwendet wird.

Die FFT entfernt redundante Berechnungen in der diskreten Fourier-Transformation, indem die Periodizität von Wn ^ R ausgenutzt wird. Die spektrale Rekonstruktion wird in log2 (N) -Stufen der Schmetterlingsberechnungen mit X [K] abgeschlossen. die realen und imaginären Frequenzbereichsdaten in rechteckiger Form. Um nach Betrag und Phase (Polarkoordinaten) zu konvertieren, müssen der absolute Wert √ (Re2 + Im2) und das Argument tan-1 (Im / Re) ermittelt werden.

Das vollständige Butterfly-Flussdiagramm für eine Radix 2-FFT mit acht Punkten ist unten dargestellt. Beachten Sie, dass die Eingangssignale zuvor entsprechend dem zuvor beschriebenen Verfahren für die Dezimierung in der Zeit umgeordnet wurden.

Die FFT arbeitet normalerweise mit komplexen Eingaben und erzeugt eine komplexe Ausgabe. Für reelle Signale kann der Imaginärteil auf Null gesetzt werden, und der Realteil kann auf das Eingangssignal x [n] gesetzt werden, es sind jedoch viele Optimierungen möglich, die die Transformation von reinen Daten beinhalten. Die während der Rekonstruktion verwendeten Werte von Wn ^ R können unter Verwendung der Exponentialgewichtungsgleichung bestimmt werden.

Der Wert von R (die exponentielle Gewichtungsstärke) wird als die aktuelle Stufe der Spektralrekonstruktion und die Stromberechnung innerhalb eines bestimmten Schmetterlings bestimmt.

Codebeispiel (C / C ++)

Ein AC / C ++ - Codebeispiel zur Berechnung der Radix 2-FFT finden Sie unten. Dies ist eine einfache Implementierung, die für jede Größe N geeignet ist, bei der N eine Potenz von 2 ist. Sie ist etwa dreimal langsamer als die schnellste FFTw-Implementierung, aber immer noch eine sehr gute Basis für zukünftige Optimierungen oder für das Lernen, wie dieser Algorithmus funktioniert.

#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

Aufgrund der starken Dualität der Fourier-Transformation kann das Anpassen der Ausgabe einer Vorwärtstransformation die inverse FFT erzeugen. Daten im Frequenzbereich können auf folgende Weise in den Zeitbereich konvertiert werden:

1. Finden Sie das komplexe Konjugat der Frequenzbereichsdaten, indem Sie die Imaginärkomponente für alle Instanzen von K invertieren.
2. Führen Sie die Vorwärts-FFT für die Daten der konjugierten Frequenzdomäne durch.
3. Dividieren Sie jede Ausgabe des Ergebnisses dieser FFT durch N, um den tatsächlichen Zeitbereichswert zu erhalten.
4. Finden Sie das komplexe Konjugat der Ausgabe, indem Sie die Imaginärkomponente der Zeitbereichsdaten für alle Instanzen von n invertieren.

Hinweis : Sowohl Frequenz- als auch Zeitbereichsdaten sind komplexe Variablen. Typischerweise ist die imaginäre Komponente des Zeitbereichssignals nach einer inversen FFT entweder Null oder wird als Rundungsfehler ignoriert. Durch die Erhöhung der Genauigkeit von Variablen von 32-Bit-Gleitkommazahlen auf 64-Bit-Doppeln oder 128-Bit-Doppeln werden Rundungsfehler, die durch mehrere aufeinander folgende FFT-Operationen verursacht werden, erheblich reduziert.

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