FFTW3学习笔记2：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），最后释放资源，还是⽐较简单的。
⼆、⼀维复数据的DFT
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⼆进制兼容，即顺序存储，实部在前，虚部在后（见报告），该解决⽅案在所有主流标准库实现中都能正确⼯作。所以实际上可以⽤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模式⾃动选最优⽅案，速度较慢；然后使⽤该模式变换数据就会较快。⽰例代码为：
sizeof是什么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的内容（⼤⼩不变）以⽤同⼀个⽅案执⾏新的变换。
三、多维复数据的DFT
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。
四、⼀维实数据的DFT
函数接⼝
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接⼝中到，它利⽤了半复数输出类型（即实部和虚部分开放在不通的区域⾥），使输出数组具
有和输⼊数组同样的长度和类型。该接⼝在多维变换中⽤处不⼤，但有时可能会有⼀些性能的提升。
五、多维实数据的DFT
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×…×n d-1的实输⼊数组（⾏优先），经过r2c正变换后，输出为⼀个n0×n1×n1×…×(n d-1/2+1)的复数（fftw_complex）数组（⾏优先），其中除
法向下取整。由于复数数据的总长度要⼤于实数据，所以如果需要进⾏原位运算则输⼊实数组必须扩展以能够容纳所有复数据，即实数组的最后⼀维必须包含2(floor(n d-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, ..., r n/2, i(n+1)/2-1, ..., i2, i1
上述(n+1)/2向下取整。r k是第k个输出的实部，i k是第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）。REDFT ab和RODFT ab中的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)；对于 R
ODFT00，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,
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等。
发表评论