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，输出格式如下：

r₀, r₁, r₂, ..., r_n/2, i_(n+1)/2-1, ..., i₂, i₁

上述(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）。REDFTab和RODFTab中的a,b为数据移位标志（1表示移位），这些构成了DCT和DST的I－IV类，比较常用的为DCT-II，FFTW支持所有这些类型的变换。

FFTW_REDFT00 (DCT-I): even around j=0 and even around j=n-1.
FFTW_REDFT10 (DCT-II, the DCT): even around j=-0.5 and even around j=n-0.5.
FFTW_REDFT01 (DCT-III, the IDCT): even around j=0 and odd around j=n.
FFTW_REDFT11 (DCT-IV): even around j=-0.5 and odd around j=n-0.5.
FFTW_RODFT00 (DST-I): odd around j=-1 and odd around j=n.
FFTW_RODFT10 (DST-II): odd around j=-0.5 and odd around j=n-0.5.
FFTW_RODFT01 (DST-III): odd around j=-1 and even around j=n-1.
FFTW_RODFT11 (DST-IV): odd around j=-0.5 and even around j=n-0.5.

对称性只是逻辑意义上的，对物理输入数据没有任何限制。比如对于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)；对于 RODFT00，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,
                             int sign, unsigned flags);

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及逆变换：

FFTW_R2HC computes a real-input DFT with output in halfcomplex format, i.e. real and imaginary parts
for a transform of size n stored as:r₀, r₁, r₂, ..., r_n/2, i_(n+1)/2-1, ..., i₂, i₁
(Logical N=n, inverse is FFTW_HC2R.)
FFTW_HC2R computes the reverse of FFTW_R2HC, above. (Logical N=n, inverse is FFTW_R2HC.)
FFTW_DHT computes a discrete Hartley transform. (Logical N=n, inverse is FFTW_DHT.)
FFTW_REDFT00 computes an REDFT00 transform, i.e. a DCT-I. (Logical N=2*(n-1), inverse is FFTW_REDFT00.)
FFTW_REDFT10 computes an REDFT10 transform, i.e. a DCT-II (sometimes called the DCT).
(Logical N=2*n, inverse is FFTW_REDFT01.)
FFTW_REDFT01 computes an REDFT01 transform, i.e. a DCT-III (sometimes called the IDCT, being the inverse of DCT-II).
(Logical N=2*n, inverse is FFTW_REDFT=10.)
FFTW_REDFT11 computes an REDFT11 transform, i.e. a DCT-IV. (Logical N=2*n, inverse is FFTW_REDFT11.)
FFTW_RODFT00 computes an RODFT00 transform, i.e. a DST-I. (Logical N=2*(n+1), inverse is FFTW_RODFT00.)
FFTW_RODFT10 computes an RODFT10 transform, i.e. a DST-II. (Logical N=2*n, inverse is FFTW_RODFT01.)
FFTW_RODFT01 computes an RODFT01 transform, i.e. a DST-III. (Logical N=2*n, inverse is FFTW_RODFT=10.)
FFTW_RODFT11 computes an RODFT11 transform, i.e. a DST-IV. (Logical N=2*n, inverse is FFTW_RODFT11.)

八、其它

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等。
所有以大写"FFTW_"开头的名字不变
将函数参数中的double替换为float或long double

最后，虽然long double是C99的标准，但你的编译器可能根本不支持该类型，或它并不能提供比double更高的精度。

2. 用同一个fftw_plan执行多个数据的变换

前面说过可以利用同一个fftw_plan通过对输入数据赋不同值来实现不同的变换，实际上还可以利用同一个fftw_plan实现对不同输入输出数据的变换，也就是说可以有多个输入输出数据数组，各自进行变换，互不影响。当然这要满足一定的条件：

输入/输出数据大小相等。
变换类型、是否原位运算不变。
对split数组（指实虚部分开），实部和虚部的分割方式与方案创建时相同。
数组的对齐方式相同。如果都是用fftw_malloc分配的则此项条件满足，除非使用 FFTW_UNALIGNED标志。

如果想对新数组，比如大小相等的一批数组执行变换，可以使用以下接口：

     void fftw_execute_dft(
          const fftw_plan p,
          fftw_complex *in, fftw_complex *out);
     
     void fftw_execute_split_dft(
          const fftw_plan p,
          double *ri, double *ii, double *ro, double *io);
     
     void fftw_execute_dft_r2c(
          const fftw_plan p,
          double *in, fftw_complex *out);
     
     void fftw_execute_split_dft_r2c(
          const fftw_plan p,
          double *in, double *ro, double *io);
     
     void fftw_execute_dft_c2r(
          const fftw_plan p,
          fftw_complex *in, double *out);
     
     void fftw_execute_split_dft_c2r(
          const fftw_plan p,
          double *ri, double *ii, double *out);
     
     void fftw_execute_r2r(
          const fftw_plan p,
          double *in, double *out);

这些函数的执行不会修改原始plan，并且可以和fftw_execute混合使用。

3. 多线程FFTW

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

4. FFTW变换公式

一维复DFT正变换
一维复DFT逆变换
一维实DFT正变换（Y_k = Y_n-k_^*）
一维实DFT逆变换
REDFT00 (DCT-I)
REDFT10 (DCT-II)
REDFT01 (DCT-III)
REDFT11 (DCT-IV)
RODFT00 (DST-I)
RODFT10 (DST-II)
RODFT01 (DST-III)
RODFT11 (DST-IV)
1d Discrete Hartley Transforms (DHTs)

六、 更多实数据的DFT

七、 函数参考

八、 其它

六、更多实数据的DFT

七、函数参考

八、其它