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>−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 √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 √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
|