使用python手寫FFT算法

FFT(Fast Fourier Transform) 是 DFT(Discrete Fourier Transform)的快讀實作，它在機理上沒有改變DFT的算法，隻是在實作上采用的巧妙的實作。 使 O ( N 2 ) O(N^2) O(N2)的實作變成了 O ( N l o g 2 N ) O(Nlog_2N) O(Nlog2​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表示

X k = ∑ n = 0 N − 1 x n e − j 2 π k n N = ∑ n = 0 N 2 − 1 x 2 r + 1 e − j 2 π k ( 2 r + 1 ) N + ∑ n = 0 N 2 − 1 x 2 r e − j 2 π k ( 2 r ) N = ∑ n = 0 N 2 − 1 x 2 r + 1 e − j 4 π k r N e − j 2 π k N + ∑ n = 0 N 2 − 1 x 2 r e − j 4 π k r N = ∑ n = 0 N 2 − 1 x 2 r + 1 e − j 2 π k r N / 2 e − j 2 π k N + ∑ n = 0 N 2 − 1 x 2 r e − j 2 π k r N / 2 = e − j 2 π k N X o d d + X e v e n X_k = \sum_{n=0}^{N-1} x_n e^{\frac{-j2{\pi}kn}{N}} \newline = \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}} \newline = \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}} \newline = \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{-j{2\pi}kr}{N/2}} \newline = e^{\frac{-j2{\pi}k}{N}} X_{odd} + X_{even} Xk​=n=0∑N−1​xn​eN−j2πkn​=n=0∑2N​−1​x2r+1​eN−j2πk(2r+1)​+n=0∑2N​−1​x2r​eN−j2πk(2r)​=n=0∑2N​−1​x2r+1​eN−j4πkr​eN−j2πk​+n=0∑2N​−1​x2r​eN−j4πkr​=n=0∑2N​−1​x2r+1​eN/2−j2πkr​eN−j2πk​+n=0∑2N​−1​x2r​eN/2−j2πkr​=eN−j2πk​Xodd​+Xeven​

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