I was trying to interface Python with C, and at some point I had to pass some numpy arrarys after being Fourier transformed. I was in pure despair, literally, for a whole week, facing a bug that was making no sense. In the end, I found out that the bug was introduced by a change in data-layout (strides) of the numpy array caused by the np.fft.fft operation.
Let me demonstrate it below:
See the following example with 1d FFTs over a 2d numpy array.
This is just a copy-paste of the python interpreter on my Linux PC:
Python 3.8.8 (default, Apr 13 2021, 19:58:26)
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> input = np.zeros((101,50), dtype=np.complex64)
>>> input.real = np.random.uniform(size=(101,50)).astype(np.float32)
>>> input.imag = np.random.uniform(size=(101,50)).astype(np.float32)
>>> input.shape
(101, 50)
>>> input.strides
(400, 8)
>>> output = np.fft.fft(input, axis=0)
>>> output.shape
(101, 50)
>>> output.strides
(16, 1616)
Here I basically create a 2d numpy array of single-precision complex floats, and then apply FFTs over axis=0. Look at the strides of the input and the output array!!! The shapes are the same but the strides differ.
I have no idea what kind of FFT algorithm is this!
Had I used scipy.fftpack.fft for this operation, the strides would be as expected:
>>> import scipy.fftpack
>>> output = scipy.fftpack.fft(input, axis=0)
>>> output.shape
(101, 50)
>>> output.strides
(400, 8)
So, in the numpy case, it seems like a virtual FFT operation... where instead of FLOPS, some elements reordering is being performed.
I solved my problem using scipy, but I'm curious to learn more. Does anyone have some comments?
