implementing Prime-factor FFT algorithm

Question

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):
    """
    radix-2,3,5 FFT algorithm
    """
    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)))

Answer 1

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):
    """
    radix-2,3 FFT algorithm
    """
    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)))

```