fftw中文参考手册.pdf

发布时间：2022-06-08 发布人：admin 分类：说明书资料大小：0.62M 资料格式：pdf 举报版权申诉

ghxandsky-3052538-4744302542970564198.pdf-第1页.png

第1页 / 共13页

ghxandsky-3052538-4744302542970564198.pdf-第2页.png

第2页 / 共13页

ghxandsky-3052538-4744302542970564198.pdf-第3页.png

第3页 / 共13页

ghxandsky-3052538-4744302542970564198.pdf-第4页.png

第4页 / 共13页

ghxandsky-3052538-4744302542970564198.pdf-第5页.png

第5页 / 共13页

ghxandsky-3052538-4744302542970564198.pdf-第6页.png

第6页 / 共13页

ghxandsky-3052538-4744302542970564198.pdf-第7页.png

第7页 / 共13页

ghxandsky-3052538-4744302542970564198.pdf-第8页.png

第8页 / 共13页

文本预览

据说 FFTW（Fastest Fourier Transform in the West）是世界上最快的 FFT。为了详细了解 FFTW 以及为编程方便，特将用户手册看了一下，并结合手册制作了以下 FFTW 中文参考。其中大部分是原文重点内容的翻译，并加入了一些注解。一、简介先看一下使用 FFTW 编程的方法： #include ... { 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 前加入#include ，这样一来 fftw_complex 就被定义为本机复数类型，而且与上述 typedef 二进制兼容（指内存排列），经测试不能用在 Windows 下。C++有一个复数模板类 complex，在头文件下定义。C++标准委员会最近同意该类的存储方式与 C99 二进制兼容，即顺序存储，实部在前，虚部在后（见报告 WG21/N1388），该解决方案在所有主流标准库实现中都能正确工作。所以实际上可以用 complex 来代替 fftw_complex，

比如有一个复数数组 complex *x，则可以将其类型转换后作为参数传递给 fftw： reinterpret_cast(x)。测试如下：开两个数组 fftw_complex x1[2] 和 complex x2[2]，然后赋相同值，在调试模式下可以看到它们的内存排列是相同的。complex类数据赋值的方式不是很直接，必须采用无名对象方式 x2[i] = complex (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 *类型使用（比如赋值或读取等），变换时再将其转换为 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=2a 3b 5c 7d 11e 13f 的变换处理得最好，其中 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 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); 多维复数据的 DFT 同一维复数据的 DFT 用法类似，数组 in/out 为行优先方式存储。 fftw_plan_dft 是一个通用的复 DFT 函数，可以执行一维、二维或多维复 DFT。比如对于

图像（2 维数据），其变换为 fftw_plan_dft_2d（height,width,85）,因为原始图像数据为 height×width 的矩阵，即第一维(n0)为行数 height。四、一维实数据的 DFT 函数接口 fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out, unsigned flags); fftw_plan fftw_plan_dft_c2r_1d(int n, fftw_complex *in, double *out, unsigned flags); r2c 版本：实输入数据，复 Hermitian 输出，正变换。 c2r 版本：复 Hermitian 输入数据，实输出数据，逆变换。 n：逻辑长度，不必为物理长度。由于实数据的 DFT 具有 Hermitian 对称性，所以只需要计算 n/2+1（向下取整）个输出就可以了。比如对于 r2c，输入 in 有 n 个数据，输出 out 有 floor(n /2)＋1 个数据。对于原位运算，in 和 out 为同一数组（out 须强制类型转换），所以其必须足够大以容纳所有数据，长度为 2*(n/2+1)，in 数组的前 n 个数据为输入数据，后面的数据用来保存输出。 c2r 逆变换在任何情况下（不管是否为原位运算）都破坏输入数组 in，如果有必要可以通过在 flags 中加入 FFTW_PRESERVE_INPUT 来阻止，但这会损失一些性能，而其这个标志位目前在多维实 DFT 中已不被支持。比如对于 n=4，in=[1 2 3 4]，out=[10 -2+2i -2 -2-2i]，out 具有共轭对称性， out 的实际内存为 10 0 -2 2 -2 0，共 3 个复数数据。对于 n=5，in=[1 2 3 4 5]，out=[15 -2.5+3.44i -2.5+0.81i -2.5-0.81i -2.5-3.44i] ，out 的实际内存为 15 0 -2.5 3.44 -2.5 0.81，共 3 个复数数据。实数据 DFT 中，首个变换数据为直流分量，为实数；如果 n 为偶数，由 Nyquist 采样定理，第 n/2 个变换数据也为实数；所以可以把这两个数据组合在一起，即将第 n/2 个变换数据作为第 0 个变换数据的虚部，这样一来输入数组就和输出数组相等（n=n/2*2）。一些 FFT 的实现就是这么做的，但 FFTW 没有这么做，因为它并不能很好地推广到多维 DFT 中，而且存储空间的节省也是非常小以至于可以忽略不计。一个一维 c2r 和 r2c DFT 的替代接口可以在 r2r 接口中找到，它利用了半复数输出类型（即实部和虚部分开放在不通的区域里），使输出数组具有和输入数组同样的长度和类型。该接口在多维变换中用处不大，但有时可能会有一些性能的提升。五、多维实数据的 DFT 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); 这是 r2c 接口（正变换），c2r 接口（逆变换）只是简单的将输入/输出类型交换一下。用法大致同 1d 情况一样，需要特别注意的是复数据的存放方式。对于 n0×n1×n1×…×nd-1 的实输入数组（行优先），经过 r2c 正变换后，输出为一个 n0×n1×n1×…×(nd-1/2+1) 的复数（fftw_complex）数组（行优先），其中除法向下取整。由于复数数据的总长度要大于实数据，所以如果需要进行原位运算则输入实数组必须扩展以能够容纳所有复数据，即实数组的最后一维必须包含 2(floor(nd-1/2)+1)个 double 元素。比如对于一个 2 维实正 DFT，输入数组为 n0×n1 大小，输出复数组大小为 n0×floor(n1/2+1)（由对称性），其长度大于实输入数组。所以对于原位运算，输入数组要扩展到 n0×2floor(n1/2+1)。如果 n1 为偶数，扩展为 n0×（n1＋2）；如果 n1 为奇数，扩展为 n0×（n1＋1）；这些扩展的内存不需要赋初值，它们只用来存放输出数据。比如对于 3×3 的实正 DFT，in=[0 2 4;6 1 3;5 7 4]，out=[32 0.5+0.86i 0.5-0.86i;-7+5.2i -1-1.73i -8.5-6.06i;-7-5.2i -8.5+6.06i -1+1.73i]；out 的实际内存为 32,0,0.5,0.86,-7,5.2,-1,-1.73,-7,-5.19,-8.5,6.06；即为 3×2 的复数组，换算成 double 类型为 3×4，所以如果要进行原位运算，in 数组大小必须为 3×4，即最后一维（每行）扩展一个 double 元素。另外，如果用这个 out 数组进行 3×3 的 c2r 逆变换，将得到实数据[0 18 36;54 9 27;45 63 36]，即原始数据的 9（n0×n1）倍，这是因为 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，输出格式如下： r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1 上述(n+1)/2 向下取整。rk 是第 k 个输出的实部，ik 是第 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);

分享到：

赞收藏

资料库

fftw中文参考手册.pdf

相关推荐

课程资源

热门标签

最新资料