FFT using recursive python function

Question

I am trying to use the following code for finding FFT of a given list.

After a lot of trials I have found that this code runs only for an input list having 2^m or 2^m+1 elements.

Can you please clarify why this is so and whether it can be modified to use an input list containing some other number of elements. (P.S. I have an input list of 16001 elements)

    from cmath import exp, pi

    def fft(x):
        N = len(x)
        if N <= 1: return x
        even = fft(x[0::2])
        odd =  fft(x[1::2])
        T= [exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)]
        return [even[k] + T[k] for k in xrange(N/2)] + \
        [even[k] - T[k] for k in xrange(N/2)]

    print( ' '.join("%5.3f" % abs(f) 
            for f in fft([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])) )

Edit 1 Could You please tell the difference between the previous and the following function definition:

def fft(x):
    N = len(x)
    T = exp(-2*pi*1j/N)
    if N > 1:
        x = fft(x[::2]) + fft(x[1::2])
        for k in xrange(N/2):
            xk = x[k]
            x[k] = xk + T**k*x[k+N/2]
            x[k+N/2] = xk - T**k*x[k+N/2]
    return x

Edit 2: In fact this code(under Edit 1) does work, (sorry for the fault in indentation and variable naming earlier) which is why I want to understand the difference between the two.(this works for 16001 elements too!)

Answer 1

Algorithm

This algorithm works only for a power of two number of input data. Have a look at the theory.

This is a improved version that checks for this pre-condition:

from __future__ import print_function

import sys

if sys.version_info.major < 3:
    range = xrange

from cmath import exp, pi

def fft(x):
    N = len(x)
    if N <= 1: 
        return x
    if N % 2 > 0:
        raise ValueError("size of x must be a power of 2")
    even = fft(x[::2])
    odd =  fft(x[1::2])
    r = range(N//2)
    T = [exp(-2j * pi * k / N) * odd[k] for k in r]
    [even[k] for k in r]
    res = ([even[k] + T[k] for k in r] +
           [even[k] - T[k] for k in r])
    return res

input_data = [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]
print(' '.join("%5.3f" % abs(f) for f in fft(input_data)))

NumPy comes with its optimized version of a FFT function. It works for sizes different than the power of two. But it will be much slower because it needs to apply a different algorithm.

For example:

from numpy.fft import fft as np_fft
​
n = 1000
a = np.arange(n ** 2)
b = np.arange(n ** 2 - 10)

The runtime for the power-of-two case:

%timeit np_fft(a)
10 loops, best of 3: 59.3 ms per loop

is much less than that where this not the case:

%timeit np_fft(b)
1 loops, best of 3: 511 ms per loop

Recursion Limit

Python has built-in recursion limit of 1000:

>>> import sys
>>> sys.getrecursionlimit()
1000

But you can increase the recursion limit:

sys.setrecursion(50000)

The docs tell you why:

getrecursionlimit()

Return the current value of the recursion limit, the maximum depth of the Python interpreter stack. This limit prevents infinite recursion from causing an overflow of the C stack and crashing Python. It can be set by setrecursionlimit().

setrecursionlimit()

Set the maximum depth of the Python interpreter stack to limit. This limit prevents infinite recursion from causing an overflow of the C stack and crashing Python.

The highest possible limit is platform-dependent. A user may need to set the limit higher when they have a program that requires deep recursion and a platform that supports a higher limit. This should be done with care, because a too-high limit can lead to a crash.

Answer 2

Answer to the Edited Version

While this version:

from __future__ import print_function

import sys

if sys.version_info.major < 3:
    range = xrange

from cmath import exp, pi

def fft2(x):
    N = len(x)
    T = exp(-2*pi*1j/N)
    if N > 1:
        x = fft2(x[::2]) + fft2(x[1::2])
        for k in range(N//2):
            xk = x[k]
            x[k] = xk + T**k*x[k+N//2]
            x[k+N//2] = xk - T**k*x[k+N//2]
    return x

seems to work for a power-of-two number of inputs:

import numpy as np
from numpy.fft import fft as np_fft

data = [1, 2, 3, 4]
np.allclose(fft2(data), np_fft(data))

is True.

It does not give correct results for a different number of inputs.

data2 = [1, 2, 3, 4, 5]
np.allclose(fft2(data2), np_fft(data2))

is False.

It is still based on the assumption that the number of inputs is a power of two, even though it does not throw an exception.

Answer 3

in the above code in this line:

T = [exp(-2j * pi * k / N) * odd[k] for k in r]

if you observe closely at this part: exp(-2j * pi * k / N) would lead all your output vectors in clockwise direction, for the right anticlockwise answer replace exp(-2j * pi * k / N) with exp(2j * pi * k / N) as omega = exp(2j * pi * k / N)