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
Let me demonstrate it below:
See the following example with 1d FFTs over a 2d
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?