FFT 알고리즘의 재귀적 구현(a recursive FFT implementation)

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 (-i\frac{2\pi nk}{N})=\sum _{n=0}^{N-1}{F}_n(W^n{)}^k,W=\exp (-i\frac{2\pi }{N})$$ 로 표현된다. $N$이 2의 지수승으로 주어지면 이 합은 $n$이 짝수인 항과 홀수인 항으로 아래처럼 분리할 수 있다:

$$\begin{array}{rl}{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{array}$$

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

$$a_{k+N/2}=e_k-W^k{o}_k,k=0,1,...,N/2-1.$$

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

N=8일 때 전개:

$$\begin{array}{rl}& 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\\ & =[\sum _{n=0}^1{F}_{4n}(W^{4n}{)}^k+W^{2k}\sum _{n=0}^1{F}_{4n+2}(W^{4n}{)}^k]+W^k[\sum _{n=0}^1{F}_{4n+1}(W^{4n}{)}^k+W^{2k}\sum _{n=0}^1{F}_{4n+3}(W^{4n}{)}^k]\\ & =[(F_0+W^{4k}{F}_4)+W^{2k}(F_2+W^{4k}{F}_6)]+W^k[(F_1+W^{4k}{F}_5)+W^{2k}(F_3+W^{4k}{F}_7)]\\ & =[(F_{000}+W^{4k}{F}_{100})+W^{2k}(F_{010}+W^{4k}{F}_{110})]\\ & +W^k[(F_{001}+W^{4k}{F}_{101})+W^{2k}(F_{011}+W^{4k}{F}_{111})]\end{array}$$

$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 결과:

real part (imaginary part = 0)

FFT2D

FFT2D