FFTW3详解：快速傅里叶变换指南

需积分: 31 145 浏览量更新于2024-07-23 收藏 645KB PDF 举报

"FFTW3是著名的开源快速傅里叶变换（FFT）库，用于执行一维、多维以及各种实数和复数数据类型的离散傅里叶变换。本手册详细介绍了版本3.3.3的相关内容，由Matteo Frigo和Steven G. Johnson编写，并且得到了麻省理工学院的版权许可。用户可以自由复制和分发手册，只要保留原始的版权信息。对于修改版和翻译版，也需遵循相同的规定，并保持许可通知的一致性。" FFTW3是一个高效的库，其核心功能是执行傅里叶变换，这是许多科学计算和信号处理任务中的关键步骤。傅里叶变换是一种将信号或函数从时域（或空间域）转换到频域的数学工具，它在图像处理、音频分析、物理模拟等领域有广泛的应用。本手册分为多个部分，首先从基础介绍开始： 1. **介绍**：这部分可能涵盖FFTW3的基本概念、设计原则和使用场景，帮助读者理解FFTW3的核心价值和适用范围。 2. **教程**：教程部分深入浅出地展示了如何使用FFTW3进行不同类型的傅里叶变换： - **复数一维DFTs**：描述了对一维复数序列进行离散傅里叶变换的步骤。 - **复数多维DFTs**：扩展到更高维度，适用于处理二维和三维数据，如图像和多维信号。 - **实数一维DFTs**：讨论了一种优化的算法，针对仅包含实数的数据，可以更高效地计算一维变换。 - **实数多维DFTs**：在多维实数数据上的应用，通常用于处理网格或阵列数据。 - **更多实数DFTs**：可能包括其他特定的实数数据处理技巧，如半复共轭对称性，以进一步提高效率。这些章节不仅会讲解基本的API调用，还可能包括计划、执行和清理变换的流程，以及如何利用FFTW3的智能内存管理和其他优化特性。 FFTW3支持多种编译选项和接口，包括C、C++、Fortran等，手册中可能会详细解释如何在不同环境中配置和使用FFTW3。此外，性能调优、并行计算支持（例如OpenMP和MPI）以及与其他库的集成也是手册的重要部分。通过深入学习FFTW3手册，用户将能够熟练掌握这个强大的工具，从而在自己的项目中高效地实现傅里叶变换。无论你是初次接触FFT，还是寻求提高现有代码性能的方法，这份详尽的指南都将提供必要的指导。

10 FFTW 3.3.3

2.5 More DFTs of Real Data

FFTW supports several other transform types via a uniﬁed r2r (real-to-real) interface, so

called because it takes a real (double) array and outputs a real array of the same size.

These r2r transforms currently fall into three categories: DFTs of real input and complex-

Hermitian output in halfcomplex format, DFTs of real input with even/odd symmetry

(a.k.a. discrete cosine/sine transforms, DCTs/DSTs), and discrete Hartley transforms

(DHTs), all described in more detail by the following sections.

The r2r transforms follow the by now familiar interface of creating an fftw_plan, exe-

cuting it with fftw_execute(plan), and destroying it with fftw_destroy_plan(plan).

Furthermore, all r2r transforms share the same planner interface:

fftw_plan fftw_plan_r2r_1d(int n, double *in, double *out,

fftw_r2r_kind kind, unsigned flags);

fftw_plan fftw_plan_r2r_2d(int n0, int n1, double *in, double *out,

fftw_r2r_kind kind0, fftw_r2r_kind kind1,

unsigned flags);

fftw_plan fftw_plan_r2r_3d(int n0, int n1, int n2,

double *in, double *out,

fftw_r2r_kind kind0,

fftw_r2r_kind kind1,

fftw_r2r_kind kind2,

unsigned flags);

fftw_plan fftw_plan_r2r(int rank, const int *n, double *in, double *out,

const fftw_r2r_kind *kind, unsigned flags);

Just as for the complex DFT, these plan 1d/2d/3d/multi-dimensional transforms for con-

tiguous arrays in row-major order, transforming (real) input to output of the same size,

where n speciﬁes the physical dimensions of the arrays. All positive n are supported (with

the exception of n=1 for the FFTW_REDFT00 kind, noted in the real-even subsection below);

products of small factors are most eﬃcient (factorizing n-1 and n+1 for FFTW_REDFT00 and

FFTW_RODFT00 kinds, described below), but an O(n log n) algorithm is used even for prime

sizes.

Each dimension has a kind parameter, of type fftw_r2r_kind, specifying the kind of r2r

transform to be used for that dimension. (In the case of fftw_plan_r2r, this is an array

kind[rank] where kind[i] is the transform kind for the dimension n[i].) The kind can

be one of a set of predeﬁned constants, deﬁned in the following subsections.

In other words, FFTW computes the separable product of the speciﬁed r2r transforms

over each dimension, which can be used e.g. for partial diﬀerential equations with mixed

boundary conditions. (For some r2r kinds, notably the halfcomplex DFT and the DHT,

such a separable product is somewhat problematic in more than one dimension, however,

as is described below.)

In the current version of FFTW, all r2r transforms except for the halfcomplex type are com-

puted via pre- or post-processing of halfcomplex transforms, and they are therefore not as

fast as they could be. Since most other general DCT/DST codes employ a similar algorithm,

however, FFTW’s implementation should provide at least competitive performance.

Chapter 2: Tutorial 11

2.5.1 The Halfcomplex-format DFT

An r2r kind of FFTW_R2HC (r2hc) corresponds to an r2c DFT (see Section 2.3

[One-Dimensional DFTs of Real Data], page 6) but with “halfcomplex” format output,

and may sometimes be faster and/or more convenient than the latter. The inverse hc2r

transform is of kind FFTW_HC2R. This consists of the non-redundant half of the complex

output for a 1d real-input DFT of size n, stored as a sequence of n real numbers (double)

in the format:

, r

, . . . , r

n/2

, i

(n+1)/2−1

, . . . , i

, i

Here, r

is the real part of the kth output, and i

is the imaginary part. (Division by 2

is rounded down.) For a halfcomplex array hc[n], the kth component thus has its real

part in hc[k] and its imaginary part in hc[n-k], with the exception of k == 0 or n/2 (the

latter only if n is even)—in these two cases, the imaginary part is zero due to symmetries

of the real-input DFT, and is not stored. Thus, the r2hc transform of n real values is a

halfcomplex array of length n, and vice versa for hc2r.

Aside from the diﬀering format, the output of FFTW_R2HC/FFTW_HC2R is otherwise exactly

the same as for the corresponding 1d r2c/c2r transform (i.e. FFTW_FORWARD/FFTW_BACKWARD

transforms, respectively). Recall that these transforms are unnormalized, so r2hc followed

by hc2r will result in the original data multiplied by n. Furthermore, like the c2r transform,

an out-of-place hc2r transform will destroy its input array.

Although these halfcomplex transforms can be used with the multi-dimensional r2r interface,

the interpretation of such a separable product of transforms along each dimension is prob-

lematic. For example, consider a two-dimensional n0 by n1, r2hc by r2hc transform planned

by fftw_plan_r2r_2d(n0, n1, in, out, FFTW_R2HC, FFTW_R2HC, FFTW_MEASURE). Con-

ceptually, FFTW ﬁrst transforms the rows (of size n1) to produce halfcomplex rows, and

then transforms the columns (of size n0). Half of these column transforms, however, are

of imaginary parts, and should therefore be multiplied by i and combined with the r2hc

transforms of the real columns to produce the 2d DFT amplitudes; FFTW’s r2r transform

does not perform this combination for you. Thus, if a multi-dimensional real-input/output

DFT is required, we recommend using the ordinary r2c/c2r interface (see Section 2.4 [Multi-

Dimensional DFTs of Real Data], page 7).

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

The Fourier transform of a real-even function f(−x) = f(x) is real-even, and i times

the Fourier transform of a real-odd function f(−x) = −f(x) is real-odd. Similar results

hold for a discrete Fourier transform, and thus for these symmetries the need for complex

inputs/outputs is entirely eliminated. Moreover, one gains a factor of two in speed/space

from the fact that the data are real, and an additional factor of two from the even/odd

symmetry: only the non-redundant (ﬁrst) half of the array need be stored. The result is

the real-even DFT (REDFT) and the real-odd DFT (RODFT), also known as the discrete

cosine and sine transforms (DCT and DST), respectively.

(In this section, we describe the 1d transforms; multi-dimensional transforms are just a

separable product of these transforms operating along each dimension.)

12 FFTW 3.3.3

Because of the discrete sampling, one has an additional choice: is the data even/odd around

a sampling point, or around the point halfway between two samples? The latter corresponds

to shifting the samples by half an interval, and gives rise to several transform variants

denoted by REDFTab and RODFTab: a and b are 0 or 1, and indicate whether the input

(a) and/or output (b) are shifted by half a sample (1 means it is shifted). These are also

known as types I-IV of the DCT and DST, and all four types are supported by FFTW’s

r2r interface.

The r2r kinds for the various REDFT and RODFT types supported by FFTW, along with

the boundary conditions at both ends of the input array (n real numbers in[j=0..n-1]),

are:

• FFTW_REDFT00 (DCT-I): even around j = 0 and even around j = n − 1.

• FFTW_REDFT10 (DCT-II, “the” DCT): even around j = −0.5 and even around j =

n − 0.5.

• FFTW_REDFT01 (DCT-III, “the” IDCT): even around j = 0 and odd around j = n.

• FFTW_REDFT11 (DCT-IV): even around j = −0.5 and odd around j = n − 0.5.

• FFTW_RODFT00 (DST-I): odd around j = −1 and odd around j = n.

• FFTW_RODFT10 (DST-II): odd around j = −0.5 and odd around j = n − 0.5.

• FFTW_RODFT01 (DST-III): odd around j = −1 and even around j = n − 1.

• FFTW_RODFT11 (DST-IV): odd around j = −0.5 and even around j = n − 0.5.

Note that these symmetries apply to the “logical” array being transformed; there are no

constraints on your physical input data. So, for example, if you specify a size-5 REDFT00

(DCT-I) of the data abcde, it corresponds to the DFT of the logical even array abcdedcb of

size 8. A size-4 REDFT10 (DCT-II) of the data abcd corresponds to the size-8 logical DFT

of the even array abcddcba, shifted by half a sample.

All of these transforms are invertible. The inverse of R*DFT00 is R*DFT00; of R*DFT10 is

R*DFT01 and vice versa (these are often called simply “the” DCT and IDCT, respectively);

and of R*DFT11 is R*DFT11. However, the transforms computed by FFTW are unnor-

malized, exactly like the corresponding real and complex DFTs, so computing a transform

followed by its inverse yields the original array scaled by N , where N is the logical DFT

size. For REDFT00, N = 2(n − 1); for RODFT00, N = 2(n + 1); otherwise, N = 2n.

Note that the boundary conditions of the transform output array are given by the input

boundary conditions of the inverse transform. Thus, the above transforms are all inequiva-

lent in terms of input/output boundary conditions, even neglecting the 0.5 shift diﬀerence.

FFTW is most eﬃcient when N is a product of small factors; note that this diﬀers from the

factorization of the physical size n for REDFT00 and RODFT00! There is another oddity:

n=1 REDFT00 transforms correspond to N = 0, and so are not deﬁned (the planner will

return NULL). Otherwise, any positive n is supported.

For the precise mathematical deﬁnitions of these transforms as used by FFTW, see

Section 4.8 [What FFTW Really Computes], page 42. (For people accustomed to the

There are also type V-VIII transforms, which correspond to a logical DFT of odd size N , independent of

whether the physical size n is odd, but we do not support these variants.

Chapter 2: Tutorial 13

DCT/DST, FFTW’s deﬁnitions have a coeﬃcient of 2 in front of the cos/sin functions so

that they correspond precisely to an even/odd DFT of size N. Some authors also include

additional multiplicative factors of

√

2 for selected inputs and outputs; this makes the

transform orthogonal, but sacriﬁces the direct equivalence to a symmetric DFT.)

Which type do you need?

Since the required ﬂavor of even/odd DFT depends upon your problem, you are the best

judge of this choice, but we can make a few comments on relative eﬃciency to help you

in your selection. In particular, R*DFT01 and R*DFT10 tend to be slightly faster than

R*DFT11 (especially for odd sizes), while the R*DFT00 transforms are sometimes signiﬁ-

cantly slower (especially for even sizes).

Thus, if only the boundary conditions on the transform inputs are speciﬁed, we generally

recommend R*DFT10 over R*DFT00 and R*DFT01 over R*DFT11 (unless the half-sample

shift or the self-inverse property is signiﬁcant for your problem).

If performance is important to you and you are using only small sizes (say n < 200), e.g. for

multi-dimensional transforms, then you might consider generating hard-coded transforms

of those sizes and types that you are interested in (see Section 10.4 [Generating your own

code], page 101).

We are interested in hearing what types of symmetric transforms you ﬁnd most useful.

2.5.3 The Discrete Hartley Transform

If you are planning to use the DHT because you’ve heard that it is “faster” than the DFT

(FFT), stop here. The DHT is not faster than the DFT. That story is an old but enduring

misconception that was debunked in 1987.

The discrete Hartley transform (DHT) is an invertible linear transform closely related to

the DFT. In the DFT, one multiplies each input by cos − i ∗ sin (a complex exponential),

whereas in the DHT each input is multiplied by simply cos+sin. Thus, the DHT transforms

n real numbers to n real numbers, and has the convenient property of being its own inverse.

In FFTW, a DHT (of any positive n) can be speciﬁed by an r2r kind of FFTW_DHT.

Like the DFT, in FFTW the DHT is unnormalized, so computing a DHT of size n followed

by another DHT of the same size will result in the original array multiplied by n.

The DHT was originally proposed as a more eﬃcient alternative to the DFT for real data,

but it was subsequently shown that a specialized DFT (such as FFTW’s r2hc or r2c trans-

forms) could be just as fast. In FFTW, the DHT is actually computed by post-processing

an r2hc transform, so there is ordinarily no reason to prefer it from a performance per-

R*DFT00 is sometimes slower in FFTW because we discovered that the standard algorithm for computing

this by a pre/post-processed real DFT—the algorithm used in FFTPACK, Numerical Recipes, and other

sources for decades now—has serious numerical problems: it already loses several decimal places of accuracy

for 16k sizes. There seem to be only two alternatives in the literature that do not suﬀer similarly: a recursive

decomposition into smaller DCTs, which would require a large set of codelets for eﬃciency and generality,

or sacriﬁcing a factor of ∼ 2 in speed to use a real DFT of twice the size. We currently employ the latter

technique for general n, as well as a limited form of the former method: a split-radix decomposition when

n is odd (N a multiple of 4). For N containing many factors of 2, the split-radix method seems to recover

most of the speed of the standard algorithm without the accuracy tradeoﬀ.

剩余115页未读，继续阅读

wanghailian528

粉丝: 0
资源: 2

FFTW3详解：快速傅里叶变换指南

FFTW中文参考手册

开源傅里叶变换运算库fftw文档fftw3.pdf

快速傅里叶变换（FFTW3)

fftw说明书

fftw3-delphi_fftw3_fftw3delphi_Tested_sevenxw7_

fftw3工具库使用说明

fftw3 manual

FFTW3的API函数说明

fftw3ubantu

fftw3 使用

最新资源