快速傅里叶变换-FFTW库的使用-参考和翻译官方文档

Complex One-Dimensional DFTs

Plan：为实现意外结果的最佳方法而烦恼。 [Ambrose Bierce, The Enlarged Devil’s Dictionary。] FFTW 计算大小为 N 的一维 DFT 的基本用法很简单，它通常看起来像这样的代码：

#include <fftw3.h>
...
{fftw_complex *in, *out;fftw_plan p;
...in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
...fftw_execute(p); /* repeat as needed */
...fftw_destroy_plan(p);fftw_free(in); fftw_free(out);
}

您必须将此代码与 fftw3 库链接。在 Unix 系统上，使用 -lfftw3 -lm 链接。示例代码首先分配输入和输出数组。您可以以任何您喜欢的方式分配它们，但我们建议使用 fftw_malloc，它的行为类似于 malloc，只是它在 SIMD 指令（如 SSE 和 Altivec）可用时正确对齐数组（参见第 3.1 节 [SIMD 对齐和 fftw malloc ]，第 15 页）。 [或者，我们提供了一个方便的包装函数 fftw_alloc_complex(N)，它具有相同的效果。] 数据是 fftw_complex 类型的数组，默认情况下是由实数 (in[i] [0]) 组成的 double[2] 和复数的虚数 (in[i] [1]) 部分。下一步是创建一个计划，这是一个包含 FFTW 计算 FFT 所需的所有数据的对象。此函数创建计划：

fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);

第一个参数 n 是您尝试计算的变换的大小。大小 n 可以是任何正整数，但作为小因子乘积的大小转换效率最高（尽管素数大小仍使用 O(n log n) 算法）。接下来的两个参数是指向转换的输入和输出数组的指针。这些指针可以相等，表示就地转换。第四个参数 sign 可以是 FFTW_FORWARD (-1) 或 FFTW_BACKWARD (+1)，表示你感兴趣的变换方向；从技术上讲，它是变换中指数的符号。 flags 参数通常是 FFTW_MEASURE 或 FFTW_ESTIMATE。 FFTW_MEASURE 指示 FFTW 运行并测量几个 FFT 的执行时间，以便找到计算大小为 n 的变换的最佳方法。此过程需要一些时间（通常为几秒钟），具体取决于您的机器和转换的大小。相反，FFTW_ESTIMATE 不运行任何计算，只是构建一个可能不是最佳的合理计划。简而言之，如果您的程序执行许多相同大小的转换并且初始化时间并不重要，请使用 FFTW_MEASURE；否则使用估计。您必须在初始化输入之前创建计划，因为 FFTW_MEASURE 会覆盖输入/输出数组。（从技术上讲，FFTW_ESTIMATE 不会触及您的数组，但您应该始终先创建计划以确保确定。）一旦创建了计划，您可以根据需要多次使用它来对指定的输入/输出数组进行转换，通过 fftw_execute(plan) 计算实际转换：

void fftw_execute(const fftw_plan plan);

DFT 结果按顺序存储在数组 out 中，零频率 (DC) 分量存储在 out[0] 中。如果 in != out，则变换是不合适的，并且输入数组 in 不会被修改。否则，输入数组将被转换覆盖。如果要转换相同大小的不同数组，可以使用 fftw_plan_dft_1d 创建一个新计划，如果可能，FFTW 会自动重用先前计划中的信息。或者，如果您小心的话，您可以使用“guru”界面将给定计划应用于不同的阵列。请参阅第 4 章 [FFTW 参考]，第 21 页。完成计划后，您可以通过调用 fftw_destroy_plan(plan) 来释放它:

void fftw_destroy_plan(fftw_plan plan);

FFTW 计算未归一化的 DFT。因此，计算前向变换和后向变换（反之亦然）会导致原始数组按 n 缩放。有关 DFT 的定义，请参见第 42 页第 4.8 节 [FFTW 真正计算的内容]。如果您有支持 C99 标准的 C 编译器，例如 gcc，并且您在 < fftw3> 之前#include <complex.h>，则 fftw_complex 是本机双精度复数类型，您可以使用普通算术对其进行操作。否则，FFTW 定义自己的复数类型，与 C99 复数类型位兼容。请参阅第 21 页的第 4.1.1 节 [复数]。（C++ < complex> 模板类也可以通过类型转换使用。）要使用 FFTW 的单精度或长双精度版本，请将 fftw_ 前缀替换为 fftwf_ 或 fftwl_并与 -lfftw3f 或 -lfftw3l 链接，但使用相同的 <fftw3.h> 头文件。除了 FFTW_MEASURE 和 FFTW_ESTIMATE 之外，还有更多的标志。例如，如果您愿意等待更长的时间以获得可能更快的计划，请使用 FFTW_PATIENT（参见第 4 章 [FFTW 参考]，第 21 页）。您还可以保存计划以备将来使用，如第 3.3 节 [Words of Wisdom-Saving Plans]，第 18 页所述.

Complex Multi-Dimensional DFTs

多维变换的工作方式与一维变换大致相同：分配 fftw_complex 数组（最好使用 fftw_malloc），创建 fftw_plan，使用 fftw_execute(plan) 执行任意多次，然后使用 fftw_destroy_plan 进行清理（计划）（和 fftw_free）。 FFTW 提供了两个用于创建 2d 和 3d 变换计划的例程，以及一个用于创建任意维度计划的例程。 2d 和 3d 例程具有以下签名：

//二维
fftw_plan fftw_plan_dft_2d(int n0, int n1, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);//三维
fftw_plan fftw_plan_dft_3d(int n0, int n1, int n2, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);

这些例程分别为 n 0 × n 1 n0 \times n1 n0×n1 二维 (2d) 变换和 n 0 × n 1 × n 2 n0 \times n1 \times n2 n0×n1×n2 3d变换创建计划。所有这些转换都以 C 标准的行优先顺序对连续数组进行操作，因此最后一个维度具有数组中变化最快的索引。这种布局在第 15 页第 3.2 节 [多维数组格式] 中有进一步描述。FFTW 还可以计算更高维的变换。为了避免混淆“维度”一词的各种含义，我们使用术语秩来表示数组中独立索引的数量。例如，我们说 2d 变换具有 2 秩，3d 变换秩为 3，依此类推。您可以通过以下函数计划任意秩的变换：

fftw_plan fftw_plan_dft(int rank, const int *n，fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);

这里，n 是指向数组 n[rank] 的指针，表示 n [ 0 ] × n [ 1 ] × ⋯ × n [ r a n k − 1 ] n[0] \times n[1] \times \cdots \times n[rank-1] n[0]×n[1]×⋯×n[rank−1]变换。因此，例如，调用

fftw_plan_dft_2d(n0, n1, in, out, sign, flags);

相当于下面的代码片段：

int n[2];
n[0] = n0;
n[1] = n1;
fftw_plan_dft(2, n, in, out, sign, flags);

然而，fftw_plan_dft 不限于 2d 和 3d 变换，但它可以规划任意秩的变换。您可能已经注意到，到目前为止描述的所有计划程序例程都具有重叠的功能。例如，您可以通过使用秩为 1 或 2 的 fftw_plan_dft 或甚至通过调用 n0 和/或 n1 等于 1 的 fftw_plan_dft_3d 来规划 1d 或 2d 变换（不会降低效率）。这种模式继续存在，FFTW 的规划例程通常形成“部分顺序”，接口序列具有严格增加的通用性，但相应地更复数。

fftw_plan_dft 是我们在本教程中描述的最通用的复数 DFT 例程，但也有高级和大师接口，它们允许将多个/跨步转换有效地组合成单个 FFTW 计划，转换更大的多维数组的子集 , 和/或处理更通用的复数格式。有关详细信息，请参阅第 4 章 [FFTW 参考]，第 21 页。

One-Dimensional DFTs of Real Data

在许多实际应用中，输入数据 in[i] 是纯实数，在这种情况下，DFT 输出满足“Hermitian”冗余：out[i] 是 out[n-i] 的共轭。可以利用这些情况来在速度和内存使用方面实现大约两倍的改进。

为了换取这些速度和空间优势，用户牺牲了 FFTW 复数变换的一些简单性。首先，输入和输出数组的大小和类型不同：输入是n个实数，而输出是n/2+1个复数（非冗余输出）；这也需要对输入数组进行轻微的“填充”以进行就地转换。其次，默认情况下，逆变换（复数到实数）具有覆盖其输入数组的副作用。这些不便都不应该给用户带来严重的问题，但重要的是要意识到它们。执行实数据转换的例程与复数转换的例程几乎相同：分配 double 和/或 fftw_complex 数组（最好使用 fftw_malloc 或 fftw_alloc_complex），创建 fftw_plan，使用 fftw_execute (plan)执行任意多次，并使用 fftw_destroy_plan(plan) （和 fftw_ free）进行清理。唯一的区别是输入（或输出）是 double 类型，并且有新的例程来创建计划。在一维中：

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);

用于complex-Hermitian输出 (r2c) 的实数输入和complex-Hermitian输入到实数输出 (c2r) 的变换。与复数的 DFT 规划器不同，没有sign参数。相反，r2c DFTs 始终为 FFTW_FORWARD，而 c2r DFT 始终为 FFTW_BACKWARD。（对于单/长双精度 fftwf 和 fftwl，double 应分别替换为 float 和 long double。）这里，n 是 DFT 的“逻辑”大小，不一定是数组的物理大小。特别是，实数 (double) 数组有 n 个元素，而复数 (fftw_complex) 数组有 n/2+1 个元素（除法向下舍入）。对于就地变换，输入和输出别名为同一个数组，该数组必须足够大以容纳两者；因此，真正的数组实际上将有 2*(n/2+1) 个元素，其中第一个 n 之外的元素是未使用的填充。（请注意，这与将变换“零填充”到更大长度的概念非常不同，后者通过实际添加新输入数据来改变 DFT 的逻辑大小。）复数数组的第 k 个元素完全相同作为相应复数 DFT 的第 k 个元素。支持所有正数 n；小因子的乘积是最有效的，但即使是素数大小也使用 O(n log n) 算法。如上所述，c2r 变换会破坏其输入数组，即使是异地变换也是如此。如有必要，可以通过在标志中包含 FFTW_PRESERVE_INPUT 来防止这种情况，但不幸的是会牺牲一些性能。多维实 DFT（下一节）目前也不支持此标志。熟悉真实数据的 DFT 的读者会记得，复数输出的第 0（“DC”）和第 n/2（“奈奎斯特”频率，当 n 为偶数时）元素是纯实数。因此，一些实现将奈奎斯特元素存储在 DC 虚部所在的位置，以使输入和输出数组大小相同。然而，这种打包并不能很好地推广到多维变换，而且在任何情况下节省的空间都是微不足道的； FFTW 不支持。一维 r2c 和 c2r DFT 的替代接口可以在“r2r”接口中找到（参见第 2.5.1 节 [半复数格式 DFT]，第 11 页），“半复数”格式输出大小相同（和类型）作为输入数组。该接口虽然对多维变换不是很有用，但有时可能会产生更好的性能。

Multi-Dimensional DFTs of Real Data

真实数据的多维 DFT 使用以下planner例程：

fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1, double *in, fftw_complex *out, unsigned flags);fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, int n2, double *in, fftw_complex *out, unsigned flags);fftw_plan fftw_plan_dft_r2c(int rank, const int *n, double *in, fftw_complex *out, unsigned flags);

以及交换了输入/输出类型的相应 c2r 例程。这些例程的工作方式与它们的复数类似物相似，除了这里的复数输出数组被大致切成两半并且实际数组需要填充以进行就地变换（如上面的 1d 中所示）。

和之前一样，n 是数组的逻辑大小，这对复数数组格式的影响值得特别注意。假设真实数据的维度为 n 0 × n 1 × n 2 × ⋯ × n d − 1 n0 \times n1 \times n2 \times \cdots \times n_{d−1} n0×n1×n2×⋯×nd−1（按行优先顺序）。然后，在 r2c 变换之后，输出是一个 n 0 × n 1 × n 2 × ⋯ × ( n d − 1 / 2 + 1 ) n0 \times n1 \times n2 \times \cdots \times (n_{d−1}/2 + 1) n0×n1×n2×⋯×(nd−1/2+1)个 fftw_complex 值的数组，以行优先顺序排列，对应于略多于一半的输出对应的复数 DFT。（除法四舍五入。）数据的排序在其他方面与复数 DFT 情况完全相同。

对于异地变换(out-of-place)，这就是故事的结尾：真实数据存储为大小为 n 0 × n 1 × n 2 × ⋯ × n d − 1 n0 \times n1 \times n2 \times \cdots \times n_{d−1} n0×n1×n2×⋯×nd−1 的行主数组，复数数据存储为行- 大小为 n 0 × n 1 × n 2 × ⋯ × ( n d − 1 / 2 + 1 ) n0 \times n1 \times n2 \times \cdots \times (n_{d−1}/2 + 1) n0×n1×n2×⋯×(nd−1/2+1)的主数组。

然而，对于就地转换(in-place)，实数数组的额外填充是必要的，因为复数数组比实数数组大，并且两个数组共享相同的内存位置。因此，对于就地转换，实际数据数组的最后一维必须用额外的值填充以适应复数数据的大小——如果最后一维是偶数则两个值，如果最后一维是奇数则一个值。也就是说，实际数据的最后一维物理上必须包含 2 ( n d − 1 / 2 + 1 ) 2(n_{d−1}/2 + 1) 2(nd−1/2+1) 个双精度值（正好足以容纳复数数据）。然而，这个物理数组大小不会改变逻辑数组大小——只有 n d − 1 n_{d-1} nd−1 个值实际存储在最后一维中，而 n d − 1 n_{d-1} nd−1 是传递给计划创建例程的最后一维。

例如，考虑一个大小 n 0 × n 1 n0 \times n1 n0×n1 的二维实数数组的变换。 r2c 变换的输出是一个大小为 n 0 × ( n 1 / 2 + 1 ) n0 \times (n1/2+1) n0×(n1/2+1) 的二维复数数组，由于输出中的冗余，其中 y 维几乎被削减了一半。因为 fftw_complex 的大小是 double 的两倍，所以输出数组略大于输入数组。因此，如果我们想就地计算变换，我们必须填充输入数组，使其大小为 n 0 × 2 ∗ ( n 1 / 2 + 1 ) n0 \times 2*(n1/2+1) n0×2∗(n1/2+1)。如果 n1 是偶数，则在每行的末尾有两个填充元素（不需要初始化，因为它们仅用于输出）。

2d n x × n y nx \times ny nx×ny real-to-complex 变换的数据布局图示:

图 2.1 描述了刚刚描述的输入和输出数组，用于异地和就地变换（箭头表示连续的内存位置）：

这些变换是非标准化的，因此 r2c 后跟 c2r 变换（反之亦然）将导致原始数据按真实数据元素的数量进行缩放，即真实数据的（逻辑）维度的乘积。

（因为最后一个维度被特殊处理，如果它等于 1，则变换不等价于低维 r2c/c2r 变换。在这种情况下，最后一个复数维度的大小也为 1（=1/2+1） , 并没有比复数的变换有任何优势。)

More DFTs of Real Data

FFTW 通过一个统一的 r2r（real-to-real）接口支持其他几种变换类型，之所以这么称呼是因为它采用一个实数（双）数组并输出一个相同大小的实数数组。这些 r2r 变换目前分为三类：半复数格式的实数输入和复数Hermitian输出的 DFT，具有偶数/奇数对称性的实数输入的 DFT（也称为离散余弦/正弦变换，DCT/DST）和离散Hartley变换 (DHT)，

以下各节将更详细地描述所有内容。 r2r 转换遵循现在熟悉的创建 fftw_plan ，使用 fftw_execute(plan) 执行它，并使用 fftw_destroy_plan(plan) 销毁它。此外，所有 r2r 转换共享相同的规划器接口：

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);

就像复杂的 DFT 一样，这些计划以行优先顺序对连续数组进行 1d/2d/3d/多维变换，将（实际）输入转换为相同大小的输出，其中 n 指定数组的物理尺寸。支持所有正数 n（FFTW_REDFT00 类型的 n=1 除外，在下面的实数小节中注明）；小因子的乘积是最有效的（对 FFTW_REDFT00 和 FFTW_RODFT00 种类进行因子分解 n-1 和 n+1，如下所述），但即使对于素数大小也使用 O(n log n) 算法。

每个维度都有一个类型为 fftw_r2r_kind 的参数，指定用于该维度的 r2r 变换的种类。（在 fftw_plan_r2r 的情况下，这是一个数组 kind[rank]，其中 kind[i] 是维度 n[i] 的变换种类。）种类可以是一组预定义常量之一，在以下小节中定义.

换句话说，FFTW 计算指定 r2r 变换在每个维度上的可分离乘积，可用于例如对于具有混合边界条件的偏微分方程。（对于某些 r2r 类型，特别是半复数 DFT 和 DHT，这样的可分离乘积在多维方面有些问题，但是，如下所述。）

在当前版本的 FFTW 中，除了半复数类型之外的所有 r2r 变换是通过半复数变换的预处理或后处理来计算的，因此它们的速度不如它们可能的快。然而，由于大多数其他通用 DCT/DST 代码采用类似的算法，因此 FFTW 的实现至少应该提供具有竞争力的性能。

参考和翻译：

FFTW Home Page