Geometry & Recognition

Forward fft: 

$$ {\tt fft}(x+iy)=X + i Y = \sum(x +i y) (W_r + i W_i)\\=\sum  (x W_r - y W_i ) + i (x W_i + y W_r)$$

Inverse fft:

$$ {\tt ifft}(X+iY)= \frac{1}{N}\sum(X +i Y) (W_r - i W_i)\\= \frac{1}{N} \sum  (X W_r + Y W_i ) + i (-X W_i + Y W_r)\\ =\frac{1}{N}\sum \left(XW_r - (-Y) W_i\right) +(-i) \left(XW_i + (-Y)W_r\right) \\ = \frac{1}{N}\left[{\tt fft}(X-iY) \right]^* = \frac{1}{N}\left[{\tt fft}((X+iY)^*)\right]^*$$

또는 

$$ \frac{1}{N} {\tt fft}(Y+iX)= \frac{1}{N} \sum(Y +i X) (W_r + i W_i)\\=\frac{1}{N} \sum  (Y W_r - X W_i ) + i (Y W_i + X W_r)\\= y + i x$$

$$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)]] 
}