FFTW 3.3.8: Real even/odd DFTs (cosine/sine transforms)

The Fourier transform of a real-even function f(-x) = f(x) is Chris@82: real-even, and i times the Fourier transform of a real-odd Chris@82: function f(-x) = -f(x) is real-odd. Similar results hold for a Chris@82: discrete Fourier transform, and thus for these symmetries the need for Chris@82: complex inputs/outputs is entirely eliminated. Moreover, one gains a Chris@82: factor of two in speed/space from the fact that the data are real, and Chris@82: an additional factor of two from the even/odd symmetry: only the Chris@82: non-redundant (first) half of the array need be stored. The result is Chris@82: the real-even DFT (REDFT) and the real-odd DFT (RODFT), also Chris@82: known as the discrete cosine and sine transforms (DCT and Chris@82: DST), respectively. Chris@82: Chris@82: Chris@82: Chris@82: Chris@82: Chris@82: Chris@82: Chris@82: Chris@82:

(In this section, we describe the 1d transforms; multi-dimensional Chris@82: transforms are just a separable product of these transforms operating Chris@82: along each dimension.) Chris@82:

Because of the discrete sampling, one has an additional choice: is the Chris@82: data even/odd around a sampling point, or around the point halfway Chris@82: between two samples? The latter corresponds to shifting the Chris@82: samples by half an interval, and gives rise to several transform Chris@82: variants denoted by REDFTab and RODFTab: a and Chris@82: b are 0 or 1, and indicate whether the input Chris@82: (a) and/or output (b) are shifted by half a sample Chris@82: (1 means it is shifted). These are also known as types I-IV of Chris@82: the DCT and DST, and all four types are supported by FFTW’s r2r Chris@82: interface.³ Chris@82:

The r2r kinds for the various REDFT and RODFT types supported by FFTW, Chris@82: along with the boundary conditions at both ends of the input Chris@82: array (n real numbers in[j=0..n-1]), are: Chris@82:

Note that these symmetries apply to the “logical” array being Chris@82: transformed; there are no constraints on your physical input Chris@82: data. So, for example, if you specify a size-5 REDFT00 (DCT-I) of the Chris@82: data abcde, it corresponds to the DFT of the logical even array Chris@82: abcdedcb of size 8. A size-4 REDFT10 (DCT-II) of the data Chris@82: abcd corresponds to the size-8 logical DFT of the even array Chris@82: abcddcba, shifted by half a sample. Chris@82:

All of these transforms are invertible. The inverse of R*DFT00 is Chris@82: R*DFT00; of R*DFT10 is R*DFT01 and vice versa (these are often called Chris@82: simply “the” DCT and IDCT, respectively); and of R*DFT11 is R*DFT11. Chris@82: However, the transforms computed by FFTW are unnormalized, exactly Chris@82: like the corresponding real and complex DFTs, so computing a transform Chris@82: followed by its inverse yields the original array scaled by N, Chris@82: where N is the logical DFT size. For REDFT00, Chris@82: N=2(n-1); for RODFT00, N=2(n+1); otherwise, N=2n. Chris@82: Chris@82: Chris@82:

Note that the boundary conditions of the transform output array are Chris@82: given by the input boundary conditions of the inverse transform. Chris@82: Thus, the above transforms are all inequivalent in terms of Chris@82: input/output boundary conditions, even neglecting the 0.5 shift Chris@82: difference. Chris@82:

FFTW is most efficient when N is a product of small factors; note Chris@82: that this differs from the factorization of the physical size Chris@82: n for REDFT00 and RODFT00! There is another oddity: n=1 Chris@82: REDFT00 transforms correspond to N=0, and so are not Chris@82: defined (the planner will return NULL). Otherwise, any positive Chris@82: n is supported. Chris@82:

For the precise mathematical definitions of these transforms as used by Chris@82: FFTW, see What FFTW Really Computes. (For people accustomed to Chris@82: the DCT/DST, FFTW’s definitions have a coefficient of 2 in front Chris@82: of the cos/sin functions so that they correspond precisely to an Chris@82: even/odd DFT of size N. Some authors also include additional Chris@82: multiplicative factors of Chris@82: √2 Chris@82: for selected inputs and outputs; this makes Chris@82: the transform orthogonal, but sacrifices the direct equivalence to a Chris@82: symmetric DFT.) Chris@82:

Which type do you need?

Since the required flavor of even/odd DFT depends upon your problem, Chris@82: you are the best judge of this choice, but we can make a few comments Chris@82: on relative efficiency to help you in your selection. In particular, Chris@82: R*DFT01 and R*DFT10 tend to be slightly faster than R*DFT11 Chris@82: (especially for odd sizes), while the R*DFT00 transforms are sometimes Chris@82: significantly slower (especially for even sizes).⁴ Chris@82:

Thus, if only the boundary conditions on the transform inputs are Chris@82: specified, we generally recommend R*DFT10 over R*DFT00 and R*DFT01 over Chris@82: R*DFT11 (unless the half-sample shift or the self-inverse property is Chris@82: significant for your problem). Chris@82:

If performance is important to you and you are using only small sizes Chris@82: (say n<200), e.g. for multi-dimensional transforms, then you Chris@82: might consider generating hard-coded transforms of those sizes and types Chris@82: that you are interested in (see Generating your own code). Chris@82:

We are interested in hearing what types of symmetric transforms you find Chris@82: most useful. Chris@82:

Chris@82:

Footnotes

Chris@82: Chris@82:

(3)

Chris@82:

There are also type V-VIII transforms, which Chris@82: correspond to a logical DFT of odd size N, independent of Chris@82: whether the physical size n is odd, but we do not support these Chris@82: variants.

Chris@82:

(4)

Chris@82:

R*DFT00 is Chris@82: sometimes slower in FFTW because we discovered that the standard Chris@82: algorithm for computing this by a pre/post-processed real DFT—the Chris@82: algorithm used in FFTPACK, Numerical Recipes, and other sources for Chris@82: decades now—has serious numerical problems: it already loses several Chris@82: decimal places of accuracy for 16k sizes. There seem to be only two Chris@82: alternatives in the literature that do not suffer similarly: a Chris@82: recursive decomposition into smaller DCTs, which would require a large Chris@82: set of codelets for efficiency and generality, or sacrificing a factor of Chris@82: 2 Chris@82: in speed to use a real DFT of twice the size. We currently Chris@82: employ the latter technique for general n, as well as a limited Chris@82: form of the former method: a split-radix decomposition when n Chris@82: is odd (N a multiple of 4). For N containing many Chris@82: factors of 2, the split-radix method seems to recover most of the Chris@82: speed of the standard algorithm without the accuracy tradeoff.

Chris@82:

2.5.2 Real even/odd DFTs (cosine/sine transforms)

Which type do you need?

Footnotes

(3)

(4)