Node:One-Dimensional DFTs of Real Data, Next:, Previous:Complex Multi-Dimensional DFTs, Up:Tutorial



2.3 One-Dimensional DFTs of Real Data

In many practical applications, the input data in[i] are purely real numbers, in which case the DFT output satisfies the "Hermitian" redundancy: out[i] is the conjugate of out[n-i]. It is possible to take advantage of these circumstances in order to achieve roughly a factor of two improvement in both speed and memory usage.

In exchange for these speed and space advantages, the user sacrifices some of the simplicity of FFTW's complex transforms. First of all, the input and output arrays are of different sizes and types: the input is n real numbers, while the output is n/2+1 complex numbers (the non-redundant outputs); this also requires slight "padding" of the input array for in-place transforms. Second, the inverse transform (complex to real) has the side-effect of destroying its input array, by default. Neither of these inconveniences should pose a serious problem for users, but it is important to be aware of them.

The routines to perform real-data transforms are almost the same as those for complex transforms: you allocate arrays of double and/or fftw_complex (preferably using fftw_malloc), create an fftw_plan, execute it as many times as you want with fftw_execute(plan), and clean up with fftw_destroy_plan(plan) (and fftw_free). The only differences are that the input (or output) is of type double and there are new routines to create the plan. In one dimension:

fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out,
                               unsigned flags);
fftw_plan fftw_plan_dft_c2r_1d(int n, fftw_complex *in, double *out,
                               unsigned flags);

for the real input to complex-Hermitian output (r2c) and complex-Hermitian input to real output (c2r) transforms. Unlike the complex DFT planner, there is no sign argument. Instead, r2c DFTs are always FFTW_FORWARD and c2r DFTs are always FFTW_BACKWARD. (For single/long-double precision fftwf and fftwl, double should be replaced by float and long double, respectively.)

Here, n is the "logical" size of the DFT, not necessarily the physical size of the array. In particular, the real (double) array has n elements, while the complex (fftw_complex) array has n/2+1 elements (where the division is rounded down). For an in-place transform, in and out are aliased to the same array, which must be big enough to hold both; so, the real array would actually have 2*(n/2+1) elements, where the elements beyond the first n are unused padding. The kth element of the complex array is exactly the same as the kth element of the corresponding complex DFT. All positive n are supported; products of small factors are most efficient, but an O(n log n) algorithm is used even for prime sizes.

As noted above, the c2r transform destroys its input array even for out-of-place transforms. This can be prevented, if necessary, by including FFTW_PRESERVE_INPUT in the flags, with unfortunately some sacrifice in performance. This flag is also not currently supported for multi-dimensional real DFTs (next section).

Readers familiar with DFTs of real data will recall that the 0th (the "DC") and n/2-th (the "Nyquist" frequency, when n is even) elements of the complex output are purely real. Some implementations therefore store the Nyquist element where the DC imaginary part would go, in order to make the input and output arrays the same size. Such packing, however, does not generalize well to multi-dimensional transforms, and the space savings are miniscule in any case; FFTW does not support it.

An alternate interface for one-dimensional r2c and c2r DFTs can be found in the r2r interface (see The Halfcomplex-format DFT), with "halfcomplex"-format output that is the same size (and type) as the input array. That interface, although it is not very useful for multi-dimensional transforms, may sometimes yield better performance.