FFT 구현

$$X_m = \sum_{n=0}^{N/2-1} x_{2n} W_N^{mn} +W_N^m \sum _{n=0}^{N/2-1} x_{2n+1} W_N^{mn}$$

$$X_{m+N/2} = \sum_{n=0}^{N/2-1} x_{2n} W_N^{mn} -W_N^m \sum _{n=0}^{N/2-1} x_{2n+1} W_N^{mn}$$

$$m=0,1,...,N/2-1$$

 

Butterfly Structure:

FFT for N=8:

/* sin(), cos() 함수를 $\log_2(N/4)$번만 호출하면 된다. 이마저도 2배각 공식을 쓰면 sqrt() 호출로 대체할 수 있다 */

#define PI	3.1415926535897932384
#define SWAP(a, b) {double temp = (a); (a) = (b); (b)=temp;}
int fft(int n, double x[], double y[]) {
    if ((n !=0) && ((n & (n-1)) != 0) ) return 0;
    // bit-reversal;
    for (int i = 0, j = 0; i < n-1; i++) {
        if (i < j) {
            SWAP(x[i], x[j]); SWAP(y[i], y[j]);
        }
        int k = n >> 1;
        while (k <= j) {
            j -= k; k >>= 1;
        }
        j += k;
    }
    // Danileson-Lanczos section;
    double rotwr = -1, rotwi = 0;
    for (int loop = 1; loop < n; loop <<= 1) {	
        int block = loop << 1;				// block-size;2->4->8->...
        if (loop > 2) {
            // rotation factor of twiddle factor;
            double theta = PI / loop;
            rotwr = cos(theta);
            rotwi = -sin(theta);
        } else if (loop == 2) {
            rotwr = 0;
            rotwi = -1;
        }	
        // starting twiddle factor;
        double wr = 1;
        double wi = 0;
        for (int j = 0; j < loop; ++j) { 
            // 각 block의 같은 위치에서 동일한 twiddle factor가 곱해진다는 
            // 사실을 이용함;
            for (int i = j; i < n; i += block) {
                int ip = i + loop;
                // step-1; X[ip]*W
                double xwre = x[ip] * wr - y[ip] * wi;
                double xwim = x[ip] * wi + y[ip] * wr;
                // step-2; X[ip] = X[i] - X[ip]*W;
                x[ip] = x[i] - xwre;
                y[ip] = y[i] - xwim;
                // step-3; X[i] = X[i] + X[ip]*W;
                x[i] += xwre;
                y[i] += xwim;
            };
            // 각 block의 다음 차례에 곱해지는 twiddle factor 계산; T = W * rotW;
            double tre = wr * rotwr - wi * rotwi;
            double tim = wi * rotwr + wr * rotwi;
            wr = tre;
            wi = tim;
        }
    }
    return 1;
    // 역변환: (1/N)*conjugate[FFT[conjugate(X)]] 
}