FFTW中文参考
据说FFTW（Fastest Fourier Transform in the West）是世界上最快的FFT。为了详细了解FFTW以及为编程方便，特将用户手册看了一下，并结合手册制作了以下FFTW中文参考。其中大部分是原文重点内容的翻译，并加入了一些注解。
先看一下使用FFTW编程的方法：
#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); // 输入数据in赋值 p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); fftw_execute(p); // 执行变换 ... fftw_destroy_plan(p); fftw_free(in); fftw_free(out); } |
大致是先用fftw_malloc分配输入输出内存，然后输入数据赋值，然后创建变换方案（fftw_plan），然后执行变换（fftw_execute），最后释放资源，还是比较简单的。
1. 数据类型
fftw_complex默认由两个double组成，在内存中顺序排列，实部在 前，虚部在后，即typedef double fftw_complex[2]。FFTW文档指出如果有一个支持C99标准的C编译器（如gcc），可以在#include <fftw3.h>前加入#include <complex.h>，这样一来fftw_complex就被定义为本机复数类型，而且与上述typedef二进制兼容（指内存排列），经 测试不能用在Windows下。C++有一个复数模板类complex<T>，在头文件<complex>下定义。C++标准委
员会最近同意该类的存储方式与C99二进制兼容，即顺序存储，实部在前，虚部在后（见报告WG21/N1388），该解决方案在所有主流标准库实现中都能正确工作。所以实际上可以用complex <double>
来代替fftw_complex，比如有一个复数数组complex<double> *x，则可以将其类型转换后作为参数传递给fftw：reinterpret_cast<fftw_complex*>(x)。测试如下：开 两个数组fftw_complex x1[2]和complex<double> x2[2]，然后赋相同值，在调试模式下可以看到它们的内存排列是相同的。complex<T>类数据赋值的方式不是很直接，必须采用无名对象方式x2[i]
= complex <double>(1,2) 或成员函数方式x2[i].real(1);x2[i].imag(2);不能直接写x2[i].real=1;x2[i].imag=2。 fftw_complex赋值方式比较直接：x1[i][0]=1;x1[i][1]=2。最后，考虑到数据对齐（见后），最好使用 fftw_malloc分配内存，所以可以将其返回的指针转换为complex <double>
*类型使用（比如赋值或读取等），变换时再将其转换为fftw_complex*。
2. 函数接口
fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); |
n为数据个数，可以为任意正整数，但如果为一些小因子的乘积计算起来可以更有效，不过即使n为素数算法仍然能够达到O(nlogn)的复杂度。FFTW对N=2a 3b 5c 7d 11e 13f的变换处理得最好，其中e+f=0/1，其它幂指数可以为任意值。
如果in和out指针相同为原位运算，否则为非原位运算。
sign可以为正变换FFTW_FORWARD(-1)，也可以为逆变换FFTW_BACKWORD(+1)，实际上就是变换公式中指数项的符号。需注意FFTW的逆变换没有除以N，即数据正变换再反变换后是原始数据的N倍。
flags参数一般情况下为FFTW_MEASURE 或 FFTW_ESTIMATE。FFTW_MEASURE表示FFTW会先计算一些FFT并测量所用的时间，以便为大小为n的变换寻找最优的计算方法。依据 机器配置和变换的大小（n），这个过程耗费约数秒（时钟clock精度）。FFTW_ESTIMATE则相反，它直接构造一个合理的但可能是次最优的方 案。总体来说，如果你的程序需要进行大量相同大小的FFT，并且初始化时间不重要，可以使用FFTW_MEASURE，否则应使用
FFTW_ESTIMATE。FFTW_MEASURE模式下in和out数组中的值会被覆盖，所以该方式应该在用户初始化输入数据in之前完成。
不知道上述说法是不是这个意思：先用FFTW_MEASURE模式自动选最优方案，速度较慢；然后使用该模式变换数据就会较快。示例代码为：
int length = 50000; fftw_complex* din = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2); fftw_complex* dout = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2); fftw_plan p = fftw_plan_dft_1d(length, din, din, FFTW_FORWARD, FFTW_MEASURE); fftw_execute(p); // 输入数据din赋值 // ... fftw_execute(p); // 读取变换结果 // ... fftw_destroy_plan(p); fftw_free(din); fftw_free(dout); |
实验发现第一个fftw_execute耗费了数秒，而第二个fftw_execute则瞬间完成，说明上述猜想可能是对的。
创建完方案（fftw_plan）后，就可以用fftw_execute对指定的 数据in/out做任意次变换。如果想变换一个相同大小（N相等）但数据不同的另外一个数组in，可以创建一个新方案，FFTW会自动重用上次方案的信 息。这一点其实是非常好的，比如你首先用FFTW_MEASURE模式创建了一个最优的变换方案，只要变换数据的大小不变，你可以用 fftw_plan_dft_1d创建新的方案以对新数据执行变换，同时新变换仍然是最优的。一个fftw_plan只能对固定的in/out进行变换，
但可以在变换后改变in的内容（大小不变）以用同一个方案执行新的变换。
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); fftw_plan fftw_plan_dft(int rank, const int *n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); |
多维复数据的DFT同一维复数据的DFT用法类似，数组in/out为行优先方式 存储。fftw_plan_dft是一个通用的复DFT函数，可以执行一维、二维或多维复DFT。比如对于图像（2维数据），其变换为 fftw_plan_dft_2d（height,width,85）,因为原始图像数据为height×width的矩阵，即第一维(n0)为行数 height。
函数接口
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); |
r2c版本：实输入数据，复Hermitian输出，正变换。
c2r版本：复Hermitian输入数据，实输出数据，逆变换。
n：逻辑长度，不必为物理长度。由于实数据的DFT具有 Hermitian对称性，所以只需要计算n/2+1（向下取整）个输出就可以了。比如对于r2c，输入in有n个数据，输出out有floor(n /2)＋1个数据。对于原位运算，in和out为同一数组（out须强制类型转换），所以其必须足够大以容纳所有数据，长度为2*(n/2+1)，in数 组的前n个数据为输入数据，后面的数据用来保存输出。
c2r逆变换在任何情况下（不管是否为原位运算）都破坏输入数组in，如果有必要可以通过在flags中加入FFTW_PRESERVE_INPUT来阻止，但这会损失一些性能，而其这个标志位目前在多维实DFT中已不被支持。
比如对于n=4，in=[1 2 3 4]，out=[10 -2+2i -2 -2-2i]，out具有共轭对称性，out的实际内存为10 0 -2 2 -2 0，共3个复数数据。对于n=5，in=[1 2 3 4 5]，out=[15 -2.5+3.44i -2.5+0.81i -2.5-0.81i -2.5-3.44i] ，out的实际内存为15 0 -2.5 3.44 -2.5 0.81，共3个复数数据。
实数据DFT中，首个变换数据为直流分量，为实数；如果n为偶数，由 Nyquist采样定理，第n/2个变换数据也为实数；所以可以把这两个数据组合在一起，即将第n/2个变换数据作为第0个变换数据的虚部，这样一来输入 数组就和输出数组相等（n=n/2*2）。一些FFT的实现就是这么做的，但FFTW没有这么做，因为它并不能很好地推广到多维DFT中，而且存储空间的 节省也是非常小以至于可以忽略不计。
一个一维c2r和r2c DFT的替代接口可以在r2r接口中找到，它利用了半复数输出类型（即实部和虚部分开放在不通的区域里），使输出数组具有和输入数组同样的长度和类型。该接口在多维变换中用处不大，但有时可能会有一些性能的提升。
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); |
这是r2c接口（正变换），c2r接口（逆变换）只是简单的将输入/输出类型交换一下。用法大致同1d情况一样，需要特别注意的是复数据的存放方式。对于n0×n1×n1×…×nd-1的实输入数组（行优先），经过r2c正变换后，输出为一个n0×n1×n1×…×(nd-1/2+1)的复数（fftw_complex）数组（行优先），其中除法向下取整。由于复数数据的总长度要大于实数据，所以如果需要进行原位运算则输入实数组必须扩展以能够容纳所有复数据，即实数组的最后一维必须包含2(floor(nd-1/2)+1)个double元素。比如对于一个2维实正DFT，输入数组为n0×n1大小，输出复数组大小为n0×floor(n1/2+1)（由对称性），其长度大于实输入数组。所以对于原位运算，输入数组要扩展到n0×2floor(n1/2+1)。如果n1为偶数，扩展为n0×（n1＋2）；如果n1为奇数，扩展为n0×（n1＋1）；这些扩展的内存不需要赋初值，它们只用来存放输出数据。
比如对于3×3的实正DFT，in=[0 2 4;6 1 3;5 7 4]，out=[32 0.5+0.86i 0.5-0.86i;-7+5.2i -1-1.73i -8.5-6.06i;-7-5.2i -8.5+6.06i -1+1.73i]；out的实际内存为32,0,0.5,0.86,-7,5.2,-1,-1.73,-7,-5.19,-8.5,6.06；即为 3×2的复数组，换算成double类型为3×4，所以如果要进行原位运算，in数组大小必须为3×4，即最后一维（每行）扩展一个double元素。另
外，如果用这个out数组进行3×3的c2r逆变换，将得到实数据[0 18 36;54 9 27;45 63 36]，即原始数据的9（n0×n1）倍，这是因为FFTW的逆变换没有归一化。