FFTW使用小结

简介
================
FFTW—Fastest Fourier Transform in the West,是由 MIT 的 Matteo Frigo 博士和 Steven G. Johnson 博士开发的一个完全免费的软件包。FFTW 最初的 release 版本于 1997 年发布,最新的 release 版本 fftw-3.3.4。git路径:
https://github.com/FFTW/fftw3.git
它是一个 C 语言开发的库,支持任意大小的、任意维数的数据的离散傅里叶变换(DFT),并且还支持离散余弦变换(DCT)、离散正弦变换(DST)和离散哈特莱变换(DHT)

数据类型
================
FFTW 有三个版本的数据类型:double、float 和 long double,使用方法如下:
1.链接对应的库(比如 libfftw3-3、libfftw3f-3、或 ibfftw3l-3)
2.包含同样的头文件 fftw3.h
将所有以小写"fftw_"开头的名字替换为"fftwf_"(float 版本)或"fftwl_"(long double 版本)。比如将 fftw_complex 替换为 fftwf_complex,将 fftw_execute
替换为 fftwf_execute 等。
3.所有以大写"FFTW_"开头的名字不变
4.将函数参数中的 double 替换为 float 或 long double
5.最后,虽然 long double 是 C99 的标准,但你的编译器可能根本不支持该类型,或它并不能提供比 double 更高的精度。
6.fftw_malloc 考虑了数据对齐,以便使用 SIMD 指令加速,所以最好不要用 C 函数malloc 替代,而且不要将 fftw_malloc、fftw_free 和 malloc、free、 delete 等混用。
尽量使用 fftw_malloc 分配空间,而不是使用的静态数组,因为静态数组是在栈上分配的,而栈空间较小;还因为这种方式没有考虑数据对齐,不便应用SIMD 指令。


编译构造
=================
1.默认FFTW编译生成double类型，加入参数“--enable-single”或“--enable-float”编译单精度(float)，加入参数“--enable-long-double”支持长双进度类型
2.另外,在ARM平台上可增加“--enable-neon”,使能ARM NEON流媒体加速核心

多线程
=================
FFTW 可以多线程执行,但是多线程存在线程同步问题,这可能会降低性能。所以除非问题规模非常大,一般并不能从多线程中获益。

plan复用
=================
用同一个 fftw_plan 执行多个数据的变换
同一个 fftw_plan 通过对输入数据赋不同值来实现不同的变换,实际上还可以利用同一个 fftw_plan 实现对不同输入输出数据的变换,也就是说可以有多个
输入输出数据数组,各自进行变换,互不影响。当然这要满足一定的条件:
1.输入/输出数据大小相等。
2.变换类型、是否原位运算不变。
3.对 split 数组(指实虚部分开),实部和虚部的分割方式与方案创建时相同。
4.数组的对齐方式相同。如果都是用 fftw_malloc 分配的则此项条件满足,除非使用
FFTW_UNALIGNED 标志

总结：是用现有PLAN必须满足所有参数条件和申请PLAN时一致。

一维数据
=================
目前使用最广泛，项目也只使用1D，并且输入数据为实数，因此变化有两种方法
1.将输入部分扩展成虚数，即所有虚数为0，然后使用一维复数变化接口

fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
参数说明:
n -- 复数数据点数
in/out -- 输入数据和输出数据，可以相同(原位相同)
sign -- FFTW_FORWARD(-1)正变换 FFTW_BACKWORD(+1) 逆变换
flags -- 主要有两个参数 FFTW_MEASURE/FFTW_ESTIMATE
    FFTW_MEASURE 先进行预处理，对大数据(数据长度不变)连续FFT变化效果明显
    FFTW_ESTIMATE 直接构造一个次最优的方案，非连续 实时选择这种方案

2.使用1D实数DFT变化，由于实数DFT变化具有Hmitian对称性，所以只需计算n/2 + 1个输出数据即可
需要注意：
a.如果采用原位运算则out空间必须>= 2*（n/2 + 1)
b.1D实数DFT逆变换任何情况都会破坏IN数据
c.当n为偶数，由Nyquust采样定理，第n/2个变化结果也是实数，所以可以把第0个和n/2个数据合并(n/2实部做0虚部)，但实际上FFTW没有这样实现，因为在多维数据时不兼容

几种方案测试

=============

使用C2C MEASURE
===============================
FFT-DFT :309us
FFT-DFT :345us
FFT-DFT :320us
FFT-DFT :349us
FFT-DFT :319us
FFT-DFT :320us
FFT-DFT :321us
FFT-DFT :324us
FFT-DFT :320us
FFT-DFT :319us
FFT-DFT :757919us
幅值:9.465116    相位:-169.040115

real    0m2.064s
user    0m2.044s
sys    0m0.020s
===============================

===============================
使用R2C
FFT-DFT :208us
FFT-DFT :195us
FFT-DFT :198us
FFT-DFT :214us
FFT-DFT :212us
FFT-DFT :212us
FFT-DFT :221us
FFT-DFT :215us
FFT-DFT :214us
FFT-DFT :216us
FFT-DFT :564592us
幅值:9.465116    相位:-169.040115

real    0m0.584s
user    0m0.552s
sys    0m0.032
================================
================================
使用C2C 次最优方案
FFT-DFT :213us
FFT-DFT :204us
FFT-DFT :195us
FFT-DFT :195us
FFT-DFT :195us
FFT-DFT :196us
FFT-DFT :197us
FFT-DFT :196us
FFT-DFT :198us
FFT-DFT :197us
FFT-DFT :554788us
幅值:9.465116    相位:-169.040115

real    0m0.571s
user    0m0.544s
sys    0m0.028s
================================

采样点数为24000，进行1维复数/实数变化，参数选择有FFTW_MEASURE/FFTW_ESTIMATE

1.经测试，发现使用MEASURE方案时，在第一构造plan时花费约1.5秒，而一次DFT变化在300us-,构造的时间能做普通(FFTW_ESTIMATE)5000次左右，故在一般场合FFTW_ESTIMATE完全满足需求。

2.实数FFTW_ESTIMATE,时间上并没有优于复数DFT变化,只是内存开辟只有复数DFT的一半。实数DFT在构造plan时间上比复数慢，但变换比复数快，所以整体持平

3.最终选择复数和实数 ESTIMATE方案

FFTW_WISDOM
============
关于WISDOM,翻译自FFTW官网
FFTW -wisdom是一种实用工具来生成wisdom文件，其中包含有关如何以最佳方式计算（傅立叶）变换不同大小的保存的信息。
wisdom 的大体思路就是把生成好的策略相关的配置信息存储在磁盘里,然后在下次重新运行程序的时候,把策略相关的配置信息重新载入到内存中,这样在重新生成 plan 的时候就可以节约大量的时间。

FFTW 提供了多种方式来生成wisdom文件,使用时不必关心其中具体格式是怎么样，可以导入/出到文件/字符串等

322 FFTW_EXTERN int X(export_wisdom_to_filename)(const char *filename);        \
323 FFTW_EXTERN void X(export_wisdom_to_file)(FILE *output_file);              \
324 FFTW_EXTERN char *X(export_wisdom_to_string)(void);                        \
325 FFTW_EXTERN void X(export_wisdom)(X(write_char_func) write_char,           \
326                                   void *data);                             \
327 FFTW_EXTERN int X(import_system_wisdom)(void);                             \
328 FFTW_EXTERN int X(import_wisdom_from_filename)(const char *filename);      \
329 FFTW_EXTERN int X(import_wisdom_from_file)(FILE *input_file);              \
330 FFTW_EXTERN int X(import_wisdom_from_string)(const char *input_string);    \
331 FFTW_EXTERN int X(import_wisdom)(X(read_char_func) read_char, void *data); 

wisdom 存储起来的不是 plan 本身,而是和 plan 相关的配置信息,例如内存、寄存器等。故在导入wisdom后还是需要执行plan初始化。
example:
    fftw_export_wisdom_to_filename("wisdom.conf");
    p_str = fftw_export_wisdom_to_string();
    
    fftw_export_wisdom_from_filename("wisdom.conf");
    fftw_export_wisdom_from_string(p_str);