algorithm
고속 푸리에 변환

소개

기수 2 FFT

FFT를 계산하는 가장 간단하고 아마도 가장 잘 알려진 방법은 Radix-2 Decimation in Time 알고리즘입니다. 기수 -2 FFT는 N 포인트 시간 도메인 신호를 각각 단일 포인트로 구성된 N 시간 도메인 신호들로 분해함으로써 작용한다 .

신호 분해 또는 '시간의 데시 메이션 (decimation in time)'은 시간 도메인 데이터의 배열에 대한 인덱스를 역전시키는 비트에 의해 달성된다. 따라서 16 점 신호의 경우 샘플 1 (Binary 0001)은 샘플 8 (1000)로 스왑되고 샘플 2 (0010)는 4 (0100)로 스왑됩니다. 비트 리버스 기법을 사용한 샘플 스와핑은 소프트웨어에서 간단하게 구현할 수 있지만 Radix 2 FFT를 길이 N = 2 ^ M의 신호로 사용하는 것을 제한합니다.

시간 영역에서 1- 점 신호의 값은 주파수 영역에서의 값과 동일하므로 분해 된 단일 시간 영역 점의 배열은 변환이 필요 없으므로 주파수 영역 점의 배열이됩니다. N 개의 단일 포인트; 그러나, 하나의 N- 포인트 주파수 스펙트럼으로 재구성 될 필요가있다. 완벽한 주파수 스펙트럼의 최적 재구성은 버터 플라이 계산을 사용하여 수행됩니다. Radix-2 FFT의 각 재구성 단계는 유사한 지수 가중 함수 세트 인 Wn ^ R을 사용하여 여러 개의 2 포인트 나비를 수행합니다.

FFT는 Wn ^ R의 주기성을 이용하여 Discrete Fourier Transform에서 중복 계산을 제거합니다. 분광 재구성은 X (K)를주는 나비 계산의 log2 (N) 단계에서 완료됩니다. 직각 형태의 실수 및 허수의 주파수 도메인 데이터. 진폭 및 위상 (극좌표)으로 변환하려면 절대 값 √ (Re2 + Im2) 및 인수 tan-1 (Im / Re)을 찾아야합니다.

8 포인트 기수 2 FFT에 대한 완전한 버터 플라이 흐름 다이어그램은 아래와 같습니다. 앞에서 설명한 데시 메이션 절차에 따라 입력 신호가 이전에 재정렬되었습니다.

FFT는 일반적으로 복잡한 입력에서 작동하고 복잡한 출력을 생성합니다. 실제 신호의 경우, 허수 부는 0으로 설정되고 실수 부는 입력 신호 x [n]로 설정 될 수 있지만, 실제 전용 데이터의 변환을 포함하여 많은 최적화가 가능합니다. 재구성 전체에 걸쳐 사용 된 Wn ^ R의 값은 지수 가중 방정식을 사용하여 결정될 수있다.

R의 값 (지수 가중력)은 스펙트럼 재구성의 현재 단계와 특정 나비 내 현재 계산에서 결정됩니다.

코드 예제 (C / C ++)

기수 2 FFT를 계산하기위한 AC / C ++ 코드 샘플은 아래에서 찾을 수 있습니다. 이는 N이 2의 거듭 제곱 인 모든 크기 N에서 작동하는 간단한 구현입니다. 가장 빠른 FFTw 구현보다 약 3 배 느리지 만 향후 최적화 또는이 알고리즘의 작동 방식을 배우는 데는 매우 좋은 기초가됩니다.

#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 역 FFT

푸리에 변환의 강한 이중성 때문에, 순방향 변환의 출력을 조정하면 역 FFT가 생성 될 수 있습니다. 주파수 도메인의 데이터는 다음 방법으로 시간 도메인으로 변환 할 수 있습니다.

1. K의 모든 경우에 대한 허수 성분을 반전시킴으로써 주파수 영역 데이터의 복소 공액을 찾는다.
2. 공액 주파수 영역 데이터에 대해 순방향 FFT를 수행합니다.
3. 이 FFT 결과의 각 출력을 N으로 나누어 실제 시간 도메인 값을 제공하십시오.
4. n의 모든 인스턴스에 대해 시간 도메인 데이터의 허수 성분을 반전시켜 출력의 복소 공액을 찾습니다.

참고 : 빈도 및 시간 도메인 데이터는 모두 복잡한 변수입니다. 전형적으로, 역 FFT 이후의 시간 영역 신호의 허수 성분은 제로이거나 반올림 에러로서 무시된다. 변수의 정밀도를 32 비트 float에서 64 비트 double 또는 128 비트 long double로 늘리면 여러 번의 연속 FFT 연산에서 발생하는 반올림 오류가 크게 줄어 듭니다.

코드 예제 (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
    }    
}