Chris@19: Chris@19: Chris@19: Chris@19: FFTW FAQ - Section 3 Chris@19: Chris@19: Chris@19: Chris@19: Chris@19: Chris@19:

Chris@19: FFTW FAQ - Section 3
Chris@19: Using FFTW Chris@19:

Chris@19: Chris@19:
Chris@19: Chris@19:

Chris@19: Question 3.1. Why not support the FFTW 2 interface in FFTW Chris@19: 3? Chris@19:

Chris@19: Chris@19: FFTW 3 has semantics incompatible with earlier versions: its plans can Chris@19: only be used for a given stride, multiplicity, and other Chris@19: characteristics of the input and output arrays; these stronger Chris@19: semantics are necessary for performance reasons. Thus, it is Chris@19: impossible to efficiently emulate the older interface (whose plans can Chris@19: be used for any transform of the same size). We believe that it Chris@19: should be possible to upgrade most programs without any difficulty, Chris@19: however. Chris@19:

Chris@19: Question 3.2. Why do FFTW 3 plans encapsulate the input/output arrays Chris@19: and not just the algorithm? Chris@19:

Chris@19: Chris@19: There are several reasons: Chris@19: Chris@19: Chris@19:

Chris@19: Question 3.3. FFTW seems really slow. Chris@19:

Chris@19: Chris@19: You are probably recreating the plan before every transform, rather Chris@19: than creating it once and reusing it for all transforms of the same Chris@19: size. FFTW is designed to be used in the following way: Chris@19: Chris@19: Chris@19: If you don't need to compute many transforms and the time for the Chris@19: planner is significant, you have two options. First, you can use the Chris@19: FFTW_ESTIMATE option in the planner, which uses heuristics Chris@19: instead of runtime measurements and produces a good plan in a short Chris@19: time. Second, you can use the wisdom feature to precompute the plan; Chris@19: see Q3.9 `Can I save FFTW's plans?' Chris@19:

Chris@19: Question 3.4. FFTW slows down after repeated Chris@19: calls. Chris@19:

Chris@19: Chris@19: Probably, NaNs or similar are creeping into your data, and the Chris@19: slowdown is due to the resulting floating-point exceptions. For Chris@19: example, be aware that repeatedly FFTing the same array is a diverging Chris@19: process (because FFTW computes the unnormalized transform). Chris@19: Chris@19:

Chris@19: Question 3.5. An FFTW routine is crashing when I call Chris@19: it. Chris@19:

Chris@19: Chris@19: Did the FFTW test programs pass (make check, or cd tests; make bigcheck if you want to be paranoid)? If so, you almost Chris@19: certainly have a bug in your own code. For example, you could be Chris@19: passing invalid arguments (such as wrongly-sized arrays) to FFTW, or Chris@19: you could simply have memory corruption elsewhere in your program that Chris@19: causes random crashes later on. Please don't complain to us unless Chris@19: you can come up with a minimal self-contained program (preferably Chris@19: under 30 lines) that illustrates the problem. Chris@19: Chris@19:

Chris@19: Question 3.6. My Fortran program crashes when calling Chris@19: FFTW. Chris@19:

Chris@19: Chris@19: As described in the manual, on 64-bit machines you must store the Chris@19: plans in variables large enough to hold a pointer, for example Chris@19: integer*8. We recommend using integer*8 on 32-bit machines as well, to simplify porting. Chris@19: Chris@19:

Chris@19: Question 3.7. FFTW gives results different from my old Chris@19: FFT. Chris@19:

Chris@19: Chris@19: People follow many different conventions for the DFT, and you should Chris@19: be sure to know the ones that we use (described in the FFTW manual). Chris@19: In particular, you should be aware that the Chris@19: FFTW_FORWARD/FFTW_BACKWARD directions correspond to signs of -1/+1 in the exponent of the DFT definition. Chris@19: (Numerical Recipes uses the opposite convention.) Chris@19:

Chris@19: You should also know that we compute an unnormalized transform. In Chris@19: contrast, Matlab is an example of program that computes a normalized Chris@19: transform. See Q3.10 `Why does your inverse transform return a scaled Chris@19: result?'. Chris@19:

Chris@19: Finally, note that floating-point arithmetic is not exact, so Chris@19: different FFT algorithms will give slightly different results (on the Chris@19: order of the numerical accuracy; typically a fractional difference of Chris@19: 1e-15 or so in double precision). Chris@19:

Chris@19: Question 3.8. FFTW gives different results between Chris@19: runs Chris@19:

Chris@19: Chris@19: If you use FFTW_MEASURE or FFTW_PATIENT mode, then the algorithm FFTW employs is not deterministic: it depends on Chris@19: runtime performance measurements. This will cause the results to vary Chris@19: slightly from run to run. However, the differences should be slight, Chris@19: on the order of the floating-point precision, and therefore should Chris@19: have no practical impact on most applications. Chris@19: Chris@19:

Chris@19: If you use saved plans (wisdom) or FFTW_ESTIMATE mode, however, then the algorithm is deterministic and the results should be Chris@19: identical between runs. Chris@19:

Chris@19: Question 3.9. Can I save FFTW's plans? Chris@19:

Chris@19: Chris@19: Yes. Starting with version 1.2, FFTW provides the Chris@19: wisdom mechanism for saving plans; see the FFTW manual. Chris@19: Chris@19:

Chris@19: Question 3.10. Why does your inverse transform return a scaled Chris@19: result? Chris@19:

Chris@19: Chris@19: Computing the forward transform followed by the backward transform (or Chris@19: vice versa) yields the original array scaled by the size of the array. Chris@19: (For multi-dimensional transforms, the size of the array is the Chris@19: product of the dimensions.) We could, instead, have chosen a Chris@19: normalization that would have returned the unscaled array. Or, to Chris@19: accomodate the many conventions in this matter, the transform routines Chris@19: could have accepted a "scale factor" parameter. We did not Chris@19: do this, however, for two reasons. First, we didn't want to sacrifice Chris@19: performance in the common case where the scale factor is 1. Second, in Chris@19: real applications the FFT is followed or preceded by some computation Chris@19: on the data, into which the scale factor can typically be absorbed at Chris@19: little or no cost. Chris@19:

Chris@19: Question 3.11. How can I make FFTW put the origin (zero frequency) at Chris@19: the center of its output? Chris@19:

Chris@19: Chris@19: For human viewing of a spectrum, it is often convenient to put the Chris@19: origin in frequency space at the center of the output array, rather Chris@19: than in the zero-th element (the default in FFTW). If all of the Chris@19: dimensions of your array are even, you can accomplish this by simply Chris@19: multiplying each element of the input array by (-1)^(i + j + ...), Chris@19: where i, j, etcetera are the indices of the element. (This trick is a Chris@19: general property of the DFT, and is not specific to FFTW.) Chris@19: Chris@19:

Chris@19: Question 3.12. How do I FFT an image/audio file in Chris@19: foobar format? Chris@19:

Chris@19: Chris@19: FFTW performs an FFT on an array of floating-point values. You can Chris@19: certainly use it to compute the transform of an image or audio stream, Chris@19: but you are responsible for figuring out your data format and Chris@19: converting it to the form FFTW requires. Chris@19: Chris@19:

Chris@19: Question 3.13. My program does not link (on Chris@19: Unix). Chris@19:

Chris@19: Chris@19: The libraries must be listed in the correct order Chris@19: (-lfftw3 -lm for FFTW 3.x) and after your program sources/objects. (The general rule is that if A uses B, then A must be listed before B in the link command.). Chris@19:

Chris@19: Question 3.14. I included your header, but linking still Chris@19: fails. Chris@19:

Chris@19: Chris@19: You're a C++ programmer, aren't you? You have to compile the FFTW Chris@19: library and link it into your program, not just Chris@19: #include <fftw3.h>. (Yes, this is really a FAQ.) Chris@19:

Chris@19: Question 3.15. My program crashes, complaining about stack Chris@19: space. Chris@19:

Chris@19: Chris@19: You cannot declare large arrays with automatic storage (e.g. via Chris@19: fftw_complex array[N]); you should use fftw_malloc (or equivalent) to allocate the arrays you want Chris@19: to transform if they are larger than a few hundred elements. Chris@19: Chris@19:

Chris@19: Question 3.16. FFTW seems to have a memory Chris@19: leak. Chris@19:

Chris@19: Chris@19: After you create a plan, FFTW caches the information required to Chris@19: quickly recreate the plan. (See Q3.9 `Can I save FFTW's plans?') It also maintains a small amount of other persistent memory. You can deallocate all of Chris@19: FFTW's internally allocated memory, if you wish, by calling Chris@19: fftw_cleanup(), as documented in the manual. Chris@19:

Chris@19: Question 3.17. The output of FFTW's transform is all Chris@19: zeros. Chris@19:

Chris@19: Chris@19: You should initialize your input array after creating the plan, unless you use FFTW_ESTIMATE: planning with FFTW_MEASURE or FFTW_PATIENT overwrites the input/output arrays, as described in the manual. Chris@19: Chris@19:

Chris@19: Question 3.18. How do I call FFTW from the Microsoft language du Chris@19: jour? Chris@19:

Chris@19: Chris@19: Please do not ask us Windows-specific questions. We do not Chris@19: use Windows. We know nothing about Visual Basic, Visual C++, or .NET. Chris@19: Please find the appropriate Usenet discussion group and ask your Chris@19: question there. See also Q2.2 `Does FFTW run on Windows?'. Chris@19:

Chris@19: Question 3.19. Can I compute only a subset of the DFT Chris@19: outputs? Chris@19:

Chris@19: Chris@19: In general, no, an FFT intrinsically computes all outputs from all Chris@19: inputs. In principle, there is something called a Chris@19: pruned FFT that can do what you want, but to compute K outputs out of N the Chris@19: complexity is in general O(N log K) instead of O(N log N), thus saving Chris@19: only a small additive factor in the log. (The same argument holds if Chris@19: you instead have only K nonzero inputs.) Chris@19: Chris@19:

Chris@19: There are some specific cases in which you can get the O(N log K) Chris@19: performance benefits easily, however, by combining a few ordinary Chris@19: FFTs. In particular, the case where you want the first K outputs, Chris@19: where K divides N, can be handled by performing N/K transforms of size Chris@19: K and then summing the outputs multiplied by appropriate phase Chris@19: factors. For more details, see pruned FFTs with FFTW. Chris@19:

Chris@19: There are also some algorithms that compute pruned transforms Chris@19: approximately, but they are beyond the scope of this FAQ. Chris@19: Chris@19:

Chris@19: Question 3.20. Can I use FFTW's routines for in-place and Chris@19: out-of-place matrix transposition? Chris@19:

Chris@19: Chris@19: You can use the FFTW guru interface to create a rank-0 transform of Chris@19: vector rank 2 where the vector strides are transposed. (A rank-0 Chris@19: transform is equivalent to a 1D transform of size 1, which. just Chris@19: copies the input into the output.) Specifying the same location for Chris@19: the input and output makes the transpose in-place. Chris@19: Chris@19:

Chris@19: For double-valued data stored in row-major format, plan creation looks Chris@19: like this:

Chris@19: fftw_plan plan_transpose(int rows, int cols, double *in, double *out)
Chris@19: {
Chris@19:     const unsigned flags = FFTW_ESTIMATE; /* other flags are possible */
Chris@19:     fftw_iodim howmany_dims[2];
Chris@19: 
Chris@19:     howmany_dims[0].n  = rows;
Chris@19:     howmany_dims[0].is = cols;
Chris@19:     howmany_dims[0].os = 1;
Chris@19: 
Chris@19:     howmany_dims[1].n  = cols;
Chris@19:     howmany_dims[1].is = 1;
Chris@19:     howmany_dims[1].os = rows;
Chris@19: 
Chris@19:     return fftw_plan_guru_r2r(/*rank=*/ 0, /*dims=*/ NULL,
Chris@19:                               /*howmany_rank=*/ 2, howmany_dims,
Chris@19:                               in, out, /*kind=*/ NULL, flags);
Chris@19: }
Chris@19: 
Chris@19: (This entry was written by Rhys Ulerich.) Chris@19:
Chris@19: Next: Internals of FFTW.
Chris@19: Back: Installing FFTW.
Chris@19: Return to contents.

Chris@19:

Chris@19: Matteo Frigo and Steven G. Johnson / fftw@fftw.org Chris@19: - 04 March 2014 Chris@19:

Chris@19: Extracted from FFTW Frequently Asked Questions with Answers, Chris@19: Copyright © 2014 Matteo Frigo and Massachusetts Institute of Technology. Chris@19: