algorithm
Быстрое преобразование Фурье

Вступление

Radix 2 FFT

Простейшим и, возможно, самым известным методом вычисления БПФ является алгоритм Decix-2 Decimation in Time. Radix-2 FFT работает, разлагая N-точечный сигнал во временной области в N сигналов временной области, каждый из которых состоит из одной точки ,

Разложение сигнала или «прореживание во времени» достигается путем бина, изменяющего индексы для массива данных временной области. Таким образом, для шестнадцатиточечного сигнала образец 1 (двоичный 0001) заменяется образцом 8 (1000), образец 2 (0010) заменяется на 4 (0100) и так далее. Образцовая свопинг с использованием метода обратного бита может быть достигнута просто в программном обеспечении, но ограничивает использование Radix 2 FFT сигналами длиной N = 2 ^ M.

Значение одноточечного сигнала во временной области равно его значению в частотной области, поэтому этот массив разложенных единичных точек временной области не требует преобразования, чтобы стать массивом точек частотной области. N отдельных точек; однако необходимо перестроить в один N-точечный частотный спектр. Оптимальная реконструкция полного частотного спектра выполняется с использованием расчетов бабочки. Каждый этап реконструкции в Fix Radix-2 выполняет множество двухточечных бабочек, используя аналогичный набор экспоненциальных весовых функций Wn ^ R.

БПФ удаляет избыточные вычисления в Дискретном преобразовании Фурье, используя периодичность Wn ^ R. Спектральная реконструкция завершена в log2 (N) этапах расчетов бабочек, дающих X [K]; реальные и мнимые данные частотной области в прямоугольной форме. Для преобразования в амплитуду и фазу (полярные координаты) требуется найти абсолютное значение, √ (Re2 + Im2) и аргумент tan-1 (Im / Re).

Ниже показана полная диаграмма потока бабочек для восьмиступенчатого Fix Radix 2. Обратите внимание, что входные сигналы ранее были переупорядочены в соответствии с процедурой прореживания во времени, описанной ранее.

БПФ обычно работает на сложных входах и производит сложный вывод. Для реальных сигналов мнимая часть может быть установлена ​​равной нулю, а действительная часть установлена ​​на входной сигнал x [n], однако возможны многие оптимизации, связанные с преобразованием данных реального времени. Значения Wn ^ R, используемые во время реконструкции, могут быть определены с использованием экспоненциального весового уравнения.

Значение R (экспоненциальная весовая мощность) определяется текущим этапом в спектральной реконструкции и вычислением тока в конкретной бабочке.

Пример кода (C / C ++)

Пример кода AC / C ++ для вычисления Radix 2 FFT можно найти ниже. Это простая реализация, которая работает для любого размера N, где N - мощность 2. Она примерно на 3 раза медленнее, чем самая быстрая реализация FFTw, но все же является хорошей основой для будущей оптимизации или для изучения того, как работает этот алгоритм.

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

Радикс 2 Обратный БПФ

Из-за сильной двойственности преобразования Фурье, настройка выхода прямого преобразования может привести к обратному БПФ. Данные в частотной области могут быть преобразованы во временную область следующим образом:

1. Найти комплексное сопряжение данных частотной области путем инвертирования мнимой составляющей для всех экземпляров K.
2. Выполнение прямого БПФ на данных сопряженной частотной области.
3. Разделите каждый вывод результата этого FFT на N, чтобы указать значение истинной временной области.
4. Найти комплексное сопряжение вывода путем инвертирования мнимой составляющей данных во временной области для всех экземпляров n.

Примечание : данные частоты и временной области являются комплексными переменными. Обычно мнимая составляющая сигнала временной области после обратного БПФ либо равна нулю, либо игнорируется как ошибка округления. Повышение точности переменных с 32-битного поплавка до 64-битного двойного или 128-битного удвоения значительно снижает ошибки округления, возникающие в результате нескольких последовательных операций БПФ.

Пример кода (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
    }    
}