FFTW3学习笔记2：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），最后释放资源，还是比较简单的。

二、 一维复数据的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二进制兼容，即顺序存储，实部在前，虚部在后（见报告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=2^a 3^b 5^c 7^d 11^e 13^f的变换处理得最好，其中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的内容（大小不变）以用同一个方案执行新的变换。

三、 多维复数据的DFT

fftw_plan fftw_plan_dft_2d(int n0, int n1, fftw_complex *in, fftw_complex *out, int sign, unsigned f