使用python手写FFT算法

FFT(Fast Fourier Transform) 是 DFT(Discrete Fourier Transform)的快读实现，它在机理上没有改变DFT的算法，只是在实现上采用的巧妙的实现。 使$O(N^2)$的实现变成了$O(Nlo{g}_{2}N)$的实现，优化算法的复杂度。

DFT 公式

DFT的公式我们先简单回忆一下：

DFT的python实现

def dft(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

x = np.random.random(1024)
np.allclose(dft(x), np.fft.fft(x))

---

FFT的思想其实分而治之，将整个DFT安装奇偶来分开计算

FFT之奇偶分治

首先根据DFT公式，将x分为odd(奇数),下标用2r+1表示；和even(偶数)下标用2r表示
$$\begin{gathered} X_k=\sum _{n=0}^{N-1}{x}_n{e}^{\frac{-j2\pi kn}{N}} \\ =\sum _{n=0}^{\frac{N}{2}-1}{x}_{2r+1}{e}^{\frac{-j2\pi k(2r+1)}{N}}+\sum _{n=0}^{\frac{N}{2}-1}{x}_{2r}{e}^{\frac{-j2\pi k(2r)}{N}} \\ =\sum _{n=0}^{\frac{N}{2}-1}{x}_{2r+1}{e}^{\frac{-j4\pi kr}{N}}{e}^{\frac{-j2\pi k}{N}}+\sum _{n=0}^{\frac{N}{2}-1}{x}_{2r}{e}^{\frac{-j4\pi kr}{N}} \\ =\sum _{n=0}^{\frac{N}{2}-1}{x}_{2r+1}{e}^{\frac{-j2\pi kr}{N/2}}{e}^{\frac{-j2\pi k}{N}}+\sum _{n=0}^{\frac{N}{2}-1}{x}_{2r}{e}^{\frac{-j2\pi kr}{N/2}} \\ =e^{\frac{-j2\pi k}{N}}{X}_{odd}+X_{even} \end{gathered}$$

python实现

1. 根据奇偶分治和迭代的方式实现

def fft(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    if N % 2 > 0:
        raise ValueError("must be a power of 2")
    elif N <= 2:
        return dft(x)
    else:
        X_even = fft(x[::2])
        X_odd = fft(x[1::2])
        terms = np.exp(-2j * np.pi * np.arange(N) / N)
        return np.concatenate([X_even + terms[:int(N/2)] * X_odd,
                               X_even + terms[int(N/2):] * X_odd])
x = np.random.random(1024)
np.allclose(fft(x), np.fft.fft(x))

2. 使用矩阵运算加速

def fft_v(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    if np.log2(N) % 1 > 0:
        raise ValueError("must be a power of 2")
        
    N_min = min(N, 2)
    
    n = np.arange(N_min)
    k = n[:, None]
    M = np.exp(-2j * np.pi * n * k / N_min)
    X = np.dot(M, x.reshape((N_min, -1)))
	while X.shape[0] < N:
        X_even = X[:, :int(X.shape[1] / 2)]
        X_odd = X[:, int(X.shape[1] / 2):]
        terms = np.exp(-1j * np.pi * np.arange(X.shape[0])
                        / X.shape[0])[:, None]
        X = np.vstack([X_even + terms * X_odd,
                       X_even - terms * X_odd])
return X.ravel()

参考文档 How to implement the Fast Fourier Transform algorithm in Python from scratch