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的逆变换没有归一化。
六、 更多实数据的DFT
通过一个统一的r2r（real-to-real，实数－实数）接口，FFTW支持其它的一些变换类型，这些变换的输入和输出数组大小相同。这些r2r变换可以分为3个类型：DFT的实 数据输入，complex-Hermitian（指复Hermitian对称）以半复数格式的输出；DCT/DST（离散正余弦变换）；DHT（离散 Hartley变换）。接口如下：
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);
|
这里n为数组的物理尺寸。对于多维变换，数组按行优先方式存储（与C++标准相同，与Fortran不同）。由于DFT是可分离变换，所以2维/3维/多维的变换是在每个维度上分别进行变换得到的，每个维度都可指定一个kind参数，指定该维的变换类型。
1. 半复数格式DFT（HalfComplex-format）
对于大小为n的1维DFT，输出格式如下：
r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1
上述(n+1)/2向下取整。rk是第k个输出的实部，ik是 第k个输出的虚部。对于一个半复数格式的数组hc[n]，第k个元素的实部为hc[k]，虚部为[n-k]；k==0或n/2（n为偶数）情况除外，这两 种情况下虚部为0，不存储。所以对于r2hc(real-half complex，正变换)变换，输入输出数组大小都为n，hc2r（half complex- real，逆变换）变换也是如此。除了数据格式的差异，r2hc和hc2r变换的结果与前述r2c和c2r的变换结果是相同的。
对于多维比如2维变换，由可分离性，可以先对行 变换，再对列变换，FFTW_R2HC方式行变换的结果是半复数行，如果采用FFTW_R2HC方式进行列变换，需要进行一些数据处理，r2r变换不会自 动处理，需要手动进行，所以对于多维实数据变换，推荐使用普通的r2c/c2r接口。
2. DCT/DST
DCT可以认为是实偶对称数据DFT（REDFT,Real-Even DFT）, DST可以认为是实奇对称数据DFT（RODFT,Real-Odd DFT）。REDFTab和RODFTab中的a,b为数据移位标志（1表示移位），这些构成了DCT和DST的I－IV类，比较常用的为DCT-II，FFTW支持所有这些类型的变换。
|
对称性只是逻辑意义上的，对物理输入数据没有任何限制。比如对于n＝5的REDFT00 (DCT-I)，输入数据为abcde，它对应n=8的abcdedcb的常规DFT；对于n＝4的REDFT10 (DCT-II)，输入数据为abcd，它对应n=8的abcddcba的常规DFT。
所有这些变换都是可逆的。R*DFT00的逆变 换是R*DFT00，R*DFT10的逆变换是R*DFT01（即DCT和IDCT），R*DFT11的逆变换是R*DFT11。如同DFT一样，这些变 换的结果都没有归一化，所以正变换再逆变换后数据变为原始数据的N倍，N为逻辑DFT大小。比如对于REDFT00变换，N=2(n-1)；对于 RODFT00，N=2n。
注意n=1的REDFT00对应与N=0的DFT，所以它是未定义的（返回值为NULL的fftw_plan）。
R*DFT01和R*DFT10要比R*DFT11稍微快一些，尤其对于奇数长度数据；而R*DFT00则要慢一些，尤其对于奇数长度数据。
比如对于in=[1 2 3 4]，n=4，N=2n=8。Matlab下dct变换的结果为[5 -2.2304 0 -0.15851]；FFTW的结果为（FFTW_REDFT10）out=[20 -6.3086 0 -0.4483415]，为matlab结果的√8（√N）倍；用out进行逆dct变换（FFTW_REDFT01）的结果为[8 16 24 32]，为原始数据的8（N）倍。
再比如对于in=[0 2 4;6 1 3;5 7 4]的二维DCT变换，n=3，N=2n=6。Matlab下dct2的变换结果为out_matlab=[10.667 0 0.4714;-4.0825 -2.5 1.4434;0.4714 -2.5981 -3.1667]；FFTW的结果为（fftw_plan_r2r_2d(3, 3, in, out, FFTW_REDFT10, FFTW_REDFT10, FFTW_ESTIMATE）out_fftw=[128 0 4;-34.641 -15 8.66;4 -15.588 -19]，这与Matlab的结果同样是有差别的。将Matlab的结果变换到FFTW结果的步骤为：
1. 将out_matlab乘以√6×√6（即√N×√N）；
2. 再将上一步得到的out_matlab的第一行和第一列都乘以√2，因此第一个元素（即左上角的元素）要乘以2。
第一个是归一化的原因，第二个是FFTW为了将DCT变换与实偶对称FFT相对应的结果。这些对于DCT变换的应用都影响不大。（见此）
最后对out_fftw进行IDCT变换 （fftw_plan_r2r_2d(3, 3, in, out, FFTW_REDFT01, FFTW_REDFT01, FFTW_ESTIMATE），将得到[0 72 144;216 36 108;180 252 144]；它是原始数据in的36（N×N，N=6）倍。
3. 其它
fftw_malloc考虑了数据对齐，以便使 用SIMD指令加速，所以最好不要用C函数malloc替代，而且不要将fftw_malloc、fftw_free和malloc、free、 delete等混用。尽量使用fftw_malloc分配数组，而不是下面的固定数组，因为固定数组是在栈上分配的，而栈空间较小；还因为这种方式没有考 虑数据对齐，不便应用SIMD指令。
fftw_complex data[N0][N1][N2]; fftw_plan plan; ... plan = fftw_plan_dft_3d(N0, N1, N2, &data[0][0][0], &data[0][0][0], FFTW_FORWARD, FFTW_ESTIMATE); ... |
对于多维数组也尽量使用fftw_malloc(n0*n1*n285*sizeof(double))，不要使用C的malloc。
fftw_complex *a_good_array; a_good_array = (fftw_complex*) fftw_malloc(5*12*27* sizeof(fftw_complex)); fftw_complex ***a_bad_array; /* another way to make a 5x12x27 array */ a_bad_array = (fftw_complex ***) malloc(5 * sizeof(fftw_complex **)); |
1. 复数DFT
fftw_plan fftw_plan_dft_1d(int n,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
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);
|
2. 实数DFT
fftw_plan fftw_plan_dft_r2c_1d(int n,
double *in, fftw_complex *out,
unsigned flags);
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);
|
3. 实数－实数变换
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);
|
4. 实数－实数变换类型
对于大小为n的下列变换，对应的逻辑DFT大小为N，N用来进行归一化。FFTW的变换没有归一化，正变换后再逆变换为原数据的N倍（不是n倍），对于多维变换，为N的乘积(N0*N1*N285)。下表列出了变换类型及其对应的逻辑大小N及逆变换：
|
1. 数据类型
FFTW有三个版本的数据类型：double、float和long double，使用方法如下：
- 链接对应的库（比如libfftw3-3、libfftw3f-3、或ibfftw3l-3）
- 包含同样的头文件fftw3.h
- 将所有以小写"fftw_"开头的名字替换为"fftwf_"（float版本）或"fftwl_"（long double版本）。比如将fftw_complex替换为fftwf_complex，将fftw_execute替换为fftwf_execute等。
- 所有以大写"FFTW_"开头的名字不变
- 将函数参数中的double替换为float或long double
最后，虽然long double是C99的标准，但你的编译器可能根本不支持该类型，或它并不能提供比double更高的精度。
2. 用同一个fftw_plan执行多个数据的变换
前面说过可以利用同一个fftw_plan通过对输入数据赋不同值来实现不同的变换，实际上还可以利用同一个fftw_plan实现对不同输入输出数据的变换，也就是说可以有多个输入输出数据数组，各自进行变换，互不影响。当然这要满足一定的条件：
- 输入/输出数据大小相等。
- 变换类型、是否原位运算不变。
- 对split数组（指实虚部分开），实部和虚部的分割方式与方案创建时相同。
- 数组的对齐方式相同。如果都是用fftw_malloc分配的则此项条件满足，除非使用 FFTW_UNALIGNED标志。
如果想对新数组，比如大小相等的一批数组执行变换，可以使用以下接口：
void fftw_execute_dft(
const fftw_plan p,
fftw_complex *in, fftw_complex *out);
void fftw_execute_split_dft(
const fftw_plan p,
double *ri, double *ii, double *ro, double *io);
void fftw_execute_dft_r2c(
const fftw_plan p,
double *in, fftw_complex *out);
void fftw_execute_split_dft_r2c(
const fftw_plan p,
double *in, double *ro, double *io);
void fftw_execute_dft_c2r(
const fftw_plan p,
fftw_complex *in, double *out);
void fftw_execute_split_dft_c2r(
const fftw_plan p,
double *ri, double *ii, double *out);
void fftw_execute_r2r(
const fftw_plan p,
double *in, double *out);
|
这些函数的执行不会修改原始plan，并且可以和fftw_execute混合使用。
3. 多线程FFTW
FFTW可以多线程执行，但是多线程存在线程同步问题，这可能会降低性能。所以除非问题规模非常大，一般并不能从多线程中获益。
4. FFTW变换公式
一维复DFT正变换 |
|
一维复DFT逆变换 |
|
一维实DFT正变换（Yk = Yn-k* ） |
|
一维实DFT逆变换 |
|
REDFT00 (DCT-I) |
|
REDFT10 (DCT-II) |
|
REDFT01 (DCT-III) |
|
REDFT11 (DCT-IV) |
|
RODFT00 (DST-I) |
|
RODFT10 (DST-II) |
|
RODFT01 (DST-III) |
|
RODFT11 (DST-IV) |
|
1d Discrete Hartley Transforms (DHTs) |
|
转：http://blog.csdn.net/congwulong/article/details/7576012
