Subversion Repositories shark

\input texinfo   @c -*- texinfo -*-

@c % $Id: fftw.tex,v 1.1.1.1 2002-03-29 14:12:55 pj Exp $

@c %**start of header

@setfilename fftw.info

@settitle FFTW

@c %**end of header

@include version.texi

@setchapternewpage odd

@c define constant index (ct)

@defcodeindex ct

@syncodeindex ct fn

@syncodeindex vr fn

@syncodeindex pg fn

@syncodeindex tp fn

@c define foreign function index (ff)

@defcodeindex ff

@syncodeindex ff cp

@c define foreign constant index (fc)

@defcodeindex fc

@syncodeindex fc cp

@c define foreign program index (fp)

@defcodeindex fp

@syncodeindex fp cp

@ifinfo

This is the FFTW User's manual.

Copyright @copyright{} 1997--1999 Massachusetts Institute of Technology

Permission is granted to make and distribute verbatim copies of this

manual provided the copyright notice and this permission notice are

preserved on all copies.

Permission is granted to copy and distribute modified versions of this

manual under the conditions for verbatim copying, provided that the

entire resulting derived work is distributed under the terms of a

permission notice identical to this one.

Permission is granted to copy and distribute translations of this manual

into another language, under the above conditions for modified versions,

except that this permission notice may be stated in a translation

approved by the Free Software Foundation.

@end ifinfo

@titlepage

@sp 10

@comment The title is printed in a large font.

@title{FFTW User's Manual}

@subtitle For version @value{VERSION}, @value{UPDATED}

@author{Matteo Frigo}

@author{Steven G. Johnson}

@c The following two commands start the copyright page.

@page

@vskip 0pt plus 1filll

Copyright @copyright{} 1997--1999 Massachusetts Institute of Technology.

Permission is granted to make and distribute verbatim copies of this

manual provided the copyright notice and this permission notice are

preserved on all copies.

Permission is granted to copy and distribute modified versions of this

manual under the conditions for verbatim copying, provided that the

entire resulting derived work is distributed under the terms of a

permission notice identical to this one.

Permission is granted to copy and distribute translations of this manual

into another language, under the above conditions for modified versions,

except that this permission notice may be stated in a translation

approved by the Free Software Foundation.

@end titlepage

@node    Top, Introduction, (dir), (dir)

@ifinfo

@top FFTW User Manual

Welcome to FFTW, the Fastest Fourier Transform in the West.  FFTW is a

collection of fast C routines to compute the discrete Fourier transform.

This manual documents FFTW version @value{VERSION}.

@end ifinfo

@menu

* Introduction::                

* Tutorial::                    

* FFTW Reference::              

* Parallel FFTW::              

* Calling FFTW from Fortran::  

* Installation and Customization::  

* Acknowledgments::            

* License and Copyright::      

* Concept Index::              

* Library Index::              

@detailmenu --- The Detailed Node Listing ---

Tutorial

* Complex One-dimensional Transforms Tutorial::  

* Complex Multi-dimensional Transforms Tutorial::  

* Real One-dimensional Transforms Tutorial::  

* Real Multi-dimensional Transforms Tutorial::  

* Multi-dimensional Array Format::  

* Words of Wisdom::            

Multi-dimensional Array Format

* Row-major Format::            

* Column-major Format::        

* Static Arrays in C::          

* Dynamic Arrays in C::        

* Dynamic Arrays in C-The Wrong Way::  

Words of Wisdom

* Caveats in Using Wisdom::     What you should worry about in using wisdom

* Importing and Exporting Wisdom::  I/O of wisdom to disk and other media

FFTW Reference

* Data Types::                  real, complex, and halfcomplex numbers

* One-dimensional Transforms Reference::  

* Multi-dimensional Transforms Reference::  

* Real One-dimensional Transforms Reference::  

* Real Multi-dimensional Transforms Reference::  

* Wisdom Reference::            

* Memory Allocator Reference::  

* Thread safety::              

One-dimensional Transforms Reference

* fftw_create_plan::            Plan Creation

* Discussion on Specific Plans::  

* fftw::                        Plan Execution

* fftw_destroy_plan::           Plan Destruction

* What FFTW Really Computes::   Definition of the DFT.

Multi-dimensional Transforms Reference

* fftwnd_create_plan::          Plan Creation

* fftwnd::                      Plan Execution

* fftwnd_destroy_plan::         Plan Destruction

* What FFTWND Really Computes::  

Real One-dimensional Transforms Reference

* rfftw_create_plan::           Plan Creation  

* rfftw::                       Plan Execution  

* rfftw_destroy_plan::          Plan Destruction

* What RFFTW Really Computes::  

Real Multi-dimensional Transforms Reference

* rfftwnd_create_plan::         Plan Creation

* rfftwnd::                     Plan Execution

* Array Dimensions for Real Multi-dimensional Transforms::  

* Strides in In-place RFFTWND::  

* rfftwnd_destroy_plan::        Plan Destruction

* What RFFTWND Really Computes::  

Wisdom Reference

* fftw_export_wisdom::          

* fftw_import_wisdom::          

* fftw_forget_wisdom::          

Parallel FFTW

* Multi-threaded FFTW::        

* MPI FFTW::                    

Multi-threaded FFTW

* Installation and Supported Hardware/Software::  

* Usage of Multi-threaded FFTW::  

* How Many Threads to Use?::    

* Using Multi-threaded FFTW in a Multi-threaded Program::  

* Tips for Optimal Threading::  

MPI FFTW

* MPI FFTW Installation::      

* Usage of MPI FFTW for Complex Multi-dimensional Transforms::  

* MPI Data Layout::            

* Usage of MPI FFTW for Real Multi-dimensional Transforms::  

* Usage of MPI FFTW for Complex One-dimensional Transforms::  

* MPI Tips::                    

Calling FFTW from Fortran

* Wrapper Routines::            

* FFTW Constants in Fortran::  

* Fortran Examples::            

Installation and Customization

* Installation on Unix::        

* Installation on non-Unix Systems::  

* Installing FFTW in both single and double precision::  

* gcc and Pentium/PentiumPro hacks::  

* Customizing the timer::      

* Generating your own code::    

@end detailmenu

@end menu

@c ************************************************************

@node    Introduction, Tutorial, Top, Top

@chapter Introduction

This manual documents version @value{VERSION} of FFTW, the @emph{Fastest

Fourier Transform in the West}.  FFTW is a comprehensive collection of

fast C routines for computing the discrete Fourier transform (DFT) in

one or more dimensions, of both real and complex data, and of arbitrary

input size.  FFTW also includes parallel transforms for both shared- and

distributed-memory systems.  We assume herein that the reader is already

familiar with the properties and uses of the DFT that are relevant to

her application.  Otherwise, see e.g. @cite{The Fast Fourier Transform}

by E. O. Brigham (Prentice-Hall, Englewood Cliffs, NJ, 1974).

@uref{http://theory.lcs.mit.edu/~fftw, Our web page} also has links to

FFT-related information online.

@cindex FFTW

FFTW is usually faster (and sometimes much faster) than all other

freely-available Fourier transform programs found on the Net.  For

transforms whose size is a power of two, it compares favorably with the

FFT codes in Sun's Performance Library and IBM's ESSL library, which are

targeted at specific machines.  Moreover, FFTW's performance is

@emph{portable}.  Indeed, FFTW is unique in that it automatically adapts

itself to your machine, your cache, the size of your memory, the number

of registers, and all the other factors that normally make it impossible

to optimize a program for more than one machine.  An extensive

comparison of FFTW's performance with that of other Fourier transform

codes has been made. The results are available on the Web at

@uref{http://theory.lcs.mit.edu/~benchfft, the benchFFT home page}.

@cindex benchmark

@fpindex benchfft

In order to use FFTW effectively, you need to understand one basic

concept of FFTW's internal structure.  FFTW does not used a fixed

algorithm for computing the transform, but it can adapt the DFT

algorithm to details of the underlying hardware in order to achieve best

performance.  Hence, the computation of the transform is split into two

phases.  First, FFTW's @dfn{planner} is called, which ``learns'' the

@cindex plan

fastest way to compute the transform on your machine.  The planner

@cindex planner

produces a data structure called a @dfn{plan} that contains this

information.  Subsequently, the plan is passed to FFTW's @dfn{executor},

@cindex executor

along with an array of input data.  The executor computes the actual

transform, as dictated by the plan.  The plan can be reused as many

times as needed.  In typical high-performance applications, many

transforms of the same size are computed, and consequently a

relatively-expensive initialization of this sort is acceptable.  On the

other hand, if you need a single transform of a given size, the one-time

cost of the planner becomes significant.  For this case, FFTW provides

fast planners based on heuristics or on previously computed plans.

The pattern of planning/execution applies to all four operation modes of

FFTW, that is, @w{I) one-dimensional} complex transforms (FFTW), @w{II)

multi-dimensional} complex transforms (FFTWND), @w{III) one-dimensional}

transforms of real data (RFFTW), @w{IV) multi-dimensional} transforms of

real data (RFFTWND).  Each mode comes with its own planner and executor.

Besides the automatic performance adaptation performed by the planner,

it is also possible for advanced users to customize FFTW for their

special needs.  As distributed, FFTW works most efficiently for arrays

whose size can be factored into small primes (@math{2}, @math{3},

@math{5}, and @math{7}), and uses a slower general-purpose routine for

other factors.  FFTW, however, comes with a code generator that can

produce fast C programs for any particular array size you may care

about.

@cindex code generator

For example, if you need transforms of size

@ifinfo

@math{513 = 19 x 3^3},

@end ifinfo

@iftex

@tex

$513 = 19 \cdot 3^3$,

@end tex

@end iftex

@ifhtml

513&nbsp;=&nbsp;19*3<sup>3</sup>,

@end ifhtml

you can customize FFTW to support the factor @math{19} efficiently.

FFTW can exploit multiple processors if you have them.  FFTW comes with

a shared-memory implementation on top of POSIX (and similar) threads, as

well as a distributed-memory implementation based on MPI.

@cindex parallel transform

@cindex threads

@cindex MPI

We also provide an experimental parallel implementation written in Cilk,

@emph{the superior programming tool of choice for discriminating

hackers} (Olin Shivers).  (See @uref{http://supertech.lcs.mit.edu/cilk,

the Cilk home page}.)

@cindex Cilk

For more information regarding FFTW, see the paper, ``The Fastest

Fourier Transform in the West,'' by M. Frigo and S. G. Johnson, which is

the technical report MIT-LCS-TR-728 (Sep. '97).  See also, ``FFTW: An

Adaptive Software Architecture for the FFT,'' by M. Frigo and

S. G. Johnson, which appeared in the 23rd International Conference on

Acoustics, Speech, and Signal Processing (@cite{Proc. ICASSP 1998}

@b{3}, p. 1381).  The code generator is described in the paper ``A Fast

Fourier Transform Compiler'',

@cindex compiler

by M. Frigo, to appear in the @cite{Proceedings of the 1999 ACM SIGPLAN

Conference on Programming Language Design and Implementation (PLDI),

Atlanta, Georgia, May 1999}.  These papers, along with the latest

version of FFTW, the FAQ, benchmarks, and other links, are available at

@uref{http://theory.lcs.mit.edu/~fftw, the FFTW home page}.  The current

version of FFTW incorporates many good ideas from the past thirty years

of FFT literature.  In one way or another, FFTW uses the Cooley-Tukey

algorithm, the Prime Factor algorithm, Rader's algorithm for prime

sizes, and the split-radix algorithm (with a variation due to Dan

Bernstein).  Our code generator also produces new algorithms that we do

not yet completely understand.

@cindex algorithm

The reader is referred to the cited papers for the appropriate

references.

The rest of this manual is organized as follows.  We first discuss the

sequential (one-processor) implementation.  We start by describing the

basic features of FFTW in @ref{Tutorial}.  This discussion includes the

storage scheme of multi-dimensional arrays (@ref{Multi-dimensional Array

Format}) and FFTW's mechanisms for storing plans on disk (@ref{Words of

Wisdom}).  Next, @ref{FFTW Reference} provides comprehensive

documentation of all FFTW's features.  Parallel transforms are discussed

in their own chapter @ref{Parallel FFTW}.  Fortran programmers can also

use FFTW, as described in @ref{Calling FFTW from Fortran}.

@ref{Installation and Customization} explains how to install FFTW in

your computer system and how to adapt FFTW to your needs.  License and

copyright information is given in @ref{License and Copyright}.  Finally,

we thank all the people who helped us in @ref{Acknowledgments}.

@c ************************************************************

@node  Tutorial, FFTW Reference, Introduction, Top

@chapter Tutorial

@cindex Tutorial

This chapter describes the basic usage of FFTW, i.e., how to compute the

Fourier transform of a single array.  This chapter tells the truth, but

not the @emph{whole} truth. Specifically, FFTW implements additional

routines and flags, providing extra functionality, that are not

documented here.  @xref{FFTW Reference}, for more complete information.

(Note that you need to compile and install FFTW before you can use it in

a program.  @xref{Installation and Customization}, for the details of

the installation.)

This tutorial chapter is structured as follows.  @ref{Complex

One-dimensional Transforms Tutorial} describes the basic usage of the

one-dimensional transform of complex data.  @ref{Complex

Multi-dimensional Transforms Tutorial} describes the basic usage of the

multi-dimensional transform of complex data.  @ref{Real One-dimensional

Transforms Tutorial} describes the one-dimensional transform of real

data and its inverse.  Finally, @ref{Real Multi-dimensional Transforms

Tutorial} describes the multi-dimensional transform of real data and its

inverse.  We recommend that you read these sections in the order that

they are presented.  We then discuss two topics in detail.  In

@ref{Multi-dimensional Array Format}, we discuss the various

alternatives for storing multi-dimensional arrays in memory.  @ref{Words

of Wisdom} shows how you can save FFTW's plans for future use.

@menu

* Complex One-dimensional Transforms Tutorial::  

* Complex Multi-dimensional Transforms Tutorial::  

* Real One-dimensional Transforms Tutorial::  

* Real Multi-dimensional Transforms Tutorial::  

* Multi-dimensional Array Format::  

* Words of Wisdom::            

@end menu

@node  Complex One-dimensional Transforms Tutorial, Complex Multi-dimensional Transforms Tutorial, Tutorial, Tutorial

@section Complex One-dimensional Transforms Tutorial

@cindex complex one-dimensional transform

@cindex complex transform

The basic usage of FFTW is simple.  A typical call to FFTW looks like:

@example

#include <fftw.h>

...

@{

     fftw_complex in[N], out[N];

     fftw_plan p;

     ...

     p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);

     ...

     fftw_one(p, in, out);

     ...

     fftw_destroy_plan(p);  

@}

@end example

The first thing we do is to create a @dfn{plan}, which is an object

@cindex plan

that contains all the data that FFTW needs to compute the FFT, using the

following function:

@example

fftw_plan fftw_create_plan(int n, fftw_direction dir, int flags);

@end example

@findex fftw_create_plan

@findex fftw_direction

@tindex fftw_plan

The first argument, @code{n}, is the size of the transform you are

trying to compute.  The size @code{n} can be any positive integer, but

sizes that are products of small factors are transformed most

efficiently.  The second argument, @code{dir}, can be either

@code{FFTW_FORWARD} or @code{FFTW_BACKWARD}, and indicates the direction

of the transform you

@ctindex FFTW_FORWARD

@ctindex FFTW_BACKWARD

are interested in.  Alternatively, you can use the sign of the exponent

in the transform, @math{-1} or @math{+1}, which corresponds to

@code{FFTW_FORWARD} or @code{FFTW_BACKWARD} respectively.  The

@code{flags} argument is either @code{FFTW_MEASURE} or

@cindex flags

@code{FFTW_ESTIMATE}.  @code{FFTW_MEASURE} means that FFTW actually runs

@ctindex FFTW_MEASURE

and measures the execution time of several FFTs in order to find the

best way to compute the transform of size @code{n}.  This may take some

time, depending on your installation and on the precision of the timer

in your machine.  @code{FFTW_ESTIMATE}, on the contrary, does not run

any computation, and just builds a

@ctindex FFTW_ESTIMATE

reasonable plan, which may be sub-optimal.  In other words, if your

program performs many transforms of the same size and initialization

time is not important, use @code{FFTW_MEASURE}; otherwise use the

estimate.  (A compromise between these two extremes exists.  @xref{Words

of Wisdom}.)

Once the plan has been created, you can use it as many times as you like

for transforms on arrays of the same size.  When you are done with the

plan, you deallocate it by calling @code{fftw_destroy_plan(plan)}.

@findex fftw_destroy_plan

The transform itself is computed by passing the plan along with the

input and output arrays to @code{fftw_one}:

@example

void fftw_one(fftw_plan plan, fftw_complex *in, fftw_complex *out);

@end example

@findex fftw_one

Note that the transform is out of place: @code{in} and @code{out} must

point to distinct arrays. It operates on data of type

@code{fftw_complex}, a data structure with real (@code{in[i].re}) and

imaginary (@code{in[i].im}) floating-point components.  The @code{in}

and @code{out} arrays should have the length specified when the plan was

created.  An alternative function, @code{fftw}, allows you to

efficiently perform multiple and/or strided transforms (@pxref{FFTW

Reference}).

@tindex fftw_complex

The DFT results are stored in-order in the array @code{out}, with the

zero-frequency (DC) component in @code{out[0]}.

@cindex frequency

The array @code{in} is not modified.  Users should note that FFTW

computes an unnormalized DFT, the sign of whose exponent is given by the

@code{dir} parameter of @code{fftw_create_plan}.  Thus, computing a

forward followed by a backward transform (or vice versa) results in the

original array scaled by @code{n}.  @xref{What FFTW Really Computes},

for the definition of DFT.

@cindex normalization

A program using FFTW should be linked with @code{-lfftw -lm} on Unix

systems, or with the FFTW and standard math libraries in general.

@cindex linking on Unix

@node Complex Multi-dimensional Transforms Tutorial, Real One-dimensional Transforms Tutorial, Complex One-dimensional Transforms Tutorial, Tutorial

@section Complex Multi-dimensional Transforms Tutorial

@cindex complex multi-dimensional transform

@cindex multi-dimensional transform

FFTW can also compute transforms of any number of dimensions

(@dfn{rank}).  The syntax is similar to that for the one-dimensional

@cindex rank

transforms, with @samp{fftw_} replaced by @samp{fftwnd_} (which stands

for ``@code{fftw} in @code{N} dimensions'').

As before, we @code{#include <fftw.h>} and create a plan for the

transforms, this time of type @code{fftwnd_plan}:

@example

fftwnd_plan fftwnd_create_plan(int rank, const int *n,

                               fftw_direction dir, int flags);

@end example

@tindex fftwnd_plan

@tindex fftw_direction

@findex fftwnd_create_plan

@code{rank} is the dimensionality of the array, and can be any

non-negative integer.  The next argument, @code{n}, is a pointer to an

integer array of length @code{rank} containing the (positive) sizes of

each dimension of the array.  (Note that the array will be stored in

row-major order. @xref{Multi-dimensional Array Format}, for information

on row-major order.)  The last two parameters are the same as in

@code{fftw_create_plan}.  We now, however, have an additional possible

flag, @code{FFTW_IN_PLACE}, since @code{fftwnd} supports true in-place

@cindex flags

@ctindex FFTW_IN_PLACE

@findex fftwnd

transforms.  Multiple flags are combined using a bitwise @dfn{or}

(@samp{|}).  (An @dfn{in-place} transform is one in which the output

data overwrite the input data.  It thus requires half as much memory

as---and is often faster than---its opposite, an @dfn{out-of-place}

transform.)

@cindex in-place transform

@cindex out-of-place transform

For two- and three-dimensional transforms, FFTWND provides alternative

routines that accept the sizes of each dimension directly, rather than

indirectly through a rank and an array of sizes.  These are otherwise

identical to @code{fftwnd_create_plan}, and are sometimes more

convenient:

@example

fftwnd_plan fftw2d_create_plan(int nx, int ny,

                               fftw_direction dir, int flags);

fftwnd_plan fftw3d_create_plan(int nx, int ny, int nz,

                               fftw_direction dir, int flags);

@end example

@findex fftw2d_create_plan

@findex fftw3d_create_plan

Once the plan has been created, you can use it any number of times for

transforms of the same size.  When you do not need a plan anymore, you

can deallocate the plan by calling @code{fftwnd_destroy_plan(plan)}.

@findex fftwnd_destroy_plan

Given a plan, you can compute the transform of an array of data by

calling:

@example

void fftwnd_one(fftwnd_plan plan, fftw_complex *in, fftw_complex *out);

@end example

@findex fftwnd_one

Here, @code{in} and @code{out} point to multi-dimensional arrays in

row-major order, of the size specified when the plan was created.  In

the case of an in-place transform, the @code{out} parameter is ignored

and the output data are stored in the @code{in} array.  The results are

stored in-order, unnormalized, with the zero-frequency component in

@code{out[0]}.

@cindex frequency

A forward followed by a backward transform (or vice-versa) yields the

original data multiplied by the size of the array (i.e. the product of

the dimensions).  @xref{What FFTWND Really Computes}, for a discussion

of what FFTWND computes.

@cindex normalization

For example, code to perform an in-place FFT of a three-dimensional

array might look like:

@example

#include <fftw.h>

...

@{

     fftw_complex in[L][M][N];

     fftwnd_plan p;

     ...

     p = fftw3d_create_plan(L, M, N, FFTW_FORWARD,

                            FFTW_MEASURE | FFTW_IN_PLACE);

     ...

     fftwnd_one(p, &in[0][0][0], NULL);

     ...

     fftwnd_destroy_plan(p);  

@}

@end example

Note that @code{in} is a statically-declared array, which is

automatically in row-major order, but we must take the address of the

first element in order to fit the type expected by @code{fftwnd_one}.

(@xref{Multi-dimensional Array Format}.)

@node Real One-dimensional Transforms Tutorial, Real Multi-dimensional Transforms Tutorial, Complex Multi-dimensional Transforms Tutorial, Tutorial

@section Real One-dimensional Transforms Tutorial

@cindex real transform

@cindex complex to real transform

@cindex RFFTW

If the input data are purely real, you can save roughly a factor of two

in both time and storage by using the @dfn{rfftw} transforms, which are

FFTs specialized for real data.  The output of a such a transform is a

@dfn{halfcomplex} array, which consists of only half of the complex DFT

amplitudes (since the negative-frequency amplitudes for real data are

the complex conjugate of the positive-frequency amplitudes).

@cindex halfcomplex array

In exchange for these speed and space advantages, the user sacrifices

some of the simplicity of FFTW's complex transforms.  First of all, to

allow maximum performance, the output format of the one-dimensional real

transforms is different from that used by the multi-dimensional

transforms.  Second, the inverse transform (halfcomplex to real) has the

side-effect of destroying its input array.  Neither of these

inconveniences should pose a serious problem for users, but it is

important to be aware of them.  (Both the inconvenient output format

and the side-effect of the inverse transform can be ameliorated for

one-dimensional transforms, at the expense of some performance, by using

instead the multi-dimensional transform routines with a rank of one.)

The computation of the plan is similar to that for the complex

transforms.  First, you @code{#include <rfftw.h>}.  Then, you create a

plan (of type @code{rfftw_plan}) by calling:

@example

rfftw_plan rfftw_create_plan(int n, fftw_direction dir, int flags);

@end example

@tindex rfftw_plan

@tindex fftw_direction

@findex rfftw_create_plan

@code{n} is the length of the @emph{real} array in the transform (even

for halfcomplex-to-real transforms), and can be any positive integer

(although sizes with small factors are transformed more efficiently).

@code{dir} is either @code{FFTW_REAL_TO_COMPLEX} or

@code{FFTW_COMPLEX_TO_REAL}.

@ctindex FFTW_REAL_TO_COMPLEX

@ctindex FFTW_COMPLEX_TO_REAL

The @code{flags} parameter is the same as in @code{fftw_create_plan}.

Once created, a plan can be used for any number of transforms, and is

deallocated when you are done with it by calling

@code{rfftw_destroy_plan(plan)}.

@findex rfftw_destroy_plan

Given a plan, a real-to-complex or complex-to-real transform is computed

by calling:

@example

void rfftw_one(rfftw_plan plan, fftw_real *in, fftw_real *out);

@end example

@findex rfftw_one

(Note that @code{fftw_real} is an alias for the floating-point type for

which FFTW was compiled.)  Depending upon the direction of the plan,

either the input or the output array is halfcomplex, and is stored in

the following format:

@cindex halfcomplex array

@iftex

@tex

$$

r_0, r_1, r_2, \ldots, r_{n/2}, i_{(n+1)/2-1}, \ldots, i_2, i_1

$$

@end tex

@end iftex

@ifinfo

r0, r1, r2, r(n/2), i((n+1)/2-1), ..., i2, i1

@end ifinfo

@ifhtml

<p align=center>

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>

</p>

@end ifhtml

Here,

@ifinfo

rk

@end ifinfo

@iftex

@tex

$r_k$

@end tex

@end iftex

@ifhtml

r<sub>k</sub>

@end ifhtml

is the real part of the @math{k}th output, and

@ifinfo

ik

@end ifinfo

@iftex

@tex

$i_k$

@end tex

@end iftex

@ifhtml

i<sub>k</sub>

@end ifhtml

is the imaginary part.  (We follow here the C convention that integer

division is rounded down, e.g. @math{7 / 2 = 3}.) For a halfcomplex

array @code{hc[]}, the @math{k}th component has its real part in

@code{hc[k]} and its imaginary part in @code{hc[n-k]}, with the

exception of @code{k} @code{==} @code{0} or @code{n/2} (the latter only

if n is even)---in these two cases, the imaginary part is zero due to

symmetries of the real-complex transform, and is not stored.  Thus, the

transform of @code{n} real values is a halfcomplex array of length

@code{n}, and vice versa.  @footnote{The output for the

multi-dimensional rfftw is a more-conventional array of

@code{fftw_complex} values, but the format here permitted us greater

efficiency in one dimension.}  This is actually only half of the DFT

spectrum of the data.  Although the other half can be obtained by

complex conjugation, it is not required by many applications such as

convolution and filtering.

Like the complex transforms, the RFFTW transforms are unnormalized, so a

forward followed by a backward transform (or vice-versa) yields the

original data scaled by the length of the array, @code{n}.

@cindex normalization

Let us reiterate here our warning that an @code{FFTW_COMPLEX_TO_REAL}

transform has the side-effect of destroying its (halfcomplex) input.

The @code{FFTW_REAL_TO_COMPLEX} transform, however, leaves its (real)

input untouched, just as you would hope.

As an example, here is an outline of how you might use RFFTW to compute

the power spectrum of a real array (i.e. the squares of the absolute

values of the DFT amplitudes):

@cindex power spectrum

@example

#include <rfftw.h>

...

@{

     fftw_real in[N], out[N], power_spectrum[N/2+1];

     rfftw_plan p;

     int k;

     ...

     p = rfftw_create_plan(N, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);

     ...

     rfftw_one(p, in, out);

     power_spectrum[0] = out[0]*out[0];  /* DC component */

     for (k = 1; k < (N+1)/2; ++k)  /* (k < N/2 rounded up) */

          power_spectrum[k] = out[k]*out[k] + out[N-k]*out[N-k];

     if (N % 2 == 0) /* N is even */

          power_spectrum[N/2] = out[N/2]*out[N/2];  /* Nyquist freq. */

     ...

     rfftw_destroy_plan(p);

@}

@end example

Programs using RFFTW should link with @code{-lrfftw -lfftw -lm} on Unix,

or with the FFTW, RFFTW, and math libraries in general.

@cindex linking on Unix

@node Real Multi-dimensional Transforms Tutorial, Multi-dimensional Array Format, Real One-dimensional Transforms Tutorial, Tutorial

@section Real Multi-dimensional Transforms Tutorial

@cindex real multi-dimensional transform

FFTW includes multi-dimensional transforms for real data of any rank.

As with the one-dimensional real transforms, they save roughly a factor

of two in time and storage over complex transforms of the same size.

Also as in one dimension, these gains come at the expense of some

increase in complexity---the output format is different from the

one-dimensional RFFTW (and is more similar to that of the complex FFTW)

and the inverse (complex to real) transforms have the side-effect of

overwriting their input data.  

To use the real multi-dimensional transforms, you first @code{#include

<rfftw.h>} and then create a plan for the size and direction of

transform that you are interested in:

@example

rfftwnd_plan rfftwnd_create_plan(int rank, const int *n,

                                 fftw_direction dir, int flags);

@end example

@tindex rfftwnd_plan

@findex rfftwnd_create_plan

The first two parameters describe the size of the real data (not the

halfcomplex data, which will have different dimensions).  The last two

parameters are the same as those for @code{rfftw_create_plan}.  Just as

for fftwnd, there are two alternate versions of this routine,

@code{rfftw2d_create_plan} and @code{rfftw3d_create_plan}, that are

sometimes more convenient for two- and three-dimensional transforms.

@findex rfftw2d_create_plan

@findex rfftw3d_create_plan

Also as in fftwnd, rfftwnd supports true in-place transforms, specified

by including @code{FFTW_IN_PLACE} in the flags.

Once created, a plan can be used for any number of transforms, and is

deallocated by calling @code{rfftwnd_destroy_plan(plan)}.

Given a plan, the transform is computed by calling one of the following

two routines:

@example

void rfftwnd_one_real_to_complex(rfftwnd_plan plan,

                                 fftw_real *in, fftw_complex *out);

void rfftwnd_one_complex_to_real(rfftwnd_plan plan,

                                 fftw_complex *in, fftw_real *out);

@end example

@findex rfftwnd_one_real_to_complex

@findex rfftwnd_one_complex_to_real

As is clear from their names and parameter types, the former function is

for @code{FFTW_REAL_TO_COMPLEX} transforms and the latter is for

@code{FFTW_COMPLEX_TO_REAL} transforms.  (We could have used only a

single routine, since the direction of the transform is encoded in the

plan, but we wanted to correctly express the datatypes of the

parameters.)  The latter routine, as we discuss elsewhere, has the

side-effect of overwriting its input (except when the rank of the array

is one).  In both cases, the @code{out} parameter is ignored for

in-place transforms.

The format of the complex arrays deserves careful attention.

@cindex rfftwnd array format

Suppose that the real data has dimensions

@iftex

@tex

$n_1 \times n_2 \times \cdots \times n_d$

@end tex

@end iftex

@ifinfo

n1 x n2 x ... x nd

@end ifinfo

@ifhtml

n<sub>1</sub> x n<sub>2</sub> x ... x n<sub>d</sub>

@end ifhtml

(in row-major order).  Then, after a real-to-complex transform, the

output is an

@iftex

@tex

$n_1 \times n_2 \times \cdots \times (n_d/2+1)$

@end tex

@end iftex

@ifinfo

n1 x n2 x ... x (nd/2+1)

@end ifinfo

@ifhtml

n<sub>1</sub> x n<sub>2</sub> x ... x (n<sub>d</sub>/2+1)

@end ifhtml

array of @code{fftw_complex} values in row-major order, corresponding to

slightly over half of the output of the corresponding complex transform.

(Note that the division is rounded down.)  The ordering of the data is

otherwise exactly the same as in the complex case.  (In principle, the

output could be exactly half the size of the complex transform output,

but in more than one dimension this requires too complicated a format to

be practical.)  Note that, unlike the one-dimensional RFFTW, the real

and imaginary parts of the DFT amplitudes are here stored together in

the natural way.

Since the complex data is slightly larger than the real data, some

complications arise for in-place transforms.  In this case, the final

dimension of the real data must be padded with extra values to

accommodate the size of the complex data---two extra if the last

dimension is even and one if it is odd.

@cindex padding

That is, the last dimension of the real data must physically contain

@iftex

@tex

$2 (n_d/2+1)$

@end tex

@end iftex

@ifinfo

2 * (nd/2+1)

@end ifinfo

@ifhtml

2 * (n<sub>d</sub>/2+1)

@end ifhtml

@code{fftw_real} values (exactly enough to hold the complex data).

This physical array size does not, however, change the @emph{logical}

array size---only

@iftex

@tex

$n_d$

@end tex

@end iftex

@ifinfo

nd

@end ifinfo

@ifhtml

n<sub>d</sub>

@end ifhtml

values are actually stored in the last dimension, and

@iftex

@tex

$n_d$

@end tex

@end iftex

@ifinfo

nd

@end ifinfo

@ifhtml

n<sub>d</sub>

@end ifhtml

is the last dimension passed to @code{rfftwnd_create_plan}.

For example, consider the transform of a two-dimensional real array of

size @code{nx} by @code{ny}.  The output of the @code{rfftwnd} transform

is a two-dimensional real array of size @code{nx} by @code{ny/2+1},

where the @code{y} dimension has been cut nearly in half because of

redundancies in the output.  Because @code{fftw_complex} is twice the

size of @code{fftw_real}, the output array is slightly bigger than the

input array.  Thus, if we want to compute the transform in place, we

must @emph{pad} the input array so that it is of size @code{nx} by

@code{2*(ny/2+1)}.  If @code{ny} is even, then there are two padding

elements at the end of each row (which need not be initialized, as they

are only used for output).

@ifhtml

The following illustration depicts the input and output arrays just

described, for both the out-of-place and in-place transforms (with the

arrows indicating consecutive memory locations):

<p align=center><img src="rfftwnd.gif" width=389 height=583>

@end ifhtml

@iftex

Figure 1 depicts the input and output arrays just

described, for both the out-of-place and in-place transforms (with the

arrows indicating consecutive memory locations).

@tex

{

\pageinsert

\vfill

\vskip405pt

\hskip40pt

\special{psfile="rfftwnd.eps"

}

\vskip 24pt

Figure 1: Illustration of the data layout for real to complex

transforms.  

\vfill

\endinsert}

@end tex

@end iftex

The RFFTWND transforms are unnormalized, so a forward followed by a

backward transform will result in the original data scaled by the number

of real data elements---that is, the product of the (logical) dimensions

of the real data.

@cindex normalization

Below, we illustrate the use of RFFTWND by showing how you might use it

to compute the (cyclic) convolution of two-dimensional real arrays

@code{a} and @code{b} (using the identity that a convolution corresponds

to a pointwise product of the Fourier transforms).  For variety,

in-place transforms are used for the forward FFTs and an out-of-place

transform is used for the inverse transform.

@cindex convolution

@cindex cyclic convolution

@example

#include <rfftw.h>

...

@{

     fftw_real a[M][2*(N/2+1)], b[M][2*(N/2+1)], c[M][N];

     fftw_complex *A, *B, C[M][N/2+1];

     rfftwnd_plan p, pinv;

     fftw_real scale = 1.0 / (M * N);

     int i, j;

     ...

     p    = rfftw2d_create_plan(M, N, FFTW_REAL_TO_COMPLEX,