annotate src/fftw-3.3.3/doc/reference.texi @ 23:619f715526df sv_v2.1

Update Vamp plugin SDK to 2.5
author Chris Cannam
date Thu, 09 May 2013 10:52:46 +0100
parents 37bf6b4a2645
children
rev   line source
Chris@10 1 @node FFTW Reference, Multi-threaded FFTW, Other Important Topics, Top
Chris@10 2 @chapter FFTW Reference
Chris@10 3
Chris@10 4 This chapter provides a complete reference for all sequential (i.e.,
Chris@10 5 one-processor) FFTW functions. Parallel transforms are described in
Chris@10 6 later chapters.
Chris@10 7
Chris@10 8 @menu
Chris@10 9 * Data Types and Files::
Chris@10 10 * Using Plans::
Chris@10 11 * Basic Interface::
Chris@10 12 * Advanced Interface::
Chris@10 13 * Guru Interface::
Chris@10 14 * New-array Execute Functions::
Chris@10 15 * Wisdom::
Chris@10 16 * What FFTW Really Computes::
Chris@10 17 @end menu
Chris@10 18
Chris@10 19 @c ------------------------------------------------------------
Chris@10 20 @node Data Types and Files, Using Plans, FFTW Reference, FFTW Reference
Chris@10 21 @section Data Types and Files
Chris@10 22
Chris@10 23 All programs using FFTW should include its header file:
Chris@10 24
Chris@10 25 @example
Chris@10 26 #include <fftw3.h>
Chris@10 27 @end example
Chris@10 28
Chris@10 29 You must also link to the FFTW library. On Unix, this
Chris@10 30 means adding @code{-lfftw3 -lm} at the @emph{end} of the link command.
Chris@10 31
Chris@10 32 @menu
Chris@10 33 * Complex numbers::
Chris@10 34 * Precision::
Chris@10 35 * Memory Allocation::
Chris@10 36 @end menu
Chris@10 37
Chris@10 38 @c =========>
Chris@10 39 @node Complex numbers, Precision, Data Types and Files, Data Types and Files
Chris@10 40 @subsection Complex numbers
Chris@10 41
Chris@10 42 The default FFTW interface uses @code{double} precision for all
Chris@10 43 floating-point numbers, and defines a @code{fftw_complex} type to hold
Chris@10 44 complex numbers as:
Chris@10 45
Chris@10 46 @example
Chris@10 47 typedef double fftw_complex[2];
Chris@10 48 @end example
Chris@10 49 @tindex fftw_complex
Chris@10 50
Chris@10 51 Here, the @code{[0]} element holds the real part and the @code{[1]}
Chris@10 52 element holds the imaginary part.
Chris@10 53
Chris@10 54 Alternatively, if you have a C compiler (such as @code{gcc}) that
Chris@10 55 supports the C99 revision of the ANSI C standard, you can use C's new
Chris@10 56 native complex type (which is binary-compatible with the typedef above).
Chris@10 57 In particular, if you @code{#include <complex.h>} @emph{before}
Chris@10 58 @code{<fftw3.h>}, then @code{fftw_complex} is defined to be the native
Chris@10 59 complex type and you can manipulate it with ordinary arithmetic
Chris@10 60 (e.g. @code{x = y * (3+4*I)}, where @code{x} and @code{y} are
Chris@10 61 @code{fftw_complex} and @code{I} is the standard symbol for the
Chris@10 62 imaginary unit);
Chris@10 63 @cindex C99
Chris@10 64
Chris@10 65
Chris@10 66 C++ has its own @code{complex<T>} template class, defined in the
Chris@10 67 standard @code{<complex>} header file. Reportedly, the C++ standards
Chris@10 68 committee has recently agreed to mandate that the storage format used
Chris@10 69 for this type be binary-compatible with the C99 type, i.e. an array
Chris@10 70 @code{T[2]} with consecutive real @code{[0]} and imaginary @code{[1]}
Chris@10 71 parts. (See report
Chris@10 72 @uref{http://www.open-std.org/jtc1/sc22/WG21/docs/papers/2002/n1388.pdf
Chris@10 73 WG21/N1388}.) Although not part of the official standard as of this
Chris@10 74 writing, the proposal stated that: ``This solution has been tested with
Chris@10 75 all current major implementations of the standard library and shown to
Chris@10 76 be working.'' To the extent that this is true, if you have a variable
Chris@10 77 @code{complex<double> *x}, you can pass it directly to FFTW via
Chris@10 78 @code{reinterpret_cast<fftw_complex*>(x)}.
Chris@10 79 @cindex C++
Chris@10 80 @cindex portability
Chris@10 81
Chris@10 82 @c =========>
Chris@10 83 @node Precision, Memory Allocation, Complex numbers, Data Types and Files
Chris@10 84 @subsection Precision
Chris@10 85 @cindex precision
Chris@10 86
Chris@10 87 You can install single and long-double precision versions of FFTW,
Chris@10 88 which replace @code{double} with @code{float} and @code{long double},
Chris@10 89 respectively (@pxref{Installation and Customization}). To use these
Chris@10 90 interfaces, you:
Chris@10 91
Chris@10 92 @itemize @bullet
Chris@10 93
Chris@10 94 @item
Chris@10 95 Link to the single/long-double libraries; on Unix, @code{-lfftw3f} or
Chris@10 96 @code{-lfftw3l} instead of (or in addition to) @code{-lfftw3}. (You
Chris@10 97 can link to the different-precision libraries simultaneously.)
Chris@10 98
Chris@10 99 @item
Chris@10 100 Include the @emph{same} @code{<fftw3.h>} header file.
Chris@10 101
Chris@10 102 @item
Chris@10 103 Replace all lowercase instances of @samp{fftw_} with @samp{fftwf_} or
Chris@10 104 @samp{fftwl_} for single or long-double precision, respectively.
Chris@10 105 (@code{fftw_complex} becomes @code{fftwf_complex}, @code{fftw_execute}
Chris@10 106 becomes @code{fftwf_execute}, etcetera.)
Chris@10 107
Chris@10 108 @item
Chris@10 109 Uppercase names, i.e. names beginning with @samp{FFTW_}, remain the
Chris@10 110 same.
Chris@10 111
Chris@10 112 @item
Chris@10 113 Replace @code{double} with @code{float} or @code{long double} for
Chris@10 114 subroutine parameters.
Chris@10 115
Chris@10 116 @end itemize
Chris@10 117
Chris@10 118 Depending upon your compiler and/or hardware, @code{long double} may not
Chris@10 119 be any more precise than @code{double} (or may not be supported at all,
Chris@10 120 although it is standard in C99).
Chris@10 121 @cindex C99
Chris@10 122
Chris@10 123
Chris@10 124 We also support using the nonstandard @code{__float128}
Chris@10 125 quadruple-precision type provided by recent versions of @code{gcc} on
Chris@10 126 32- and 64-bit x86 hardware (@pxref{Installation and Customization}).
Chris@10 127 To use this type, link with @code{-lfftw3q -lquadmath -lm} (the
Chris@10 128 @code{libquadmath} library provided by @code{gcc} is needed for
Chris@10 129 quadruple-precision trigonometric functions) and use @samp{fftwq_}
Chris@10 130 identifiers.
Chris@10 131
Chris@10 132 @c =========>
Chris@10 133 @node Memory Allocation, , Precision, Data Types and Files
Chris@10 134 @subsection Memory Allocation
Chris@10 135
Chris@10 136 @example
Chris@10 137 void *fftw_malloc(size_t n);
Chris@10 138 void fftw_free(void *p);
Chris@10 139 @end example
Chris@10 140 @findex fftw_malloc
Chris@10 141 @findex fftw_free
Chris@10 142
Chris@10 143 These are functions that behave identically to @code{malloc} and
Chris@10 144 @code{free}, except that they guarantee that the returned pointer obeys
Chris@10 145 any special alignment restrictions imposed by any algorithm in FFTW
Chris@10 146 (e.g. for SIMD acceleration). @xref{SIMD alignment and fftw_malloc}.
Chris@10 147 @cindex alignment
Chris@10 148
Chris@10 149
Chris@10 150 Data allocated by @code{fftw_malloc} @emph{must} be deallocated by
Chris@10 151 @code{fftw_free} and not by the ordinary @code{free}.
Chris@10 152
Chris@10 153 These routines simply call through to your operating system's
Chris@10 154 @code{malloc} or, if necessary, its aligned equivalent
Chris@10 155 (e.g. @code{memalign}), so you normally need not worry about any
Chris@10 156 significant time or space overhead. You are @emph{not required} to use
Chris@10 157 them to allocate your data, but we strongly recommend it.
Chris@10 158
Chris@10 159 Note: in C++, just as with ordinary @code{malloc}, you must typecast
Chris@10 160 the output of @code{fftw_malloc} to whatever pointer type you are
Chris@10 161 allocating.
Chris@10 162 @cindex C++
Chris@10 163
Chris@10 164
Chris@10 165 We also provide the following two convenience functions to allocate
Chris@10 166 real and complex arrays with @code{n} elements, which are equivalent
Chris@10 167 to @code{(double *) fftw_malloc(sizeof(double) * n)} and
Chris@10 168 @code{(fftw_complex *) fftw_malloc(sizeof(fftw_complex) * n)},
Chris@10 169 respectively:
Chris@10 170
Chris@10 171 @example
Chris@10 172 double *fftw_alloc_real(size_t n);
Chris@10 173 fftw_complex *fftw_alloc_complex(size_t n);
Chris@10 174 @end example
Chris@10 175 @findex fftw_alloc_real
Chris@10 176 @findex fftw_alloc_complex
Chris@10 177
Chris@10 178 The equivalent functions in other precisions allocate arrays of @code{n}
Chris@10 179 elements in that precision. e.g. @code{fftwf_alloc_real(n)} is
Chris@10 180 equivalent to @code{(float *) fftwf_malloc(sizeof(float) * n)}.
Chris@10 181 @cindex precision
Chris@10 182
Chris@10 183 @c ------------------------------------------------------------
Chris@10 184 @node Using Plans, Basic Interface, Data Types and Files, FFTW Reference
Chris@10 185 @section Using Plans
Chris@10 186
Chris@10 187 Plans for all transform types in FFTW are stored as type
Chris@10 188 @code{fftw_plan} (an opaque pointer type), and are created by one of the
Chris@10 189 various planning routines described in the following sections.
Chris@10 190 @tindex fftw_plan
Chris@10 191 An @code{fftw_plan} contains all information necessary to compute the
Chris@10 192 transform, including the pointers to the input and output arrays.
Chris@10 193
Chris@10 194 @example
Chris@10 195 void fftw_execute(const fftw_plan plan);
Chris@10 196 @end example
Chris@10 197 @findex fftw_execute
Chris@10 198
Chris@10 199 This executes the @code{plan}, to compute the corresponding transform on
Chris@10 200 the arrays for which it was planned (which must still exist). The plan
Chris@10 201 is not modified, and @code{fftw_execute} can be called as many times as
Chris@10 202 desired.
Chris@10 203
Chris@10 204 To apply a given plan to a different array, you can use the new-array execute
Chris@10 205 interface. @xref{New-array Execute Functions}.
Chris@10 206
Chris@10 207 @code{fftw_execute} (and equivalents) is the only function in FFTW
Chris@10 208 guaranteed to be thread-safe; see @ref{Thread safety}.
Chris@10 209
Chris@10 210 This function:
Chris@10 211 @example
Chris@10 212 void fftw_destroy_plan(fftw_plan plan);
Chris@10 213 @end example
Chris@10 214 @findex fftw_destroy_plan
Chris@10 215 deallocates the @code{plan} and all its associated data.
Chris@10 216
Chris@10 217 FFTW's planner saves some other persistent data, such as the
Chris@10 218 accumulated wisdom and a list of algorithms available in the current
Chris@10 219 configuration. If you want to deallocate all of that and reset FFTW
Chris@10 220 to the pristine state it was in when you started your program, you can
Chris@10 221 call:
Chris@10 222
Chris@10 223 @example
Chris@10 224 void fftw_cleanup(void);
Chris@10 225 @end example
Chris@10 226 @findex fftw_cleanup
Chris@10 227
Chris@10 228 After calling @code{fftw_cleanup}, all existing plans become undefined,
Chris@10 229 and you should not attempt to execute them nor to destroy them. You can
Chris@10 230 however create and execute/destroy new plans, in which case FFTW starts
Chris@10 231 accumulating wisdom information again.
Chris@10 232
Chris@10 233 @code{fftw_cleanup} does not deallocate your plans, however. To prevent
Chris@10 234 memory leaks, you must still call @code{fftw_destroy_plan} before
Chris@10 235 executing @code{fftw_cleanup}.
Chris@10 236
Chris@10 237 Occasionally, it may useful to know FFTW's internal ``cost'' metric
Chris@10 238 that it uses to compare plans to one another; this cost is
Chris@10 239 proportional to an execution time of the plan, in undocumented units,
Chris@10 240 if the plan was created with the @code{FFTW_MEASURE} or other
Chris@10 241 timing-based options, or alternatively is a heuristic cost function
Chris@10 242 for @code{FFTW_ESTIMATE} plans. (The cost values of measured and
Chris@10 243 estimated plans are not comparable, being in different units. Also,
Chris@10 244 costs from different FFTW versions or the same version compiled
Chris@10 245 differently may not be in the same units. Plans created from wisdom
Chris@10 246 have a cost of 0 since no timing measurement is performed for them.
Chris@10 247 Finally, certain problems for which only one top-level algorithm was
Chris@10 248 possible may have required no measurements of the cost of the whole
Chris@10 249 plan, in which case @code{fftw_cost} will also return 0.) The cost
Chris@10 250 metric for a given plan is returned by:
Chris@10 251
Chris@10 252 @example
Chris@10 253 double fftw_cost(const fftw_plan plan);
Chris@10 254 @end example
Chris@10 255 @findex fftw_cost
Chris@10 256
Chris@10 257 The following two routines are provided purely for academic purposes
Chris@10 258 (that is, for entertainment).
Chris@10 259
Chris@10 260 @example
Chris@10 261 void fftw_flops(const fftw_plan plan,
Chris@10 262 double *add, double *mul, double *fma);
Chris@10 263 @end example
Chris@10 264 @findex fftw_flops
Chris@10 265
Chris@10 266 Given a @code{plan}, set @code{add}, @code{mul}, and @code{fma} to an
Chris@10 267 exact count of the number of floating-point additions, multiplications,
Chris@10 268 and fused multiply-add operations involved in the plan's execution. The
Chris@10 269 total number of floating-point operations (flops) is @code{add + mul +
Chris@10 270 2*fma}, or @code{add + mul + fma} if the hardware supports fused
Chris@10 271 multiply-add instructions (although the number of FMA operations is only
Chris@10 272 approximate because of compiler voodoo). (The number of operations
Chris@10 273 should be an integer, but we use @code{double} to avoid overflowing
Chris@10 274 @code{int} for large transforms; the arguments are of type @code{double}
Chris@10 275 even for single and long-double precision versions of FFTW.)
Chris@10 276
Chris@10 277 @example
Chris@10 278 void fftw_fprint_plan(const fftw_plan plan, FILE *output_file);
Chris@10 279 void fftw_print_plan(const fftw_plan plan);
Chris@10 280 @end example
Chris@10 281 @findex fftw_fprint_plan
Chris@10 282 @findex fftw_print_plan
Chris@10 283
Chris@10 284 This outputs a ``nerd-readable'' representation of the @code{plan} to
Chris@10 285 the given file or to @code{stdout}, respectively.
Chris@10 286
Chris@10 287 @c ------------------------------------------------------------
Chris@10 288 @node Basic Interface, Advanced Interface, Using Plans, FFTW Reference
Chris@10 289 @section Basic Interface
Chris@10 290 @cindex basic interface
Chris@10 291
Chris@10 292 Recall that the FFTW API is divided into three parts@footnote{@i{Gallia est
Chris@10 293 omnis divisa in partes tres} (Julius Caesar).}: the @dfn{basic interface}
Chris@10 294 computes a single transform of contiguous data, the @dfn{advanced
Chris@10 295 interface} computes transforms of multiple or strided arrays, and the
Chris@10 296 @dfn{guru interface} supports the most general data layouts,
Chris@10 297 multiplicities, and strides. This section describes the the basic
Chris@10 298 interface, which we expect to satisfy the needs of most users.
Chris@10 299
Chris@10 300 @menu
Chris@10 301 * Complex DFTs::
Chris@10 302 * Planner Flags::
Chris@10 303 * Real-data DFTs::
Chris@10 304 * Real-data DFT Array Format::
Chris@10 305 * Real-to-Real Transforms::
Chris@10 306 * Real-to-Real Transform Kinds::
Chris@10 307 @end menu
Chris@10 308
Chris@10 309 @c =========>
Chris@10 310 @node Complex DFTs, Planner Flags, Basic Interface, Basic Interface
Chris@10 311 @subsection Complex DFTs
Chris@10 312
Chris@10 313 @example
Chris@10 314 fftw_plan fftw_plan_dft_1d(int n0,
Chris@10 315 fftw_complex *in, fftw_complex *out,
Chris@10 316 int sign, unsigned flags);
Chris@10 317 fftw_plan fftw_plan_dft_2d(int n0, int n1,
Chris@10 318 fftw_complex *in, fftw_complex *out,
Chris@10 319 int sign, unsigned flags);
Chris@10 320 fftw_plan fftw_plan_dft_3d(int n0, int n1, int n2,
Chris@10 321 fftw_complex *in, fftw_complex *out,
Chris@10 322 int sign, unsigned flags);
Chris@10 323 fftw_plan fftw_plan_dft(int rank, const int *n,
Chris@10 324 fftw_complex *in, fftw_complex *out,
Chris@10 325 int sign, unsigned flags);
Chris@10 326 @end example
Chris@10 327 @findex fftw_plan_dft_1d
Chris@10 328 @findex fftw_plan_dft_2d
Chris@10 329 @findex fftw_plan_dft_3d
Chris@10 330 @findex fftw_plan_dft
Chris@10 331
Chris@10 332 Plan a complex input/output discrete Fourier transform (DFT) in zero or
Chris@10 333 more dimensions, returning an @code{fftw_plan} (@pxref{Using Plans}).
Chris@10 334
Chris@10 335 Once you have created a plan for a certain transform type and
Chris@10 336 parameters, then creating another plan of the same type and parameters,
Chris@10 337 but for different arrays, is fast and shares constant data with the
Chris@10 338 first plan (if it still exists).
Chris@10 339
Chris@10 340 The planner returns @code{NULL} if the plan cannot be created. In the
Chris@10 341 standard FFTW distribution, the basic interface is guaranteed to return
Chris@10 342 a non-@code{NULL} plan. A plan may be @code{NULL}, however, if you are
Chris@10 343 using a customized FFTW configuration supporting a restricted set of
Chris@10 344 transforms.
Chris@10 345
Chris@10 346 @subsubheading Arguments
Chris@10 347 @itemize @bullet
Chris@10 348
Chris@10 349 @item
Chris@10 350 @code{rank} is the rank of the transform (it should be the size of the
Chris@10 351 array @code{*n}), and can be any non-negative integer. (@xref{Complex
Chris@10 352 Multi-Dimensional DFTs}, for the definition of ``rank''.) The
Chris@10 353 @samp{_1d}, @samp{_2d}, and @samp{_3d} planners correspond to a
Chris@10 354 @code{rank} of @code{1}, @code{2}, and @code{3}, respectively. The rank
Chris@10 355 may be zero, which is equivalent to a rank-1 transform of size 1, i.e. a
Chris@10 356 copy of one number from input to output.
Chris@10 357
Chris@10 358 @item
Chris@10 359 @code{n0}, @code{n1}, @code{n2}, or @code{n[0..rank-1]} (as appropriate
Chris@10 360 for each routine) specify the size of the transform dimensions. They
Chris@10 361 can be any positive integer.
Chris@10 362
Chris@10 363 @itemize @minus
Chris@10 364 @item
Chris@10 365 @cindex row-major
Chris@10 366 Multi-dimensional arrays are stored in row-major order with dimensions:
Chris@10 367 @code{n0} x @code{n1}; or @code{n0} x @code{n1} x @code{n2}; or
Chris@10 368 @code{n[0]} x @code{n[1]} x ... x @code{n[rank-1]}.
Chris@10 369 @xref{Multi-dimensional Array Format}.
Chris@10 370 @item
Chris@10 371 FFTW is best at handling sizes of the form
Chris@10 372 @ifinfo
Chris@10 373 @math{2^a 3^b 5^c 7^d 11^e 13^f},
Chris@10 374 @end ifinfo
Chris@10 375 @tex
Chris@10 376 $2^a 3^b 5^c 7^d 11^e 13^f$,
Chris@10 377 @end tex
Chris@10 378 @html
Chris@10 379 2<sup>a</sup> 3<sup>b</sup> 5<sup>c</sup> 7<sup>d</sup>
Chris@10 380 11<sup>e</sup> 13<sup>f</sup>,
Chris@10 381 @end html
Chris@10 382 where @math{e+f} is either @math{0} or @math{1}, and the other exponents
Chris@10 383 are arbitrary. Other sizes are computed by means of a slow,
Chris@10 384 general-purpose algorithm (which nevertheless retains @Onlogn{} performance even for prime sizes). It is possible to customize FFTW
Chris@10 385 for different array sizes; see @ref{Installation and Customization}.
Chris@10 386 Transforms whose sizes are powers of @math{2} are especially fast.
Chris@10 387 @end itemize
Chris@10 388
Chris@10 389 @item
Chris@10 390 @code{in} and @code{out} point to the input and output arrays of the
Chris@10 391 transform, which may be the same (yielding an in-place transform).
Chris@10 392 @cindex in-place
Chris@10 393 These arrays are overwritten during planning, unless
Chris@10 394 @code{FFTW_ESTIMATE} is used in the flags. (The arrays need not be
Chris@10 395 initialized, but they must be allocated.)
Chris@10 396
Chris@10 397 If @code{in == out}, the transform is @dfn{in-place} and the input
Chris@10 398 array is overwritten. If @code{in != out}, the two arrays must
Chris@10 399 not overlap (but FFTW does not check for this condition).
Chris@10 400
Chris@10 401 @item
Chris@10 402 @ctindex FFTW_FORWARD
Chris@10 403 @ctindex FFTW_BACKWARD
Chris@10 404 @code{sign} is the sign of the exponent in the formula that defines the
Chris@10 405 Fourier transform. It can be @math{-1} (= @code{FFTW_FORWARD}) or
Chris@10 406 @math{+1} (= @code{FFTW_BACKWARD}).
Chris@10 407
Chris@10 408 @item
Chris@10 409 @cindex flags
Chris@10 410 @code{flags} is a bitwise OR (@samp{|}) of zero or more planner flags,
Chris@10 411 as defined in @ref{Planner Flags}.
Chris@10 412
Chris@10 413 @end itemize
Chris@10 414
Chris@10 415 FFTW computes an unnormalized transform: computing a forward followed by
Chris@10 416 a backward transform (or vice versa) will result in the original data
Chris@10 417 multiplied by the size of the transform (the product of the dimensions).
Chris@10 418 @cindex normalization
Chris@10 419 For more information, see @ref{What FFTW Really Computes}.
Chris@10 420
Chris@10 421 @c =========>
Chris@10 422 @node Planner Flags, Real-data DFTs, Complex DFTs, Basic Interface
Chris@10 423 @subsection Planner Flags
Chris@10 424
Chris@10 425 All of the planner routines in FFTW accept an integer @code{flags}
Chris@10 426 argument, which is a bitwise OR (@samp{|}) of zero or more of the flag
Chris@10 427 constants defined below. These flags control the rigor (and time) of
Chris@10 428 the planning process, and can also impose (or lift) restrictions on the
Chris@10 429 type of transform algorithm that is employed.
Chris@10 430
Chris@10 431 @emph{Important:} the planner overwrites the input array during
Chris@10 432 planning unless a saved plan (@pxref{Wisdom}) is available for that
Chris@10 433 problem, so you should initialize your input data after creating the
Chris@10 434 plan. The only exceptions to this are the @code{FFTW_ESTIMATE} and
Chris@10 435 @code{FFTW_WISDOM_ONLY} flags, as mentioned below.
Chris@10 436
Chris@10 437 In all cases, if wisdom is available for the given problem that was
Chris@10 438 created with equal-or-greater planning rigor, then the more rigorous
Chris@10 439 wisdom is used. For example, in @code{FFTW_ESTIMATE} mode any available
Chris@10 440 wisdom is used, whereas in @code{FFTW_PATIENT} mode only wisdom created
Chris@10 441 in patient or exhaustive mode can be used. @xref{Words of Wisdom-Saving
Chris@10 442 Plans}.
Chris@10 443
Chris@10 444 @subsubheading Planning-rigor flags
Chris@10 445 @itemize @bullet
Chris@10 446
Chris@10 447 @item
Chris@10 448 @ctindex FFTW_ESTIMATE
Chris@10 449 @code{FFTW_ESTIMATE} specifies that, instead of actual measurements of
Chris@10 450 different algorithms, a simple heuristic is used to pick a (probably
Chris@10 451 sub-optimal) plan quickly. With this flag, the input/output arrays are
Chris@10 452 not overwritten during planning.
Chris@10 453
Chris@10 454 @item
Chris@10 455 @ctindex FFTW_MEASURE
Chris@10 456 @code{FFTW_MEASURE} tells FFTW to find an optimized plan by actually
Chris@10 457 @emph{computing} several FFTs and measuring their execution time.
Chris@10 458 Depending on your machine, this can take some time (often a few
Chris@10 459 seconds). @code{FFTW_MEASURE} is the default planning option.
Chris@10 460
Chris@10 461 @item
Chris@10 462 @ctindex FFTW_PATIENT
Chris@10 463 @code{FFTW_PATIENT} is like @code{FFTW_MEASURE}, but considers a wider
Chris@10 464 range of algorithms and often produces a ``more optimal'' plan
Chris@10 465 (especially for large transforms), but at the expense of several times
Chris@10 466 longer planning time (especially for large transforms).
Chris@10 467
Chris@10 468 @item
Chris@10 469 @ctindex FFTW_EXHAUSTIVE
Chris@10 470 @code{FFTW_EXHAUSTIVE} is like @code{FFTW_PATIENT}, but considers an
Chris@10 471 even wider range of algorithms, including many that we think are
Chris@10 472 unlikely to be fast, to produce the most optimal plan but with a
Chris@10 473 substantially increased planning time.
Chris@10 474
Chris@10 475 @item
Chris@10 476 @ctindex FFTW_WISDOM_ONLY
Chris@10 477 @code{FFTW_WISDOM_ONLY} is a special planning mode in which the plan
Chris@10 478 is only created if wisdom is available for the given problem, and
Chris@10 479 otherwise a @code{NULL} plan is returned. This can be combined with
Chris@10 480 other flags, e.g. @samp{FFTW_WISDOM_ONLY | FFTW_PATIENT} creates a
Chris@10 481 plan only if wisdom is available that was created in
Chris@10 482 @code{FFTW_PATIENT} or @code{FFTW_EXHAUSTIVE} mode. The
Chris@10 483 @code{FFTW_WISDOM_ONLY} flag is intended for users who need to detect
Chris@10 484 whether wisdom is available; for example, if wisdom is not available
Chris@10 485 one may wish to allocate new arrays for planning so that user data is
Chris@10 486 not overwritten.
Chris@10 487
Chris@10 488 @end itemize
Chris@10 489
Chris@10 490 @subsubheading Algorithm-restriction flags
Chris@10 491 @itemize @bullet
Chris@10 492
Chris@10 493 @item
Chris@10 494 @ctindex FFTW_DESTROY_INPUT
Chris@10 495 @code{FFTW_DESTROY_INPUT} specifies that an out-of-place transform is
Chris@10 496 allowed to @emph{overwrite its input} array with arbitrary data; this
Chris@10 497 can sometimes allow more efficient algorithms to be employed.
Chris@10 498 @cindex out-of-place
Chris@10 499
Chris@10 500 @item
Chris@10 501 @ctindex FFTW_PRESERVE_INPUT
Chris@10 502 @code{FFTW_PRESERVE_INPUT} specifies that an out-of-place transform must
Chris@10 503 @emph{not change its input} array. This is ordinarily the
Chris@10 504 @emph{default}, except for c2r and hc2r (i.e. complex-to-real)
Chris@10 505 transforms for which @code{FFTW_DESTROY_INPUT} is the default. In the
Chris@10 506 latter cases, passing @code{FFTW_PRESERVE_INPUT} will attempt to use
Chris@10 507 algorithms that do not destroy the input, at the expense of worse
Chris@10 508 performance; for multi-dimensional c2r transforms, however, no
Chris@10 509 input-preserving algorithms are implemented and the planner will return
Chris@10 510 @code{NULL} if one is requested.
Chris@10 511 @cindex c2r
Chris@10 512 @cindex hc2r
Chris@10 513
Chris@10 514 @item
Chris@10 515 @ctindex FFTW_UNALIGNED
Chris@10 516 @cindex alignment
Chris@10 517 @code{FFTW_UNALIGNED} specifies that the algorithm may not impose any
Chris@10 518 unusual alignment requirements on the input/output arrays (i.e. no
Chris@10 519 SIMD may be used). This flag is normally @emph{not necessary}, since
Chris@10 520 the planner automatically detects misaligned arrays. The only use for
Chris@10 521 this flag is if you want to use the new-array execute interface to
Chris@10 522 execute a given plan on a different array that may not be aligned like
Chris@10 523 the original. (Using @code{fftw_malloc} makes this flag unnecessary
Chris@10 524 even then.)
Chris@10 525
Chris@10 526 @end itemize
Chris@10 527
Chris@10 528 @subsubheading Limiting planning time
Chris@10 529
Chris@10 530 @example
Chris@10 531 extern void fftw_set_timelimit(double seconds);
Chris@10 532 @end example
Chris@10 533 @findex fftw_set_timelimit
Chris@10 534
Chris@10 535 This function instructs FFTW to spend at most @code{seconds} seconds
Chris@10 536 (approximately) in the planner. If @code{seconds ==
Chris@10 537 FFTW_NO_TIMELIMIT} (the default value, which is negative), then
Chris@10 538 planning time is unbounded. Otherwise, FFTW plans with a
Chris@10 539 progressively wider range of algorithms until the the given time limit
Chris@10 540 is reached or the given range of algorithms is explored, returning the
Chris@10 541 best available plan.
Chris@10 542 @ctindex FFTW_NO_TIMELIMIT
Chris@10 543
Chris@10 544
Chris@10 545 For example, specifying @code{FFTW_PATIENT} first plans in
Chris@10 546 @code{FFTW_ESTIMATE} mode, then in @code{FFTW_MEASURE} mode, then
Chris@10 547 finally (time permitting) in @code{FFTW_PATIENT}. If
Chris@10 548 @code{FFTW_EXHAUSTIVE} is specified instead, the planner will further
Chris@10 549 progress to @code{FFTW_EXHAUSTIVE} mode.
Chris@10 550
Chris@10 551 Note that the @code{seconds} argument specifies only a rough limit; in
Chris@10 552 practice, the planner may use somewhat more time if the time limit is
Chris@10 553 reached when the planner is in the middle of an operation that cannot
Chris@10 554 be interrupted. At the very least, the planner will complete planning
Chris@10 555 in @code{FFTW_ESTIMATE} mode (which is thus equivalent to a time limit
Chris@10 556 of 0).
Chris@10 557
Chris@10 558
Chris@10 559 @c =========>
Chris@10 560 @node Real-data DFTs, Real-data DFT Array Format, Planner Flags, Basic Interface
Chris@10 561 @subsection Real-data DFTs
Chris@10 562
Chris@10 563 @example
Chris@10 564 fftw_plan fftw_plan_dft_r2c_1d(int n0,
Chris@10 565 double *in, fftw_complex *out,
Chris@10 566 unsigned flags);
Chris@10 567 fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1,
Chris@10 568 double *in, fftw_complex *out,
Chris@10 569 unsigned flags);
Chris@10 570 fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, int n2,
Chris@10 571 double *in, fftw_complex *out,
Chris@10 572 unsigned flags);
Chris@10 573 fftw_plan fftw_plan_dft_r2c(int rank, const int *n,
Chris@10 574 double *in, fftw_complex *out,
Chris@10 575 unsigned flags);
Chris@10 576 @end example
Chris@10 577 @findex fftw_plan_dft_r2c_1d
Chris@10 578 @findex fftw_plan_dft_r2c_2d
Chris@10 579 @findex fftw_plan_dft_r2c_3d
Chris@10 580 @findex fftw_plan_dft_r2c
Chris@10 581 @cindex r2c
Chris@10 582
Chris@10 583 Plan a real-input/complex-output discrete Fourier transform (DFT) in
Chris@10 584 zero or more dimensions, returning an @code{fftw_plan} (@pxref{Using
Chris@10 585 Plans}).
Chris@10 586
Chris@10 587 Once you have created a plan for a certain transform type and
Chris@10 588 parameters, then creating another plan of the same type and parameters,
Chris@10 589 but for different arrays, is fast and shares constant data with the
Chris@10 590 first plan (if it still exists).
Chris@10 591
Chris@10 592 The planner returns @code{NULL} if the plan cannot be created. A
Chris@10 593 non-@code{NULL} plan is always returned by the basic interface unless
Chris@10 594 you are using a customized FFTW configuration supporting a restricted
Chris@10 595 set of transforms, or if you use the @code{FFTW_PRESERVE_INPUT} flag
Chris@10 596 with a multi-dimensional out-of-place c2r transform (see below).
Chris@10 597
Chris@10 598 @subsubheading Arguments
Chris@10 599 @itemize @bullet
Chris@10 600
Chris@10 601 @item
Chris@10 602 @code{rank} is the rank of the transform (it should be the size of the
Chris@10 603 array @code{*n}), and can be any non-negative integer. (@xref{Complex
Chris@10 604 Multi-Dimensional DFTs}, for the definition of ``rank''.) The
Chris@10 605 @samp{_1d}, @samp{_2d}, and @samp{_3d} planners correspond to a
Chris@10 606 @code{rank} of @code{1}, @code{2}, and @code{3}, respectively. The rank
Chris@10 607 may be zero, which is equivalent to a rank-1 transform of size 1, i.e. a
Chris@10 608 copy of one real number (with zero imaginary part) from input to output.
Chris@10 609
Chris@10 610 @item
Chris@10 611 @code{n0}, @code{n1}, @code{n2}, or @code{n[0..rank-1]}, (as appropriate
Chris@10 612 for each routine) specify the size of the transform dimensions. They
Chris@10 613 can be any positive integer. This is different in general from the
Chris@10 614 @emph{physical} array dimensions, which are described in @ref{Real-data
Chris@10 615 DFT Array Format}.
Chris@10 616
Chris@10 617 @itemize @minus
Chris@10 618 @item
Chris@10 619 FFTW is best at handling sizes of the form
Chris@10 620 @ifinfo
Chris@10 621 @math{2^a 3^b 5^c 7^d 11^e 13^f},
Chris@10 622 @end ifinfo
Chris@10 623 @tex
Chris@10 624 $2^a 3^b 5^c 7^d 11^e 13^f$,
Chris@10 625 @end tex
Chris@10 626 @html
Chris@10 627 2<sup>a</sup> 3<sup>b</sup> 5<sup>c</sup> 7<sup>d</sup>
Chris@10 628 11<sup>e</sup> 13<sup>f</sup>,
Chris@10 629 @end html
Chris@10 630 where @math{e+f} is either @math{0} or @math{1}, and the other exponents
Chris@10 631 are arbitrary. Other sizes are computed by means of a slow,
Chris@10 632 general-purpose algorithm (which nevertheless retains @Onlogn{} performance even for prime sizes). (It is possible to customize FFTW
Chris@10 633 for different array sizes; see @ref{Installation and Customization}.)
Chris@10 634 Transforms whose sizes are powers of @math{2} are especially fast, and
Chris@10 635 it is generally beneficial for the @emph{last} dimension of an r2c/c2r
Chris@10 636 transform to be @emph{even}.
Chris@10 637 @end itemize
Chris@10 638
Chris@10 639 @item
Chris@10 640 @code{in} and @code{out} point to the input and output arrays of the
Chris@10 641 transform, which may be the same (yielding an in-place transform).
Chris@10 642 @cindex in-place
Chris@10 643 These arrays are overwritten during planning, unless
Chris@10 644 @code{FFTW_ESTIMATE} is used in the flags. (The arrays need not be
Chris@10 645 initialized, but they must be allocated.) For an in-place transform, it
Chris@10 646 is important to remember that the real array will require padding,
Chris@10 647 described in @ref{Real-data DFT Array Format}.
Chris@10 648 @cindex padding
Chris@10 649
Chris@10 650 @item
Chris@10 651 @cindex flags
Chris@10 652 @code{flags} is a bitwise OR (@samp{|}) of zero or more planner flags,
Chris@10 653 as defined in @ref{Planner Flags}.
Chris@10 654
Chris@10 655 @end itemize
Chris@10 656
Chris@10 657 The inverse transforms, taking complex input (storing the non-redundant
Chris@10 658 half of a logically Hermitian array) to real output, are given by:
Chris@10 659
Chris@10 660 @example
Chris@10 661 fftw_plan fftw_plan_dft_c2r_1d(int n0,
Chris@10 662 fftw_complex *in, double *out,
Chris@10 663 unsigned flags);
Chris@10 664 fftw_plan fftw_plan_dft_c2r_2d(int n0, int n1,
Chris@10 665 fftw_complex *in, double *out,
Chris@10 666 unsigned flags);
Chris@10 667 fftw_plan fftw_plan_dft_c2r_3d(int n0, int n1, int n2,
Chris@10 668 fftw_complex *in, double *out,
Chris@10 669 unsigned flags);
Chris@10 670 fftw_plan fftw_plan_dft_c2r(int rank, const int *n,
Chris@10 671 fftw_complex *in, double *out,
Chris@10 672 unsigned flags);
Chris@10 673 @end example
Chris@10 674 @findex fftw_plan_dft_c2r_1d
Chris@10 675 @findex fftw_plan_dft_c2r_2d
Chris@10 676 @findex fftw_plan_dft_c2r_3d
Chris@10 677 @findex fftw_plan_dft_c2r
Chris@10 678 @cindex c2r
Chris@10 679
Chris@10 680 The arguments are the same as for the r2c transforms, except that the
Chris@10 681 input and output data formats are reversed.
Chris@10 682
Chris@10 683 FFTW computes an unnormalized transform: computing an r2c followed by a
Chris@10 684 c2r transform (or vice versa) will result in the original data
Chris@10 685 multiplied by the size of the transform (the product of the logical
Chris@10 686 dimensions).
Chris@10 687 @cindex normalization
Chris@10 688 An r2c transform produces the same output as a @code{FFTW_FORWARD}
Chris@10 689 complex DFT of the same input, and a c2r transform is correspondingly
Chris@10 690 equivalent to @code{FFTW_BACKWARD}. For more information, see @ref{What
Chris@10 691 FFTW Really Computes}.
Chris@10 692
Chris@10 693 @c =========>
Chris@10 694 @node Real-data DFT Array Format, Real-to-Real Transforms, Real-data DFTs, Basic Interface
Chris@10 695 @subsection Real-data DFT Array Format
Chris@10 696 @cindex r2c/c2r multi-dimensional array format
Chris@10 697
Chris@10 698 The output of a DFT of real data (r2c) contains symmetries that, in
Chris@10 699 principle, make half of the outputs redundant (@pxref{What FFTW Really
Chris@10 700 Computes}). (Similarly for the input of an inverse c2r transform.) In
Chris@10 701 practice, it is not possible to entirely realize these savings in an
Chris@10 702 efficient and understandable format that generalizes to
Chris@10 703 multi-dimensional transforms. Instead, the output of the r2c
Chris@10 704 transforms is @emph{slightly} over half of the output of the
Chris@10 705 corresponding complex transform. We do not ``pack'' the data in any
Chris@10 706 way, but store it as an ordinary array of @code{fftw_complex} values.
Chris@10 707 In fact, this data is simply a subsection of what would be the array in
Chris@10 708 the corresponding complex transform.
Chris@10 709
Chris@10 710 Specifically, for a real transform of @math{d} (= @code{rank})
Chris@10 711 dimensions @ndims{}, the complex data is an @ndimshalf array of
Chris@10 712 @code{fftw_complex} values in row-major order (with the division rounded
Chris@10 713 down). That is, we only store the @emph{lower} half (non-negative
Chris@10 714 frequencies), plus one element, of the last dimension of the data from
Chris@10 715 the ordinary complex transform. (We could have instead taken half of
Chris@10 716 any other dimension, but implementation turns out to be simpler if the
Chris@10 717 last, contiguous, dimension is used.)
Chris@10 718
Chris@10 719 @cindex out-of-place
Chris@10 720 For an out-of-place transform, the real data is simply an array with
Chris@10 721 physical dimensions @ndims in row-major order.
Chris@10 722
Chris@10 723 @cindex in-place
Chris@10 724 @cindex padding
Chris@10 725 For an in-place transform, some complications arise since the complex data
Chris@10 726 is slightly larger than the real data. In this case, the final
Chris@10 727 dimension of the real data must be @emph{padded} with extra values to
Chris@10 728 accommodate the size of the complex data---two extra if the last
Chris@10 729 dimension is even and one if it is odd. That is, the last dimension of
Chris@10 730 the real data must physically contain
Chris@10 731 @tex
Chris@10 732 $2 (n_{d-1}/2+1)$
Chris@10 733 @end tex
Chris@10 734 @ifinfo
Chris@10 735 2 * (n[d-1]/2+1)
Chris@10 736 @end ifinfo
Chris@10 737 @html
Chris@10 738 2 * (n<sub>d-1</sub>/2+1)
Chris@10 739 @end html
Chris@10 740 @code{double} values (exactly enough to hold the complex data). This
Chris@10 741 physical array size does not, however, change the @emph{logical} array
Chris@10 742 size---only
Chris@10 743 @tex
Chris@10 744 $n_{d-1}$
Chris@10 745 @end tex
Chris@10 746 @ifinfo
Chris@10 747 n[d-1]
Chris@10 748 @end ifinfo
Chris@10 749 @html
Chris@10 750 n<sub>d-1</sub>
Chris@10 751 @end html
Chris@10 752 values are actually stored in the last dimension, and
Chris@10 753 @tex
Chris@10 754 $n_{d-1}$
Chris@10 755 @end tex
Chris@10 756 @ifinfo
Chris@10 757 n[d-1]
Chris@10 758 @end ifinfo
Chris@10 759 @html
Chris@10 760 n<sub>d-1</sub>
Chris@10 761 @end html
Chris@10 762 is the last dimension passed to the planner.
Chris@10 763
Chris@10 764 @c =========>
Chris@10 765 @node Real-to-Real Transforms, Real-to-Real Transform Kinds, Real-data DFT Array Format, Basic Interface
Chris@10 766 @subsection Real-to-Real Transforms
Chris@10 767 @cindex r2r
Chris@10 768
Chris@10 769 @example
Chris@10 770 fftw_plan fftw_plan_r2r_1d(int n, double *in, double *out,
Chris@10 771 fftw_r2r_kind kind, unsigned flags);
Chris@10 772 fftw_plan fftw_plan_r2r_2d(int n0, int n1, double *in, double *out,
Chris@10 773 fftw_r2r_kind kind0, fftw_r2r_kind kind1,
Chris@10 774 unsigned flags);
Chris@10 775 fftw_plan fftw_plan_r2r_3d(int n0, int n1, int n2,
Chris@10 776 double *in, double *out,
Chris@10 777 fftw_r2r_kind kind0,
Chris@10 778 fftw_r2r_kind kind1,
Chris@10 779 fftw_r2r_kind kind2,
Chris@10 780 unsigned flags);
Chris@10 781 fftw_plan fftw_plan_r2r(int rank, const int *n, double *in, double *out,
Chris@10 782 const fftw_r2r_kind *kind, unsigned flags);
Chris@10 783 @end example
Chris@10 784 @findex fftw_plan_r2r_1d
Chris@10 785 @findex fftw_plan_r2r_2d
Chris@10 786 @findex fftw_plan_r2r_3d
Chris@10 787 @findex fftw_plan_r2r
Chris@10 788
Chris@10 789 Plan a real input/output (r2r) transform of various kinds in zero or
Chris@10 790 more dimensions, returning an @code{fftw_plan} (@pxref{Using Plans}).
Chris@10 791
Chris@10 792 Once you have created a plan for a certain transform type and
Chris@10 793 parameters, then creating another plan of the same type and parameters,
Chris@10 794 but for different arrays, is fast and shares constant data with the
Chris@10 795 first plan (if it still exists).
Chris@10 796
Chris@10 797 The planner returns @code{NULL} if the plan cannot be created. A
Chris@10 798 non-@code{NULL} plan is always returned by the basic interface unless
Chris@10 799 you are using a customized FFTW configuration supporting a restricted
Chris@10 800 set of transforms, or for size-1 @code{FFTW_REDFT00} kinds (which are
Chris@10 801 not defined).
Chris@10 802 @ctindex FFTW_REDFT00
Chris@10 803
Chris@10 804 @subsubheading Arguments
Chris@10 805 @itemize @bullet
Chris@10 806
Chris@10 807 @item
Chris@10 808 @code{rank} is the dimensionality of the transform (it should be the
Chris@10 809 size of the arrays @code{*n} and @code{*kind}), and can be any
Chris@10 810 non-negative integer. The @samp{_1d}, @samp{_2d}, and @samp{_3d}
Chris@10 811 planners correspond to a @code{rank} of @code{1}, @code{2}, and
Chris@10 812 @code{3}, respectively. A @code{rank} of zero is equivalent to a copy
Chris@10 813 of one number from input to output.
Chris@10 814
Chris@10 815 @item
Chris@10 816 @code{n}, or @code{n0}/@code{n1}/@code{n2}, or @code{n[rank]},
Chris@10 817 respectively, gives the (physical) size of the transform dimensions.
Chris@10 818 They can be any positive integer.
Chris@10 819
Chris@10 820 @itemize @minus
Chris@10 821 @item
Chris@10 822 @cindex row-major
Chris@10 823 Multi-dimensional arrays are stored in row-major order with dimensions:
Chris@10 824 @code{n0} x @code{n1}; or @code{n0} x @code{n1} x @code{n2}; or
Chris@10 825 @code{n[0]} x @code{n[1]} x ... x @code{n[rank-1]}.
Chris@10 826 @xref{Multi-dimensional Array Format}.
Chris@10 827 @item
Chris@10 828 FFTW is generally best at handling sizes of the form
Chris@10 829 @ifinfo
Chris@10 830 @math{2^a 3^b 5^c 7^d 11^e 13^f},
Chris@10 831 @end ifinfo
Chris@10 832 @tex
Chris@10 833 $2^a 3^b 5^c 7^d 11^e 13^f$,
Chris@10 834 @end tex
Chris@10 835 @html
Chris@10 836 2<sup>a</sup> 3<sup>b</sup> 5<sup>c</sup> 7<sup>d</sup>
Chris@10 837 11<sup>e</sup> 13<sup>f</sup>,
Chris@10 838 @end html
Chris@10 839 where @math{e+f} is either @math{0} or @math{1}, and the other exponents
Chris@10 840 are arbitrary. Other sizes are computed by means of a slow,
Chris@10 841 general-purpose algorithm (which nevertheless retains @Onlogn{} performance even for prime sizes). (It is possible to customize FFTW
Chris@10 842 for different array sizes; see @ref{Installation and Customization}.)
Chris@10 843 Transforms whose sizes are powers of @math{2} are especially fast.
Chris@10 844 @item
Chris@10 845 For a @code{REDFT00} or @code{RODFT00} transform kind in a dimension of
Chris@10 846 size @math{n}, it is @math{n-1} or @math{n+1}, respectively, that
Chris@10 847 should be factorizable in the above form.
Chris@10 848 @end itemize
Chris@10 849
Chris@10 850 @item
Chris@10 851 @code{in} and @code{out} point to the input and output arrays of the
Chris@10 852 transform, which may be the same (yielding an in-place transform).
Chris@10 853 @cindex in-place
Chris@10 854 These arrays are overwritten during planning, unless
Chris@10 855 @code{FFTW_ESTIMATE} is used in the flags. (The arrays need not be
Chris@10 856 initialized, but they must be allocated.)
Chris@10 857
Chris@10 858 @item
Chris@10 859 @code{kind}, or @code{kind0}/@code{kind1}/@code{kind2}, or
Chris@10 860 @code{kind[rank]}, is the kind of r2r transform used for the
Chris@10 861 corresponding dimension. The valid kind constants are described in
Chris@10 862 @ref{Real-to-Real Transform Kinds}. In a multi-dimensional transform,
Chris@10 863 what is computed is the separable product formed by taking each
Chris@10 864 transform kind along the corresponding dimension, one dimension after
Chris@10 865 another.
Chris@10 866
Chris@10 867 @item
Chris@10 868 @cindex flags
Chris@10 869 @code{flags} is a bitwise OR (@samp{|}) of zero or more planner flags,
Chris@10 870 as defined in @ref{Planner Flags}.
Chris@10 871
Chris@10 872 @end itemize
Chris@10 873
Chris@10 874 @c =========>
Chris@10 875 @node Real-to-Real Transform Kinds, , Real-to-Real Transforms, Basic Interface
Chris@10 876 @subsection Real-to-Real Transform Kinds
Chris@10 877 @cindex kind (r2r)
Chris@10 878
Chris@10 879 FFTW currently supports 11 different r2r transform kinds, specified by
Chris@10 880 one of the constants below. For the precise definitions of these
Chris@10 881 transforms, see @ref{What FFTW Really Computes}. For a more colloquial
Chris@10 882 introduction to these transform kinds, see @ref{More DFTs of Real Data}.
Chris@10 883
Chris@10 884 For dimension of size @code{n}, there is a corresponding ``logical''
Chris@10 885 dimension @code{N} that determines the normalization (and the optimal
Chris@10 886 factorization); the formula for @code{N} is given for each kind below.
Chris@10 887 Also, with each transform kind is listed its corrsponding inverse
Chris@10 888 transform. FFTW computes unnormalized transforms: a transform followed
Chris@10 889 by its inverse will result in the original data multiplied by @code{N}
Chris@10 890 (or the product of the @code{N}'s for each dimension, in
Chris@10 891 multi-dimensions).
Chris@10 892 @cindex normalization
Chris@10 893
Chris@10 894 @itemize @bullet
Chris@10 895
Chris@10 896 @item
Chris@10 897 @ctindex FFTW_R2HC
Chris@10 898 @code{FFTW_R2HC} computes a real-input DFT with output in
Chris@10 899 ``halfcomplex'' format, i.e. real and imaginary parts for a transform of
Chris@10 900 size @code{n} stored as:
Chris@10 901 @tex
Chris@10 902 $$
Chris@10 903 r_0, r_1, r_2, \ldots, r_{n/2}, i_{(n+1)/2-1}, \ldots, i_2, i_1
Chris@10 904 $$
Chris@10 905 @end tex
Chris@10 906 @ifinfo
Chris@10 907 r0, r1, r2, r(n/2), i((n+1)/2-1), ..., i2, i1
Chris@10 908 @end ifinfo
Chris@10 909 @html
Chris@10 910 <p align=center>
Chris@10 911 r<sub>0</sub>, r<sub>1</sub>, r<sub>2</sub>, ..., r<sub>n/2</sub>, i<sub>(n+1)/2-1</sub>, ..., i<sub>2</sub>, i<sub>1</sub>
Chris@10 912 </p>
Chris@10 913 @end html
Chris@10 914 (Logical @code{N=n}, inverse is @code{FFTW_HC2R}.)
Chris@10 915
Chris@10 916 @item
Chris@10 917 @ctindex FFTW_HC2R
Chris@10 918 @code{FFTW_HC2R} computes the reverse of @code{FFTW_R2HC}, above.
Chris@10 919 (Logical @code{N=n}, inverse is @code{FFTW_R2HC}.)
Chris@10 920
Chris@10 921 @item
Chris@10 922 @ctindex FFTW_DHT
Chris@10 923 @code{FFTW_DHT} computes a discrete Hartley transform.
Chris@10 924 (Logical @code{N=n}, inverse is @code{FFTW_DHT}.)
Chris@10 925 @cindex discrete Hartley transform
Chris@10 926
Chris@10 927 @item
Chris@10 928 @ctindex FFTW_REDFT00
Chris@10 929 @code{FFTW_REDFT00} computes an REDFT00 transform, i.e. a DCT-I.
Chris@10 930 (Logical @code{N=2*(n-1)}, inverse is @code{FFTW_REDFT00}.)
Chris@10 931 @cindex discrete cosine transform
Chris@10 932 @cindex DCT
Chris@10 933
Chris@10 934 @item
Chris@10 935 @ctindex FFTW_REDFT10
Chris@10 936 @code{FFTW_REDFT10} computes an REDFT10 transform, i.e. a DCT-II (sometimes called ``the'' DCT).
Chris@10 937 (Logical @code{N=2*n}, inverse is @code{FFTW_REDFT01}.)
Chris@10 938
Chris@10 939 @item
Chris@10 940 @ctindex FFTW_REDFT01
Chris@10 941 @code{FFTW_REDFT01} computes an REDFT01 transform, i.e. a DCT-III (sometimes called ``the'' IDCT, being the inverse of DCT-II).
Chris@10 942 (Logical @code{N=2*n}, inverse is @code{FFTW_REDFT=10}.)
Chris@10 943 @cindex IDCT
Chris@10 944
Chris@10 945 @item
Chris@10 946 @ctindex FFTW_REDFT11
Chris@10 947 @code{FFTW_REDFT11} computes an REDFT11 transform, i.e. a DCT-IV.
Chris@10 948 (Logical @code{N=2*n}, inverse is @code{FFTW_REDFT11}.)
Chris@10 949
Chris@10 950 @item
Chris@10 951 @ctindex FFTW_RODFT00
Chris@10 952 @code{FFTW_RODFT00} computes an RODFT00 transform, i.e. a DST-I.
Chris@10 953 (Logical @code{N=2*(n+1)}, inverse is @code{FFTW_RODFT00}.)
Chris@10 954 @cindex discrete sine transform
Chris@10 955 @cindex DST
Chris@10 956
Chris@10 957 @item
Chris@10 958 @ctindex FFTW_RODFT10
Chris@10 959 @code{FFTW_RODFT10} computes an RODFT10 transform, i.e. a DST-II.
Chris@10 960 (Logical @code{N=2*n}, inverse is @code{FFTW_RODFT01}.)
Chris@10 961
Chris@10 962 @item
Chris@10 963 @ctindex FFTW_RODFT01
Chris@10 964 @code{FFTW_RODFT01} computes an RODFT01 transform, i.e. a DST-III.
Chris@10 965 (Logical @code{N=2*n}, inverse is @code{FFTW_RODFT=10}.)
Chris@10 966
Chris@10 967 @item
Chris@10 968 @ctindex FFTW_RODFT11
Chris@10 969 @code{FFTW_RODFT11} computes an RODFT11 transform, i.e. a DST-IV.
Chris@10 970 (Logical @code{N=2*n}, inverse is @code{FFTW_RODFT11}.)
Chris@10 971
Chris@10 972 @end itemize
Chris@10 973
Chris@10 974 @c ------------------------------------------------------------
Chris@10 975 @node Advanced Interface, Guru Interface, Basic Interface, FFTW Reference
Chris@10 976 @section Advanced Interface
Chris@10 977 @cindex advanced interface
Chris@10 978
Chris@10 979 FFTW's ``advanced'' interface supplements the basic interface with four
Chris@10 980 new planner routines, providing a new level of flexibility: you can plan
Chris@10 981 a transform of multiple arrays simultaneously, operate on non-contiguous
Chris@10 982 (strided) data, and transform a subset of a larger multi-dimensional
Chris@10 983 array. Other than these additional features, the planner operates in
Chris@10 984 the same fashion as in the basic interface, and the resulting
Chris@10 985 @code{fftw_plan} is used in the same way (@pxref{Using Plans}).
Chris@10 986
Chris@10 987 @menu
Chris@10 988 * Advanced Complex DFTs::
Chris@10 989 * Advanced Real-data DFTs::
Chris@10 990 * Advanced Real-to-real Transforms::
Chris@10 991 @end menu
Chris@10 992
Chris@10 993 @c =========>
Chris@10 994 @node Advanced Complex DFTs, Advanced Real-data DFTs, Advanced Interface, Advanced Interface
Chris@10 995 @subsection Advanced Complex DFTs
Chris@10 996
Chris@10 997 @example
Chris@10 998 fftw_plan fftw_plan_many_dft(int rank, const int *n, int howmany,
Chris@10 999 fftw_complex *in, const int *inembed,
Chris@10 1000 int istride, int idist,
Chris@10 1001 fftw_complex *out, const int *onembed,
Chris@10 1002 int ostride, int odist,
Chris@10 1003 int sign, unsigned flags);
Chris@10 1004 @end example
Chris@10 1005 @findex fftw_plan_many_dft
Chris@10 1006
Chris@10 1007 This routine plans multiple multidimensional complex DFTs, and it
Chris@10 1008 extends the @code{fftw_plan_dft} routine (@pxref{Complex DFTs}) to
Chris@10 1009 compute @code{howmany} transforms, each having rank @code{rank} and size
Chris@10 1010 @code{n}. In addition, the transform data need not be contiguous, but
Chris@10 1011 it may be laid out in memory with an arbitrary stride. To account for
Chris@10 1012 these possibilities, @code{fftw_plan_many_dft} adds the new parameters
Chris@10 1013 @code{howmany}, @{@code{i},@code{o}@}@code{nembed},
Chris@10 1014 @{@code{i},@code{o}@}@code{stride}, and
Chris@10 1015 @{@code{i},@code{o}@}@code{dist}. The FFTW basic interface
Chris@10 1016 (@pxref{Complex DFTs}) provides routines specialized for ranks 1, 2,
Chris@10 1017 and@tie{}3, but the advanced interface handles only the general-rank
Chris@10 1018 case.
Chris@10 1019
Chris@10 1020 @code{howmany} is the number of transforms to compute. The resulting
Chris@10 1021 plan computes @code{howmany} transforms, where the input of the
Chris@10 1022 @code{k}-th transform is at location @code{in+k*idist} (in C pointer
Chris@10 1023 arithmetic), and its output is at location @code{out+k*odist}. Plans
Chris@10 1024 obtained in this way can often be faster than calling FFTW multiple
Chris@10 1025 times for the individual transforms. The basic @code{fftw_plan_dft}
Chris@10 1026 interface corresponds to @code{howmany=1} (in which case the @code{dist}
Chris@10 1027 parameters are ignored).
Chris@10 1028 @cindex howmany parameter
Chris@10 1029 @cindex dist
Chris@10 1030
Chris@10 1031
Chris@10 1032 Each of the @code{howmany} transforms has rank @code{rank} and size
Chris@10 1033 @code{n}, as in the basic interface. In addition, the advanced
Chris@10 1034 interface allows the input and output arrays of each transform to be
Chris@10 1035 row-major subarrays of larger rank-@code{rank} arrays, described by
Chris@10 1036 @code{inembed} and @code{onembed} parameters, respectively.
Chris@10 1037 @{@code{i},@code{o}@}@code{nembed} must be arrays of length @code{rank},
Chris@10 1038 and @code{n} should be elementwise less than or equal to
Chris@10 1039 @{@code{i},@code{o}@}@code{nembed}. Passing @code{NULL} for an
Chris@10 1040 @code{nembed} parameter is equivalent to passing @code{n} (i.e. same
Chris@10 1041 physical and logical dimensions, as in the basic interface.)
Chris@10 1042
Chris@10 1043 The @code{stride} parameters indicate that the @code{j}-th element of
Chris@10 1044 the input or output arrays is located at @code{j*istride} or
Chris@10 1045 @code{j*ostride}, respectively. (For a multi-dimensional array,
Chris@10 1046 @code{j} is the ordinary row-major index.) When combined with the
Chris@10 1047 @code{k}-th transform in a @code{howmany} loop, from above, this means
Chris@10 1048 that the (@code{j},@code{k})-th element is at @code{j*stride+k*dist}.
Chris@10 1049 (The basic @code{fftw_plan_dft} interface corresponds to a stride of 1.)
Chris@10 1050 @cindex stride
Chris@10 1051
Chris@10 1052
Chris@10 1053 For in-place transforms, the input and output @code{stride} and
Chris@10 1054 @code{dist} parameters should be the same; otherwise, the planner may
Chris@10 1055 return @code{NULL}.
Chris@10 1056
Chris@10 1057 Arrays @code{n}, @code{inembed}, and @code{onembed} are not used after
Chris@10 1058 this function returns. You can safely free or reuse them.
Chris@10 1059
Chris@10 1060 @strong{Examples}:
Chris@10 1061 One transform of one 5 by 6 array contiguous in memory:
Chris@10 1062 @example
Chris@10 1063 int rank = 2;
Chris@10 1064 int n[] = @{5, 6@};
Chris@10 1065 int howmany = 1;
Chris@10 1066 int idist = odist = 0; /* unused because howmany = 1 */
Chris@10 1067 int istride = ostride = 1; /* array is contiguous in memory */
Chris@10 1068 int *inembed = n, *onembed = n;
Chris@10 1069 @end example
Chris@10 1070
Chris@10 1071 Transform of three 5 by 6 arrays, each contiguous in memory,
Chris@10 1072 stored in memory one after another:
Chris@10 1073 @example
Chris@10 1074 int rank = 2;
Chris@10 1075 int n[] = @{5, 6@};
Chris@10 1076 int howmany = 3;
Chris@10 1077 int idist = odist = n[0]*n[1]; /* = 30, the distance in memory
Chris@10 1078 between the first element
Chris@10 1079 of the first array and the
Chris@10 1080 first element of the second array */
Chris@10 1081 int istride = ostride = 1; /* array is contiguous in memory */
Chris@10 1082 int *inembed = n, *onembed = n;
Chris@10 1083 @end example
Chris@10 1084
Chris@10 1085 Transform each column of a 2d array with 10 rows and 3 columns:
Chris@10 1086 @example
Chris@10 1087 int rank = 1; /* not 2: we are computing 1d transforms */
Chris@10 1088 int n[] = @{10@}; /* 1d transforms of length 10 */
Chris@10 1089 int howmany = 3;
Chris@10 1090 int idist = odist = 1;
Chris@10 1091 int istride = ostride = 3; /* distance between two elements in
Chris@10 1092 the same column */
Chris@10 1093 int *inembed = n, *onembed = n;
Chris@10 1094 @end example
Chris@10 1095
Chris@10 1096 @c =========>
Chris@10 1097 @node Advanced Real-data DFTs, Advanced Real-to-real Transforms, Advanced Complex DFTs, Advanced Interface
Chris@10 1098 @subsection Advanced Real-data DFTs
Chris@10 1099
Chris@10 1100 @example
Chris@10 1101 fftw_plan fftw_plan_many_dft_r2c(int rank, const int *n, int howmany,
Chris@10 1102 double *in, const int *inembed,
Chris@10 1103 int istride, int idist,
Chris@10 1104 fftw_complex *out, const int *onembed,
Chris@10 1105 int ostride, int odist,
Chris@10 1106 unsigned flags);
Chris@10 1107 fftw_plan fftw_plan_many_dft_c2r(int rank, const int *n, int howmany,
Chris@10 1108 fftw_complex *in, const int *inembed,
Chris@10 1109 int istride, int idist,
Chris@10 1110 double *out, const int *onembed,
Chris@10 1111 int ostride, int odist,
Chris@10 1112 unsigned flags);
Chris@10 1113 @end example
Chris@10 1114 @findex fftw_plan_many_dft_r2c
Chris@10 1115 @findex fftw_plan_many_dft_c2r
Chris@10 1116
Chris@10 1117 Like @code{fftw_plan_many_dft}, these two functions add @code{howmany},
Chris@10 1118 @code{nembed}, @code{stride}, and @code{dist} parameters to the
Chris@10 1119 @code{fftw_plan_dft_r2c} and @code{fftw_plan_dft_c2r} functions, but
Chris@10 1120 otherwise behave the same as the basic interface.
Chris@10 1121
Chris@10 1122 The interpretation of @code{howmany}, @code{stride}, and @code{dist} are
Chris@10 1123 the same as for @code{fftw_plan_many_dft}, above. Note that the
Chris@10 1124 @code{stride} and @code{dist} for the real array are in units of
Chris@10 1125 @code{double}, and for the complex array are in units of
Chris@10 1126 @code{fftw_complex}.
Chris@10 1127
Chris@10 1128 If an @code{nembed} parameter is @code{NULL}, it is interpreted as what
Chris@10 1129 it would be in the basic interface, as described in @ref{Real-data DFT
Chris@10 1130 Array Format}. That is, for the complex array the size is assumed to be
Chris@10 1131 the same as @code{n}, but with the last dimension cut roughly in half.
Chris@10 1132 For the real array, the size is assumed to be @code{n} if the transform
Chris@10 1133 is out-of-place, or @code{n} with the last dimension ``padded'' if the
Chris@10 1134 transform is in-place.
Chris@10 1135
Chris@10 1136 If an @code{nembed} parameter is non-@code{NULL}, it is interpreted as
Chris@10 1137 the physical size of the corresponding array, in row-major order, just
Chris@10 1138 as for @code{fftw_plan_many_dft}. In this case, each dimension of
Chris@10 1139 @code{nembed} should be @code{>=} what it would be in the basic
Chris@10 1140 interface (e.g. the halved or padded @code{n}).
Chris@10 1141
Chris@10 1142 Arrays @code{n}, @code{inembed}, and @code{onembed} are not used after
Chris@10 1143 this function returns. You can safely free or reuse them.
Chris@10 1144
Chris@10 1145 @c =========>
Chris@10 1146 @node Advanced Real-to-real Transforms, , Advanced Real-data DFTs, Advanced Interface
Chris@10 1147 @subsection Advanced Real-to-real Transforms
Chris@10 1148
Chris@10 1149 @example
Chris@10 1150 fftw_plan fftw_plan_many_r2r(int rank, const int *n, int howmany,
Chris@10 1151 double *in, const int *inembed,
Chris@10 1152 int istride, int idist,
Chris@10 1153 double *out, const int *onembed,
Chris@10 1154 int ostride, int odist,
Chris@10 1155 const fftw_r2r_kind *kind, unsigned flags);
Chris@10 1156 @end example
Chris@10 1157 @findex fftw_plan_many_r2r
Chris@10 1158
Chris@10 1159 Like @code{fftw_plan_many_dft}, this functions adds @code{howmany},
Chris@10 1160 @code{nembed}, @code{stride}, and @code{dist} parameters to the
Chris@10 1161 @code{fftw_plan_r2r} function, but otherwise behave the same as the
Chris@10 1162 basic interface. The interpretation of those additional parameters are
Chris@10 1163 the same as for @code{fftw_plan_many_dft}. (Of course, the
Chris@10 1164 @code{stride} and @code{dist} parameters are now in units of
Chris@10 1165 @code{double}, not @code{fftw_complex}.)
Chris@10 1166
Chris@10 1167 Arrays @code{n}, @code{inembed}, @code{onembed}, and @code{kind} are not
Chris@10 1168 used after this function returns. You can safely free or reuse them.
Chris@10 1169
Chris@10 1170 @c ------------------------------------------------------------
Chris@10 1171 @node Guru Interface, New-array Execute Functions, Advanced Interface, FFTW Reference
Chris@10 1172 @section Guru Interface
Chris@10 1173 @cindex guru interface
Chris@10 1174
Chris@10 1175 The ``guru'' interface to FFTW is intended to expose as much as possible
Chris@10 1176 of the flexibility in the underlying FFTW architecture. It allows one
Chris@10 1177 to compute multi-dimensional ``vectors'' (loops) of multi-dimensional
Chris@10 1178 transforms, where each vector/transform dimension has an independent
Chris@10 1179 size and stride.
Chris@10 1180 @cindex vector
Chris@10 1181 One can also use more general complex-number formats, e.g. separate real
Chris@10 1182 and imaginary arrays.
Chris@10 1183
Chris@10 1184 For those users who require the flexibility of the guru interface, it is
Chris@10 1185 important that they pay special attention to the documentation lest they
Chris@10 1186 shoot themselves in the foot.
Chris@10 1187
Chris@10 1188 @menu
Chris@10 1189 * Interleaved and split arrays::
Chris@10 1190 * Guru vector and transform sizes::
Chris@10 1191 * Guru Complex DFTs::
Chris@10 1192 * Guru Real-data DFTs::
Chris@10 1193 * Guru Real-to-real Transforms::
Chris@10 1194 * 64-bit Guru Interface::
Chris@10 1195 @end menu
Chris@10 1196
Chris@10 1197 @c =========>
Chris@10 1198 @node Interleaved and split arrays, Guru vector and transform sizes, Guru Interface, Guru Interface
Chris@10 1199 @subsection Interleaved and split arrays
Chris@10 1200
Chris@10 1201 The guru interface supports two representations of complex numbers,
Chris@10 1202 which we call the interleaved and the split format.
Chris@10 1203
Chris@10 1204 The @dfn{interleaved} format is the same one used by the basic and
Chris@10 1205 advanced interfaces, and it is documented in @ref{Complex numbers}.
Chris@10 1206 In the interleaved format, you provide pointers to the real part of a
Chris@10 1207 complex number, and the imaginary part understood to be stored in the
Chris@10 1208 next memory location.
Chris@10 1209 @cindex interleaved format
Chris@10 1210
Chris@10 1211
Chris@10 1212 The @dfn{split} format allows separate pointers to the real and
Chris@10 1213 imaginary parts of a complex array.
Chris@10 1214 @cindex split format
Chris@10 1215
Chris@10 1216
Chris@10 1217 Technically, the interleaved format is redundant, because you can
Chris@10 1218 always express an interleaved array in terms of a split array with
Chris@10 1219 appropriate pointers and strides. On the other hand, the interleaved
Chris@10 1220 format is simpler to use, and it is common in practice. Hence, FFTW
Chris@10 1221 supports it as a special case.
Chris@10 1222
Chris@10 1223 @c =========>
Chris@10 1224 @node Guru vector and transform sizes, Guru Complex DFTs, Interleaved and split arrays, Guru Interface
Chris@10 1225 @subsection Guru vector and transform sizes
Chris@10 1226
Chris@10 1227 The guru interface introduces one basic new data structure,
Chris@10 1228 @code{fftw_iodim}, that is used to specify sizes and strides for
Chris@10 1229 multi-dimensional transforms and vectors:
Chris@10 1230
Chris@10 1231 @example
Chris@10 1232 typedef struct @{
Chris@10 1233 int n;
Chris@10 1234 int is;
Chris@10 1235 int os;
Chris@10 1236 @} fftw_iodim;
Chris@10 1237 @end example
Chris@10 1238 @tindex fftw_iodim
Chris@10 1239
Chris@10 1240 Here, @code{n} is the size of the dimension, and @code{is} and @code{os}
Chris@10 1241 are the strides of that dimension for the input and output arrays. (The
Chris@10 1242 stride is the separation of consecutive elements along this dimension.)
Chris@10 1243
Chris@10 1244 The meaning of the stride parameter depends on the type of the array
Chris@10 1245 that the stride refers to. @emph{If the array is interleaved complex,
Chris@10 1246 strides are expressed in units of complex numbers
Chris@10 1247 (@code{fftw_complex}). If the array is split complex or real, strides
Chris@10 1248 are expressed in units of real numbers (@code{double}).} This
Chris@10 1249 convention is consistent with the usual pointer arithmetic in the C
Chris@10 1250 language. An interleaved array is denoted by a pointer @code{p} to
Chris@10 1251 @code{fftw_complex}, so that @code{p+1} points to the next complex
Chris@10 1252 number. Split arrays are denoted by pointers to @code{double}, in
Chris@10 1253 which case pointer arithmetic operates in units of
Chris@10 1254 @code{sizeof(double)}.
Chris@10 1255 @cindex stride
Chris@10 1256
Chris@10 1257
Chris@10 1258 The guru planner interfaces all take a (@code{rank}, @code{dims[rank]})
Chris@10 1259 pair describing the transform size, and a (@code{howmany_rank},
Chris@10 1260 @code{howmany_dims[howmany_rank]}) pair describing the ``vector'' size (a
Chris@10 1261 multi-dimensional loop of transforms to perform), where @code{dims} and
Chris@10 1262 @code{howmany_dims} are arrays of @code{fftw_iodim}.
Chris@10 1263
Chris@10 1264 For example, the @code{howmany} parameter in the advanced complex-DFT
Chris@10 1265 interface corresponds to @code{howmany_rank} = 1,
Chris@10 1266 @code{howmany_dims[0].n} = @code{howmany}, @code{howmany_dims[0].is} =
Chris@10 1267 @code{idist}, and @code{howmany_dims[0].os} = @code{odist}.
Chris@10 1268 @cindex howmany loop
Chris@10 1269 @cindex dist
Chris@10 1270 (To compute a single transform, you can just use @code{howmany_rank} = 0.)
Chris@10 1271
Chris@10 1272
Chris@10 1273 A row-major multidimensional array with dimensions @code{n[rank]}
Chris@10 1274 (@pxref{Row-major Format}) corresponds to @code{dims[i].n} =
Chris@10 1275 @code{n[i]} and the recurrence @code{dims[i].is} = @code{n[i+1] *
Chris@10 1276 dims[i+1].is} (similarly for @code{os}). The stride of the last
Chris@10 1277 (@code{i=rank-1}) dimension is the overall stride of the array.
Chris@10 1278 e.g. to be equivalent to the advanced complex-DFT interface, you would
Chris@10 1279 have @code{dims[rank-1].is} = @code{istride} and
Chris@10 1280 @code{dims[rank-1].os} = @code{ostride}.
Chris@10 1281 @cindex row-major
Chris@10 1282
Chris@10 1283
Chris@10 1284 In general, we only guarantee FFTW to return a non-@code{NULL} plan if
Chris@10 1285 the vector and transform dimensions correspond to a set of distinct
Chris@10 1286 indices, and for in-place transforms the input/output strides should
Chris@10 1287 be the same.
Chris@10 1288
Chris@10 1289 @c =========>
Chris@10 1290 @node Guru Complex DFTs, Guru Real-data DFTs, Guru vector and transform sizes, Guru Interface
Chris@10 1291 @subsection Guru Complex DFTs
Chris@10 1292
Chris@10 1293 @example
Chris@10 1294 fftw_plan fftw_plan_guru_dft(
Chris@10 1295 int rank, const fftw_iodim *dims,
Chris@10 1296 int howmany_rank, const fftw_iodim *howmany_dims,
Chris@10 1297 fftw_complex *in, fftw_complex *out,
Chris@10 1298 int sign, unsigned flags);
Chris@10 1299
Chris@10 1300 fftw_plan fftw_plan_guru_split_dft(
Chris@10 1301 int rank, const fftw_iodim *dims,
Chris@10 1302 int howmany_rank, const fftw_iodim *howmany_dims,
Chris@10 1303 double *ri, double *ii, double *ro, double *io,
Chris@10 1304 unsigned flags);
Chris@10 1305 @end example
Chris@10 1306 @findex fftw_plan_guru_dft
Chris@10 1307 @findex fftw_plan_guru_split_dft
Chris@10 1308
Chris@10 1309 These two functions plan a complex-data, multi-dimensional DFT
Chris@10 1310 for the interleaved and split format, respectively.
Chris@10 1311 Transform dimensions are given by (@code{rank}, @code{dims}) over a
Chris@10 1312 multi-dimensional vector (loop) of dimensions (@code{howmany_rank},
Chris@10 1313 @code{howmany_dims}). @code{dims} and @code{howmany_dims} should point
Chris@10 1314 to @code{fftw_iodim} arrays of length @code{rank} and
Chris@10 1315 @code{howmany_rank}, respectively.
Chris@10 1316
Chris@10 1317 @cindex flags
Chris@10 1318 @code{flags} is a bitwise OR (@samp{|}) of zero or more planner flags,
Chris@10 1319 as defined in @ref{Planner Flags}.
Chris@10 1320
Chris@10 1321 In the @code{fftw_plan_guru_dft} function, the pointers @code{in} and
Chris@10 1322 @code{out} point to the interleaved input and output arrays,
Chris@10 1323 respectively. The sign can be either @math{-1} (=
Chris@10 1324 @code{FFTW_FORWARD}) or @math{+1} (= @code{FFTW_BACKWARD}). If the
Chris@10 1325 pointers are equal, the transform is in-place.
Chris@10 1326
Chris@10 1327 In the @code{fftw_plan_guru_split_dft} function,
Chris@10 1328 @code{ri} and @code{ii} point to the real and imaginary input arrays,
Chris@10 1329 and @code{ro} and @code{io} point to the real and imaginary output
Chris@10 1330 arrays. The input and output pointers may be the same, indicating an
Chris@10 1331 in-place transform. For example, for @code{fftw_complex} pointers
Chris@10 1332 @code{in} and @code{out}, the corresponding parameters are:
Chris@10 1333
Chris@10 1334 @example
Chris@10 1335 ri = (double *) in;
Chris@10 1336 ii = (double *) in + 1;
Chris@10 1337 ro = (double *) out;
Chris@10 1338 io = (double *) out + 1;
Chris@10 1339 @end example
Chris@10 1340
Chris@10 1341 Because @code{fftw_plan_guru_split_dft} accepts split arrays, strides
Chris@10 1342 are expressed in units of @code{double}. For a contiguous
Chris@10 1343 @code{fftw_complex} array, the overall stride of the transform should
Chris@10 1344 be 2, the distance between consecutive real parts or between
Chris@10 1345 consecutive imaginary parts; see @ref{Guru vector and transform
Chris@10 1346 sizes}. Note that the dimension strides are applied equally to the
Chris@10 1347 real and imaginary parts; real and imaginary arrays with different
Chris@10 1348 strides are not supported.
Chris@10 1349
Chris@10 1350 There is no @code{sign} parameter in @code{fftw_plan_guru_split_dft}.
Chris@10 1351 This function always plans for an @code{FFTW_FORWARD} transform. To
Chris@10 1352 plan for an @code{FFTW_BACKWARD} transform, you can exploit the
Chris@10 1353 identity that the backwards DFT is equal to the forwards DFT with the
Chris@10 1354 real and imaginary parts swapped. For example, in the case of the
Chris@10 1355 @code{fftw_complex} arrays above, the @code{FFTW_BACKWARD} transform
Chris@10 1356 is computed by the parameters:
Chris@10 1357
Chris@10 1358 @example
Chris@10 1359 ri = (double *) in + 1;
Chris@10 1360 ii = (double *) in;
Chris@10 1361 ro = (double *) out + 1;
Chris@10 1362 io = (double *) out;
Chris@10 1363 @end example
Chris@10 1364
Chris@10 1365 @c =========>
Chris@10 1366 @node Guru Real-data DFTs, Guru Real-to-real Transforms, Guru Complex DFTs, Guru Interface
Chris@10 1367 @subsection Guru Real-data DFTs
Chris@10 1368
Chris@10 1369 @example
Chris@10 1370 fftw_plan fftw_plan_guru_dft_r2c(
Chris@10 1371 int rank, const fftw_iodim *dims,
Chris@10 1372 int howmany_rank, const fftw_iodim *howmany_dims,
Chris@10 1373 double *in, fftw_complex *out,
Chris@10 1374 unsigned flags);
Chris@10 1375
Chris@10 1376 fftw_plan fftw_plan_guru_split_dft_r2c(
Chris@10 1377 int rank, const fftw_iodim *dims,
Chris@10 1378 int howmany_rank, const fftw_iodim *howmany_dims,
Chris@10 1379 double *in, double *ro, double *io,
Chris@10 1380 unsigned flags);
Chris@10 1381
Chris@10 1382 fftw_plan fftw_plan_guru_dft_c2r(
Chris@10 1383 int rank, const fftw_iodim *dims,
Chris@10 1384 int howmany_rank, const fftw_iodim *howmany_dims,
Chris@10 1385 fftw_complex *in, double *out,
Chris@10 1386 unsigned flags);
Chris@10 1387
Chris@10 1388 fftw_plan fftw_plan_guru_split_dft_c2r(
Chris@10 1389 int rank, const fftw_iodim *dims,
Chris@10 1390 int howmany_rank, const fftw_iodim *howmany_dims,
Chris@10 1391 double *ri, double *ii, double *out,
Chris@10 1392 unsigned flags);
Chris@10 1393 @end example
Chris@10 1394 @findex fftw_plan_guru_dft_r2c
Chris@10 1395 @findex fftw_plan_guru_split_dft_r2c
Chris@10 1396 @findex fftw_plan_guru_dft_c2r
Chris@10 1397 @findex fftw_plan_guru_split_dft_c2r
Chris@10 1398
Chris@10 1399 Plan a real-input (r2c) or real-output (c2r), multi-dimensional DFT with
Chris@10 1400 transform dimensions given by (@code{rank}, @code{dims}) over a
Chris@10 1401 multi-dimensional vector (loop) of dimensions (@code{howmany_rank},
Chris@10 1402 @code{howmany_dims}). @code{dims} and @code{howmany_dims} should point
Chris@10 1403 to @code{fftw_iodim} arrays of length @code{rank} and
Chris@10 1404 @code{howmany_rank}, respectively. As for the basic and advanced
Chris@10 1405 interfaces, an r2c transform is @code{FFTW_FORWARD} and a c2r transform
Chris@10 1406 is @code{FFTW_BACKWARD}.
Chris@10 1407
Chris@10 1408 The @emph{last} dimension of @code{dims} is interpreted specially:
Chris@10 1409 that dimension of the real array has size @code{dims[rank-1].n}, but
Chris@10 1410 that dimension of the complex array has size @code{dims[rank-1].n/2+1}
Chris@10 1411 (division rounded down). The strides, on the other hand, are taken to
Chris@10 1412 be exactly as specified. It is up to the user to specify the strides
Chris@10 1413 appropriately for the peculiar dimensions of the data, and we do not
Chris@10 1414 guarantee that the planner will succeed (return non-@code{NULL}) for
Chris@10 1415 any dimensions other than those described in @ref{Real-data DFT Array
Chris@10 1416 Format} and generalized in @ref{Advanced Real-data DFTs}. (That is,
Chris@10 1417 for an in-place transform, each individual dimension should be able to
Chris@10 1418 operate in place.)
Chris@10 1419 @cindex in-place
Chris@10 1420
Chris@10 1421
Chris@10 1422 @code{in} and @code{out} point to the input and output arrays for r2c
Chris@10 1423 and c2r transforms, respectively. For split arrays, @code{ri} and
Chris@10 1424 @code{ii} point to the real and imaginary input arrays for a c2r
Chris@10 1425 transform, and @code{ro} and @code{io} point to the real and imaginary
Chris@10 1426 output arrays for an r2c transform. @code{in} and @code{ro} or
Chris@10 1427 @code{ri} and @code{out} may be the same, indicating an in-place
Chris@10 1428 transform. (In-place transforms where @code{in} and @code{io} or
Chris@10 1429 @code{ii} and @code{out} are the same are not currently supported.)
Chris@10 1430
Chris@10 1431 @cindex flags
Chris@10 1432 @code{flags} is a bitwise OR (@samp{|}) of zero or more planner flags,
Chris@10 1433 as defined in @ref{Planner Flags}.
Chris@10 1434
Chris@10 1435 In-place transforms of rank greater than 1 are currently only
Chris@10 1436 supported for interleaved arrays. For split arrays, the planner will
Chris@10 1437 return @code{NULL}.
Chris@10 1438 @cindex in-place
Chris@10 1439
Chris@10 1440 @c =========>
Chris@10 1441 @node Guru Real-to-real Transforms, 64-bit Guru Interface, Guru Real-data DFTs, Guru Interface
Chris@10 1442 @subsection Guru Real-to-real Transforms
Chris@10 1443
Chris@10 1444 @example
Chris@10 1445 fftw_plan fftw_plan_guru_r2r(int rank, const fftw_iodim *dims,
Chris@10 1446 int howmany_rank,
Chris@10 1447 const fftw_iodim *howmany_dims,
Chris@10 1448 double *in, double *out,
Chris@10 1449 const fftw_r2r_kind *kind,
Chris@10 1450 unsigned flags);
Chris@10 1451 @end example
Chris@10 1452 @findex fftw_plan_guru_r2r
Chris@10 1453
Chris@10 1454 Plan a real-to-real (r2r) multi-dimensional @code{FFTW_FORWARD}
Chris@10 1455 transform with transform dimensions given by (@code{rank}, @code{dims})
Chris@10 1456 over a multi-dimensional vector (loop) of dimensions
Chris@10 1457 (@code{howmany_rank}, @code{howmany_dims}). @code{dims} and
Chris@10 1458 @code{howmany_dims} should point to @code{fftw_iodim} arrays of length
Chris@10 1459 @code{rank} and @code{howmany_rank}, respectively.
Chris@10 1460
Chris@10 1461 The transform kind of each dimension is given by the @code{kind}
Chris@10 1462 parameter, which should point to an array of length @code{rank}. Valid
Chris@10 1463 @code{fftw_r2r_kind} constants are given in @ref{Real-to-Real Transform
Chris@10 1464 Kinds}.
Chris@10 1465
Chris@10 1466 @code{in} and @code{out} point to the real input and output arrays; they
Chris@10 1467 may be the same, indicating an in-place transform.
Chris@10 1468
Chris@10 1469 @cindex flags
Chris@10 1470 @code{flags} is a bitwise OR (@samp{|}) of zero or more planner flags,
Chris@10 1471 as defined in @ref{Planner Flags}.
Chris@10 1472
Chris@10 1473 @c =========>
Chris@10 1474 @node 64-bit Guru Interface, , Guru Real-to-real Transforms, Guru Interface
Chris@10 1475 @subsection 64-bit Guru Interface
Chris@10 1476 @cindex 64-bit architecture
Chris@10 1477
Chris@10 1478 When compiled in 64-bit mode on a 64-bit architecture (where addresses
Chris@10 1479 are 64 bits wide), FFTW uses 64-bit quantities internally for all
Chris@10 1480 transform sizes, strides, and so on---you don't have to do anything
Chris@10 1481 special to exploit this. However, in the ordinary FFTW interfaces,
Chris@10 1482 you specify the transform size by an @code{int} quantity, which is
Chris@10 1483 normally only 32 bits wide. This means that, even though FFTW is
Chris@10 1484 using 64-bit sizes internally, you cannot specify a single transform
Chris@10 1485 dimension larger than
Chris@10 1486 @ifinfo
Chris@10 1487 2^31-1
Chris@10 1488 @end ifinfo
Chris@10 1489 @html
Chris@10 1490 2<sup><small>31</small></sup>&minus;1
Chris@10 1491 @end html
Chris@10 1492 @tex
Chris@10 1493 $2^31-1$
Chris@10 1494 @end tex
Chris@10 1495 numbers.
Chris@10 1496
Chris@10 1497 We expect that few users will require transforms larger than this, but,
Chris@10 1498 for those who do, we provide a 64-bit version of the guru interface in
Chris@10 1499 which all sizes are specified as integers of type @code{ptrdiff_t}
Chris@10 1500 instead of @code{int}. (@code{ptrdiff_t} is a signed integer type
Chris@10 1501 defined by the C standard to be wide enough to represent address
Chris@10 1502 differences, and thus must be at least 64 bits wide on a 64-bit
Chris@10 1503 machine.) We stress that there is @emph{no performance advantage} to
Chris@10 1504 using this interface---the same internal FFTW code is employed
Chris@10 1505 regardless---and it is only necessary if you want to specify very
Chris@10 1506 large transform sizes.
Chris@10 1507 @tindex ptrdiff_t
Chris@10 1508
Chris@10 1509
Chris@10 1510 In particular, the 64-bit guru interface is a set of planner routines
Chris@10 1511 that are exactly the same as the guru planner routines, except that
Chris@10 1512 they are named with @samp{guru64} instead of @samp{guru} and they take
Chris@10 1513 arguments of type @code{fftw_iodim64} instead of @code{fftw_iodim}.
Chris@10 1514 For example, instead of @code{fftw_plan_guru_dft}, we have
Chris@10 1515 @code{fftw_plan_guru64_dft}.
Chris@10 1516
Chris@10 1517 @example
Chris@10 1518 fftw_plan fftw_plan_guru64_dft(
Chris@10 1519 int rank, const fftw_iodim64 *dims,
Chris@10 1520 int howmany_rank, const fftw_iodim64 *howmany_dims,
Chris@10 1521 fftw_complex *in, fftw_complex *out,
Chris@10 1522 int sign, unsigned flags);
Chris@10 1523 @end example
Chris@10 1524 @findex fftw_plan_guru64_dft
Chris@10 1525
Chris@10 1526 The @code{fftw_iodim64} type is similar to @code{fftw_iodim}, with the
Chris@10 1527 same interpretation, except that it uses type @code{ptrdiff_t} instead
Chris@10 1528 of type @code{int}.
Chris@10 1529
Chris@10 1530 @example
Chris@10 1531 typedef struct @{
Chris@10 1532 ptrdiff_t n;
Chris@10 1533 ptrdiff_t is;
Chris@10 1534 ptrdiff_t os;
Chris@10 1535 @} fftw_iodim64;
Chris@10 1536 @end example
Chris@10 1537 @tindex fftw_iodim64
Chris@10 1538
Chris@10 1539 Every other @samp{fftw_plan_guru} function also has a
Chris@10 1540 @samp{fftw_plan_guru64} equivalent, but we do not repeat their
Chris@10 1541 documentation here since they are identical to the 32-bit versions
Chris@10 1542 except as noted above.
Chris@10 1543
Chris@10 1544 @c -----------------------------------------------------------
Chris@10 1545 @node New-array Execute Functions, Wisdom, Guru Interface, FFTW Reference
Chris@10 1546 @section New-array Execute Functions
Chris@10 1547 @cindex execute
Chris@10 1548 @cindex new-array execution
Chris@10 1549
Chris@10 1550 Normally, one executes a plan for the arrays with which the plan was
Chris@10 1551 created, by calling @code{fftw_execute(plan)} as described in @ref{Using
Chris@10 1552 Plans}.
Chris@10 1553 @findex fftw_execute
Chris@10 1554 However, it is possible for sophisticated users to apply a given plan
Chris@10 1555 to a @emph{different} array using the ``new-array execute'' functions
Chris@10 1556 detailed below, provided that the following conditions are met:
Chris@10 1557
Chris@10 1558 @itemize @bullet
Chris@10 1559
Chris@10 1560 @item
Chris@10 1561 The array size, strides, etcetera are the same (since those are set by
Chris@10 1562 the plan).
Chris@10 1563
Chris@10 1564 @item
Chris@10 1565 The input and output arrays are the same (in-place) or different
Chris@10 1566 (out-of-place) if the plan was originally created to be in-place or
Chris@10 1567 out-of-place, respectively.
Chris@10 1568
Chris@10 1569 @item
Chris@10 1570 For split arrays, the separations between the real and imaginary
Chris@10 1571 parts, @code{ii-ri} and @code{io-ro}, are the same as they were for
Chris@10 1572 the input and output arrays when the plan was created. (This
Chris@10 1573 condition is automatically satisfied for interleaved arrays.)
Chris@10 1574
Chris@10 1575 @item
Chris@10 1576 The @dfn{alignment} of the new input/output arrays is the same as that
Chris@10 1577 of the input/output arrays when the plan was created, unless the plan
Chris@10 1578 was created with the @code{FFTW_UNALIGNED} flag.
Chris@10 1579 @ctindex FFTW_UNALIGNED
Chris@10 1580 Here, the alignment is a platform-dependent quantity (for example, it is
Chris@10 1581 the address modulo 16 if SSE SIMD instructions are used, but the address
Chris@10 1582 modulo 4 for non-SIMD single-precision FFTW on the same machine). In
Chris@10 1583 general, only arrays allocated with @code{fftw_malloc} are guaranteed to
Chris@10 1584 be equally aligned (@pxref{SIMD alignment and fftw_malloc}).
Chris@10 1585
Chris@10 1586 @end itemize
Chris@10 1587
Chris@10 1588 @cindex alignment
Chris@10 1589 The alignment issue is especially critical, because if you don't use
Chris@10 1590 @code{fftw_malloc} then you may have little control over the alignment
Chris@10 1591 of arrays in memory. For example, neither the C++ @code{new} function
Chris@10 1592 nor the Fortran @code{allocate} statement provide strong enough
Chris@10 1593 guarantees about data alignment. If you don't use @code{fftw_malloc},
Chris@10 1594 therefore, you probably have to use @code{FFTW_UNALIGNED} (which
Chris@10 1595 disables most SIMD support). If possible, it is probably better for
Chris@10 1596 you to simply create multiple plans (creating a new plan is quick once
Chris@10 1597 one exists for a given size), or better yet re-use the same array for
Chris@10 1598 your transforms.
Chris@10 1599
Chris@10 1600 If you are tempted to use the new-array execute interface because you
Chris@10 1601 want to transform a known bunch of arrays of the same size, you should
Chris@10 1602 probably go use the advanced interface instead (@pxref{Advanced
Chris@10 1603 Interface})).
Chris@10 1604
Chris@10 1605 The new-array execute functions are:
Chris@10 1606
Chris@10 1607 @example
Chris@10 1608 void fftw_execute_dft(
Chris@10 1609 const fftw_plan p,
Chris@10 1610 fftw_complex *in, fftw_complex *out);
Chris@10 1611
Chris@10 1612 void fftw_execute_split_dft(
Chris@10 1613 const fftw_plan p,
Chris@10 1614 double *ri, double *ii, double *ro, double *io);
Chris@10 1615
Chris@10 1616 void fftw_execute_dft_r2c(
Chris@10 1617 const fftw_plan p,
Chris@10 1618 double *in, fftw_complex *out);
Chris@10 1619
Chris@10 1620 void fftw_execute_split_dft_r2c(
Chris@10 1621 const fftw_plan p,
Chris@10 1622 double *in, double *ro, double *io);
Chris@10 1623
Chris@10 1624 void fftw_execute_dft_c2r(
Chris@10 1625 const fftw_plan p,
Chris@10 1626 fftw_complex *in, double *out);
Chris@10 1627
Chris@10 1628 void fftw_execute_split_dft_c2r(
Chris@10 1629 const fftw_plan p,
Chris@10 1630 double *ri, double *ii, double *out);
Chris@10 1631
Chris@10 1632 void fftw_execute_r2r(
Chris@10 1633 const fftw_plan p,
Chris@10 1634 double *in, double *out);
Chris@10 1635 @end example
Chris@10 1636 @findex fftw_execute_dft
Chris@10 1637 @findex fftw_execute_split_dft
Chris@10 1638 @findex fftw_execute_dft_r2c
Chris@10 1639 @findex fftw_execute_split_dft_r2c
Chris@10 1640 @findex fftw_execute_dft_c2r
Chris@10 1641 @findex fftw_execute_split_dft_c2r
Chris@10 1642 @findex fftw_execute_r2r
Chris@10 1643
Chris@10 1644 These execute the @code{plan} to compute the corresponding transform on
Chris@10 1645 the input/output arrays specified by the subsequent arguments. The
Chris@10 1646 input/output array arguments have the same meanings as the ones passed
Chris@10 1647 to the guru planner routines in the preceding sections. The @code{plan}
Chris@10 1648 is not modified, and these routines can be called as many times as
Chris@10 1649 desired, or intermixed with calls to the ordinary @code{fftw_execute}.
Chris@10 1650
Chris@10 1651 The @code{plan} @emph{must} have been created for the transform type
Chris@10 1652 corresponding to the execute function, e.g. it must be a complex-DFT
Chris@10 1653 plan for @code{fftw_execute_dft}. Any of the planner routines for that
Chris@10 1654 transform type, from the basic to the guru interface, could have been
Chris@10 1655 used to create the plan, however.
Chris@10 1656
Chris@10 1657 @c ------------------------------------------------------------
Chris@10 1658 @node Wisdom, What FFTW Really Computes, New-array Execute Functions, FFTW Reference
Chris@10 1659 @section Wisdom
Chris@10 1660 @cindex wisdom
Chris@10 1661 @cindex saving plans to disk
Chris@10 1662
Chris@10 1663 This section documents the FFTW mechanism for saving and restoring
Chris@10 1664 plans from disk. This mechanism is called @dfn{wisdom}.
Chris@10 1665
Chris@10 1666 @menu
Chris@10 1667 * Wisdom Export::
Chris@10 1668 * Wisdom Import::
Chris@10 1669 * Forgetting Wisdom::
Chris@10 1670 * Wisdom Utilities::
Chris@10 1671 @end menu
Chris@10 1672
Chris@10 1673 @c =========>
Chris@10 1674 @node Wisdom Export, Wisdom Import, Wisdom, Wisdom
Chris@10 1675 @subsection Wisdom Export
Chris@10 1676
Chris@10 1677 @example
Chris@10 1678 int fftw_export_wisdom_to_filename(const char *filename);
Chris@10 1679 void fftw_export_wisdom_to_file(FILE *output_file);
Chris@10 1680 char *fftw_export_wisdom_to_string(void);
Chris@10 1681 void fftw_export_wisdom(void (*write_char)(char c, void *), void *data);
Chris@10 1682 @end example
Chris@10 1683 @findex fftw_export_wisdom
Chris@10 1684 @findex fftw_export_wisdom_to_filename
Chris@10 1685 @findex fftw_export_wisdom_to_file
Chris@10 1686 @findex fftw_export_wisdom_to_string
Chris@10 1687
Chris@10 1688 These functions allow you to export all currently accumulated wisdom
Chris@10 1689 in a form from which it can be later imported and restored, even
Chris@10 1690 during a separate run of the program. (@xref{Words of Wisdom-Saving
Chris@10 1691 Plans}.) The current store of wisdom is not affected by calling any
Chris@10 1692 of these routines.
Chris@10 1693
Chris@10 1694 @code{fftw_export_wisdom} exports the wisdom to any output
Chris@10 1695 medium, as specified by the callback function
Chris@10 1696 @code{write_char}. @code{write_char} is a @code{putc}-like function that
Chris@10 1697 writes the character @code{c} to some output; its second parameter is
Chris@10 1698 the @code{data} pointer passed to @code{fftw_export_wisdom}. For
Chris@10 1699 convenience, the following three ``wrapper'' routines are provided:
Chris@10 1700
Chris@10 1701 @code{fftw_export_wisdom_to_filename} writes wisdom to a file named
Chris@10 1702 @code{filename} (which is created or overwritten), returning @code{1}
Chris@10 1703 on success and @code{0} on failure. A lower-level function, which
Chris@10 1704 requires you to open and close the file yourself (e.g. if you want to
Chris@10 1705 write wisdom to a portion of a larger file) is
Chris@10 1706 @code{fftw_export_wisdom_to_file}. This writes the wisdom to the
Chris@10 1707 current position in @code{output_file}, which should be open with
Chris@10 1708 write permission; upon exit, the file remains open and is positioned
Chris@10 1709 at the end of the wisdom data.
Chris@10 1710
Chris@10 1711 @code{fftw_export_wisdom_to_string} returns a pointer to a
Chris@10 1712 @code{NULL}-terminated string holding the wisdom data. This string is
Chris@10 1713 dynamically allocated, and it is the responsibility of the caller to
Chris@10 1714 deallocate it with @code{free} when it is no longer needed.
Chris@10 1715
Chris@10 1716 All of these routines export the wisdom in the same format, which we
Chris@10 1717 will not document here except to say that it is LISP-like ASCII text
Chris@10 1718 that is insensitive to white space.
Chris@10 1719
Chris@10 1720 @c =========>
Chris@10 1721 @node Wisdom Import, Forgetting Wisdom, Wisdom Export, Wisdom
Chris@10 1722 @subsection Wisdom Import
Chris@10 1723
Chris@10 1724 @example
Chris@10 1725 int fftw_import_system_wisdom(void);
Chris@10 1726 int fftw_import_wisdom_from_filename(const char *filename);
Chris@10 1727 int fftw_import_wisdom_from_string(const char *input_string);
Chris@10 1728 int fftw_import_wisdom(int (*read_char)(void *), void *data);
Chris@10 1729 @end example
Chris@10 1730 @findex fftw_import_wisdom
Chris@10 1731 @findex fftw_import_system_wisdom
Chris@10 1732 @findex fftw_import_wisdom_from_filename
Chris@10 1733 @findex fftw_import_wisdom_from_file
Chris@10 1734 @findex fftw_import_wisdom_from_string
Chris@10 1735
Chris@10 1736 These functions import wisdom into a program from data stored by the
Chris@10 1737 @code{fftw_export_wisdom} functions above. (@xref{Words of
Chris@10 1738 Wisdom-Saving Plans}.) The imported wisdom replaces any wisdom
Chris@10 1739 already accumulated by the running program.
Chris@10 1740
Chris@10 1741 @code{fftw_import_wisdom} imports wisdom from any input medium, as
Chris@10 1742 specified by the callback function @code{read_char}. @code{read_char} is
Chris@10 1743 a @code{getc}-like function that returns the next character in the
Chris@10 1744 input; its parameter is the @code{data} pointer passed to
Chris@10 1745 @code{fftw_import_wisdom}. If the end of the input data is reached
Chris@10 1746 (which should never happen for valid data), @code{read_char} should
Chris@10 1747 return @code{EOF} (as defined in @code{<stdio.h>}). For convenience,
Chris@10 1748 the following three ``wrapper'' routines are provided:
Chris@10 1749
Chris@10 1750 @code{fftw_import_wisdom_from_filename} reads wisdom from a file named
Chris@10 1751 @code{filename}. A lower-level function, which requires you to open
Chris@10 1752 and close the file yourself (e.g. if you want to read wisdom from a
Chris@10 1753 portion of a larger file) is @code{fftw_import_wisdom_from_file}. This
Chris@10 1754 reads wisdom from the current position in @code{input_file} (which
Chris@10 1755 should be open with read permission); upon exit, the file remains
Chris@10 1756 open, but the position of the read pointer is unspecified.
Chris@10 1757
Chris@10 1758 @code{fftw_import_wisdom_from_string} reads wisdom from the
Chris@10 1759 @code{NULL}-terminated string @code{input_string}.
Chris@10 1760
Chris@10 1761 @code{fftw_import_system_wisdom} reads wisdom from an
Chris@10 1762 implementation-defined standard file (@code{/etc/fftw/wisdom} on Unix
Chris@10 1763 and GNU systems).
Chris@10 1764 @cindex wisdom, system-wide
Chris@10 1765
Chris@10 1766
Chris@10 1767 The return value of these import routines is @code{1} if the wisdom was
Chris@10 1768 read successfully and @code{0} otherwise. Note that, in all of these
Chris@10 1769 functions, any data in the input stream past the end of the wisdom data
Chris@10 1770 is simply ignored.
Chris@10 1771
Chris@10 1772 @c =========>
Chris@10 1773 @node Forgetting Wisdom, Wisdom Utilities, Wisdom Import, Wisdom
Chris@10 1774 @subsection Forgetting Wisdom
Chris@10 1775
Chris@10 1776 @example
Chris@10 1777 void fftw_forget_wisdom(void);
Chris@10 1778 @end example
Chris@10 1779 @findex fftw_forget_wisdom
Chris@10 1780
Chris@10 1781 Calling @code{fftw_forget_wisdom} causes all accumulated @code{wisdom}
Chris@10 1782 to be discarded and its associated memory to be freed. (New
Chris@10 1783 @code{wisdom} can still be gathered subsequently, however.)
Chris@10 1784
Chris@10 1785 @c =========>
Chris@10 1786 @node Wisdom Utilities, , Forgetting Wisdom, Wisdom
Chris@10 1787 @subsection Wisdom Utilities
Chris@10 1788
Chris@10 1789 FFTW includes two standalone utility programs that deal with wisdom. We
Chris@10 1790 merely summarize them here, since they come with their own @code{man}
Chris@10 1791 pages for Unix and GNU systems (with HTML versions on our web site).
Chris@10 1792
Chris@10 1793 The first program is @code{fftw-wisdom} (or @code{fftwf-wisdom} in
Chris@10 1794 single precision, etcetera), which can be used to create a wisdom file
Chris@10 1795 containing plans for any of the transform sizes and types supported by
Chris@10 1796 FFTW. It is preferable to create wisdom directly from your executable
Chris@10 1797 (@pxref{Caveats in Using Wisdom}), but this program is useful for
Chris@10 1798 creating global wisdom files for @code{fftw_import_system_wisdom}.
Chris@10 1799 @cindex fftw-wisdom utility
Chris@10 1800
Chris@10 1801
Chris@10 1802 The second program is @code{fftw-wisdom-to-conf}, which takes a wisdom
Chris@10 1803 file as input and produces a @dfn{configuration routine} as output. The
Chris@10 1804 latter is a C subroutine that you can compile and link into your
Chris@10 1805 program, replacing a routine of the same name in the FFTW library, that
Chris@10 1806 determines which parts of FFTW are callable by your program.
Chris@10 1807 @code{fftw-wisdom-to-conf} produces a configuration routine that links
Chris@10 1808 to only those parts of FFTW needed by the saved plans in the wisdom,
Chris@10 1809 greatly reducing the size of statically linked executables (which should
Chris@10 1810 only attempt to create plans corresponding to those in the wisdom,
Chris@10 1811 however).
Chris@10 1812 @cindex fftw-wisdom-to-conf utility
Chris@10 1813 @cindex configuration routines
Chris@10 1814
Chris@10 1815 @c ------------------------------------------------------------
Chris@10 1816 @node What FFTW Really Computes, , Wisdom, FFTW Reference
Chris@10 1817 @section What FFTW Really Computes
Chris@10 1818
Chris@10 1819 In this section, we provide precise mathematical definitions for the
Chris@10 1820 transforms that FFTW computes. These transform definitions are fairly
Chris@10 1821 standard, but some authors follow slightly different conventions for the
Chris@10 1822 normalization of the transform (the constant factor in front) and the
Chris@10 1823 sign of the complex exponent. We begin by presenting the
Chris@10 1824 one-dimensional (1d) transform definitions, and then give the
Chris@10 1825 straightforward extension to multi-dimensional transforms.
Chris@10 1826
Chris@10 1827 @menu
Chris@10 1828 * The 1d Discrete Fourier Transform (DFT)::
Chris@10 1829 * The 1d Real-data DFT::
Chris@10 1830 * 1d Real-even DFTs (DCTs)::
Chris@10 1831 * 1d Real-odd DFTs (DSTs)::
Chris@10 1832 * 1d Discrete Hartley Transforms (DHTs)::
Chris@10 1833 * Multi-dimensional Transforms::
Chris@10 1834 @end menu
Chris@10 1835
Chris@10 1836 @c =========>
Chris@10 1837 @node The 1d Discrete Fourier Transform (DFT), The 1d Real-data DFT, What FFTW Really Computes, What FFTW Really Computes
Chris@10 1838 @subsection The 1d Discrete Fourier Transform (DFT)
Chris@10 1839
Chris@10 1840 @cindex discrete Fourier transform
Chris@10 1841 @cindex DFT
Chris@10 1842 The forward (@code{FFTW_FORWARD}) discrete Fourier transform (DFT) of a
Chris@10 1843 1d complex array @math{X} of size @math{n} computes an array @math{Y},
Chris@10 1844 where:
Chris@10 1845 @tex
Chris@10 1846 $$
Chris@10 1847 Y_k = \sum_{j = 0}^{n - 1} X_j e^{-2\pi j k \sqrt{-1}/n} \ .
Chris@10 1848 $$
Chris@10 1849 @end tex
Chris@10 1850 @ifinfo
Chris@10 1851 @center Y[k] = sum for j = 0 to (n - 1) of X[j] * exp(-2 pi j k sqrt(-1)/n) .
Chris@10 1852 @end ifinfo
Chris@10 1853 @html
Chris@10 1854 <center><img src="equation-dft.png" align="top">.</center>
Chris@10 1855 @end html
Chris@10 1856 The backward (@code{FFTW_BACKWARD}) DFT computes:
Chris@10 1857 @tex
Chris@10 1858 $$
Chris@10 1859 Y_k = \sum_{j = 0}^{n - 1} X_j e^{2\pi j k \sqrt{-1}/n} \ .
Chris@10 1860 $$
Chris@10 1861 @end tex
Chris@10 1862 @ifinfo
Chris@10 1863 @center Y[k] = sum for j = 0 to (n - 1) of X[j] * exp(2 pi j k sqrt(-1)/n) .
Chris@10 1864 @end ifinfo
Chris@10 1865 @html
Chris@10 1866 <center><img src="equation-idft.png" align="top">.</center>
Chris@10 1867 @end html
Chris@10 1868
Chris@10 1869 @cindex normalization
Chris@10 1870 FFTW computes an unnormalized transform, in that there is no coefficient
Chris@10 1871 in front of the summation in the DFT. In other words, applying the
Chris@10 1872 forward and then the backward transform will multiply the input by
Chris@10 1873 @math{n}.
Chris@10 1874
Chris@10 1875 @cindex frequency
Chris@10 1876 From above, an @code{FFTW_FORWARD} transform corresponds to a sign of
Chris@10 1877 @math{-1} in the exponent of the DFT. Note also that we use the
Chris@10 1878 standard ``in-order'' output ordering---the @math{k}-th output
Chris@10 1879 corresponds to the frequency @math{k/n} (or @math{k/T}, where @math{T}
Chris@10 1880 is your total sampling period). For those who like to think in terms of
Chris@10 1881 positive and negative frequencies, this means that the positive
Chris@10 1882 frequencies are stored in the first half of the output and the negative
Chris@10 1883 frequencies are stored in backwards order in the second half of the
Chris@10 1884 output. (The frequency @math{-k/n} is the same as the frequency
Chris@10 1885 @math{(n-k)/n}.)
Chris@10 1886
Chris@10 1887 @c =========>
Chris@10 1888 @node The 1d Real-data DFT, 1d Real-even DFTs (DCTs), The 1d Discrete Fourier Transform (DFT), What FFTW Really Computes
Chris@10 1889 @subsection The 1d Real-data DFT
Chris@10 1890
Chris@10 1891 The real-input (r2c) DFT in FFTW computes the @emph{forward} transform
Chris@10 1892 @math{Y} of the size @code{n} real array @math{X}, exactly as defined
Chris@10 1893 above, i.e.
Chris@10 1894 @tex
Chris@10 1895 $$
Chris@10 1896 Y_k = \sum_{j = 0}^{n - 1} X_j e^{-2\pi j k \sqrt{-1}/n} \ .
Chris@10 1897 $$
Chris@10 1898 @end tex
Chris@10 1899 @ifinfo
Chris@10 1900 @center Y[k] = sum for j = 0 to (n - 1) of X[j] * exp(-2 pi j k sqrt(-1)/n) .
Chris@10 1901 @end ifinfo
Chris@10 1902 @html
Chris@10 1903 <center><img src="equation-dft.png" align="top">.</center>
Chris@10 1904 @end html
Chris@10 1905 This output array @math{Y} can easily be shown to possess the
Chris@10 1906 ``Hermitian'' symmetry
Chris@10 1907 @cindex Hermitian
Chris@10 1908 @tex
Chris@10 1909 $Y_k = Y_{n-k}^*$,
Chris@10 1910 @end tex
Chris@10 1911 @ifinfo
Chris@10 1912 Y[k] = Y[n-k]*,
Chris@10 1913 @end ifinfo
Chris@10 1914 @html
Chris@10 1915 <i>Y<sub>k</sub> = Y<sub>n-k</sub></i><sup>*</sup>,
Chris@10 1916 @end html
Chris@10 1917 where we take @math{Y} to be periodic so that
Chris@10 1918 @tex
Chris@10 1919 $Y_n = Y_0$.
Chris@10 1920 @end tex
Chris@10 1921 @ifinfo
Chris@10 1922 Y[n] = Y[0].
Chris@10 1923 @end ifinfo
Chris@10 1924 @html
Chris@10 1925 <i>Y<sub>n</sub> = Y</i><sub>0</sub>.
Chris@10 1926 @end html
Chris@10 1927
Chris@10 1928 As a result of this symmetry, half of the output @math{Y} is redundant
Chris@10 1929 (being the complex conjugate of the other half), and so the 1d r2c
Chris@10 1930 transforms only output elements @math{0}@dots{}@math{n/2} of @math{Y}
Chris@10 1931 (@math{n/2+1} complex numbers), where the division by @math{2} is
Chris@10 1932 rounded down.
Chris@10 1933
Chris@10 1934 Moreover, the Hermitian symmetry implies that
Chris@10 1935 @tex
Chris@10 1936 $Y_0$
Chris@10 1937 @end tex
Chris@10 1938 @ifinfo
Chris@10 1939 Y[0]
Chris@10 1940 @end ifinfo
Chris@10 1941 @html
Chris@10 1942 <i>Y</i><sub>0</sub>
Chris@10 1943 @end html
Chris@10 1944 and, if @math{n} is even, the
Chris@10 1945 @tex
Chris@10 1946 $Y_{n/2}$
Chris@10 1947 @end tex
Chris@10 1948 @ifinfo
Chris@10 1949 Y[n/2]
Chris@10 1950 @end ifinfo
Chris@10 1951 @html
Chris@10 1952 <i>Y</i><sub><i>n</i>/2</sub>
Chris@10 1953 @end html
Chris@10 1954 element, are purely real. So, for the @code{R2HC} r2r transform, these
Chris@10 1955 elements are not stored in the halfcomplex output format.
Chris@10 1956 @cindex r2r
Chris@10 1957 @ctindex R2HC
Chris@10 1958 @cindex halfcomplex format
Chris@10 1959
Chris@10 1960
Chris@10 1961 The c2r and @code{H2RC} r2r transforms compute the backward DFT of the
Chris@10 1962 @emph{complex} array @math{X} with Hermitian symmetry, stored in the
Chris@10 1963 r2c/@code{R2HC} output formats, respectively, where the backward
Chris@10 1964 transform is defined exactly as for the complex case:
Chris@10 1965 @tex
Chris@10 1966 $$
Chris@10 1967 Y_k = \sum_{j = 0}^{n - 1} X_j e^{2\pi j k \sqrt{-1}/n} \ .
Chris@10 1968 $$
Chris@10 1969 @end tex
Chris@10 1970 @ifinfo
Chris@10 1971 @center Y[k] = sum for j = 0 to (n - 1) of X[j] * exp(2 pi j k sqrt(-1)/n) .
Chris@10 1972 @end ifinfo
Chris@10 1973 @html
Chris@10 1974 <center><img src="equation-idft.png" align="top">.</center>
Chris@10 1975 @end html
Chris@10 1976 The outputs @code{Y} of this transform can easily be seen to be purely
Chris@10 1977 real, and are stored as an array of real numbers.
Chris@10 1978
Chris@10 1979 @cindex normalization
Chris@10 1980 Like FFTW's complex DFT, these transforms are unnormalized. In other
Chris@10 1981 words, applying the real-to-complex (forward) and then the
Chris@10 1982 complex-to-real (backward) transform will multiply the input by
Chris@10 1983 @math{n}.
Chris@10 1984
Chris@10 1985 @c =========>
Chris@10 1986 @node 1d Real-even DFTs (DCTs), 1d Real-odd DFTs (DSTs), The 1d Real-data DFT, What FFTW Really Computes
Chris@10 1987 @subsection 1d Real-even DFTs (DCTs)
Chris@10 1988
Chris@10 1989 The Real-even symmetry DFTs in FFTW are exactly equivalent to the unnormalized
Chris@10 1990 forward (and backward) DFTs as defined above, where the input array
Chris@10 1991 @math{X} of length @math{N} is purely real and is also @dfn{even} symmetry. In
Chris@10 1992 this case, the output array is likewise real and even symmetry.
Chris@10 1993 @cindex real-even DFT
Chris@10 1994 @cindex REDFT
Chris@10 1995
Chris@10 1996
Chris@10 1997 @ctindex REDFT00
Chris@10 1998 For the case of @code{REDFT00}, this even symmetry means that
Chris@10 1999 @tex
Chris@10 2000 $X_j = X_{N-j}$,
Chris@10 2001 @end tex
Chris@10 2002 @ifinfo
Chris@10 2003 X[j] = X[N-j],
Chris@10 2004 @end ifinfo
Chris@10 2005 @html
Chris@10 2006 <i>X<sub>j</sub> = X<sub>N-j</sub></i>,
Chris@10 2007 @end html
Chris@10 2008 where we take @math{X} to be periodic so that
Chris@10 2009 @tex
Chris@10 2010 $X_N = X_0$.
Chris@10 2011 @end tex
Chris@10 2012 @ifinfo
Chris@10 2013 X[N] = X[0].
Chris@10 2014 @end ifinfo
Chris@10 2015 @html
Chris@10 2016 <i>X<sub>N</sub> = X</i><sub>0</sub>.
Chris@10 2017 @end html
Chris@10 2018 Because of this redundancy, only the first @math{n} real numbers are
Chris@10 2019 actually stored, where @math{N = 2(n-1)}.
Chris@10 2020
Chris@10 2021 The proper definition of even symmetry for @code{REDFT10},
Chris@10 2022 @code{REDFT01}, and @code{REDFT11} transforms is somewhat more intricate
Chris@10 2023 because of the shifts by @math{1/2} of the input and/or output, although
Chris@10 2024 the corresponding boundary conditions are given in @ref{Real even/odd
Chris@10 2025 DFTs (cosine/sine transforms)}. Because of the even symmetry, however,
Chris@10 2026 the sine terms in the DFT all cancel and the remaining cosine terms are
Chris@10 2027 written explicitly below. This formulation often leads people to call
Chris@10 2028 such a transform a @dfn{discrete cosine transform} (DCT), although it is
Chris@10 2029 really just a special case of the DFT.
Chris@10 2030 @cindex discrete cosine transform
Chris@10 2031 @cindex DCT
Chris@10 2032
Chris@10 2033
Chris@10 2034 In each of the definitions below, we transform a real array @math{X} of
Chris@10 2035 length @math{n} to a real array @math{Y} of length @math{n}:
Chris@10 2036
Chris@10 2037 @subsubheading REDFT00 (DCT-I)
Chris@10 2038 @ctindex REDFT00
Chris@10 2039 An @code{REDFT00} transform (type-I DCT) in FFTW is defined by:
Chris@10 2040 @tex
Chris@10 2041 $$
Chris@10 2042 Y_k = X_0 + (-1)^k X_{n-1}
Chris@10 2043 + 2 \sum_{j=1}^{n-2} X_j \cos [ \pi j k / (n-1)].
Chris@10 2044 $$
Chris@10 2045 @end tex
Chris@10 2046 @ifinfo
Chris@10 2047 Y[k] = X[0] + (-1)^k X[n-1] + 2 (sum for j = 1 to n-2 of X[j] cos(pi jk /(n-1))).
Chris@10 2048 @end ifinfo
Chris@10 2049 @html
Chris@10 2050 <center><img src="equation-redft00.png" align="top">.</center>
Chris@10 2051 @end html
Chris@10 2052 Note that this transform is not defined for @math{n=1}. For @math{n=2},
Chris@10 2053 the summation term above is dropped as you might expect.
Chris@10 2054
Chris@10 2055 @subsubheading REDFT10 (DCT-II)
Chris@10 2056 @ctindex REDFT10
Chris@10 2057 An @code{REDFT10} transform (type-II DCT, sometimes called ``the'' DCT) in FFTW is defined by:
Chris@10 2058 @tex
Chris@10 2059 $$
Chris@10 2060 Y_k = 2 \sum_{j=0}^{n-1} X_j \cos [\pi (j+1/2) k / n].
Chris@10 2061 $$
Chris@10 2062 @end tex
Chris@10 2063 @ifinfo
Chris@10 2064 Y[k] = 2 (sum for j = 0 to n-1 of X[j] cos(pi (j+1/2) k / n)).
Chris@10 2065 @end ifinfo
Chris@10 2066 @html
Chris@10 2067 <center><img src="equation-redft10.png" align="top">.</center>
Chris@10 2068 @end html
Chris@10 2069
Chris@10 2070 @subsubheading REDFT01 (DCT-III)
Chris@10 2071 @ctindex REDFT01
Chris@10 2072 An @code{REDFT01} transform (type-III DCT) in FFTW is defined by:
Chris@10 2073 @tex
Chris@10 2074 $$
Chris@10 2075 Y_k = X_0 + 2 \sum_{j=1}^{n-1} X_j \cos [\pi j (k+1/2) / n].
Chris@10 2076 $$
Chris@10 2077 @end tex
Chris@10 2078 @ifinfo
Chris@10 2079 Y[k] = X[0] + 2 (sum for j = 1 to n-1 of X[j] cos(pi j (k+1/2) / n)).
Chris@10 2080 @end ifinfo
Chris@10 2081 @html
Chris@10 2082 <center><img src="equation-redft01.png" align="top">.</center>
Chris@10 2083 @end html
Chris@10 2084 In the case of @math{n=1}, this reduces to
Chris@10 2085 @tex
Chris@10 2086 $Y_0 = X_0$.
Chris@10 2087 @end tex
Chris@10 2088 @ifinfo
Chris@10 2089 Y[0] = X[0].
Chris@10 2090 @end ifinfo
Chris@10 2091 @html
Chris@10 2092 <i>Y</i><sub>0</sub> = <i>X</i><sub>0</sub>.
Chris@10 2093 @end html
Chris@10 2094 Up to a scale factor (see below), this is the inverse of @code{REDFT10} (``the'' DCT), and so the @code{REDFT01} (DCT-III) is sometimes called the ``IDCT''.
Chris@10 2095 @cindex IDCT
Chris@10 2096
Chris@10 2097 @subsubheading REDFT11 (DCT-IV)
Chris@10 2098 @ctindex REDFT11
Chris@10 2099 An @code{REDFT11} transform (type-IV DCT) in FFTW is defined by:
Chris@10 2100 @tex
Chris@10 2101 $$
Chris@10 2102 Y_k = 2 \sum_{j=0}^{n-1} X_j \cos [\pi (j+1/2) (k+1/2) / n].
Chris@10 2103 $$
Chris@10 2104 @end tex
Chris@10 2105 @ifinfo
Chris@10 2106 Y[k] = 2 (sum for j = 0 to n-1 of X[j] cos(pi (j+1/2) (k+1/2) / n)).
Chris@10 2107 @end ifinfo
Chris@10 2108 @html
Chris@10 2109 <center><img src="equation-redft11.png" align="top">.</center>
Chris@10 2110 @end html
Chris@10 2111
Chris@10 2112 @subsubheading Inverses and Normalization
Chris@10 2113
Chris@10 2114 These definitions correspond directly to the unnormalized DFTs used
Chris@10 2115 elsewhere in FFTW (hence the factors of @math{2} in front of the
Chris@10 2116 summations). The unnormalized inverse of @code{REDFT00} is
Chris@10 2117 @code{REDFT00}, of @code{REDFT10} is @code{REDFT01} and vice versa, and
Chris@10 2118 of @code{REDFT11} is @code{REDFT11}. Each unnormalized inverse results
Chris@10 2119 in the original array multiplied by @math{N}, where @math{N} is the
Chris@10 2120 @emph{logical} DFT size. For @code{REDFT00}, @math{N=2(n-1)} (note that
Chris@10 2121 @math{n=1} is not defined); otherwise, @math{N=2n}.
Chris@10 2122 @cindex normalization
Chris@10 2123
Chris@10 2124
Chris@10 2125 In defining the discrete cosine transform, some authors also include
Chris@10 2126 additional factors of
Chris@10 2127 @ifinfo
Chris@10 2128 sqrt(2)
Chris@10 2129 @end ifinfo
Chris@10 2130 @html
Chris@10 2131 &radic;2
Chris@10 2132 @end html
Chris@10 2133 @tex
Chris@10 2134 $\sqrt{2}$
Chris@10 2135 @end tex
Chris@10 2136 (or its inverse) multiplying selected inputs and/or outputs. This is a
Chris@10 2137 mostly cosmetic change that makes the transform orthogonal, but
Chris@10 2138 sacrifices the direct equivalence to a symmetric DFT.
Chris@10 2139
Chris@10 2140 @c =========>
Chris@10 2141 @node 1d Real-odd DFTs (DSTs), 1d Discrete Hartley Transforms (DHTs), 1d Real-even DFTs (DCTs), What FFTW Really Computes
Chris@10 2142 @subsection 1d Real-odd DFTs (DSTs)
Chris@10 2143
Chris@10 2144 The Real-odd symmetry DFTs in FFTW are exactly equivalent to the unnormalized
Chris@10 2145 forward (and backward) DFTs as defined above, where the input array
Chris@10 2146 @math{X} of length @math{N} is purely real and is also @dfn{odd} symmetry. In
Chris@10 2147 this case, the output is odd symmetry and purely imaginary.
Chris@10 2148 @cindex real-odd DFT
Chris@10 2149 @cindex RODFT
Chris@10 2150
Chris@10 2151
Chris@10 2152 @ctindex RODFT00
Chris@10 2153 For the case of @code{RODFT00}, this odd symmetry means that
Chris@10 2154 @tex
Chris@10 2155 $X_j = -X_{N-j}$,
Chris@10 2156 @end tex
Chris@10 2157 @ifinfo
Chris@10 2158 X[j] = -X[N-j],
Chris@10 2159 @end ifinfo
Chris@10 2160 @html
Chris@10 2161 <i>X<sub>j</sub> = -X<sub>N-j</sub></i>,
Chris@10 2162 @end html
Chris@10 2163 where we take @math{X} to be periodic so that
Chris@10 2164 @tex
Chris@10 2165 $X_N = X_0$.
Chris@10 2166 @end tex
Chris@10 2167 @ifinfo
Chris@10 2168 X[N] = X[0].
Chris@10 2169 @end ifinfo
Chris@10 2170 @html
Chris@10 2171 <i>X<sub>N</sub> = X</i><sub>0</sub>.
Chris@10 2172 @end html
Chris@10 2173 Because of this redundancy, only the first @math{n} real numbers
Chris@10 2174 starting at @math{j=1} are actually stored (the @math{j=0} element is
Chris@10 2175 zero), where @math{N = 2(n+1)}.
Chris@10 2176
Chris@10 2177 The proper definition of odd symmetry for @code{RODFT10},
Chris@10 2178 @code{RODFT01}, and @code{RODFT11} transforms is somewhat more intricate
Chris@10 2179 because of the shifts by @math{1/2} of the input and/or output, although
Chris@10 2180 the corresponding boundary conditions are given in @ref{Real even/odd
Chris@10 2181 DFTs (cosine/sine transforms)}. Because of the odd symmetry, however,
Chris@10 2182 the cosine terms in the DFT all cancel and the remaining sine terms are
Chris@10 2183 written explicitly below. This formulation often leads people to call
Chris@10 2184 such a transform a @dfn{discrete sine transform} (DST), although it is
Chris@10 2185 really just a special case of the DFT.
Chris@10 2186 @cindex discrete sine transform
Chris@10 2187 @cindex DST
Chris@10 2188
Chris@10 2189
Chris@10 2190 In each of the definitions below, we transform a real array @math{X} of
Chris@10 2191 length @math{n} to a real array @math{Y} of length @math{n}:
Chris@10 2192
Chris@10 2193 @subsubheading RODFT00 (DST-I)
Chris@10 2194 @ctindex RODFT00
Chris@10 2195 An @code{RODFT00} transform (type-I DST) in FFTW is defined by:
Chris@10 2196 @tex
Chris@10 2197 $$
Chris@10 2198 Y_k = 2 \sum_{j=0}^{n-1} X_j \sin [ \pi (j+1) (k+1) / (n+1)].
Chris@10 2199 $$
Chris@10 2200 @end tex
Chris@10 2201 @ifinfo
Chris@10 2202 Y[k] = 2 (sum for j = 0 to n-1 of X[j] sin(pi (j+1)(k+1) / (n+1))).
Chris@10 2203 @end ifinfo
Chris@10 2204 @html
Chris@10 2205 <center><img src="equation-rodft00.png" align="top">.</center>
Chris@10 2206 @end html
Chris@10 2207
Chris@10 2208 @subsubheading RODFT10 (DST-II)
Chris@10 2209 @ctindex RODFT10
Chris@10 2210 An @code{RODFT10} transform (type-II DST) in FFTW is defined by:
Chris@10 2211 @tex
Chris@10 2212 $$
Chris@10 2213 Y_k = 2 \sum_{j=0}^{n-1} X_j \sin [\pi (j+1/2) (k+1) / n].
Chris@10 2214 $$
Chris@10 2215 @end tex
Chris@10 2216 @ifinfo
Chris@10 2217 Y[k] = 2 (sum for j = 0 to n-1 of X[j] sin(pi (j+1/2) (k+1) / n)).
Chris@10 2218 @end ifinfo
Chris@10 2219 @html
Chris@10 2220 <center><img src="equation-rodft10.png" align="top">.</center>
Chris@10 2221 @end html
Chris@10 2222
Chris@10 2223 @subsubheading RODFT01 (DST-III)
Chris@10 2224 @ctindex RODFT01
Chris@10 2225 An @code{RODFT01} transform (type-III DST) in FFTW is defined by:
Chris@10 2226 @tex
Chris@10 2227 $$
Chris@10 2228 Y_k = (-1)^k X_{n-1} + 2 \sum_{j=0}^{n-2} X_j \sin [\pi (j+1) (k+1/2) / n].
Chris@10 2229 $$
Chris@10 2230 @end tex
Chris@10 2231 @ifinfo
Chris@10 2232 Y[k] = (-1)^k X[n-1] + 2 (sum for j = 0 to n-2 of X[j] sin(pi (j+1) (k+1/2) / n)).
Chris@10 2233 @end ifinfo
Chris@10 2234 @html
Chris@10 2235 <center><img src="equation-rodft01.png" align="top">.</center>
Chris@10 2236 @end html
Chris@10 2237 In the case of @math{n=1}, this reduces to
Chris@10 2238 @tex
Chris@10 2239 $Y_0 = X_0$.
Chris@10 2240 @end tex
Chris@10 2241 @ifinfo
Chris@10 2242 Y[0] = X[0].
Chris@10 2243 @end ifinfo
Chris@10 2244 @html
Chris@10 2245 <i>Y</i><sub>0</sub> = <i>X</i><sub>0</sub>.
Chris@10 2246 @end html
Chris@10 2247
Chris@10 2248 @subsubheading RODFT11 (DST-IV)
Chris@10 2249 @ctindex RODFT11
Chris@10 2250 An @code{RODFT11} transform (type-IV DST) in FFTW is defined by:
Chris@10 2251 @tex
Chris@10 2252 $$
Chris@10 2253 Y_k = 2 \sum_{j=0}^{n-1} X_j \sin [\pi (j+1/2) (k+1/2) / n].
Chris@10 2254 $$
Chris@10 2255 @end tex
Chris@10 2256 @ifinfo
Chris@10 2257 Y[k] = 2 (sum for j = 0 to n-1 of X[j] sin(pi (j+1/2) (k+1/2) / n)).
Chris@10 2258 @end ifinfo
Chris@10 2259 @html
Chris@10 2260 <center><img src="equation-rodft11.png" align="top">.</center>
Chris@10 2261 @end html
Chris@10 2262
Chris@10 2263 @subsubheading Inverses and Normalization
Chris@10 2264
Chris@10 2265 These definitions correspond directly to the unnormalized DFTs used
Chris@10 2266 elsewhere in FFTW (hence the factors of @math{2} in front of the
Chris@10 2267 summations). The unnormalized inverse of @code{RODFT00} is
Chris@10 2268 @code{RODFT00}, of @code{RODFT10} is @code{RODFT01} and vice versa, and
Chris@10 2269 of @code{RODFT11} is @code{RODFT11}. Each unnormalized inverse results
Chris@10 2270 in the original array multiplied by @math{N}, where @math{N} is the
Chris@10 2271 @emph{logical} DFT size. For @code{RODFT00}, @math{N=2(n+1)};
Chris@10 2272 otherwise, @math{N=2n}.
Chris@10 2273 @cindex normalization
Chris@10 2274
Chris@10 2275
Chris@10 2276 In defining the discrete sine transform, some authors also include
Chris@10 2277 additional factors of
Chris@10 2278 @ifinfo
Chris@10 2279 sqrt(2)
Chris@10 2280 @end ifinfo
Chris@10 2281 @html
Chris@10 2282 &radic;2
Chris@10 2283 @end html
Chris@10 2284 @tex
Chris@10 2285 $\sqrt{2}$
Chris@10 2286 @end tex
Chris@10 2287 (or its inverse) multiplying selected inputs and/or outputs. This is a
Chris@10 2288 mostly cosmetic change that makes the transform orthogonal, but
Chris@10 2289 sacrifices the direct equivalence to an antisymmetric DFT.
Chris@10 2290
Chris@10 2291 @c =========>
Chris@10 2292 @node 1d Discrete Hartley Transforms (DHTs), Multi-dimensional Transforms, 1d Real-odd DFTs (DSTs), What FFTW Really Computes
Chris@10 2293 @subsection 1d Discrete Hartley Transforms (DHTs)
Chris@10 2294
Chris@10 2295 @cindex discrete Hartley transform
Chris@10 2296 @cindex DHT
Chris@10 2297 The discrete Hartley transform (DHT) of a 1d real array @math{X} of size
Chris@10 2298 @math{n} computes a real array @math{Y} of the same size, where:
Chris@10 2299 @tex
Chris@10 2300 $$
Chris@10 2301 Y_k = \sum_{j = 0}^{n - 1} X_j [ \cos(2\pi j k / n) + \sin(2\pi j k / n)].
Chris@10 2302 $$
Chris@10 2303 @end tex
Chris@10 2304 @ifinfo
Chris@10 2305 @center Y[k] = sum for j = 0 to (n - 1) of X[j] * [cos(2 pi j k / n) + sin(2 pi j k / n)].
Chris@10 2306 @end ifinfo
Chris@10 2307 @html
Chris@10 2308 <center><img src="equation-dht.png" align="top">.</center>
Chris@10 2309 @end html
Chris@10 2310
Chris@10 2311 @cindex normalization
Chris@10 2312 FFTW computes an unnormalized transform, in that there is no coefficient
Chris@10 2313 in front of the summation in the DHT. In other words, applying the
Chris@10 2314 transform twice (the DHT is its own inverse) will multiply the input by
Chris@10 2315 @math{n}.
Chris@10 2316
Chris@10 2317 @c =========>
Chris@10 2318 @node Multi-dimensional Transforms, , 1d Discrete Hartley Transforms (DHTs), What FFTW Really Computes
Chris@10 2319 @subsection Multi-dimensional Transforms
Chris@10 2320
Chris@10 2321 The multi-dimensional transforms of FFTW, in general, compute simply the
Chris@10 2322 separable product of the given 1d transform along each dimension of the
Chris@10 2323 array. Since each of these transforms is unnormalized, computing the
Chris@10 2324 forward followed by the backward/inverse multi-dimensional transform
Chris@10 2325 will result in the original array scaled by the product of the
Chris@10 2326 normalization factors for each dimension (e.g. the product of the
Chris@10 2327 dimension sizes, for a multi-dimensional DFT).
Chris@10 2328
Chris@10 2329 @tex
Chris@10 2330 As an explicit example, consider the following exact mathematical
Chris@10 2331 definition of our multi-dimensional DFT. Let $X$ be a $d$-dimensional
Chris@10 2332 complex array whose elements are $X[j_1, j_2, \ldots, j_d]$, where $0
Chris@10 2333 \leq j_s < n_s$ for all~$s \in \{ 1, 2, \ldots, d \}$. Let also
Chris@10 2334 $\omega_s = e^{2\pi \sqrt{-1}/n_s}$, for all ~$s \in \{ 1, 2, \ldots, d
Chris@10 2335 \}$.
Chris@10 2336
Chris@10 2337 The forward transform computes a complex array~$Y$, whose
Chris@10 2338 structure is the same as that of~$X$, defined by
Chris@10 2339
Chris@10 2340 $$
Chris@10 2341 Y[k_1, k_2, \ldots, k_d] =
Chris@10 2342 \sum_{j_1 = 0}^{n_1 - 1}
Chris@10 2343 \sum_{j_2 = 0}^{n_2 - 1}
Chris@10 2344 \cdots
Chris@10 2345 \sum_{j_d = 0}^{n_d - 1}
Chris@10 2346 X[j_1, j_2, \ldots, j_d]
Chris@10 2347 \omega_1^{-j_1 k_1}
Chris@10 2348 \omega_2^{-j_2 k_2}
Chris@10 2349 \cdots
Chris@10 2350 \omega_d^{-j_d k_d} \ .
Chris@10 2351 $$
Chris@10 2352
Chris@10 2353 The backward transform computes
Chris@10 2354 $$
Chris@10 2355 Y[k_1, k_2, \ldots, k_d] =
Chris@10 2356 \sum_{j_1 = 0}^{n_1 - 1}
Chris@10 2357 \sum_{j_2 = 0}^{n_2 - 1}
Chris@10 2358 \cdots
Chris@10 2359 \sum_{j_d = 0}^{n_d - 1}
Chris@10 2360 X[j_1, j_2, \ldots, j_d]
Chris@10 2361 \omega_1^{j_1 k_1}
Chris@10 2362 \omega_2^{j_2 k_2}
Chris@10 2363 \cdots
Chris@10 2364 \omega_d^{j_d k_d} \ .
Chris@10 2365 $$
Chris@10 2366
Chris@10 2367 Computing the forward transform followed by the backward transform
Chris@10 2368 will multiply the array by $\prod_{s=1}^{d} n_d$.
Chris@10 2369 @end tex
Chris@10 2370
Chris@10 2371 @cindex r2c
Chris@10 2372 The definition of FFTW's multi-dimensional DFT of real data (r2c)
Chris@10 2373 deserves special attention. In this case, we logically compute the full
Chris@10 2374 multi-dimensional DFT of the input data; since the input data are purely
Chris@10 2375 real, the output data have the Hermitian symmetry and therefore only one
Chris@10 2376 non-redundant half need be stored. More specifically, for an @ndims multi-dimensional real-input DFT, the full (logical) complex output array
Chris@10 2377 @tex
Chris@10 2378 $Y[k_0, k_1, \ldots, k_{d-1}]$
Chris@10 2379 @end tex
Chris@10 2380 @html
Chris@10 2381 <i>Y</i>[<i>k</i><sub>0</sub>, <i>k</i><sub>1</sub>, ...,
Chris@10 2382 <i>k</i><sub><i>d-1</i></sub>]
Chris@10 2383 @end html
Chris@10 2384 @ifinfo
Chris@10 2385 Y[k[0], k[1], ..., k[d-1]]
Chris@10 2386 @end ifinfo
Chris@10 2387 has the symmetry:
Chris@10 2388 @tex
Chris@10 2389 $$
Chris@10 2390 Y[k_0, k_1, \ldots, k_{d-1}] = Y[n_0 - k_0, n_1 - k_1, \ldots, n_{d-1} - k_{d-1}]^*
Chris@10 2391 $$
Chris@10 2392 @end tex
Chris@10 2393 @html
Chris@10 2394 <i>Y</i>[<i>k</i><sub>0</sub>, <i>k</i><sub>1</sub>, ...,
Chris@10 2395 <i>k</i><sub><i>d-1</i></sub>] = <i>Y</i>[<i>n</i><sub>0</sub> -
Chris@10 2396 <i>k</i><sub>0</sub>, <i>n</i><sub>1</sub> - <i>k</i><sub>1</sub>, ...,
Chris@10 2397 <i>n</i><sub><i>d-1</i></sub> - <i>k</i><sub><i>d-1</i></sub>]<sup>*</sup>
Chris@10 2398 @end html
Chris@10 2399 @ifinfo
Chris@10 2400 Y[k[0], k[1], ..., k[d-1]] = Y[n[0] - k[0], n[1] - k[1], ..., n[d-1] - k[d-1]]*
Chris@10 2401 @end ifinfo
Chris@10 2402 (where each dimension is periodic). Because of this symmetry, we only
Chris@10 2403 store the
Chris@10 2404 @tex
Chris@10 2405 $k_{d-1} = 0 \cdots n_{d-1}/2$
Chris@10 2406 @end tex
Chris@10 2407 @html
Chris@10 2408 <i>k</i><sub><i>d-1</i></sub> = 0...<i>n</i><sub><i>d-1</i></sub>/2+1
Chris@10 2409 @end html
Chris@10 2410 @ifinfo
Chris@10 2411 k[d-1] = 0...n[d-1]/2
Chris@10 2412 @end ifinfo
Chris@10 2413 elements of the @emph{last} dimension (division by @math{2} is rounded
Chris@10 2414 down). (We could instead have cut any other dimension in half, but the
Chris@10 2415 last dimension proved computationally convenient.) This results in the
Chris@10 2416 peculiar array format described in more detail by @ref{Real-data DFT
Chris@10 2417 Array Format}.
Chris@10 2418
Chris@10 2419 The multi-dimensional c2r transform is simply the unnormalized inverse
Chris@10 2420 of the r2c transform. i.e. it is the same as FFTW's complex backward
Chris@10 2421 multi-dimensional DFT, operating on a Hermitian input array in the
Chris@10 2422 peculiar format mentioned above and outputting a real array (since the
Chris@10 2423 DFT output is purely real).
Chris@10 2424
Chris@10 2425 We should remind the user that the separable product of 1d transforms
Chris@10 2426 along each dimension, as computed by FFTW, is not always the same thing
Chris@10 2427 as the usual multi-dimensional transform. A multi-dimensional
Chris@10 2428 @code{R2HC} (or @code{HC2R}) transform is not identical to the
Chris@10 2429 multi-dimensional DFT, requiring some post-processing to combine the
Chris@10 2430 requisite real and imaginary parts, as was described in @ref{The
Chris@10 2431 Halfcomplex-format DFT}. Likewise, FFTW's multidimensional
Chris@10 2432 @code{FFTW_DHT} r2r transform is not the same thing as the logical
Chris@10 2433 multi-dimensional discrete Hartley transform defined in the literature,
Chris@10 2434 as discussed in @ref{The Discrete Hartley Transform}.
Chris@10 2435