Below is a Python script which contains functions for computing the naive Discrete Fourier Transform (DFT), a recursive radix-2 Fast Fourier Transform (FFT) algorithm, and an in-place FFT radix-2 algorithm.
TODO: add timing information for FFT in-place; add timeplot output; compare outputs with FFT in-place; reorganise structure of Gist
import numpy as np
import matplotlib.pyplot as plt
from time import perf_counter
def dft(x):
""" Compute DFT naively using matrix multiplication """
N = x.size
inds = np.arange(N)
coeffs = np.exp(-2j * np.pi * np.outer(inds, inds) / N)
return np.matmul(coeffs, x)
def fft(x):
"""
Wrapper for recursive radix-2 FFT computation, which computes all of the
odd-term rotations, then calls the recursive sub-function
"""
N = x.size
if np.floor(np.log2(N)) != np.log2(N):
raise ValueError("Array size must be a power of 2")
M = np.arange(N / 2) # Indices
rotations = np.exp(-2j * np.pi * M / N) # Odd-term rotations
return fft_recurse(x, rotations)
def fft_recurse(x, rotations):
"""
Sub-function for recursive radix-2 FFT computation, which appropriately
combines FFTs of even terms and odd terms using rotations
"""
N = x.size
if N == 1: return x
if N == 2: return np.array([x[0] + x[1], x[0] - x[1]])
E = fft_recurse(x[0::2], rotations[0::2]) # Even terms
O = fft_recurse(x[1::2], rotations[0::2]) # Odd terms
O_rot = O * rotations # Odd terms rotated
return np.concatenate([E + O_rot, E - O_rot])
def fft_inplace(x):
"""
Implementation of radix-2 in-place (without recursion) on the input vector
x. The length of x must be an integer power of 2.
"""
# Initialise number of samples (N) and initial number of blocks (N / 2)
N = x.size
num_blocks = N >> 1
# Do initial 2-point FFT on all samples and convert the answer to complex
for i in range(num_blocks):
x[[i, i + num_blocks]] = (
x[i] + x[i + num_blocks],
x[i] - x[i + num_blocks]
)
x = x.astype(complex)
# Initialise odd-term rotations
M = np.arange(num_blocks)
rotations = np.exp((-2j * np.pi / N) * M)
# Update number of blocks and previous block size for next iteration
num_blocks >>= 1
prev_block_size = 2
# Iterate through each stage (each stage has a different number of blocks)
while prev_block_size < N:
# Get the rotations for each block of this stage
block_rotations = rotations[::num_blocks]
# Iterate through each block of this stage
for block_start_ind in range(num_blocks):
# Get the indices and even and rotated odd terms for this block
block_inds = np.arange(block_start_ind, N, num_blocks)
even_terms = x[block_inds[0::2]]
odd_terms = x[block_inds[1::2]]
odd_terms_rotated = odd_terms * block_rotations
# Calculate FFT of this block using blocks from previous stage
x[block_inds[:prev_block_size]] = even_terms + odd_terms_rotated
x[block_inds[prev_block_size:]] = even_terms - odd_terms_rotated
# Next iteration there are half as many blocks, each double the size
num_blocks >>= 1
prev_block_size <<= 1
return x
def ifft(X):
""" Calculate inverse FFT """
return np.conj(fft(np.conj(X))) / X.size
def plot_sig_and_response(t, x, f, X, name="Signal and response"):
fig, axes = plt.subplots(2, 1)
fig.set_size_inches([8, 6])
axes[0].plot(t, x, "b")
axes[0].grid(True)
axes[0].set(xlabel="Time (s)", title="Time domain")
axes[1].plot(f, np.abs(X), "b")
axes[1].grid(True)
axes[1].set(xlabel="Frequency (Hz)", title="Frequency domain")
fig.suptitle(name, fontsize=20)
fig.tight_layout(rect=[0, 0, 1, 0.95])
plt.savefig(name)
plt.close()
if __name__ == "__main__":
# Sampling frequency
f_s = 1e3
# Make plots for N = 256
N = 256
t = np.linspace(0, N/f_s, N, endpoint=False)
x = np.sin(2*np.pi*100*t) + np.sin(2*np.pi*200*t)
X_dft = dft(x)
X_fft = fft(x)
f = np.linspace(0, f_s, N, endpoint=False)
plot_sig_and_response(t, x, f, X_dft, name="DFT results")
plot_sig_and_response(t, x, f, X_fft, name="FFT results")
# Check error of inverse FFT
x_ifft = ifft(X_fft)
print("Inverse FFT error = {:.4e}".format(np.max(np.abs(x - x_ifft))))
# Print table headers for speed tests
table_headers = ["Input size", "Time taken for DFT (s)",
"Time taken for FFT (s)", "Factor difference", "Maximum absolute error"]
print(" | ".join(table_headers))
print(" | ".join(["---"] * len(table_headers)))
# Run speed tests (TODO: make log-log graph of speed results with repeats)
for N_pow in range(8, 14):
N = 2 ** N_pow
t = np.linspace(0, N/f_s, N, endpoint=False)
x = np.sin(2*np.pi*100*t) + np.sin(2*np.pi*200*t)
t0 = perf_counter()
X_dft = dft(x)
t1 = perf_counter()
X_fft = fft(x)
t2 = perf_counter()
print(" | ".join(["{}".format(N), "{:.4g}".format(t1 - t0),
"{:.4g}".format(t2 - t1), "{:.1f}".format((t1 - t0) / (t2 - t1)),
"{:.4e}".format(np.max(np.abs(X_dft - X_fft)))]))Inverse FFT error = 9.4717e-16
Below are time and accuracy results for different sized inputs:
| Input size | Time taken for DFT (s) | Time taken for FFT (s) | Factor difference | Maximum absolute error |
|---|---|---|---|---|
| 256 | 0.004445 | 0.00093 | 4.8 | 3.6437e-12 |
| 512 | 0.01783 | 0.001716 | 10.4 | 1.4453e-11 |
| 1024 | 0.06661 | 0.003368 | 19.8 | 5.6465e-11 |
| 2048 | 0.2659 | 0.00662 | 40.2 | 2.1411e-10 |
| 4096 | 1.056 | 0.01365 | 77.3 | 7.7747e-10 |
| 8192 | 4.172 | 0.02662 | 156.7 | 3.4760e-09 |
The following code can be used to print out the indices for each block in each stage of the in-place FFT, which is useful for gaining intuition of how the in-place FFT has been coded:
import numpy as np
def fft_inplace_inds(num_stages):
# Initialise number of samples and initial block size
N = 1 << num_stages
block_size = N
# Initialise x
x = np.empty([N, num_stages], dtype=int)
x[:, 0] = np.arange(N)
# Iterate through each column (1 column per stage)
for col in range(1, num_stages):
prev_block_size = block_size
block_size >>= 1
# Iterate through each block in the column
for prev_block_start_i in range(0, N, prev_block_size):
# Calculate even indices
x[
prev_block_start_i:(
prev_block_start_i + block_size
),
col
] = x[
prev_block_start_i:(
prev_block_start_i + prev_block_size
):2,
col - 1
]
# Calculate odd indices
x[
(prev_block_start_i + block_size):(
prev_block_start_i + prev_block_size
),
col
] = x[
(prev_block_start_i + 1):(
prev_block_start_i + prev_block_size
):2,
col - 1
]
return x
if __name__ == "__main__":
for num_stages in range(1, 6):
print(
1 << num_stages,
fft_inplace_inds(num_stages),
sep=":\n",
end="\n\n"
)Console outout:
2:
[[0]
[1]]
4:
[[0 0]
[1 2]
[2 1]
[3 3]]
8:
[[0 0 0]
[1 2 4]
[2 4 2]
[3 6 6]
[4 1 1]
[5 3 5]
[6 5 3]
[7 7 7]]
16:
[[ 0 0 0 0]
[ 1 2 4 8]
[ 2 4 8 4]
[ 3 6 12 12]
[ 4 8 2 2]
[ 5 10 6 10]
[ 6 12 10 6]
[ 7 14 14 14]
[ 8 1 1 1]
[ 9 3 5 9]
[10 5 9 5]
[11 7 13 13]
[12 9 3 3]
[13 11 7 11]
[14 13 11 7]
[15 15 15 15]]
32:
[[ 0 0 0 0 0]
[ 1 2 4 8 16]
[ 2 4 8 16 8]
[ 3 6 12 24 24]
[ 4 8 16 4 4]
[ 5 10 20 12 20]
[ 6 12 24 20 12]
[ 7 14 28 28 28]
[ 8 16 2 2 2]
[ 9 18 6 10 18]
[10 20 10 18 10]
[11 22 14 26 26]
[12 24 18 6 6]
[13 26 22 14 22]
[14 28 26 22 14]
[15 30 30 30 30]
[16 1 1 1 1]
[17 3 5 9 17]
[18 5 9 17 9]
[19 7 13 25 25]
[20 9 17 5 5]
[21 11 21 13 21]
[22 13 25 21 13]
[23 15 29 29 29]
[24 17 3 3 3]
[25 19 7 11 19]
[26 21 11 19 11]
[27 23 15 27 27]
[28 25 19 7 7]
[29 27 23 15 23]
[30 29 27 23 15]
[31 31 31 31 31]]
Below are plots of the results for an input size of 256:
