3D 阵列的 1D 实数 FFT 和 IFFT
1D real FFT and IFFT of a 3D array
我有一个维度为 (Nx, Ny, Nz) 的 3D 数组。
我想使用 FFTW3 库沿 z 轴应用真实的 FFT 和 IFFT。
这里，'z' 是变化最快的指数。
我已经使用
编写了与python相同的功能代码
numpy.fft.rfft 和 numpy.fft.irfft。它完全符合我的预期。
但是速度太慢了。于是尝试用C语言写代码
我尝试比较 IFFT(FFT(f)) 和 f 的结果，其中 f 是任意数组。
我使用 fft_plan_many_dft_r2c / fft_plan_many_dft_c2r 进行正向/反向 FFT。
这是我的代码。
（在 Ubuntu 16.04 中使用 gcc 和 -lfftw3 -lm 选项编译）
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>
void rfft_1d_of_3d(int Nx, int Ny, int Nz, double* input, double* output);
int main(){
double* input;
double* output;
int Nx=2, Ny=4, Nz=8;
int i,j,k,idx;
input = (double*) malloc(Nx*Ny*Nz*sizeof(double));
output = (double*) malloc(Nx*Ny*Nz*sizeof(double));
for(i=0; i<Nx; i++){
for(j=0; j<Ny; j++){
for(k=0; k<Nz; k++){
idx = k + j*Nz + i*Nz*Ny;
input[idx] = idx;
output[idx] = 0.;
}
}
}
rfft_1d_of_3d(Nx, Ny, Nz, input, output);
for(i=0; i<Nx; i++){
for(j=0; j<Ny; j++){
for(k=0; k<Nz; k++){
idx = k + j*Nz + i*Nz*Ny;
printf("%f, %f\n", input[idx], output[idx]);
}
}
}
return 0;
}
void rfft_1d_of_3d(int Nx, int Ny, int Nz, double* input, double* output){
int i,j,k,idx;
// Allocate memory space.
fftw_complex *FFTz = fftw_alloc_complex(Nx*Ny*(Nz/2+1));
// Set forward FFTz parameters
int rankz = 1;
int nz[] = {Nz};
const int *inembedz = NULL, *onembedz = NULL;
int istridez = 1, ostridez = 1;
int idistz = Nz, odistz= (Nz+2)/2;
int howmanyz = (Nx*Ny);
// Setup Forward plans.
fftw_plan FFTz_for_plan = fftw_plan_many_dft_r2c(rankz, nz, howmanyz, input, inembedz, istridez, idistz, FFTz, onembedz, ostridez, odistz, FFTW_ESTIMATE);
// Set backward FFTz parameters
int rankbz = 1;
int nbz[] = {(Nz+2)/2};
const int *inembedbz = NULL, *onembedbz = NULL;
int istridebz = 1, ostridebz = 1;
int idistbz = (Nz+2)/2, odistbz = Nz;
int howmanybz = (Nx*Ny);
// Setup Backward plans.
fftw_plan FFTz_bak_plan = fftw_plan_many_dft_c2r(rankbz, nbz, howmanybz, FFTz, inembedbz, istridebz, idistbz, output, onembedbz, ostridebz, odistbz, FFTW_ESTIMATE);
fftw_execute(FFTz_for_plan);
fftw_execute(FFTz_bak_plan);
fftw_free(FFTz);
return;
}
输入和输出结果应该是一样的，但实际上不是。
输入数组是一个 3D 数组（Nx=2，Ny=4，Nz=8），
[[[ 0. 1. 2. 3. 4. 5. 6. 7.]
[ 8. 9. 10. 11. 12. 13. 14. 15.]
[ 16. 17. 18. 19. 20. 21. 22. 23.]
[ 24. 25. 26. 27. 28. 29. 30. 31.]]
[[ 32. 33. 34. 35. 36. 37. 38. 39.]
[ 40. 41. 42. 43. 44. 45. 46. 47.]
[ 48. 49. 50. 51. 52. 53. 54. 55.]
[ 56. 57. 58. 59. 60. 61. 62. 63.]]]
我把它压平了，
[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63.]
我预计输出与输入相同，但实际结果是
0.000000, 12.000000
1.000000, 8.929290
2.000000, 28.256139
3.000000, 35.743861
4.000000, 55.070710
5.000000, 0.000000
6.000000, 0.000000
7.000000, 0.000000
8.000000, 76.000000
9.000000, 72.929290
10.000000, 92.256139
11.000000, 99.743861
12.000000, 119.070710
13.000000, 0.000000
14.000000, 0.000000
15.000000, 0.000000
16.000000, 140.000000
17.000000, 136.929290
18.000000, 156.256139
19.000000, 163.743861
20.000000, 183.070710
21.000000, 0.000000
22.000000, 0.000000
23.000000, 0.000000
24.000000, 204.000000
25.000000, 200.929290
26.000000, 220.256139
27.000000, 227.743861
28.000000, 247.070710
29.000000, 0.000000
30.000000, 0.000000
31.000000, 0.000000
32.000000, 268.000000
33.000000, 264.929290
34.000000, 284.256139
35.000000, 291.743861
36.000000, 311.070710
37.000000, 0.000000
38.000000, 0.000000
39.000000, 0.000000
40.000000, 332.000000
41.000000, 328.929290
42.000000, 348.256139
43.000000, 355.743861
44.000000, 375.070710
45.000000, 0.000000
46.000000, 0.000000
47.000000, 0.000000
48.000000, 396.000000
49.000000, 392.929290
50.000000, 412.256139
51.000000, 419.743861
52.000000, 439.070710
53.000000, 0.000000
54.000000, 0.000000
55.000000, 0.000000
56.000000, 460.000000
57.000000, 456.929290
58.000000, 476.256139
59.000000, 483.743861
60.000000, 503.070710
61.000000, 0.000000
62.000000, 0.000000
63.000000, 0.000000
左边是输入数组的元素，右边是输出数组的元素。
我哪里弄错了？
对于fftw_plan_many_dft_c2r()：nbz，变换的大小，要设置为Nz，实际输出的大小 大批。 fftw_plan_dft_c2r_1d().
可能类似
// Set backward FFTz parameters
int rankbz = 1;
int nbz[1] = {(Nz)}; // here !!!
const int *inembedbz = NULL, *onembedbz = NULL;
//int inembedbz[1]={Nz/2+1};
//int onembedbz[1]={Nz};
int istridebz = 1, ostridebz = 1;
int idistbz = (Nz/2+1), odistbz = (Nz);
int howmanybz = (Nx*Ny);
// Setup Backward plans.
fftw_plan FFTz_bak_plan = fftw_plan_many_dft_c2r(rankbz, nbz, howmanybz, FFTz, inembedbz, istridebz, idistbz, output, onembedbz, ostridebz, odistbz, FFTW_ESTIMATE);
由于 FFTW 变换使用 input 和 output 数组，建议使用 fftw_malloc() 或类似函数分配它们，就像您为 [=18= 所做的那样].这些函数确保满足有关内存对齐的要求。请参阅 fftw-3.3.6-pl2/kernel/kalloc.c 中的函数 *X(kernel_malloc)(size_t n)。它会调用 memalign() or _aligned_malloc() 等函数。这两个 return NULL 就像 malloc() 以防失败。
有尺度上的差异。实际上，链接长度为 Nz 的正向和反向 FFTW 变换会导致缩放 Nz。
最后，允许某些 FFTW 变换（例如 c2r）覆盖其输入数组，除非添加了 flag FFTW_PRESERVE_INPUT：
FFTW_PRESERVE_INPUT specifies that an out-of-place transform must not change its input array. This is ordinarily the default, except for c2r and hc2r (i.e. complex-to-real) transforms for which FFTW_DESTROY_INPUT is the default. In the latter cases, passing FFTW_PRESERVE_INPUT will attempt to use algorithms that do not destroy the input, at the expense of worse performance; for multi-dimensional c2r transforms, however, no input-preserving algorithms are implemented and the planner will return NULL if one is requested.
我有一个维度为 (Nx, Ny, Nz) 的 3D 数组。
我想使用 FFTW3 库沿 z 轴应用真实的 FFT 和 IFFT。
这里，'z' 是变化最快的指数。
我已经使用
编写了与python相同的功能代码numpy.fft.rfft 和 numpy.fft.irfft。它完全符合我的预期。
但是速度太慢了。于是尝试用C语言写代码
我尝试比较 IFFT(FFT(f)) 和 f 的结果，其中 f 是任意数组。
我使用 fft_plan_many_dft_r2c / fft_plan_many_dft_c2r 进行正向/反向 FFT。
这是我的代码。
（在 Ubuntu 16.04 中使用 gcc 和 -lfftw3 -lm 选项编译）
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>
void rfft_1d_of_3d(int Nx, int Ny, int Nz, double* input, double* output);
int main(){
double* input;
double* output;
int Nx=2, Ny=4, Nz=8;
int i,j,k,idx;
input = (double*) malloc(Nx*Ny*Nz*sizeof(double));
output = (double*) malloc(Nx*Ny*Nz*sizeof(double));
for(i=0; i<Nx; i++){
for(j=0; j<Ny; j++){
for(k=0; k<Nz; k++){
idx = k + j*Nz + i*Nz*Ny;
input[idx] = idx;
output[idx] = 0.;
}
}
}
rfft_1d_of_3d(Nx, Ny, Nz, input, output);
for(i=0; i<Nx; i++){
for(j=0; j<Ny; j++){
for(k=0; k<Nz; k++){
idx = k + j*Nz + i*Nz*Ny;
printf("%f, %f\n", input[idx], output[idx]);
}
}
}
return 0;
}
void rfft_1d_of_3d(int Nx, int Ny, int Nz, double* input, double* output){
int i,j,k,idx;
// Allocate memory space.
fftw_complex *FFTz = fftw_alloc_complex(Nx*Ny*(Nz/2+1));
// Set forward FFTz parameters
int rankz = 1;
int nz[] = {Nz};
const int *inembedz = NULL, *onembedz = NULL;
int istridez = 1, ostridez = 1;
int idistz = Nz, odistz= (Nz+2)/2;
int howmanyz = (Nx*Ny);
// Setup Forward plans.
fftw_plan FFTz_for_plan = fftw_plan_many_dft_r2c(rankz, nz, howmanyz, input, inembedz, istridez, idistz, FFTz, onembedz, ostridez, odistz, FFTW_ESTIMATE);
// Set backward FFTz parameters
int rankbz = 1;
int nbz[] = {(Nz+2)/2};
const int *inembedbz = NULL, *onembedbz = NULL;
int istridebz = 1, ostridebz = 1;
int idistbz = (Nz+2)/2, odistbz = Nz;
int howmanybz = (Nx*Ny);
// Setup Backward plans.
fftw_plan FFTz_bak_plan = fftw_plan_many_dft_c2r(rankbz, nbz, howmanybz, FFTz, inembedbz, istridebz, idistbz, output, onembedbz, ostridebz, odistbz, FFTW_ESTIMATE);
fftw_execute(FFTz_for_plan);
fftw_execute(FFTz_bak_plan);
fftw_free(FFTz);
return;
}
输入和输出结果应该是一样的，但实际上不是。
输入数组是一个 3D 数组（Nx=2，Ny=4，Nz=8），
[[[ 0. 1. 2. 3. 4. 5. 6. 7.]
[ 8. 9. 10. 11. 12. 13. 14. 15.]
[ 16. 17. 18. 19. 20. 21. 22. 23.]
[ 24. 25. 26. 27. 28. 29. 30. 31.]]
[[ 32. 33. 34. 35. 36. 37. 38. 39.]
[ 40. 41. 42. 43. 44. 45. 46. 47.]
[ 48. 49. 50. 51. 52. 53. 54. 55.]
[ 56. 57. 58. 59. 60. 61. 62. 63.]]]
我把它压平了，
[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63.]
我预计输出与输入相同，但实际结果是
0.000000, 12.000000
1.000000, 8.929290
2.000000, 28.256139
3.000000, 35.743861
4.000000, 55.070710
5.000000, 0.000000
6.000000, 0.000000
7.000000, 0.000000
8.000000, 76.000000
9.000000, 72.929290
10.000000, 92.256139
11.000000, 99.743861
12.000000, 119.070710
13.000000, 0.000000
14.000000, 0.000000
15.000000, 0.000000
16.000000, 140.000000
17.000000, 136.929290
18.000000, 156.256139
19.000000, 163.743861
20.000000, 183.070710
21.000000, 0.000000
22.000000, 0.000000
23.000000, 0.000000
24.000000, 204.000000
25.000000, 200.929290
26.000000, 220.256139
27.000000, 227.743861
28.000000, 247.070710
29.000000, 0.000000
30.000000, 0.000000
31.000000, 0.000000
32.000000, 268.000000
33.000000, 264.929290
34.000000, 284.256139
35.000000, 291.743861
36.000000, 311.070710
37.000000, 0.000000
38.000000, 0.000000
39.000000, 0.000000
40.000000, 332.000000
41.000000, 328.929290
42.000000, 348.256139
43.000000, 355.743861
44.000000, 375.070710
45.000000, 0.000000
46.000000, 0.000000
47.000000, 0.000000
48.000000, 396.000000
49.000000, 392.929290
50.000000, 412.256139
51.000000, 419.743861
52.000000, 439.070710
53.000000, 0.000000
54.000000, 0.000000
55.000000, 0.000000
56.000000, 460.000000
57.000000, 456.929290
58.000000, 476.256139
59.000000, 483.743861
60.000000, 503.070710
61.000000, 0.000000
62.000000, 0.000000
63.000000, 0.000000
左边是输入数组的元素，右边是输出数组的元素。
我哪里弄错了？
对于fftw_plan_many_dft_c2r()：nbz，变换的大小，要设置为Nz，实际输出的大小 大批。 fftw_plan_dft_c2r_1d().
// Set backward FFTz parameters
int rankbz = 1;
int nbz[1] = {(Nz)}; // here !!!
const int *inembedbz = NULL, *onembedbz = NULL;
//int inembedbz[1]={Nz/2+1};
//int onembedbz[1]={Nz};
int istridebz = 1, ostridebz = 1;
int idistbz = (Nz/2+1), odistbz = (Nz);
int howmanybz = (Nx*Ny);
// Setup Backward plans.
fftw_plan FFTz_bak_plan = fftw_plan_many_dft_c2r(rankbz, nbz, howmanybz, FFTz, inembedbz, istridebz, idistbz, output, onembedbz, ostridebz, odistbz, FFTW_ESTIMATE);
由于 FFTW 变换使用 input 和 output 数组，建议使用 fftw_malloc() 或类似函数分配它们，就像您为 [=18= 所做的那样].这些函数确保满足有关内存对齐的要求。请参阅 fftw-3.3.6-pl2/kernel/kalloc.c 中的函数 *X(kernel_malloc)(size_t n)。它会调用 memalign() or _aligned_malloc() 等函数。这两个 return NULL 就像 malloc() 以防失败。
有尺度上的差异。实际上，链接长度为 Nz 的正向和反向 FFTW 变换会导致缩放 Nz。
最后，允许某些 FFTW 变换（例如 c2r）覆盖其输入数组，除非添加了 flag FFTW_PRESERVE_INPUT：
FFTW_PRESERVE_INPUT specifies that an out-of-place transform must not change its input array. This is ordinarily the default, except for c2r and hc2r (i.e. complex-to-real) transforms for which FFTW_DESTROY_INPUT is the default. In the latter cases, passing FFTW_PRESERVE_INPUT will attempt to use algorithms that do not destroy the input, at the expense of worse performance; for multi-dimensional c2r transforms, however, no input-preserving algorithms are implemented and the planner will return NULL if one is requested.