FFTW算法计算结果为空/为0

FFTW算法计算结果为空

可以快速跳到个人总结看结论。

1.问题的出现

最近做音频处理有关的需要用到FFT算法，考虑到自己直接研究算法周期长效率低，所以考虑使用“最快FFT”之称的FFTW3.0 。

由于官方的ftp下载我这么都连不上，花了一天找库文件和学习使用，不过官方文档是全英文就没有仔细啃（危）。写了一个demo测试算法，核心代码如下：

fftw_plan fp,ifp;				//正变换和反变换的方案
double *data = NULL;			//数据输入和逆变换输出（偷懒了
fftw_complex *result = NULL;	//正变换结果输出
data = (double*)fftw_malloc(sizeof(double)*SIZE);
result = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*(SIZE/2+1));
/*****给in赋值，代码省略*****/
fp = fftw_plan_dft_r2c_1d(SIZE1, data, result, FFTW_MEASURE);
fftw_execute(fp);			//执行正变换
ifp = fftw_plan_dft_c2r_1d(SIZE1, result, data, FFTW_MEASURE);
fftw_execute(ifp);			//执行逆变换

fftw_destroy_plan(fp);		//清除计划和空间
fftw_destroy_plan(ifp);
fftw_free(data);
fftw_free(result);

运行后结果发现如果输入的序列长度大于17，那么正变换结果全是0，更不用说逆变换也都是0了。

2.查找原因

网上搜索许久无果，于是找了一份代码跑一遍试试，无论序列长短都能正常得到结果，于是对比我自己的代码和copy的代码，发现fftw_plan_dft_r2c_1d的flag参数不同，我使用的是plan速度慢、后续计算更快的FFTW_MEASURE，而copy的代码中使用了plan更快的FFTW_ESTIMATE，把我自己的参数换掉，结果正常~

但这并不是我想要的，之后需要对音频实时处理，后续处理速度才是第一位的，所以必不可能使用FFTW_ESTIMATE（不是真的因为慢而是任性）。

于是到官方手册中查找关键词FFTW_MEASURE，查到了下面一段文字：

The flags argument is usually either FFTW_MEASURE or FFTW_ESTIMATE. FFTW_MEASURE instructs FFTW to run and measure the execution time of several FFTs in order to find the best way to compute the transform of size n. This process takes some time (usually a few seconds), depending on your machine and on the size of the transform. FFTW_ESTIMATE, on the contrary, does not run any computation and just builds a reasonable plan that is probably sub-optimal. In short, if your program performs many transforms of the same size and initialization time is not important, use FFTW_MEASURE; otherwise use the estimate. You must create the plan before initializing the input, because FFTW_MEASURE overwrites the in/out arrays. (Technically, FFTW_ESTIMATE does not touch your arrays, but you should always create plans first just to be sure.)

大意是MEASURE比ESTIMATE的plan速度慢计算速度快，但是MEASURE会破坏输入，需要先plan后初始化。我一开始没有理解这段话， 于是到google上继续以这两个flag为关键词搜索，找到了一篇和我问题一样的问答。查看这个网页

The comment of @Paul_R. is sufficient to solve the problem. The input array can be modified as fftw_plan_dft_r2c_2d() is called. Hence, the input array must be initialized after the creation of the fftw plan.

回答者加粗了after，我才恍然大悟，于是调换顺序，把fp和ifp的plan移到input初始化之前，重新计算任何长度都可以完美得到结果。

fp = fftw_plan_dft_r2c_1d(SIZE1, data, result, FFTW_MEASURE);
ifp = fftw_plan_dft_c2r_1d(SIZE1, result, data, FFTW_MEASURE);
/*****给in赋值，代码省略*****/
fftw_execute(fp);			//执行正变换
fftw_execute(ifp);			//执行逆变换

3.个人总结

当plan参数使用MEASURE时，函数会计算各种方法比较时间，这样才知道那个方案最快，这也是为什么plan时间比较长的原因。计算过程中input参与计算，值会发生改变，所以fftw_execute(fp)的时候input已经变成了全0序列，后续计算结果也自然是0。当逆变换时，同样的原理，复数序列也会参与计算而失去本身。

官方文档中提醒，

but you shouldalways create plans first just to be sure

建立plan永远放在赋值之前