我有一个三维数组，它有一个维度（Nx，Ny，Nz）。
我想使用FFWT3库沿Z轴应用实际FFT和IFFT。
这里，'z'是变化最快的指数。
我已经用python编写了相同的函数代码
numpy.fft.rfft和numpy.fft.irfft它的工作和我预期的完全一样。
但是太慢了所以我试着用C语言写代码。
我试图比较I FFT（FFT（f））和f的结果，其中f是任意数组。
我使用fft-plan-dft-r2c/fft-plan-dft-c2r来进行前向/后向fft。
这是我的密码。
（使用gcc和-lfftw3-lm选项在Ubuntu 16.04中编译）
#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()或类似的函数来分配它们，就像对FFTz所做的那样这些功能确保满足有关内存对齐的要求见fftw-3.3.6-pl2/kernel/kalloc.c中的函数*X(kernel_malloc)(size_t n)它调用memalign()或_aligned_malloc()等函数这两个都返回NULL就像失败时的malloc()一样。规模不同实际上，链接长度
Nz的前向和后向FFTW变换会导致缩放Nz。通常，允许某些FFTW转换（如c2r）覆盖其输入数组，除非添加了flag FFTW_PRESERVE_INPUT：
FFTW_PRESERVE_INPUT指定位置不正确的转换不能更改其输入数组这通常是默认的，除了c2r和hc2r（即复杂到真实）转换，其中FFTW_DESTROY_INPUT是默认的在后一种情况下，传递FFTW_PRESERVE_输入将尝试使用不破坏输入的算法，以降低性能；然而，对于多维c2r变换，没有实现输入保持算法，如果请求，规划器将返回NULL。
关于c - 3D阵列的1D实数FFT和IFFT，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/55488831/