[PYTHON] FFT (Fast Fourier Transform): Formulas and Implementation Examples for Implementation

FFT principles and formulas

There was no article focusing on [Formula → Implementation], so I will spell it. Introducing this time is the simplest one called radix-2 of Cooley-Tukey type FFT.

First, the DFT (Discrete Fourier Transform) is expressed by the following equation. $F(t) = \displaystyle \sum_{x=0}^{N-1}f(x) e^{-i \frac{2πtx}{N}} \qquad t=0,1,...,N-1\tag{1}$ Here, $ N $ is the input size and is limited to powers of ** 2. ** ** Set $ W_N = e ^ {-i \ frac {2π} {N}} $ and rewrite as follows. $F(t) = \displaystyle \sum_{x=0}^{N-1}f(x) W_{N}^{tx}\qquad t=0,1,...,N-1\tag{2}$ The DFT performed on only the even-numbered elements of the input is expressed by the following formula. $F_{even}(t) = \displaystyle \sum_{k=0}^{\frac{N}{2}-1}f(2k) W_\frac{N} {2}^{tk}\tag{3}$ Similarly, the DFT performed on only the odd-numbered elements is expressed by the following formula. $F_{odd}(t) = \displaystyle \sum_{k=0}^{\frac{N}{2}-1}f(2k+1) W_\frac{N}{2}^{tk}\tag{4}$ By using $ (3) (4) $, $ (2) $ can be rewritten as follows. $F(t) = F_{even}(t) +W_N^t F_{odd}(t)\tag{5}$ So far, you've just transformed the DFT into odd and odd parts. Now, here is the important point of FFT. In fact, $ F_ {even} (t), F_ {odd} (t), W_N ^ t $ have the following periodicity.

F_{even}(t) = F_{even}(t+\frac{N}{2})\tag{6} F_{odd}(t) = F_{odd}(t+\frac{N}{2})\tag{7} W_N^t = -W_N^{t+\frac{N}{2}}\tag{8}

By using $ (6) (7) (8) $ for $ (5) $, the DFT can finally be rewritten as follows.

F(t) = F_{even}(t) +W_N^t F_{odd}(t)\qquad t=0,1,...,\frac{N}{2}-1\tag{9} F(t+\frac{N}{2}) = F_{even}(t) - W_N^tF_{odd}(t)\qquad t=0,1,...,\frac{N}{2}-1\tag{10}

Since $ F_ {even} (t) $ and $ F_ {odd} (t) $ are common, $ F (t) $ and $ F (t + \ frac {N} {2}) $ are calculated at the same time. That's why you want it. This is the mechanism for reducing the amount of calculation.

Based on the above, let's move on to implementation.

FFT program

Let's write an FFT function in Python. Here's the basic idea. After dividing the input data into even elements and odd elements and performing FFT (recursive) like $ (3) (4) $, 0 ~ (N / 2)-using $ (9) (10) $ Calculate the 1st and N / 2 ~ N-1th. Therefore, the FFT is a recursive function.

def FFT(f):
    #f:Input data of size N
    N = len(f)
    
    if N == 1:    #When N is 1, the input data is returned as it is
        return f[0]

    f_even = f[0:N:2]    #Even-numbered element of f
    f_odd = f[1:N:2]    #Odd element of f
    F_even = FFT(f_even)    #(3)FFT on even-numbered elements
    F_odd = FFT(f_odd)    #(4)FFT on even-numbered elements

    W_N = np.exp(-1j * (2 * np.pi * np.arange(0, N // 2)) / N)    #t is 0~N/2-Array that calculated W up to the first

    F = np.zeros(N, dtype ='complex')    #FFT output
    F[0:N//2] = F_even + W_N * F_odd    #(9)Calculate(t:0~N/2-1)
    F[N//2:N] = F_even - W_N * F_odd    #(10)Calculate(t:N/2~N-1)

    return F

Operation check

Let's check if we can get the same output as numpy FFT. Create a composite wave of two sine waves with the following parameters and use this as an input.

import numpy as np
import matplotlib.pyplot as plt

N = 64 #The number of data
n = np.arange(N)
f1 = 2 #Cycle 1
f2 = 5 #Cycle 2
A1 = 3 #Amplitude 1
A2 = 5 #Amplitude 2
INPUT_f = A1 * np.sin(f1 * 2 * np.pi * (n/N)) + A2 * np.sin(f2 * 2 * np.pi * (n/N)) #Synthesis of two sine waves

#display
plt.figure(figsize=(10, 5))
plt.style.use('seaborn')
plt.plot(INPUT_f)

Perform an FFT. Since the output of the FFT is a complex number, the amplitude component is extracted by multiplying it by the absolute value and visualized here.

OUTPUT_F_numpy = np.fft.fft(INPUT_f)    #numpy version
OUTPUT_F_original = FFT(INPUT_f)    #Self-made version

plt.figure(figsize=(10, 5))

plt.subplot(1,2,1)    #Amplitude component display(numpy version)
plt.title('numpy_ver')
plt.plot(np.abs(OUTPUT_F_numpy))

plt.subplot(1,2,2)    #Amplitude component display(Self-made version)
plt.title('original')
plt.plot(np.abs(OUTPUT_F_original))

I got the same result safely.