搜索
FFTW使用手册使用手册使用手册使用手册Ver1.0熊金水xjs.xjtu@gmail.com2011-5-11一、FFTW是什么FFTW—FastestFourierTransformintheWest，一个颇为古怪的名字，是由MIT的MatteoFrigo博士和StevenG.Johnson博士开发的一个完全免费的软件包。FFTW最初的release版本于1997年发布，最新的release版本3.2.2于2009年7月发布。它是一个C语言开发的库，支持任意大小的、任意维数的数据的离散傅里叶变换(DFT)，并且还支持离散余弦变换(DCT)、离散正弦变换(DST)和离散哈特莱变换(DHT)。根据两位博士的测试，FFTW在计算速度上远远优于目前其他计算DFT的免费的库，甚至可以和收费的DFT库相媲美。但是，和收费的DFT库相比，FFTW能够轻松的在不同的平台上移植。我们熟悉的MATLAB，就是调用了FFTW来实现DFT/IDFT变换的。另外，值得一提的是，FFTW还获得了1999年的J.H.WilkinsonPrizeforNumericalSoftware奖。J.H.WilkinsonPrizeforNumericalSoftware每四年颁发一次，用于奖励那些“最好的解决了高质量的数值计算问题的软件设计”。二、FFTW有哪些优越的性能高速高速高速高速——远优于目前其他的免费DFT库。支持任意维度任意维度任意维度任意维度的变换。支持任意大小任意大小任意大小任意大小的变换——FFTW对N=2^a*3^b*5^c*7^d*11^e*13^f的变换处理的最好，其中e+f=0或1，其他指数可以为任意值。支持快速的输入为输入为输入为输入为实数实数实数实数的DFT变换。支持DCT(I-IV)和和和和DST(I-IV)。支持多线程多线程多线程多线程。支持并行处理并行处理并行处理并行处理。可移植性可移植性可移植性可移植性——任意包含C编译器的平台都可以使用FFTW。同时包含C和和和和Fortran接口接口接口接口。完全免费免费免费免费——如果您在您的系统里使用了FFTW，请您尊重两位博士的辛勤劳作，在您的参考文献(Reference)中添加下面这篇文章。MatteoFrigoandStevenG.Johnson,TheDesignandImplementationofFFTW3,ProceedingsoftheIEEE93(2),216–231(2005).三、如何在windows上安装FFTW第一步，在上下载32位/64位版本的windows安装程序fftw-3.2.2.pl1-dll32.zip或fftw-3.2.2-dll64.zip，并解压。第二步，打开windows命令行窗口，把命令行的目录转到第一步中解压的目录下。运行一下三条指令：lib/machine:ix86/def:libfftw3-3.deflib/machine:ix86/def:libfftw3f-3.deflib/machine:ix86/def:libfftw3l-3.def这样会在您的当前目录下生成三个lib文件和三个相应的dll文件，其中文件名包含fftw3-3的为double精度，fftw3f-3的为float精度，fftw3l-3的为longdouble精度——在32的机器中，longdouble精度与double没有什么区别。第三步，对于lib文件，在您的VS开发工具里添加lib文件的目录，并在VS工程里添加您需要的精度所对应的lib文件。对于相应的dll文件，您可以把第二步中的目录设为环境变量的目录，然后注销重启就可以使用了，也可以拷贝到system32目录下，或者拷贝到您的工程目录下。对于头文件fftw3.h，您可以在您的VS工具下添加fftw3.h的目录，也可以直接把fftw3.h拷贝到您的工程目录下，我推荐使用后者，因为您在开发过程中，可能需要修改头文件，如果您再别的工程里再次使用该库时，就难以找到原始的头文件了。四、第一次使用FFTW第一次使用FFTW编程，可以使用下面的结构：总的来说就是先输入，然后构造策略plan，最后执行plan就可以了，还是很简单的。下面，我们一步一步来解析这个程序，并剖析每一步中需要注意些什么。1、内存空间的申请与释放以in的分配为例。在这里，空间分配有三种可以选择的方式：第一，直接in[N]；第二，使用ANSIC或者C++语言中的malloc，new等动态分配；第三，使用FFTW提供的fftw_malloc函数动态分配。第一种方法的问题在于，这种方法申请的空间由编译器在栈内存里分配，众所周知，栈内存是非常有限的，在windows上为2M大小，所以对于大的数据容易导致StackOverflow的致命性错误。当然，你可以在编译器设置里边它调大一些，但是换了个环境，你是不是就愣了~第二种方法虽然解决了第一种方法存在的问题，但是，由于不同的编译器内存分配策略的不同，可能使分配的数组并没有内存对齐，而内存没有对齐，对FFTW性能上的影响将是#includefftw3.h...{fftw_complex*in,*out;fftw_planp;...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提供的内存分配接口，它在堆内存上申请空间，并且能够保证内存对齐，同时可以在不同的平台之间顺利的移植，何乐而不为~空间的分配问题解决了，空间的释放问题当然也就水到渠成了，如果我们使用了fftw_malloc()来分配空间，那么使用fftw_free()释放相应的空间就行了。最后，在内存方面值得注意的问题是，fftw_plan是一种声明了一个变量就需要使用fftw_destroy_plan()来销毁的类型。2、实数FFT中内存的复用对于实数输入的序列，其DFT输出是一个有共轭对称性质的复数序列，在FFTW中，两位博士利用了这一点来减少内存空间的使用。对于一维的情况，假设输入序列为x(n),n=0~N-1，对应的DFT为X[k],k=0~N-1。我们知道，X(n)有共轭对称性，即：X[k]X[N-k]在FFTW中的X[k]的存储策略如下图所示。对于N为偶数的情况——以N=4为例：01230r0i1r1i2r2i变换以后，内存空间多了两个单元。对于N为奇数的情况——以N=5为例：012340R0I1R1I2R2I变换以后，内存空间多了一个单元。更一般的，我们可以总结一下，对于输入为N的序列，我们需要2*(floor(N/2)+1)个存储单元。N为奇数时，增加一个单元；N为偶数时，增加两个存储单元。这里有两个问题，当in和out不是一块内存时，即使我们为out申请了N个fftw_complex类型的内存单元，但是后面会有一部分的内存空间不被操作，如下图所示，FFT的output的后两个存储单元都没有操作（全是0是因为我在程序中手动赋值为0）。对于二维的情况，对称性显得稍微复杂一些。假设输入为x[m][n],m=0~M-1,n=0~N-1,对应的DFT序列为X[u][v],u=0~M-1,v=0~N-1。从二维傅里叶变换的公式中，我们可以很轻松的得到X有下面的共轭对称性质：X[u][v]X[M-u][N-v]u=0时，X[0][v]X[M][N-v]=X[0][N-v];v=0时，X[u][0]X[M-u][N]=X[M-u][0];u!=0且v!=0时，X[u][v]X[M-u][N-v]，二者的对称中心在(M/2,N/2)处。存储策略和一维的情况有些相似。012345678910111213141500112244556688991010121213131414以上是列数为偶数的情况，以下是列数为奇数的情况。01234567891011121314151617181920212223240011225566771010111112121515161617172020212122223、FFTW中提供了哪些策略，各有什么不同FFTW中，针对不同的情况，提供了下面一些生成策略(plan)的函数接口。fftw_planfftw_plan_dft_1d(intn,fftw_complex*in,fftw_complex*out,intsign,unsignedflags);//一维复数据的DFT、IDFTfftw_planfftw_plan_dft_2d(intn0,intn1,fftw_complex*in,fftw_complex*out,intsign,unsignedflags);//二维复数据DFT、IDFTfftw_planfftw_plan_dft_3d(intn0,intn1,intn2,fftw_complex*in,fftw_complex*out,intsign,unsignedflags);//三维复数据DFT、IDFTfftw_planfftw_plan_dft(intrank,constint*n,fftw_complex*in,fftw_complex*out,intsign,unsignedflags);//rank维复数据DFT、IDFTfftw_planfftw_plan_dft_r2c_1d(intn,double*in,fftw_complex*out,unsignedflags);//一维实数DFTfftw_planfftw_plan_dft_c2r_1d(intn,fftw_complex*in,double*out,unsignedflags);//一维实数IDFTfftw_planfftw_plan_dft_r2c_2d(intn0,intn1,double*in,fftw_complex*out,unsignedflags);//二维实数DFTfftw_planfftw_plan_dft_r2c_3d(intn0,intn1,intn2,double*in,fftw_complex*out,unsignedflags);//三维实数DFTfftw_planfftw_plan_dft_r2c(intrank,constint*n,double*in,fftw_complex*out,unsignedflags);//rank维实数DFT在复数据的DFT中，sign表示了是做DFT变换(FFTW_FORWARD)还是IDFT变换(FFTW_BACKWARD)。另外，在每一个函数中，我们都能找到一个变量叫flags，这个变量就是我们要重点讨论的策略生成方案。flags参数一般情况下常用的为FFTW_MEASURE或FFTW_ESTIMATE。FFTW_MEASURE表示FFTW会先计算一些FFT并测量所用的时间，以便为大小为n的变换寻找最优的计算方法。依据机器配置和变换的大小（n），这个过程耗费约数秒（时钟clock精度）。FFTW_ESTIMATE则相反，它直接构造一个合理的但可能是次最优的方案。总体来说，如果你的程序需要进行大量相同大小的FFT，并且初始化时间不重要，可以使用FFTW_MEASURE，否则应使用FFTW_ESTIMATE。FFTW_MEASURE模式下in和out数组中的值会被覆