Dear SciPy devels, Please find a scipy.fftpack replacement candidate fftpack2 from http://cens.ioc.ee/~pearu/misc/fftpack2-0.1.tar.gz for review. fftpack2 main features include: - Optional support for high-performance FFT libraries such as FFTW and DJBFFT libraries. By default, Fortran FFTPACK is used and other libraries are automatically detected by the scipy_distutils/system_info.py script. - Support for both real and complex DFT's and their inverse transforms. When used with FFTW and/or DJBFFT libraries then fftpack2 routines are considerably faster than Numeric.FFT and the current scipy.fftpack. Some comparison timings are included at the end of this message. - Implementation of differentiation and integration of periodic sequences, various pseudo-differential operators such as Tilbert transform, Hilbert transform, hyperbolic pseudo-differential operators, their inverses, etc. - Because different FFT libraries use different data storage convention, fftpack2 implements interface to these conventions so that users and developers need not to learn all these conventions, just one is enough. Namely, fftpack2 uses the same data storage specification as Fortran FFTPACK (also used by Numeric.FFT and Matlab, though with somewhat different formal definitions). - Generic convolution routines that allow quick and easy way of introducing new pseudo-differential operators with the same performance as core pseudo-differential operators. - All functions have complete documentation strings. - A rather complete testing site. Note that fftpack2 requires the latest F2PY that you can get from http://cens.ioc.ee/projects/f2py2e/2.x/F2PY-2-latest.tar.gz For testing, fftpack2 uses scipy_test from SciPy CVS repository. For building/testing instructions and for various notes, including open questions, see NOTES.txt file after unpacking fftpack2-0.1.tar.gz. As usual, comments, suggestions, and criticism are most welcome. Regards, Pearu ------------------------------- fftpack2 timing results: ------------------------------- Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.38 | 0.39 | 0.38 | 0.40 | 0.40 | 0.38 (secs for7000calls) 1000 | 0.31 | 0.63 | 0.60 | 0.54 | 0.64 | 0.61 (secs for 2000 256 | 0.60 | 0.75 | 0.62 | 0.71 | 0.78 | 0.66 (secs for 10000 512 | 0.77 | 1.30 | 1.00 | 0.95 | 1.30 | 1.01 (secs for 10000 1024 | 0.14 | 0.26 | 0.18 | 0.17 | 0.25 | 0.19 (secs for 1000 2048 | 0.21 | 0.51 | 0.36 | 0.34 | 0.65 | 0.54 (secs for 1000 4096 | 0.28 | 1.08 | 0.66 | 0.56 | 0.93 | 0.79 (secs for 500 8192 | 1.09 | 3.88 | 3.94 | 2.19 | 4.36 | 3.48 (secs for 500 Inverse Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.17 | 0.26 | 0.63 | 0.45 | 0.87 | 0.85 (secs for 7000 1000 | 0.34 | 1.24 | 1.25 | 0.72 | 1.29 | 1.29 (secs for 2000 256 | 0.62 | 1.57 | 1.44 | 0.69 | 1.62 | 1.47 (secs for 10000 512 | 0.80 | 2.69 | 2.30 | 1.00 | 2.62 | 2.20 (secs for 10000 1024 | 0.12 | 0.51 | 0.41 | 0.18 | 0.54 | 0.44 (secs for 1000 2048 | 0.23 | 1.16 | 1.09 | 0.40 | 1.33 | 1.27 (secs for 1000 4096 | 0.35 | 1.75 | 1.44 | 0.63 | 1.72 | 1.54 (secs for 500 8192 | 1.34 | 5.71 | 6.02 | 2.40 | 6.14 | 5.67 (secs for 500 Multi-dimensional Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100x100 | 0.43 | 0.87 | 0.81 | 0.48 | 0.88 | 0.82 (secs for 100 1000x100 | 0.54 | 0.78 | 0.72 | 0.56 | 0.73 | 0.75 (secs for 7 256x256 | 0.69 | 0.72 | 0.62 | 0.67 | 0.72 | 0.64 (secs for 10 512x512 | 0.81 | 0.91 | 0.84 | 0.81 | 0.87 | 0.81 (secs for 3 Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.36 | 0.47 | 0.47 (secs for 7000 calls) 1000 | 0.29 | 0.38 | 0.38 (secs for 2000 calls) 256 | 0.55 | 0.76 | 0.81 (secs for 10000 calls) 512 | 0.67 | 1.06 | 1.07 (secs for 10000 calls) 1024 | 0.10 | 0.18 | 0.17 (secs for 1000 calls) 2048 | 0.16 | 0.30 | 0.30 (secs for 1000 calls) 4096 | 0.20 | 0.37 | 0.32 (secs for 500 calls) 8192 | 0.75 | 1.29 | 1.26 (secs for 500 calls) Inverse Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.39 | 0.90 | 0.92 (secs for 7000 calls) 1000 | 0.32 | 0.63 | 0.70 (secs for 2000 calls) 256 | 0.59 | 1.37 | 1.39 (secs for 10000 calls) 512 | 0.73 | 1.70 | 1.80 (secs for 10000 calls) 1024 | 0.11 | 0.27 | 0.27 (secs for 1000 calls) 2048 | 0.20 | 0.45 | 0.44 (secs for 1000 calls) 4096 | 0.24 | 0.65 | 0.68 (secs for 500 calls) 8192 | 0.73 | 1.98 | 2.06 (secs for 500 calls) Shifting periodic functions ============================== size | optimized | naive ------------------------------ 100 | 0.09 | 0.52 (secs for 1500 calls) 1000 | 0.06 | 0.66 (secs for 300 calls) 256 | 0.10 | 0.84 (secs for 1500 calls) 512 | 0.10 | 1.02 (secs for 1000 calls) 1024 | 0.07 | 1.07 (secs for 500 calls) 2048 | 0.06 | 0.80 (secs for 200 calls) 4096 | 0.00 | 0.86 (secs for 100 calls) 8192 | 0.13 | 1.23 (secs for 50 calls) Differentiation of periodic functions ===================================== size | convolve | naive ------------------------------------- 100 | 0.09 | 0.59 (secs for 1500 calls) 1000 | 0.06 | 0.75 (secs for 300 calls) 256 | 0.11 | 0.97 (secs for 1500 calls) 512 | 0.09 | 1.21 (secs for 1000 calls) 1024 | 0.08 | 1.19 (secs for 500 calls) 2048 | 0.06 | 1.11 (secs for 200 calls) 4096 | 0.09 | 1.19 (secs for 100 calls) 8192 | 0.14 | 1.48 (secs for 50 calls) Tilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.09 | 0.54 (secs for 1500 calls) 1000 | 0.06 | 0.60 (secs for 300 calls) 256 | 0.10 | 0.77 (secs for 1500 calls) 512 | 0.09 | 0.89 (secs for 1000 calls) 1024 | 0.07 | 0.90 (secs for 500 calls) 2048 | 0.05 | 0.83 (secs for 200 calls) 4096 | 0.07 | 0.92 (secs for 100 calls) 8192 | 0.11 | 1.04 (secs for 50 calls) . Hilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.09 | 0.47 (secs for 1500 calls) 1000 | 0.06 | 0.48 (secs for 300 calls) 256 | 0.10 | 0.68 (secs for 1500 calls) 512 | 0.09 | 0.77 (secs for 1000 calls) 1024 | 0.07 | 0.74 (secs for 500 calls) 2048 | 0.06 | 0.72 (secs for 200 calls) 4096 | 0.07 | 0.79 (secs for 100 calls) 8192 | 0.11 | 0.93 (secs for 50 calls)
Hey Pearu, I'll look at this tonight. Thanks, eric
-----Original Message----- From: scipy-dev-admin@scipy.net [mailto:scipy-dev-admin@scipy.net] On Behalf Of Pearu Peterson Sent: Saturday, September 28, 2002 3:54 PM To: scipy-dev@scipy.org Subject: [SciPy-dev] fftpack2 for review
Dear SciPy devels,
Please find a scipy.fftpack replacement candidate fftpack2 from
http://cens.ioc.ee/~pearu/misc/fftpack2-0.1.tar.gz
for review. fftpack2 main features include:
- Optional support for high-performance FFT libraries such as FFTW and DJBFFT libraries. By default, Fortran FFTPACK is used and other libraries are automatically detected by the scipy_distutils/system_info.py script.
- Support for both real and complex DFT's and their inverse transforms. When used with FFTW and/or DJBFFT libraries then fftpack2 routines are considerably faster than Numeric.FFT and the current scipy.fftpack. Some comparison timings are included at the end of this message.
- Implementation of differentiation and integration of periodic sequences, various pseudo-differential operators such as Tilbert transform, Hilbert transform, hyperbolic pseudo-differential operators, their inverses, etc.
- Because different FFT libraries use different data storage convention, fftpack2 implements interface to these conventions so that users and developers need not to learn all these conventions, just one is enough. Namely, fftpack2 uses the same data storage specification as Fortran FFTPACK (also used by Numeric.FFT and Matlab, though with somewhat different formal definitions).
- Generic convolution routines that allow quick and easy way of introducing new pseudo-differential operators with the same performance as core pseudo-differential operators.
- All functions have complete documentation strings.
- A rather complete testing site.
Note that fftpack2 requires the latest F2PY that you can get from
http://cens.ioc.ee/projects/f2py2e/2.x/F2PY-2-latest.tar.gz
For testing, fftpack2 uses scipy_test from SciPy CVS repository.
For building/testing instructions and for various notes, including open questions, see NOTES.txt file after unpacking fftpack2-0.1.tar.gz.
As usual, comments, suggestions, and criticism are most welcome.
Regards, Pearu
------------------------------- fftpack2 timing results: ------------------------------- Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.38 | 0.39 | 0.38 | 0.40 | 0.40 | 0.38 (secs for7000calls) 1000 | 0.31 | 0.63 | 0.60 | 0.54 | 0.64 | 0.61 (secs for 2000 256 | 0.60 | 0.75 | 0.62 | 0.71 | 0.78 | 0.66 (secs for 10000 512 | 0.77 | 1.30 | 1.00 | 0.95 | 1.30 | 1.01 (secs for 10000 1024 | 0.14 | 0.26 | 0.18 | 0.17 | 0.25 | 0.19 (secs for 1000 2048 | 0.21 | 0.51 | 0.36 | 0.34 | 0.65 | 0.54 (secs for 1000 4096 | 0.28 | 1.08 | 0.66 | 0.56 | 0.93 | 0.79 (secs for 500 8192 | 1.09 | 3.88 | 3.94 | 2.19 | 4.36 | 3.48 (secs for 500
Inverse Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.17 | 0.26 | 0.63 | 0.45 | 0.87 | 0.85 (secs for 7000 1000 | 0.34 | 1.24 | 1.25 | 0.72 | 1.29 | 1.29 (secs for 2000 256 | 0.62 | 1.57 | 1.44 | 0.69 | 1.62 | 1.47 (secs for 10000 512 | 0.80 | 2.69 | 2.30 | 1.00 | 2.62 | 2.20 (secs for 10000 1024 | 0.12 | 0.51 | 0.41 | 0.18 | 0.54 | 0.44 (secs for 1000 2048 | 0.23 | 1.16 | 1.09 | 0.40 | 1.33 | 1.27 (secs for 1000 4096 | 0.35 | 1.75 | 1.44 | 0.63 | 1.72 | 1.54 (secs for 500 8192 | 1.34 | 5.71 | 6.02 | 2.40 | 6.14 | 5.67 (secs for 500
Multi-dimensional Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100x100 | 0.43 | 0.87 | 0.81 | 0.48 | 0.88 | 0.82 (secs for 100 1000x100 | 0.54 | 0.78 | 0.72 | 0.56 | 0.73 | 0.75 (secs for 7 256x256 | 0.69 | 0.72 | 0.62 | 0.67 | 0.72 | 0.64 (secs for 10 512x512 | 0.81 | 0.91 | 0.84 | 0.81 | 0.87 | 0.81 (secs for 3
Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.36 | 0.47 | 0.47 (secs for 7000 calls) 1000 | 0.29 | 0.38 | 0.38 (secs for 2000 calls) 256 | 0.55 | 0.76 | 0.81 (secs for 10000 calls) 512 | 0.67 | 1.06 | 1.07 (secs for 10000 calls) 1024 | 0.10 | 0.18 | 0.17 (secs for 1000 calls) 2048 | 0.16 | 0.30 | 0.30 (secs for 1000 calls) 4096 | 0.20 | 0.37 | 0.32 (secs for 500 calls) 8192 | 0.75 | 1.29 | 1.26 (secs for 500 calls)
Inverse Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.39 | 0.90 | 0.92 (secs for 7000 calls) 1000 | 0.32 | 0.63 | 0.70 (secs for 2000 calls) 256 | 0.59 | 1.37 | 1.39 (secs for 10000 calls) 512 | 0.73 | 1.70 | 1.80 (secs for 10000 calls) 1024 | 0.11 | 0.27 | 0.27 (secs for 1000 calls) 2048 | 0.20 | 0.45 | 0.44 (secs for 1000 calls) 4096 | 0.24 | 0.65 | 0.68 (secs for 500 calls) 8192 | 0.73 | 1.98 | 2.06 (secs for 500 calls)
Shifting periodic functions ============================== size | optimized | naive ------------------------------ 100 | 0.09 | 0.52 (secs for 1500 calls) 1000 | 0.06 | 0.66 (secs for 300 calls) 256 | 0.10 | 0.84 (secs for 1500 calls) 512 | 0.10 | 1.02 (secs for 1000 calls) 1024 | 0.07 | 1.07 (secs for 500 calls) 2048 | 0.06 | 0.80 (secs for 200 calls) 4096 | 0.00 | 0.86 (secs for 100 calls) 8192 | 0.13 | 1.23 (secs for 50 calls)
Differentiation of periodic functions ===================================== size | convolve | naive ------------------------------------- 100 | 0.09 | 0.59 (secs for 1500 calls) 1000 | 0.06 | 0.75 (secs for 300 calls) 256 | 0.11 | 0.97 (secs for 1500 calls) 512 | 0.09 | 1.21 (secs for 1000 calls) 1024 | 0.08 | 1.19 (secs for 500 calls) 2048 | 0.06 | 1.11 (secs for 200 calls) 4096 | 0.09 | 1.19 (secs for 100 calls) 8192 | 0.14 | 1.48 (secs for 50 calls)
Tilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.09 | 0.54 (secs for 1500 calls) 1000 | 0.06 | 0.60 (secs for 300 calls) 256 | 0.10 | 0.77 (secs for 1500 calls) 512 | 0.09 | 0.89 (secs for 1000 calls) 1024 | 0.07 | 0.90 (secs for 500 calls) 2048 | 0.05 | 0.83 (secs for 200 calls) 4096 | 0.07 | 0.92 (secs for 100 calls) 8192 | 0.11 | 1.04 (secs for 50 calls) . Hilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.09 | 0.47 (secs for 1500 calls) 1000 | 0.06 | 0.48 (secs for 300 calls) 256 | 0.10 | 0.68 (secs for 1500 calls) 512 | 0.09 | 0.77 (secs for 1000 calls) 1024 | 0.07 | 0.74 (secs for 500 calls) 2048 | 0.06 | 0.72 (secs for 200 calls) 4096 | 0.07 | 0.79 (secs for 100 calls) 8192 | 0.11 | 0.93 (secs for 50 calls)
_______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.net http://www.scipy.net/mailman/listinfo/scipy-dev
Hey Pearu, Overall, first tests went fairly smoothly. Even without djbfft and fftw, things appeared to install correctly on win32. *very* nice job on system checking and handling the lack of libraries. All tests passed. Sure am glad to finally have some unit tests for fft stuff! I do not have fftw installed. I tested once without djbfft and once with it. I didn't notice any material difference in timings (the non-djbfft set given below). My tests are on a 850 MHz PIII and they are slower than yours. This leads me to believe, even when trying to use djbfft, I'm not getting it linked correctly. Setup said it detected my djbfft libraries though (I defined the environment variable pointing to my libs to get this to work). I'll look over the tests tomorrow. Also, Travis O., could you give the package a once over? I've also copied your "proposal" here to solicit comments from the signal processing geeks: Differences between fftpack and FFT from Numeric ================================================ * Functions rfft and irfft accept and return only real sequences. So, the corresponding functions real_fft, inverse_real_fft from FFT are not equivalent with rfft and irfft. The difference is in the storage of data, see the definitions of corresponding functions for details. * Functions fftshift and ifftshift from Numeric.FFT are renamed to freqshift and ifreqshift, respectively. Rationale: there is nothing FFT algorithm specific what these functions aim to accomplish. Also freqshift puts the Nyquist component at the end of the result while fftshift puts it at the beginning of the result. This is for the consistency among varios functions in the fftpack module. * PROPOSAL: When calling ifft with forced truncation or zero-padding then I would like to propose that the corresponding action is applied to the middle of data. For example, ifft([1,2,3,4,5,6],n=8) is equivalent to ifft([1,2,3,4,0,0,5,6]), that is, the Fourier terms with higher frequencies and zero coefficients are introduced. In the Numeric.FFT case, the example above would be equivalent to ifft([1,2,3,4,5,6,0,0],n=8), which would mean that Fourier coefficients [5,6] become the coefficients of higher frequency terms and the original terms are zerod. Note that this proposal is **not implemented** because it needs to be discussed. For instance, Matlab uses the same convention as FFT and this change would be confusing for Matlab users. On the other hand, FFT or Matlab's conventions change the spectrum of the original signal and I don't see any sense in this behaviour (if you don't agree then please provide an example). Namely, one of the applications of the argument n would be to compose a new signal with a more dense or sparse grid than the original one by using new_signal = ifft(fft(signal),n) Note that the new_signal would have the same Fourier spectrum as original signal. With Matlab/FFT convention this is not true. Any thoughts? And now some installation notes about what I did to get djbfft running: Install notes: 1. win32 (cygwin and mingw) required #include "errno.h" in any djbfft file that used errno. 2. I installed the latest f2py but still got the warning error: WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING! fftpack2 requires F2PY version 2.23.190-1367 or higher but got 2.23.190-1354 WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING! I'm hitting the hay right now, so I'll check out why this happened tomorrow. 3. How hard would it be to get setup.py to build the djbfft library if one wasn't detected. The makefile builds a lot of the files on the fly using stuff like: 4c0.c: \ idea/c0.c pentium/c0.c ppro/c0.c sparc/c0.c auto_opt sed 1s/PRE/pre4.c/ `cat auto_opt`/c0.c > 4c0.c 4c0.o: \ compile 4c0.c pre4.c fftc4.h complex4.h real4.h fftr4.h real4.h 4i.c ./compile 4c0.c I'm not sure what all is going on here, but it seems like we could either do the same in python or auto-generate a generic version of this stuff and include it in SciPy. Of course, since djbfft isn't required, this isn't necessary. TEST RESULTS:
fftpack2.test() creating test suite for: fftpack2.basic creating test suite for: fftpack2.helper creating test suite for: fftpack2.pseudo_diffs ................... Fast Fourier Transform ======================================================== | real input | complex input
size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.39 | 0.36 | 0.35 | 0.40 | 0.37 | 0.36 (secs for 7000 call 1000 | 0.31 | 0.48 | 0.43 | 0.44 | 0.48 | 0.41 (secs for 2000 call 256 | 0.65 | 0.75 | 0.70 | 0.75 | 0.75 | 0.70 (secs for 10000 cal 512 | 0.89 | 1.32 | 1.11 | 1.16 | 1.33 | 1.12 (secs for 10000 cal 1024 | 0.14 | 0.35 | 0.30 | 0.21 | 0.24 | 0.20 (secs for 1000 call 2048 | 0.50 | 0.76 | 0.64 | 0.49 | 0.77 | 0.68 (secs for 1000 call 4096 | 0.60 | 0.93 | 0.82 | 0.67 | 0.79 | 0.65 (secs for 500 calls 8192 | 1.23 | 3.16 | 2.85 | 2.88 | 3.23 | 2.91 (secs for 500 calls . Inverse Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.38 | 0.88 | 0.87 | 0.45 | 0.89 | 0.88 (secs for 7000 call 1000 | 0.45 | 1.93 | 1.85 | 0.80 | 1.44 | 1.36 (secs for 2000 call 256 | 0.66 | 3.07 | 2.85 | 0.81 | 1.79 | 1.74 (secs for 10000 cal 512 | 0.91 | 2.89 | 2.78 | 1.27 | 3.22 | 3.23 (secs for 10000 cal 1024 | 0.24 | 0.80 | 0.74 | 0.37 | 0.90 | 0.85 (secs for 1000 call 2048 | 0.51 | 1.85 | 1.79 | 0.78 | 1.90 | 1.84 (secs for 1000 call 4096 | 0.35 | 1.99 | 1.86 | 0.84 | 2.18 | 2.02 (secs for 500 calls 8192 | 1.56 | 5.35 | 5.07 | 2.82 | 5.45 | 5.19 (secs for 500 calls . Multi-dimensional Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100x100 | 0.51 | 0.67 | 0.59 | 0.55 | 0.72 | 0.64 (secs for 100 c ) 1000x100 | 0.70 | 0.67 | 0.62 | 0.70 | 0.70 | 0.64 (secs for 7 cal 256x256 | 0.70 | 0.64 | 0.61 | 0.75 | 0.67 | 0.63 (secs for 10 ca 512x512 | 0.93 | 0.84 | 0.79 | 0.94 | 0.86 | 0.80 (secs for 3 cal . Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.35 | 0.43 | 0.47 (secs for 7000 calls) 1000 | 0.27 | 0.32 | 0.34 (secs for 2000 calls) 256 | 0.60 | 0.78 | 0.77 (secs for 10000 calls) 512 | 0.83 | 1.04 | 1.01 (secs for 10000 calls) 1024 | 0.13 | 0.22 | 0.20 (secs for 1000 calls) 2048 | 0.28 | 0.39 | 0.41 (secs for 1000 calls) 4096 | 0.25 | 0.47 | 0.43 (secs for 500 calls) 8192 | 0.81 | 1.15 | 1.20 (secs for 500 calls) . Inverse Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.36 | 0.92 | 0.87 (secs for 7000 calls) 1000 | 0.41 | 1.12 | 1.15 (secs for 2000 calls) 256 | 0.62 | 1.61 | 2.49 (secs for 10000 calls) 512 | 0.89 | 2.28 | 2.57 (secs for 10000 calls) 1024 | 0.17 | 0.52 | 0.24 (secs for 1000 calls) 2048 | 0.27 | 1.08 | 1.02 (secs for 1000 calls) 4096 | 0.44 | 1.10 | 1.08 (secs for 500 calls) 8192 | 1.10 | 2.22 | 2.18 (secs for 500 calls) ......................... Shifting periodic functions ============================== size | optimized | naive ------------------------------ 100 | 0.10 | 0.53 (secs for 1500 calls) 1000 | 0.07 | 0.76 (secs for 300 calls) 256 | 0.12 | 0.94 (secs for 1500 calls) 512 | 0.13 | 1.11 (secs for 1000 calls) 1024 | 0.11 | 1.08 (secs for 500 calls) 2048 | 0.08 | 1.06 (secs for 200 calls) 4096 | 0.12 | 1.16 (secs for 100 calls) 8192 | 0.17 | 1.45 (secs for 50 calls) . Differentiation of periodic functions ===================================== size | convolve | naive ------------------------------------- 100 | 0.09 | 0.78 (secs for 1500 calls) 1000 | 0.07 | 1.04 (secs for 300 calls) 256 | 0.13 | 1.44 (secs for 1500 calls) 512 | 0.13 | 1.74 (secs for 1000 calls) 1024 | 0.11 | 2.37 (secs for 500 calls) 2048 | 0.09 | 1.73 (secs for 200 calls) 4096 | 0.14 | 1.74 (secs for 100 calls) 8192 | 0.19 | 2.06 (secs for 50 calls) . Tilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.10 | 0.63 (secs for 1500 calls) 1000 | 0.07 | 0.93 (secs for 300 calls) 256 | 0.13 | 1.02 (secs for 1500 calls) 512 | 0.12 | 1.25 (secs for 1000 calls) 1024 | 0.10 | 1.62 (secs for 500 calls) 2048 | 0.08 | 1.12 (secs for 200 calls) 4096 | 0.12 | 1.25 (secs for 100 calls) 8192 | 0.17 | 1.36 (secs for 50 calls) . Hilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.09 | 0.51 (secs for 1500 calls) 1000 | 0.07 | 0.49 (secs for 300 calls) 256 | 0.12 | 0.76 (secs for 1500 calls) 512 | 0.11 | 1.27 (secs for 1000 calls) 1024 | 0.10 | 1.15 (secs for 500 calls) 2048 | 0.08 | 0.88 (secs for 200 calls) 4096 | 0.13 | 0.90 (secs for 100 calls) 8192 | 0.16 | 1.09 (secs for 50 calls) . ---------------------------------------------------------------------- Ran 52 tests in 268.065s OK
-----Original Message----- From: scipy-dev-admin@scipy.net [mailto:scipy-dev-admin@scipy.net] On Behalf Of Pearu Peterson Sent: Saturday, September 28, 2002 3:54 PM To: scipy-dev@scipy.org Subject: [SciPy-dev] fftpack2 for review
Dear SciPy devels,
Please find a scipy.fftpack replacement candidate fftpack2 from
http://cens.ioc.ee/~pearu/misc/fftpack2-0.1.tar.gz
for review. fftpack2 main features include:
- Optional support for high-performance FFT libraries such as FFTW and DJBFFT libraries. By default, Fortran FFTPACK is used and other libraries are automatically detected by the scipy_distutils/system_info.py script.
- Support for both real and complex DFT's and their inverse transforms. When used with FFTW and/or DJBFFT libraries then fftpack2 routines are considerably faster than Numeric.FFT and the current scipy.fftpack. Some comparison timings are included at the end of this message.
- Implementation of differentiation and integration of periodic sequences, various pseudo-differential operators such as Tilbert transform, Hilbert transform, hyperbolic pseudo-differential operators, their inverses, etc.
- Because different FFT libraries use different data storage convention, fftpack2 implements interface to these conventions so that users and developers need not to learn all these conventions, just one is enough. Namely, fftpack2 uses the same data storage specification as Fortran FFTPACK (also used by Numeric.FFT and Matlab, though with somewhat different formal definitions).
- Generic convolution routines that allow quick and easy way of introducing new pseudo-differential operators with the same performance as core pseudo-differential operators.
- All functions have complete documentation strings.
- A rather complete testing site.
Note that fftpack2 requires the latest F2PY that you can get from
http://cens.ioc.ee/projects/f2py2e/2.x/F2PY-2-latest.tar.gz
For testing, fftpack2 uses scipy_test from SciPy CVS repository.
For building/testing instructions and for various notes, including open questions, see NOTES.txt file after unpacking fftpack2-0.1.tar.gz.
As usual, comments, suggestions, and criticism are most welcome.
Regards, Pearu
------------------------------- fftpack2 timing results: ------------------------------- Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.38 | 0.39 | 0.38 | 0.40 | 0.40 | 0.38 (secs for7000calls) 1000 | 0.31 | 0.63 | 0.60 | 0.54 | 0.64 | 0.61 (secs for 2000 256 | 0.60 | 0.75 | 0.62 | 0.71 | 0.78 | 0.66 (secs for 10000 512 | 0.77 | 1.30 | 1.00 | 0.95 | 1.30 | 1.01 (secs for 10000 1024 | 0.14 | 0.26 | 0.18 | 0.17 | 0.25 | 0.19 (secs for 1000 2048 | 0.21 | 0.51 | 0.36 | 0.34 | 0.65 | 0.54 (secs for 1000 4096 | 0.28 | 1.08 | 0.66 | 0.56 | 0.93 | 0.79 (secs for 500 8192 | 1.09 | 3.88 | 3.94 | 2.19 | 4.36 | 3.48 (secs for 500
Inverse Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.17 | 0.26 | 0.63 | 0.45 | 0.87 | 0.85 (secs for 7000 1000 | 0.34 | 1.24 | 1.25 | 0.72 | 1.29 | 1.29 (secs for 2000 256 | 0.62 | 1.57 | 1.44 | 0.69 | 1.62 | 1.47 (secs for 10000 512 | 0.80 | 2.69 | 2.30 | 1.00 | 2.62 | 2.20 (secs for 10000 1024 | 0.12 | 0.51 | 0.41 | 0.18 | 0.54 | 0.44 (secs for 1000 2048 | 0.23 | 1.16 | 1.09 | 0.40 | 1.33 | 1.27 (secs for 1000 4096 | 0.35 | 1.75 | 1.44 | 0.63 | 1.72 | 1.54 (secs for 500 8192 | 1.34 | 5.71 | 6.02 | 2.40 | 6.14 | 5.67 (secs for 500
Multi-dimensional Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100x100 | 0.43 | 0.87 | 0.81 | 0.48 | 0.88 | 0.82 (secs for 100 1000x100 | 0.54 | 0.78 | 0.72 | 0.56 | 0.73 | 0.75 (secs for 7 256x256 | 0.69 | 0.72 | 0.62 | 0.67 | 0.72 | 0.64 (secs for 10 512x512 | 0.81 | 0.91 | 0.84 | 0.81 | 0.87 | 0.81 (secs for 3
Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.36 | 0.47 | 0.47 (secs for 7000 calls) 1000 | 0.29 | 0.38 | 0.38 (secs for 2000 calls) 256 | 0.55 | 0.76 | 0.81 (secs for 10000 calls) 512 | 0.67 | 1.06 | 1.07 (secs for 10000 calls) 1024 | 0.10 | 0.18 | 0.17 (secs for 1000 calls) 2048 | 0.16 | 0.30 | 0.30 (secs for 1000 calls) 4096 | 0.20 | 0.37 | 0.32 (secs for 500 calls) 8192 | 0.75 | 1.29 | 1.26 (secs for 500 calls)
Inverse Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.39 | 0.90 | 0.92 (secs for 7000 calls) 1000 | 0.32 | 0.63 | 0.70 (secs for 2000 calls) 256 | 0.59 | 1.37 | 1.39 (secs for 10000 calls) 512 | 0.73 | 1.70 | 1.80 (secs for 10000 calls) 1024 | 0.11 | 0.27 | 0.27 (secs for 1000 calls) 2048 | 0.20 | 0.45 | 0.44 (secs for 1000 calls) 4096 | 0.24 | 0.65 | 0.68 (secs for 500 calls) 8192 | 0.73 | 1.98 | 2.06 (secs for 500 calls)
Shifting periodic functions ============================== size | optimized | naive ------------------------------ 100 | 0.09 | 0.52 (secs for 1500 calls) 1000 | 0.06 | 0.66 (secs for 300 calls) 256 | 0.10 | 0.84 (secs for 1500 calls) 512 | 0.10 | 1.02 (secs for 1000 calls) 1024 | 0.07 | 1.07 (secs for 500 calls) 2048 | 0.06 | 0.80 (secs for 200 calls) 4096 | 0.00 | 0.86 (secs for 100 calls) 8192 | 0.13 | 1.23 (secs for 50 calls)
Differentiation of periodic functions ===================================== size | convolve | naive ------------------------------------- 100 | 0.09 | 0.59 (secs for 1500 calls) 1000 | 0.06 | 0.75 (secs for 300 calls) 256 | 0.11 | 0.97 (secs for 1500 calls) 512 | 0.09 | 1.21 (secs for 1000 calls) 1024 | 0.08 | 1.19 (secs for 500 calls) 2048 | 0.06 | 1.11 (secs for 200 calls) 4096 | 0.09 | 1.19 (secs for 100 calls) 8192 | 0.14 | 1.48 (secs for 50 calls)
Tilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.09 | 0.54 (secs for 1500 calls) 1000 | 0.06 | 0.60 (secs for 300 calls) 256 | 0.10 | 0.77 (secs for 1500 calls) 512 | 0.09 | 0.89 (secs for 1000 calls) 1024 | 0.07 | 0.90 (secs for 500 calls) 2048 | 0.05 | 0.83 (secs for 200 calls) 4096 | 0.07 | 0.92 (secs for 100 calls) 8192 | 0.11 | 1.04 (secs for 50 calls) . Hilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.09 | 0.47 (secs for 1500 calls) 1000 | 0.06 | 0.48 (secs for 300 calls) 256 | 0.10 | 0.68 (secs for 1500 calls) 512 | 0.09 | 0.77 (secs for 1000 calls) 1024 | 0.07 | 0.74 (secs for 500 calls) 2048 | 0.06 | 0.72 (secs for 200 calls) 4096 | 0.07 | 0.79 (secs for 100 calls) 8192 | 0.11 | 0.93 (secs for 50 calls)
_______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.net http://www.scipy.net/mailman/listinfo/scipy-dev
On Tue, 1 Oct 2002, eric jones wrote:
Hey Pearu,
Overall, first tests went fairly smoothly. Even without djbfft and fftw, things appeared to install correctly on win32. *very* nice job on system checking and handling the lack of libraries.
All tests passed. Sure am glad to finally have some unit tests for fft stuff!
I do not have fftw installed. I tested once without djbfft and once with it. I didn't notice any material difference in timings (the non-djbfft set given below). My tests are on a 850 MHz PIII and they are slower than yours. This leads me to believe, even when trying to use djbfft, I'm not getting it linked correctly. Setup said it detected my djbfft libraries though (I defined the environment variable pointing to my libs to get this to work).
Note that measuring time is currently reliable only on Linux platform where real jiffies [*] are used. Under windows and other unices either real jiffies are not available or I don't know how to use them, and then time.time is used instead. [*] jiffies are used to measure how much the given process uses CPU time, unit is 1/100th of a second. Under Linux one can get user jiffies and system jiffies, in scipy bench tests the last ones are ignored when measuring performance of algorithms.
Also freqshift puts the Nyquist component at the end of the result while fftshift puts it at the beginning of the result. This is for the consistency among varios functions in the fftpack module.
Let me add here that it would be good for performance not to fix the data storage order and provide helper functions such as dftfreq function that return this information, also provide functions for manipulating fft results. For example, djbfft uses a rather strange ordering (that may be the key for its good performance) and to transform the results from djbfft routines to "standard" convention (and back), two nontrivial copying of arrays is needed. And this may be the performance killer when using the djbfft library. If no one will object this, I'll put it in to my todo list. But it is fine with me to stick to one and only one convention if users really want that.
And now some installation notes about what I did to get djbfft running:
Install notes:
1. win32 (cygwin and mingw) required #include "errno.h" in any djbfft file that used errno.
So, we need to patch djbfft. Probably D.J.B. does not like this.
2. I installed the latest f2py but still got the warning error:
WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!
fftpack2 requires F2PY version 2.23.190-1367 or higher but got 2.23.190-1354
WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!
I'm hitting the hay right now, so I'll check out why this happened tomorrow.
That is strange. I would expect that Python will crash (with segfault) when running convolve tests if older F2PY version than 2.23.190-1367 was used for building. Any change that you have several f2py2e packages lying around in your system?
3. How hard would it be to get setup.py to build the djbfft library if one wasn't detected. The makefile builds a lot of the files on the fly using stuff like:
4c0.c: \ idea/c0.c pentium/c0.c ppro/c0.c sparc/c0.c auto_opt sed 1s/PRE/pre4.c/ `cat auto_opt`/c0.c > 4c0.c
4c0.o: \ compile 4c0.c pre4.c fftc4.h complex4.h real4.h fftr4.h real4.h 4i.c ./compile 4c0.c
I'm not sure what all is going on here, but it seems like we could either do the same in python or auto-generate a generic version of this stuff and include it in SciPy. Of course, since djbfft isn't required, this isn't necessary.
In principle, this could be done. However, D.J.B. certainly has made it more difficult for us to accomplish that because djbfft uses a rather unique building process, and also the conditions of altering this software are rather restrictive. I think that having setup.py to build djbfft library does not have a very high priority right now, especially if the performance boost is not so great on all platforms and arhitectures. So, I would put it into the todo list and deal with it in future. Pearu
On Tue, 1 Oct 2002, eric jones wrote:
Hey Pearu,
Overall, first tests went fairly smoothly. Even without djbfft and fftw, things appeared to install correctly on win32. *very* nice job on system checking and handling the lack of libraries.
All tests passed. Sure am glad to finally have some unit tests for fft stuff!
I do not have fftw installed. I tested once without djbfft and once with it. I didn't notice any material difference in timings (the non-djbfft set given below). My tests are on a 850 MHz PIII and they are slower than yours. This leads me to believe, even when trying to use djbfft, I'm not getting it linked correctly. Setup said it detected my djbfft libraries though (I defined the environment variable pointing to my libs to get this to work).
Note that measuring time is currently reliable only on Linux platform where real jiffies [*] are used. Under windows and other unices either real jiffies are not available or I don't know how to use them, and then time.time is used instead.
Ok.
[*] jiffies are used to measure how much the given process uses CPU
unit is 1/100th of a second. Under Linux one can get user jiffies and system jiffies, in scipy bench tests the last ones are ignored when measuring performance of algorithms.
Also freqshift puts the Nyquist component at the end of the result while fftshift puts it at the beginning of the result. This is for the consistency among varios functions in the fftpack module.
Let me add here that it would be good for performance not to fix the data storage order and provide helper functions such as dftfreq function
time, that
return this information, also provide functions for manipulating fft results. For example, djbfft uses a rather strange ordering (that may be the key for its good performance) and to transform the results from djbfft routines to "standard" convention (and back), two nontrivial copying of arrays is needed. And this may be the performance killer when using the djbfft library. If no one will object this, I'll put it in to my todo list. But it is fine with me to stick to one and only one convention if users really want that.
Hmmm. I think it is fine to stick with the current behavior for now. Having multiple output formats is fairly tricky to deal with. Well... I guess it worth asking now. Does anyone need absolute speed with fft so much that they are willing to call some obscure function name and then deal with the ordering themselves? Speak now or forever hold your peace. (Or add it on your own later...)
And now some installation notes about what I did to get djbfft
running:
Install notes:
1. win32 (cygwin and mingw) required #include "errno.h" in any djbfft
file
that used errno.
So, we need to patch djbfft. Probably D.J.B. does not like this.
Yes. I haven't ever gotten a response from him, so I think the chances are slim.
2. I installed the latest f2py but still got the warning error:
WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!
fftpack2 requires F2PY version 2.23.190-1367 or higher but got 2.23.190-1354
WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!WARNING!
I'm hitting the hay right now, so I'll check out why this happened tomorrow.
That is strange. I would expect that Python will crash (with segfault) when running convolve tests if older F2PY version than 2.23.190-1367 was used for building. Any change that you have several f2py2e packages lying around in your system?
I don't think so, but I'll look into it again tonight.
3. How hard would it be to get setup.py to build the djbfft library if
wasn't detected. The makefile builds a lot of the files on the fly using stuff like:
4c0.c: \ idea/c0.c pentium/c0.c ppro/c0.c sparc/c0.c auto_opt sed 1s/PRE/pre4.c/ `cat auto_opt`/c0.c > 4c0.c
4c0.o: \ compile 4c0.c pre4.c fftc4.h complex4.h real4.h fftr4.h real4.h 4i.c ./compile 4c0.c
I'm not sure what all is going on here, but it seems like we could either do the same in python or auto-generate a generic version of
one this
stuff and include it in SciPy. Of course, since djbfft isn't required, this isn't necessary.
In principle, this could be done. However, D.J.B. certainly has made it more difficult for us to accomplish that because djbfft uses a rather unique building process, and also the conditions of altering this software are rather restrictive.
If we do it within our own CVS tree, I think it should be Ok with him, but I'm also not sure it is worth it.
I think that having setup.py to build djbfft library does not have a
very
high priority right now, especially if the performance boost is not so great on all platforms and arhitectures. So, I would put it into the todo list and deal with it in future.
Agreed. Eric
I finally had some time to test fftpack2. It looks really great and I'm happy to support including it. I'll attach my test results at the end (with fftw and djbfft). My only comment right now, is that I don't like doing away with the word fftshift. This is a common command in MATLAB and is quite familiar. I don't necessarily mind freqshift,
Also freqshift puts the Nyquist component at the end of the result while fftshift puts it at the beginning of the result. This is for the consistency among varios functions in the fftpack module.
Is there a difference between freqshift and fftshift? If so, then I am not satisfied. I took great pains to ensure that fftshift works exactly the same way as in MATLAB. It is important that the common fft operations present a single interface to the user. If we want to offer specialized functions that do not mess with the ordering but expect the user to know what they are doing, that is fine, as long as the default functions follow the basic pattern. I think that the MATLAB convention should be followed here. -Travis Results of testing: ................... Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.31 | 0.27 | 0.27 | 0.32 | 0.28 | 0.27 (secs for 7000 calls) 1000 | 0.23 | 0.32 | 0.30 | 0.30 | 0.32 | 0.30 (secs for 2000 calls) 256 | 0.47 | 0.50 | 0.48 | 0.57 | 0.47 | 0.48 (secs for 10000 calls) 512 | 0.61 | 0.82 | 0.78 | 0.79 | 0.82 | 0.76 (secs for 10000 calls) 1024 | 0.09 | 0.16 | 0.12 | 0.13 | 0.15 | 0.12 (secs for 1000 calls) 2048 | 0.16 | 0.27 | 0.25 | 0.26 | 0.27 | 0.31 (secs for 1000 calls) 4096 | 0.15 | 0.43 | 0.39 | 0.30 | 0.44 | 0.41 (secs for 500 calls) 8192 | 0.60 | 2.53 | 2.54 | 1.45 | 2.71 | 2.57 (secs for 500 calls) . Inverse Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100 | 0.31 | 0.68 | 0.70 | 0.34 | 0.71 | 0.69 (secs for 7000 calls) 1000 | 0.25 | 0.76 | 0.74 | 0.39 | 0.78 | 0.76 (secs for 2000 calls) 256 | 0.51 | 1.28 | 1.29 | 0.58 | 1.29 | 1.32 (secs for 10000 calls) 512 | 0.67 | 2.03 | 1.98 | 0.81 | 2.09 | 2.04 (secs for 10000 calls) 1024 | 0.10 | 0.36 | 0.34 | 0.14 | 0.36 | 0.33 (secs for 1000 calls) 2048 | 0.18 | 0.69 | 0.70 | 0.25 | 0.72 | 0.73 (secs for 1000 calls) 4096 | 0.19 | 1.13 | 1.14 | 0.32 | 1.20 | 1.21 (secs for 500 calls) 8192 | 0.80 | 3.82 | 3.80 | 1.56 | 4.05 | 3.88 (secs for 500 calls) . Multi-dimensional Fast Fourier Transform ======================================================== | real input | complex input -------------------------------------------------------- size |fftpack2| FFT | scipy |fftpack2| FFT | scipy -------------------------------------------------------- 100x100 | 0.20 | 0.45 | 0.35 | 0.25 | 0.47 | 0.45 (secs for 100 calls) 1000x100 | 0.27 | 0.47 | 0.38 | 0.28 | 0.43 | 0.43 (secs for 7 calls) 256x256 | 0.61 | 0.46 | 0.49 | 0.63 | 0.47 | 0.45 (secs for 10 calls) 512x512 | 0.77 | 0.60 | 0.65 | 0.80 | 0.66 | 0.62 (secs for 3 calls) . Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.30 | 0.36 | 0.37 (secs for 7000 calls) 1000 | 0.22 | 0.26 | 0.23 (secs for 2000 calls) 256 | 0.44 | 0.58 | 0.62 (secs for 10000 calls) 512 | 0.57 | 0.79 | 0.84 (secs for 10000 calls) 1024 | 0.09 | 0.12 | 0.12 (secs for 1000 calls) 2048 | 0.13 | 0.19 | 0.21 (secs for 1000 calls) 4096 | 0.13 | 0.20 | 0.19 (secs for 500 calls) 8192 | 0.40 | 0.82 | 0.81 (secs for 500 calls) . Inverse Fast Fourier Transform (real data) ================================== size |fftpack2| FFT | scipy ---------------------------------- 100 | 0.29 | 0.80 | 0.77 (secs for 7000 calls) 1000 | 0.23 | 0.46 | 0.45 (secs for 2000 calls) 256 | 0.45 | 1.23 | 1.22 (secs for 10000 calls) 512 | 0.59 | 1.47 | 1.46 (secs for 10000 calls) 1024 | 0.08 | 0.21 | 0.20 (secs for 1000 calls) 2048 | 0.15 | 0.30 | 0.33 (secs for 1000 calls) 4096 | 0.13 | 0.36 | 0.38 (secs for 500 calls) 8192 | 0.40 | 1.28 | 1.50 (secs for 500 calls) . ---------------------------------------------------------------------- Ran 24 tests in 298.302s OK [travis@travis fftpack2-0.1]$ python tests/test_helper.py 10 .... ---------------------------------------------------------------------- Ran 4 tests in 0.013s OK [travis@travis fftpack2-0.1]$ python tests/test_pseudo_diffs.py 10 .................... Shifting periodic functions ============================== size | optimized | naive ------------------------------ 100 | 0.07 | 0.44 (secs for 1500 calls) 1000 | 0.05 | 0.50 (secs for 300 calls) 256 | 0.08 | 0.76 (secs for 1500 calls) 512 | 0.07 | 0.85 (secs for 1000 calls) 1024 | 0.06 | 0.80 (secs for 500 calls) 2048 | 0.05 | 0.74 (secs for 200 calls) 4096 | 0.04 | 0.81 (secs for 100 calls) 8192 | 0.11 | 0.90 (secs for 50 calls) . Differentiation of periodic functions ===================================== size | convolve | naive ------------------------------------- 100 | 0.14 | 0.59 (secs for 1500 calls) 1000 | 0.05 | 0.65 (secs for 300 calls) 256 | 0.09 | 0.92 (secs for 1500 calls) 512 | 0.07 | 1.12 (secs for 1000 calls) 1024 | 0.06 | 1.06 (secs for 500 calls) 2048 | 0.04 | 0.93 (secs for 200 calls) 4096 | 0.05 | 0.98 (secs for 100 calls) 8192 | 0.08 | 1.10 (secs for 50 calls) . Tilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.08 | 0.48 (secs for 1500 calls) 1000 | 0.04 | 0.41 (secs for 300 calls) 256 | 0.09 | 0.69 (secs for 1500 calls) 512 | 0.06 | 0.70 (secs for 1000 calls) 1024 | 0.06 | 0.67 (secs for 500 calls) 2048 | 0.04 | 0.64 (secs for 200 calls) 4096 | 0.04 | 0.66 (secs for 100 calls) 8192 | 0.06 | 0.72 (secs for 50 calls) . Hilbert transform of periodic functions ========================================= size | optimized | naive ----------------------------------------- 100 | 0.07 | 0.42 (secs for 1500 calls) 1000 | 0.04 | 0.36 (secs for 300 calls) 256 | 0.08 | 0.61 (secs for 1500 calls) 512 | 0.07 | 0.62 (secs for 1000 calls) 1024 | 0.05 | 0.58 (secs for 500 calls) 2048 | 0.03 | 0.53 (secs for 200 calls) 4096 | 0.05 | 0.63 (secs for 100 calls) 8192 | 0.06 | 0.72 (secs for 50 calls) . ---------------------------------------------------------------------- Ran 24 tests in 62.302s OK
On 4 Oct 2002, Travis Oliphant wrote:
I finally had some time to test fftpack2. It looks really great and I'm happy to support including it. I'll attach my test results at the end (with fftw and djbfft).
My only comment right now, is that I don't like doing away with the word fftshift. This is a common command in MATLAB and is quite familiar. I don't necessarily mind freqshift,
Also freqshift puts the Nyquist component at the end of the result while fftshift puts it at the beginning of the result. This is for the consistency among various functions in the fftpack module.
Is there a difference between freqshift and fftshift? If so, then I am not satisfied. I took great pains to ensure that fftshift works exactly the same way as in MATLAB.
To be honest, I used the docs of scipy.fftpack.fft as a basic reference and followed that while implementing fftpack2. Until fftshift, everything was consistent. Then it turned out that fftshift uses a different convention than fft. Namely, according to scipy.fftpack.fft docs, the corresponding frequencies of the fft result are given as [ 0, 1, 2, 3, 4, -3, -2, -1] for n=8 case, for example. But fftshift assumed [ 0, 1, 2, 3, -4, -3, -2, -1] Both cases are equivalent except when applying the definition of fftshift, that is, how to order the frequencies of the fft result. Now, in fftpack2 I just kept being consistent (as everything were documented already) and as a result, the results of fftpack2.freqshift and fftshift (in scipy.fftpack and Matlab) had to differ by the location of the Nyquist mode (that's also one of the hidden reasons why I introduced a new name freqshift for fftshift; about the other reason, see below). I guess it would not be a good idea to have both freqshift and fftshift at the same time (the latter behaving exactly as fftshift in MATLAB). Personally, I have no preference to either of them, it's just a matter of choosing a convention and sticking to it. But since you insist MATLAB convention, I'll change the docs in fftpack2: basically, the Nyquist mode will be corresponding to the smallest negative frequency instead of the largest positive one.
It is important that the common fft operations present a single interface to the user. If we want to offer specialized functions that do not mess with the ordering but expect the user to know what they are doing, that is fine, as long as the default functions follow the basic pattern.
Ok.
I think that the MATLAB convention should be followed here.
Ok. But I still think that "fftshift" is a misleading name, even if MATLAB uses it. Often people forget that FFT is just an efficient algorithm for performing the Fourier matrix multiplication with a vector, nothing more. And fftshift has nothing to do with this FFT algorithm. I see the following solutions: 1) continue with a bad manner learned from Matlab and use fftshift 2) use a more appropriate name, I proposed freqshift (other suggestions are welcome), and may be provide fftshift in MLab.py, for example. On the other hand, I don't want that some function name will make somebody unhappy. Right now I am just proposing this name change because it feels right to me and it is the right time to make a choice, later it would be too late. What do you think? Pearu
To be honest, I used the docs of scipy.fftpack.fft as a basic reference and followed that while implementing fftpack2. Until fftshift, everything was consistent. Then it turned out that fftshift uses a different convention than fft. Namely, according to scipy.fftpack.fft docs, the corresponding frequencies of the fft result are given as [ 0, 1, 2, 3, 4, -3, -2, -1] for n=8 case, for example. But fftshift assumed [ 0, 1, 2, 3, -4, -3, -2, -1]
O.K. so this is the issue: How to treat the Nyquist value for even length fft.
Both cases are equivalent except when applying the definition of fftshift,
Yes, this is true, the fft value for bin 4 is the same as for bin -4 so how you interpret it doesn't really matter. Except for the fact that a large group of people have been trained to calculate the FFT using MATLAB and will be used to using fftshift where the convention is to interpret the Nyquist value as the value for the negative bin number (-4). This is a sutble difference but could cause headaches for many people starting to use scipy who think they know how to shift the results of an FFT to get a plot based on their MATLAB experience. Because it really doesn't matter, I'm inclined to go with what MATLAB chose rather than have a new convention.
that is, how to order the frequencies of the fft result. Now, in fftpack2 I just kept being consistent (as everything were documented already) and as a result, the results of fftpack2.freqshift and fftshift (in scipy.fftpack and Matlab) had to differ by the location of the Nyquist mode
I don't understand how you were consistent. What is documented already? Why is interpreting the Nyquist mode as coming from a positive frequency more consistent? (that's also one of the hidden reasons why I
introduced a new name freqshift for fftshift; about the other reason, see below).
I guess it would not be a good idea to have both freqshift and fftshift at the same time (the latter behaving exactly as fftshift in MATLAB).
I suppose that is fine.
Personally, I have no preference to either of them, it's just a matter of choosing a convention and sticking to it. But since you insist MATLAB convention, I'll change the docs in fftpack2: basically, the Nyquist mode will be corresponding to the smallest negative frequency instead of the largest positive one.
That is what I would recommend.
Ok. But I still think that "fftshift" is a misleading name, even if MATLAB uses it. Often people forget that FFT is just an efficient algorithm for performing the Fourier matrix multiplication with a vector, nothing more. And fftshift has nothing to do with this FFT algorithm.
True people do forget the difference between DTFT and FFT, and fftshift doesn't have anything to do with the FFT algorithm. But, the reason for fftshift is because of the way the FFT algorithm computes the frequency leaving the positive frequencies in the first half and the negative frequencies in the last half. We usually don't draw graphs this way (or create ranges this way). So, I would submit that fftshift is still a good name because it is a shift made necessary by the use of the FFT. But, it doesn't really matter and I'm not so opposed to freqshift either. It's a fine name too. We aren't copying all of MATLAB's names for every other function. I just lean towards the name fftshift because I've gotten so used to using it the many, many, many times I've computed Fourier transforms
I see the following solutions: 1) continue with a bad manner learned from Matlab and use fftshift
I don't agree that fftshift is so bad.
2) use a more appropriate name, I proposed freqshift (other suggestions are welcome), and may be provide fftshift in MLab.py, for example.
On the other hand, I don't want that some function name will make somebody unhappy. Right now I am just proposing this name change because it feels right to me and it is the right time to make a choice, later it would be too late. What do you think?
I'm not going to be sad either way. But, I do think we should stick with the convention of interpreting the Nyquist mode as a negative frequency. I'm very happy with the great job you've done in cleaning up and testing the fft package. --Travis O. P.S. I also have a command called fftfreq that generates the frequency values for each index (bin) assumed by the FFT algorithm. I noticed you did something similar with dftfreq. This is another point we need to iron out.
Hi -- I'm not sure this is the right place for this question, please redirect me if there's a more appropriate mailing list: We're currently porting several of our boost-wrapped Numeric classes to numarray. The routines are spread across several cpp files, and **libnumarray_API is needed in each file. In Numeric, we invoke import_array() in the module file, to initialize **PyArray_API, and def NO_IMPORT_ARRAY in the other files; this produces extern void **PyArray_API in those files, and everyone's happy. Currently libnumarray.h doesn't implement this -- is there another way to compile a multiple file project which initializes the API just once? Thanks, Phil Austin
participants (5)
-
eric jones -
Pearu Peterson -
Philip Austin -
Travis Oliphant -
Travis Oliphant