FFT

fft.py

# coding: utf8
# ref:
#  http://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/
#  http://zh.wikipedia.org/zh-cn/库利－图基快速傅里叶变换算法#.E6.99.82.E5.9F.9F.E6.8A.BD.E5.8F.96.E6.B3.95
import numpy as np

def FFT(x):
    x = np.asarray(x, dtype=float)
    n = x.shape[0]

    if n&(n-1) or n==0:
        raise ValueError("Size of array must be a power of 2")

    if n <= 32:
        j = np.arange(n)
        i = j.reshape((n, 1))
        M = np.exp(-2j * np.pi * i * j / n)
        return np.dot(M, x)

    X_even = FFT(x[::2])
    X_odd = FFT(x[1::2])
    factor = np.exp(-1j * np.pi * np.arange(n/2) / (n/2)) # half
    return np.concatenate([X_even + factor * X_odd,
                           X_even - factor * X_odd])

x = np.random.random(1024)
print np.allclose(FFT(x), np.fft.fft(x))