Mercurial > hg > js-dsp-test

annotate fft/fftw/fftw-3.3.4/doc/FAQ/fftw-faq.html/section3.html @ 40:223f770b5341 kissfft-double tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Try a double-precision kissfft
author Chris Cannam
date Wed, 07 Sep 2016 10:40:32 +0100
parents 26056e866c29
children
rev   line source
Chris@19 1 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
Chris@19 2 <html>
Chris@19 3 <head><title>
Chris@19 4 FFTW FAQ - Section 3
Chris@19 5 </title>
Chris@19 6 <link rev="made" href="mailto:fftw@fftw.org">
Chris@19 7 <link rel="Contents" href="index.html">
Chris@19 8 <link rel="Start" href="index.html">
Chris@19 9 <link rel="Next" href="section4.html"><link rel="Previous" href="section2.html"><link rel="Bookmark" title="FFTW FAQ" href="index.html">
Chris@19 10 </head><body text="#000000" bgcolor="#FFFFFF"><h1>
Chris@19 11 FFTW FAQ - Section 3 <br>
Chris@19 12 Using FFTW
Chris@19 13 </h1>
Chris@19 14
Chris@19 15 <ul>
Chris@19 16 <li><a href="#fftw2to3" rel=subdocument>Q3.1. Why not support the FFTW 2 interface in FFTW
Chris@19 17 3?</a>
Chris@19 18 <li><a href="#planperarray" rel=subdocument>Q3.2. Why do FFTW 3 plans encapsulate the input/output arrays and not just
Chris@19 19 the algorithm?</a>
Chris@19 20 <li><a href="#slow" rel=subdocument>Q3.3. FFTW seems really slow.</a>
Chris@19 21 <li><a href="#slows" rel=subdocument>Q3.4. FFTW slows down after repeated calls.</a>
Chris@19 22 <li><a href="#segfault" rel=subdocument>Q3.5. An FFTW routine is crashing when I call it.</a>
Chris@19 23 <li><a href="#fortran64" rel=subdocument>Q3.6. My Fortran program crashes when calling FFTW.</a>
Chris@19 24 <li><a href="#conventions" rel=subdocument>Q3.7. FFTW gives results different from my old
Chris@19 25 FFT.</a>
Chris@19 26 <li><a href="#nondeterministic" rel=subdocument>Q3.8. FFTW gives different results between runs</a>
Chris@19 27 <li><a href="#savePlans" rel=subdocument>Q3.9. Can I save FFTW's plans?</a>
Chris@19 28 <li><a href="#whyscaled" rel=subdocument>Q3.10. Why does your inverse transform return a scaled
Chris@19 29 result?</a>
Chris@19 30 <li><a href="#centerorigin" rel=subdocument>Q3.11. How can I make FFTW put the origin (zero frequency) at the center of
Chris@19 31 its output?</a>
Chris@19 32 <li><a href="#imageaudio" rel=subdocument>Q3.12. How do I FFT an image/audio file in <i>foobar</i> format?</a>
Chris@19 33 <li><a href="#linkfails" rel=subdocument>Q3.13. My program does not link (on Unix).</a>
Chris@19 34 <li><a href="#linkheader" rel=subdocument>Q3.14. I included your header, but linking still
Chris@19 35 fails.</a>
Chris@19 36 <li><a href="#nostack" rel=subdocument>Q3.15. My program crashes, complaining about stack
Chris@19 37 space.</a>
Chris@19 38 <li><a href="#leaks" rel=subdocument>Q3.16. FFTW seems to have a memory leak.</a>
Chris@19 39 <li><a href="#allzero" rel=subdocument>Q3.17. The output of FFTW's transform is all zeros.</a>
Chris@19 40 <li><a href="#vbetalia" rel=subdocument>Q3.18. How do I call FFTW from the Microsoft language du
Chris@19 41 jour?</a>
Chris@19 42 <li><a href="#pruned" rel=subdocument>Q3.19. Can I compute only a subset of the DFT outputs?</a>
Chris@19 43 <li><a href="#transpose" rel=subdocument>Q3.20. Can I use FFTW's routines for in-place and out-of-place matrix
Chris@19 44 transposition?</a>
Chris@19 45 </ul><hr>
Chris@19 46
Chris@19 47 <h2><A name="fftw2to3">
Chris@19 48 Question 3.1. Why not support the FFTW 2 interface in FFTW
Chris@19 49 3?
Chris@19 50 </A></h2>
Chris@19 51
Chris@19 52 FFTW 3 has semantics incompatible with earlier versions: its plans can
Chris@19 53 only be used for a given stride, multiplicity, and other
Chris@19 54 characteristics of the input and output arrays; these stronger
Chris@19 55 semantics are necessary for performance reasons. Thus, it is
Chris@19 56 impossible to efficiently emulate the older interface (whose plans can
Chris@19 57 be used for any transform of the same size). We believe that it
Chris@19 58 should be possible to upgrade most programs without any difficulty,
Chris@19 59 however.
Chris@19 60 <h2><A name="planperarray">
Chris@19 61 Question 3.2. Why do FFTW 3 plans encapsulate the input/output arrays
Chris@19 62 and not just the algorithm?
Chris@19 63 </A></h2>
Chris@19 64
Chris@19 65 There are several reasons:
Chris@19 66 <ul>
Chris@19 67 <li>It was important for performance reasons that the plan be specific to
Chris@19 68 array characteristics like the stride (and alignment, for SIMD), and
Chris@19 69 requiring that the user maintain these invariants is error prone.
Chris@19 70
Chris@19 71 <li>In most high-performance applications, as far as we can tell, you are
Chris@19 72 usually transforming the same array over and over, so FFTW's semantics
Chris@19 73 should not be a burden.
Chris@19 74 <li>If you need to transform another array of the same size, creating a
Chris@19 75 new plan once the first exists is a cheap operation.
Chris@19 76
Chris@19 77 <li>If you need to transform many arrays of the same size at once, you
Chris@19 78 should really use the <code>plan_many</code> routines in FFTW's &quot;advanced&quot;
Chris@19 79 interface.
Chris@19 80 <li>If the abovementioned array characteristics are the same, you are
Chris@19 81 willing to pay close attention to the documentation, and you really
Chris@19 82 need to, we provide a &quot;new-array execution&quot; interface to
Chris@19 83 apply a plan to a new array.
Chris@19 84 </ul>
Chris@19 85
Chris@19 86 <h2><A name="slow">
Chris@19 87 Question 3.3. FFTW seems really slow.
Chris@19 88 </A></h2>
Chris@19 89
Chris@19 90 You are probably recreating the plan before every transform, rather
Chris@19 91 than creating it once and reusing it for all transforms of the same
Chris@19 92 size. FFTW is designed to be used in the following way:
Chris@19 93
Chris@19 94 <ul>
Chris@19 95 <li>First, you create a plan. This will take several seconds.
Chris@19 96
Chris@19 97 <li>Then, you reuse the plan many times to perform FFTs. These are fast.
Chris@19 98
Chris@19 99 </ul>
Chris@19 100 If you don't need to compute many transforms and the time for the
Chris@19 101 planner is significant, you have two options. First, you can use the
Chris@19 102 <code>FFTW_ESTIMATE</code> option in the planner, which uses heuristics
Chris@19 103 instead of runtime measurements and produces a good plan in a short
Chris@19 104 time. Second, you can use the wisdom feature to precompute the plan;
Chris@19 105 see <A href="#savePlans">Q3.9 `Can I save FFTW's plans?'</A>
Chris@19 106 <h2><A name="slows">
Chris@19 107 Question 3.4. FFTW slows down after repeated
Chris@19 108 calls.
Chris@19 109 </A></h2>
Chris@19 110
Chris@19 111 Probably, NaNs or similar are creeping into your data, and the
Chris@19 112 slowdown is due to the resulting floating-point exceptions. For
Chris@19 113 example, be aware that repeatedly FFTing the same array is a diverging
Chris@19 114 process (because FFTW computes the unnormalized transform).
Chris@19 115
Chris@19 116 <h2><A name="segfault">
Chris@19 117 Question 3.5. An FFTW routine is crashing when I call
Chris@19 118 it.
Chris@19 119 </A></h2>
Chris@19 120
Chris@19 121 Did the FFTW test programs pass (<code>make check</code>, or <code>cd tests; make bigcheck</code> if you want to be paranoid)? If so, you almost
Chris@19 122 certainly have a bug in your own code. For example, you could be
Chris@19 123 passing invalid arguments (such as wrongly-sized arrays) to FFTW, or
Chris@19 124 you could simply have memory corruption elsewhere in your program that
Chris@19 125 causes random crashes later on. Please don't complain to us unless
Chris@19 126 you can come up with a minimal self-contained program (preferably
Chris@19 127 under 30 lines) that illustrates the problem.
Chris@19 128
Chris@19 129 <h2><A name="fortran64">
Chris@19 130 Question 3.6. My Fortran program crashes when calling
Chris@19 131 FFTW.
Chris@19 132 </A></h2>
Chris@19 133
Chris@19 134 As described in the manual, on 64-bit machines you must store the
Chris@19 135 plans in variables large enough to hold a pointer, for example
Chris@19 136 <code>integer*8</code>. We recommend using <code>integer*8</code> on 32-bit machines as well, to simplify porting.
Chris@19 137
Chris@19 138 <h2><A name="conventions">
Chris@19 139 Question 3.7. FFTW gives results different from my old
Chris@19 140 FFT.
Chris@19 141 </A></h2>
Chris@19 142
Chris@19 143 People follow many different conventions for the DFT, and you should
Chris@19 144 be sure to know the ones that we use (described in the FFTW manual).
Chris@19 145 In particular, you should be aware that the
Chris@19 146 <code>FFTW_FORWARD</code>/<code>FFTW_BACKWARD</code> directions correspond to signs of -1/+1 in the exponent of the DFT definition.
Chris@19 147 (<i>Numerical Recipes</i> uses the opposite convention.)
Chris@19 148 <p>
Chris@19 149 You should also know that we compute an unnormalized transform. In
Chris@19 150 contrast, Matlab is an example of program that computes a normalized
Chris@19 151 transform. See <A href="#whyscaled">Q3.10 `Why does your inverse transform return a scaled
Chris@19 152 result?'</A>.
Chris@19 153 <p>
Chris@19 154 Finally, note that floating-point arithmetic is not exact, so
Chris@19 155 different FFT algorithms will give slightly different results (on the
Chris@19 156 order of the numerical accuracy; typically a fractional difference of
Chris@19 157 1e-15 or so in double precision).
Chris@19 158 <h2><A name="nondeterministic">
Chris@19 159 Question 3.8. FFTW gives different results between
Chris@19 160 runs
Chris@19 161 </A></h2>
Chris@19 162
Chris@19 163 If you use <code>FFTW_MEASURE</code> or <code>FFTW_PATIENT</code> mode, then the algorithm FFTW employs is not deterministic: it depends on
Chris@19 164 runtime performance measurements. This will cause the results to vary
Chris@19 165 slightly from run to run. However, the differences should be slight,
Chris@19 166 on the order of the floating-point precision, and therefore should
Chris@19 167 have no practical impact on most applications.
Chris@19 168
Chris@19 169 <p>
Chris@19 170 If you use saved plans (wisdom) or <code>FFTW_ESTIMATE</code> mode, however, then the algorithm is deterministic and the results should be
Chris@19 171 identical between runs.
Chris@19 172 <h2><A name="savePlans">
Chris@19 173 Question 3.9. Can I save FFTW's plans?
Chris@19 174 </A></h2>
Chris@19 175
Chris@19 176 Yes. Starting with version 1.2, FFTW provides the
Chris@19 177 <code>wisdom</code> mechanism for saving plans; see the FFTW manual.
Chris@19 178
Chris@19 179 <h2><A name="whyscaled">
Chris@19 180 Question 3.10. Why does your inverse transform return a scaled
Chris@19 181 result?
Chris@19 182 </A></h2>
Chris@19 183
Chris@19 184 Computing the forward transform followed by the backward transform (or
Chris@19 185 vice versa) yields the original array scaled by the size of the array.
Chris@19 186 (For multi-dimensional transforms, the size of the array is the
Chris@19 187 product of the dimensions.) We could, instead, have chosen a
Chris@19 188 normalization that would have returned the unscaled array. Or, to
Chris@19 189 accomodate the many conventions in this matter, the transform routines
Chris@19 190 could have accepted a &quot;scale factor&quot; parameter. We did not
Chris@19 191 do this, however, for two reasons. First, we didn't want to sacrifice
Chris@19 192 performance in the common case where the scale factor is 1. Second, in
Chris@19 193 real applications the FFT is followed or preceded by some computation
Chris@19 194 on the data, into which the scale factor can typically be absorbed at
Chris@19 195 little or no cost.
Chris@19 196 <h2><A name="centerorigin">
Chris@19 197 Question 3.11. How can I make FFTW put the origin (zero frequency) at
Chris@19 198 the center of its output?
Chris@19 199 </A></h2>
Chris@19 200
Chris@19 201 For human viewing of a spectrum, it is often convenient to put the
Chris@19 202 origin in frequency space at the center of the output array, rather
Chris@19 203 than in the zero-th element (the default in FFTW). If all of the
Chris@19 204 dimensions of your array are even, you can accomplish this by simply
Chris@19 205 multiplying each element of the input array by (-1)^(i + j + ...),
Chris@19 206 where i, j, etcetera are the indices of the element. (This trick is a
Chris@19 207 general property of the DFT, and is not specific to FFTW.)
Chris@19 208
Chris@19 209 <h2><A name="imageaudio">
Chris@19 210 Question 3.12. How do I FFT an image/audio file in
Chris@19 211 <i>foobar</i> format?
Chris@19 212 </A></h2>
Chris@19 213
Chris@19 214 FFTW performs an FFT on an array of floating-point values. You can
Chris@19 215 certainly use it to compute the transform of an image or audio stream,
Chris@19 216 but you are responsible for figuring out your data format and
Chris@19 217 converting it to the form FFTW requires.
Chris@19 218
Chris@19 219 <h2><A name="linkfails">
Chris@19 220 Question 3.13. My program does not link (on
Chris@19 221 Unix).
Chris@19 222 </A></h2>
Chris@19 223
Chris@19 224 The libraries must be listed in the correct order
Chris@19 225 (<code>-lfftw3 -lm</code> for FFTW 3.x) and <i>after</i> your program sources/objects. (The general rule is that if <i>A</i> uses <i>B</i>, then <i>A</i> must be listed before <i>B</i> in the link command.).
Chris@19 226 <h2><A name="linkheader">
Chris@19 227 Question 3.14. I included your header, but linking still
Chris@19 228 fails.
Chris@19 229 </A></h2>
Chris@19 230
Chris@19 231 You're a C++ programmer, aren't you? You have to compile the FFTW
Chris@19 232 library and link it into your program, not just
Chris@19 233 <code>#include &lt;fftw3.h&gt;</code>. (Yes, this is really a FAQ.)
Chris@19 234 <h2><A name="nostack">
Chris@19 235 Question 3.15. My program crashes, complaining about stack
Chris@19 236 space.
Chris@19 237 </A></h2>
Chris@19 238
Chris@19 239 You cannot declare large arrays with automatic storage (e.g. via
Chris@19 240 <code>fftw_complex array[N]</code>); you should use <code>fftw_malloc</code> (or equivalent) to allocate the arrays you want
Chris@19 241 to transform if they are larger than a few hundred elements.
Chris@19 242
Chris@19 243 <h2><A name="leaks">
Chris@19 244 Question 3.16. FFTW seems to have a memory
Chris@19 245 leak.
Chris@19 246 </A></h2>
Chris@19 247
Chris@19 248 After you create a plan, FFTW caches the information required to
Chris@19 249 quickly recreate the plan. (See <A href="#savePlans">Q3.9 `Can I save FFTW's plans?'</A>) It also maintains a small amount of other persistent memory. You can deallocate all of
Chris@19 250 FFTW's internally allocated memory, if you wish, by calling
Chris@19 251 <code>fftw_cleanup()</code>, as documented in the manual.
Chris@19 252 <h2><A name="allzero">
Chris@19 253 Question 3.17. The output of FFTW's transform is all
Chris@19 254 zeros.
Chris@19 255 </A></h2>
Chris@19 256
Chris@19 257 You should initialize your input array <i>after</i> creating the plan, unless you use <code>FFTW_ESTIMATE</code>: planning with <code>FFTW_MEASURE</code> or <code>FFTW_PATIENT</code> overwrites the input/output arrays, as described in the manual.
Chris@19 258
Chris@19 259 <h2><A name="vbetalia">
Chris@19 260 Question 3.18. How do I call FFTW from the Microsoft language du
Chris@19 261 jour?
Chris@19 262 </A></h2>
Chris@19 263
Chris@19 264 Please <i>do not</i> ask us Windows-specific questions. We do not
Chris@19 265 use Windows. We know nothing about Visual Basic, Visual C++, or .NET.
Chris@19 266 Please find the appropriate Usenet discussion group and ask your
Chris@19 267 question there. See also <A href="section2.html#runOnWindows">Q2.2 `Does FFTW run on Windows?'</A>.
Chris@19 268 <h2><A name="pruned">
Chris@19 269 Question 3.19. Can I compute only a subset of the DFT
Chris@19 270 outputs?
Chris@19 271 </A></h2>
Chris@19 272
Chris@19 273 In general, no, an FFT intrinsically computes all outputs from all
Chris@19 274 inputs. In principle, there is something called a
Chris@19 275 <i>pruned FFT</i> that can do what you want, but to compute K outputs out of N the
Chris@19 276 complexity is in general O(N log K) instead of O(N log N), thus saving
Chris@19 277 only a small additive factor in the log. (The same argument holds if
Chris@19 278 you instead have only K nonzero inputs.)
Chris@19 279
Chris@19 280 <p>
Chris@19 281 There are some specific cases in which you can get the O(N log K)
Chris@19 282 performance benefits easily, however, by combining a few ordinary
Chris@19 283 FFTs. In particular, the case where you want the first K outputs,
Chris@19 284 where K divides N, can be handled by performing N/K transforms of size
Chris@19 285 K and then summing the outputs multiplied by appropriate phase
Chris@19 286 factors. For more details, see <A href="http://www.fftw.org/pruned.html">pruned FFTs with FFTW</A>.
Chris@19 287 <p>
Chris@19 288 There are also some algorithms that compute pruned transforms
Chris@19 289 <i>approximately</i>, but they are beyond the scope of this FAQ.
Chris@19 290
Chris@19 291 <h2><A name="transpose">
Chris@19 292 Question 3.20. Can I use FFTW's routines for in-place and
Chris@19 293 out-of-place matrix transposition?
Chris@19 294 </A></h2>
Chris@19 295
Chris@19 296 You can use the FFTW guru interface to create a rank-0 transform of
Chris@19 297 vector rank 2 where the vector strides are transposed. (A rank-0
Chris@19 298 transform is equivalent to a 1D transform of size 1, which. just
Chris@19 299 copies the input into the output.) Specifying the same location for
Chris@19 300 the input and output makes the transpose in-place.
Chris@19 301
Chris@19 302 <p>
Chris@19 303 For double-valued data stored in row-major format, plan creation looks
Chris@19 304 like this: <pre>
Chris@19 305 fftw_plan plan_transpose(int rows, int cols, double *in, double *out)
Chris@19 306 {
Chris@19 307 const unsigned flags = FFTW_ESTIMATE; /* other flags are possible */
Chris@19 308 fftw_iodim howmany_dims[2];
Chris@19 309
Chris@19 310 howmany_dims[0].n = rows;
Chris@19 311 howmany_dims[0].is = cols;
Chris@19 312 howmany_dims[0].os = 1;
Chris@19 313
Chris@19 314 howmany_dims[1].n = cols;
Chris@19 315 howmany_dims[1].is = 1;
Chris@19 316 howmany_dims[1].os = rows;
Chris@19 317
Chris@19 318 return fftw_plan_guru_r2r(/*rank=*/ 0, /*dims=*/ NULL,
Chris@19 319 /*howmany_rank=*/ 2, howmany_dims,
Chris@19 320 in, out, /*kind=*/ NULL, flags);
Chris@19 321 }
Chris@19 322 </pre>
Chris@19 323 (This entry was written by Rhys Ulerich.)
Chris@19 324 <hr>
Chris@19 325 Next: <a href="section4.html" rel=precedes>Internals of FFTW</a>.<br>
Chris@19 326 Back: <a href="section2.html" rev=precedes>Installing FFTW</a>.<br>
Chris@19 327 <a href="index.html" rev=subdocument>Return to contents</a>.<p>
Chris@19 328 <address>
Chris@19 329 <A href="http://www.fftw.org">Matteo Frigo and Steven G. Johnson</A> / <A href="mailto:fftw@fftw.org">fftw@fftw.org</A>
Chris@19 330 - 04 March 2014
Chris@19 331 </address><br>
Chris@19 332 Extracted from FFTW Frequently Asked Questions with Answers,
Chris@19 333 Copyright &copy; 2014 Matteo Frigo and Massachusetts Institute of Technology.
Chris@19 334 </body></html>