# implementing Prime-factor FFT algorithm

I've been trying for days to implement this algorithm to work with size N samples but I can't manage to do it. my goal is to compute FFT for 100 samples, so I need factor 5 and 2, I wrote a simple FFT function and a prime factor function but I don't really understand how to go beyond that and articles with complex math aren't helping..

import numpy as np
import matplotlib.pyplot as plt

def prime_factor(n):
"""
prime factorization
"""
i = 2
factors = []
while i * i <= n:
if n % i:
i += 1
else:
n //= i
factors.append(i)
if n > 1:
factors.append(n)
return factors

def fft(x):
"""
FFT algorithm
"""
N = len(x)
if N <= 1:
return x
else:
even = fft(x[0::2])
odd =  fft(x[1::2])
T = [np.exp(-2j*np.pi*k/N)*odd[k] for k in range(N//2)]
return [even[k] + T[k] for k in range(N//2)] + \
[even[k] - T[k] for k in range(N//2)]

def main():
"""
main function
"""
x = np.random.rand(100)
y = fft(x)
plt.plot(y)

main()


UPDATE: Bob solution and my the final Radix-5 implementation:

import numpy as np

def fft(x):
"""
"""
N = len(x)
if N <= 1:
return x
elif N % 2 == 0:
# For multiples of 2 this formula works
even = fft(x[0::2])
odd =  fft(x[1::2])
T = [np.exp(-2j*np.pi*k/N)*odd[k] for k in range(N//2)]
return [even[k] + T[k] for k in range(N//2)] + \
[even[k] - T[k] for k in range(N//2)]
elif N % 3 == 0:
# Optional, implementing factor 3 decimation
p0 = fft(x[0::3])
p1 = fft(x[1::3])
p2 = fft(x[2::3])
# This will construct the output output without the simplifications
# you can do explorint symmetry
return [p0[k % (N//3)] +
p1[k % (N//3)] * np.exp(-2j*np.pi*k/N) +
p2[k % (N//3)] * np.exp(-4j*np.pi*k/N)
for k in range(N)]
elif N % 5 == 0:
# Here you must implement the factor 5 decimation
# start following the template for the factor 3 implementation given above
p0 = fft(x[0::5])
p1 = fft(x[1::5])
p2 = fft(x[2::5])
p3 = fft(x[3::5])
p4 = fft(x[4::5])

return [p0[k % (N//5)] +
p1[k % (N//5)] * np.exp(-2j*np.pi*k/N) +
p2[k % (N//5)] * np.exp(-4j*np.pi*k/N) +
p3[k % (N//5)] * np.exp(-6j*np.pi*k/N) +
p4[k % (N//5)] * np.exp(-8j*np.pi*k/N)
for k in range(N)]

x = np.random.rand(100) # 2 * 2 * 5 * 5
print(np.allclose(fft(x), np.fft.fft(x)))

• And you aren't going to use numpy.fft because why? Oct 10 at 19:18
• Because I actually need to implement this in CPP, the function itself, I would LOVE to use numpy otherwise Oct 10 at 19:20
• And you are not going to use fftw or some other existing C, C++, or Fortran-based math package because why? It's not that I'm averse to helping, but it's 2021: aside from self-learning or super-reliable systems, you need strong reasons to not use pre-made math packages for things like this. Oct 10 at 19:29
• I 100% agree with you, my professor told me it's better if I implement it myself that way I can learn.. I dived abit into this and it's super complicated, I'm going to give it few lasts try before I head over to FFTW Oct 10 at 19:35
• However, you write Articles with complex math aren't helping...: Well, guess why your professor gave you this homework: It's not about implementing a fast Fourier transform at all. You have a professor, which means you're writing some kind of thesis on signal processing, and complex math still irritates you: you need to get behind that. Being able to work on this until you solve this, on your own was the whole point of giving you this assignment. If you have problems with that, I strongly recommend making an appointment with your professor and asking them well-prepared … Oct 10 at 20:19

The problem with your implementation is that you are doing only the radix-2 decimation, that splits the input vectors (the inverse of interleaving), and then concatenate the results. Notice that they user the fact that $$e^{-1j\pi}=-1$$ to do only $$N/2$$ multiplications instead of $$N$$. Maybe this simplification made more difficult to see the link from this to the general formula.

Here I give you an implementation of the radix-2,3 FFT, the radix-3 decimation is left without without simplification. I hope from this you can write the radix-5 decimation more easily.

Code

import numpy as np

def fft(x):
"""
"""
N = len(x)
if N <= 1:
return x
elif N % 2 == 0:
# For multiples of 2 this formula works
even = fft(x[0::2])
odd =  fft(x[1::2])
T = [np.exp(-2j*np.pi*k/N)*odd[k] for k in range(N//2)]
return [even[k] + T[k] for k in range(N//2)] + \
[even[k] - T[k] for k in range(N//2)]
elif N % 3 == 0:
# Optional, implementing factor 3 decimation
p0 = fft(x[0::3])
p1 = fft(x[1::3])
p2 = fft(x[2::3])
# This will construct the output output without the simplifications
# you can do explorint symmetry
return [p0[k % (N//3)] +
p1[k % (N//3)] * np.exp(-2j*np.pi*k/N) +
p2[k % (N//3)] * np.exp(-4j*np.pi*k/N)
for k in range(N)]
elif N % 5 == 0:
# Here you must implement the factor 5 decimation
# start following the template for the factor 3 implementation given above
raise NotImplementedError("5-factor decimation not given")

x = np.random.rand(36) # 2 * 2 * 3 * 3
assert(np.allclose(fft(x), np.fft.fft(x)))

$$$$
`
• Thank you so much, that clearly showed me how exactly I need decimate, I think I managed to implement the factor 5 decimation, I added it to my post if you want to see but it seems like it works :) Oct 11 at 10:36