Geometry & Recognition

FFT는 입력 데이터의 개수$(N)$가 2의 지수승으로 주어질 때  $O(N\log N)$의 연산만으로 빠르게 DFT을 수행하는 알고리즘이다. 주어진 $N$개의 data $\{F_0, F_1, ..., F_{N-1}\}$의  DFT $\{a_0, a_1,..., a_{N-1}\}$는

$$ a_k = \sum_{n = 0}^{N-1} F_n \exp\Big( -i \frac{2\pi nk }{N} \Big)=\sum_{n = 0}^{N-1} F_n (W^n )^k , \quad W = \exp\Big( -i \frac{2\pi}{N}\Big)$$ 로 표현된다. $N$이 2의 지수승으로 주어지면 이 합은 $n$이 짝수인 항과 홀수인 항으로 아래처럼 분리할 수 있다:

\begin{align} a_k &= \sum _{n=0}^{N/2-1} F_{2n} (W^{2n})^k+W^{k} \sum_{n=0}^{N/2-1} F_{2n+1}(W^{2n})^k\\ &=\text{(N/2개 even 항에 대한 DFT)} + W^k \times \text{(N/2개 odd 항에 대한 DFT)}\\ &= e_k + W^k o_k \end{align}

더불어 주기성 때문에 $a_k$의 계산도 절반만 하면 된다: 

$$a_{k+ N/2} = e_k - W^k o_k, \quad k = 0,1,..., N/2-1.$$

이는 원래의 문제가 처음 입력 데이터의 절반만 가지고 처리하는 두 개의 동일한 문제로 줄일 수 있음을 보여준다. 전형적인 divide and conquer 기법을 적용할 수 있는 구조이므로 매우 빠르게 연산을 수행할 수 있다(Cooley-Tukey algorithm). 역변환의 경우에는 twiddle factor $W$를 $W^*$로 바꾸면 되므로(전체 normalization factor $\frac{1}{N}$가 덧붙여진다) 코드에서 쉽게 수정할 수 있다: $\tt W_{im}  \rightarrow  - { W}_{im}$. 아래는 FFT를 재귀적으로 구현한 코드이다. 비재귀적 구현은 kipl.tistory.com/22, https://kipl.tistory.com/622에서 찾을 수 있다.

N=8일 때 전개:

\begin{align}&a_k \\ &= \sum_{n=0}^7 F_n (W^{n})^k  \\ &= \sum_{n=0}^3 F_{2n} (W^{2n})^k +  W^k \sum_{n=0}^3 F_{2n+1}( W^{2n})^k   \\ &= \left[ \sum _{n=0}^1 F_{4n} (W^{4n})^k  + W^{2k} \sum _{n=0}^1 F_{4n +2} (W^{4n})^k\right]+W^k \left[\sum _{n=0}^1 F_{4n+1} (W^{4n})^k + W^{2k} \sum_{n=0}^1 F_{4n+3} (W^{4n})^k\right]\\ &= \left[ (F_0 + W^{4k} F_4) + W^{2k} (F_2 + W^{4k} F_6) \right]  + W^k \left[(F_1 + W^{4k} F_5 )+W^{2k} (F_3 + W^{4k} F_7) \right]\\ &= \left[ (F_{000} + W^{4k} F_{100}) + W^{2k} (F_{010} + W^{4k} F_{110}) \right] \\ &\qquad\qquad\qquad\qquad+ W^k \left[(F_{001} + W^{4k} F_{101} )+W^{2k} (F_{011} + W^{4k} F_{111}) \right]\end{align}

$F_n$ 인덱스의 이진비트가 twiddle factor를 어떻게 곱해야 하는지 알려주고, 비트를 역으로 읽으면 계수가 나타나는 순서를 알 수 있게 한다.

#define TWOPI 6.28318530717958647693
void split(int N, double* data) {
    if (N == 8) {
        double t1 = data[4]; data[4] = data[1];
        double t2 = data[2]; data[2] = t1; data[1] = t2;
        t1 = data[5]; data[5] = data[3];
        t2 = data[6]; data[6] = t1; data[3] = t2; 
        return;
    }
    double *tmp = new double [N / 2];
    for (int i = 0; i < N / 2; i++)   //copy odd elements to tmp buffer;
        tmp[i] = data[2 * i + 1];
    for (int i = 0; i < N / 2; i++)   // move even elements to lower half;
        data[i] = data[2 * i];            
    for (int i = 0; i < N / 2; i++)   // move odd elements to upper half;
        data[i + N / 2] = tmp[i];     
    delete [] tmp;
}
void fft4(double *re, double *im) {
    double tr1 = re[0]-re[2]+im[1]-im[3];
    double ti1 = im[0]-im[2]-re[1]+re[3];
    double tr2 = re[0]+re[2]-re[1]-re[3];
    double ti2 = im[0]+im[2]-im[1]-im[3];
    double ti3 = im[0]-im[2]+re[1]-re[3];
    double tr3 = re[0]-re[2]-im[1]+im[3];
    re[0] = re[0]+re[2]+re[1]+re[3];
    im[0] = im[0]+im[2]+im[1]+im[3];
    re[1] = tr1; im[1] = ti1;
    re[2] = tr2; im[2] = ti2;
    re[3] = tr3; im[3] = ti3;
}
int fft ( int N, double* re, double* im ) {
    if (N <= 0 || N & (N-1)) return 0; // N is not a power of 2;
    if ( N < 2 ) return 1;
    else if (N == 2) {
        double tr1 = re[0]-re[1], ti1 = im[0]-im[1];
        re[0] = re[0]+re[1]; im[0] = im[0]+im[1]; 
        re[1] = tr1; im[1] = ti1; 
        return 1;
    }  else if (N == 4) {
        fft4(re, im);
        return 1;
    }
        
    split ( N, re);
    split ( N, im);
    fft ( N / 2, &re[    0], &im[    0] );
    fft ( N / 2, &re[N / 2], &im[N / 2] );
    for ( int k = 0; k < N / 2; k++ ) {
        double Ere = re[k];
        double Eim = im[k];
        double Ore  = re[k + N / 2];
        double Oim  = im[k + N / 2];
        double Wre =  cos ( TWOPI * k / N );  //real part of twiddle factor:W
        double Wim =  -sin ( TWOPI * k / N ); //imag part of twiddlw factor:W
        double WOre = Wre * Ore - Wim * Oim;
        double WOim = Wre * Oim + Wim * Ore;
        re[k] = Ere + WOre;
        im[k] = Eim + WOim;
        re[k + N / 2] = Ere - WOre;
        im[k + N / 2] = Eim - WOim;
    }
    return 1;
};

6개의 주파수 $\{f_i = 2, 5, 9, 11, 21, 29\text{Hz}\}$가 섞인 신호 $s(t)$를 $[0,1]$ 초 구간에서 일정한 간격으로 sampling 해서 64개 data를 얻고(Nyquist frequency = 32Hz), 이에 대해 FFT를 수행한다.

$$s(t)=\sum_{i=0}^{5} (i+1) \cos(2\pi f_i t) $$

int fft_test_main() {
    const int samples = 64;
    double re[samples], im[samples];
    const int nfreqs = 6;
    double freq[nfreqs] = {2, 5, 9, 11, 21, 29};
    for (int i = 0; i < samples; i++) {
        double signal = 0;
        for (int k = 0; k < nfreqs; k++) 
            signal += (k + 1) * cos(TWOPI * freq[k] * i / samples);
        re[i] = signal;
        im[i] = 0.;
    }
    fft(samples, &re[0], &im[0]);
    for (int i = 0; i < samples; i++) TRACE("%d, %f, %f\n", i, re[i], im[i]);
    return 0;
}

FFT 결과:

FFT2D