FFT in Numpy (Python) when N is not a power of 2

4
votes

My question is about the algorithm which is used in Numpy's FFT function.

The documentation of Numpy says that it uses the Cooley-Tukey algorithm. However, as you may know, this algorithm works only if the number N of points is a power of 2.

Does numpy pad my input vector x[n] in order to calculate its FFT X[k]? (I don't think so, since the number of points I have in the output is also N). How could I actually "see" the code which is used by numpy for its FFT function?

Cheers!

2
pythonnumpyfft
stackoverflow.com/questions/6855169/… - mmgp
Well... here's one way to actually see the code. If that's not what you mean, then perhaps you should rephrase. - senderle
The easy hack is to downsample to the next lowest power of 2. Maybe even do that multiple times choosing which points to exclude at random, and average the results together. - wberry

2 Answers

11
votes

The docs say numpy's FFT is based on FFTPACK.

In the FFTPACK docs I find the following:

subroutine rffti(n,wsave)

subroutine rffti initializes the array wsave which is used in both rfftf and rfftb. the prime factorization of n together with a tabulation of the trigonometric functions are computed and stored in wsave.

The standard Cooley-Tukey algorithm is "radix-2 with decimation in time", which recursively reduces the computation of an FFT of size 2*n into 2 FFTs of size n, plus n FFTs of size 2. There is a general factorization version of the same algorithm that turns an FFT of size m*ninto n FFTs of size m plus m FFTs of size n. The fact that the preparation routines in FFTPACK compute the prime factorization of the input size, seems to indicate this is what they are doing. So unless you go for a prime number of elements, or your number of elements has a very large prime factor, you should still get a pretty good speed-up.

A few years back I blogged about both the radix-2 and general factorization versions of Cooley-Tukey's algorithmn. Reading those may help understand what seems to be going on inside NumPy. The following picture, taken from there, depicts a CT FFT:

enter image description here

1
votes

In my experience the algorithms don't do automatic padding, or at least some of them don't. For example, running the scipy.signal.hilbert method on a signal that wasn't of length == a power of two took about 45 seconds. When I padded the signal myself with zeros to such a length, it took 100ms.

YMMV but it's something to double check basically any time you run a signal processing algorithm.