在 C 中计算 3D FFT 和逆 FFT
Computing 3D FFT and Inverse FFT in C
我想计算 FFT 和反变换以检查它们是否相同。我的代码中有一个大型 3D 矩阵的应用程序,我尝试使用 4*4*4 矩阵对其进行测试,这是我的代码
`
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <time.h>
#include <math.h>
#include <fftw3.h>
int main()
{
int N = 4; //Dimension of matrix
unsigned int seed = 1;
double *in = (double*)malloc(sizeof(double)*N*N*N);
fftw_complex *out = fftw_malloc(sizeof(fftw_complex)*N*N*N);
double *out1 = (double*)malloc(sizeof(double)*N*N*N);
fftw_plan plan_backward;
fftw_plan plan_forward;
srand ( seed );
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
for(int k = 0; k < N; k++)
{
in[i*(N*N) + j*N + k] = rand ( );
}
}
}
printf(" Given matrix in\n");
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
for(int k = 0; k < N; k++)
{
printf("%f\n", in[i*(N*N) + j*N + k]);
}
}
}
printf("\n");
plan_backward = fftw_plan_dft_r2c_3d ( N, N, N, in, out, FFTW_ESTIMATE );
fftw_execute ( plan_backward );
fftw_destroy_plan ( plan_backward );
printf("out matrix\n");
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
for(int k = 0; k < N; k++)
{
printf("%f\t%f\n", creal(out[i*(N*N) + j*N + k]), cimag(out[i*(N*N) + j*N + k]));
}
}
}
printf("\n");
plan_forward = fftw_plan_dft_c2r_3d ( N, N, N, out, out1, FFTW_ESTIMATE );
fftw_execute ( plan_forward );
fftw_destroy_plan ( plan_forward );
printf("out1 matrix\n");
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
for(int k = 0; k < N; k++)
{
printf("%f\n", out1[i*(N*N) + j*N + k]);
}
}
}
fftw_free(in);
free(out);
fftw_free(out1);
return 0;
}`
显然我的转换结果不一样。我不明白这是怎么回事?
如FFTW3手册所述:
These transforms are unnormalized, so an r2c followed by a c2r transform (or vice versa) will result in the original data scaled by the number of real data elements—that is, the product of the (logical) dimensions of the real data.
在您的情况下,实际维度的数量是 4^3
,将第一个结果 115474520512
除以该数字返回第一个输入 115474520512/(4^3) = 1804289383
。
您的 FFT 未归一化。你的输入和输出之间有一个常数因子。
看看here
These transforms are unnormalized, so an r2c followed by a c2r transform (or vice versa) will result in the original data scaled by the number of real data elements—that is, the product of the (logical) dimensions of the real data.
所以系数应该是N * N * N
。只需将您的数据除以该系数,您应该会得到与输入相同的数据。
我想计算 FFT 和反变换以检查它们是否相同。我的代码中有一个大型 3D 矩阵的应用程序,我尝试使用 4*4*4 矩阵对其进行测试,这是我的代码
`
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <time.h>
#include <math.h>
#include <fftw3.h>
int main()
{
int N = 4; //Dimension of matrix
unsigned int seed = 1;
double *in = (double*)malloc(sizeof(double)*N*N*N);
fftw_complex *out = fftw_malloc(sizeof(fftw_complex)*N*N*N);
double *out1 = (double*)malloc(sizeof(double)*N*N*N);
fftw_plan plan_backward;
fftw_plan plan_forward;
srand ( seed );
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
for(int k = 0; k < N; k++)
{
in[i*(N*N) + j*N + k] = rand ( );
}
}
}
printf(" Given matrix in\n");
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
for(int k = 0; k < N; k++)
{
printf("%f\n", in[i*(N*N) + j*N + k]);
}
}
}
printf("\n");
plan_backward = fftw_plan_dft_r2c_3d ( N, N, N, in, out, FFTW_ESTIMATE );
fftw_execute ( plan_backward );
fftw_destroy_plan ( plan_backward );
printf("out matrix\n");
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
for(int k = 0; k < N; k++)
{
printf("%f\t%f\n", creal(out[i*(N*N) + j*N + k]), cimag(out[i*(N*N) + j*N + k]));
}
}
}
printf("\n");
plan_forward = fftw_plan_dft_c2r_3d ( N, N, N, out, out1, FFTW_ESTIMATE );
fftw_execute ( plan_forward );
fftw_destroy_plan ( plan_forward );
printf("out1 matrix\n");
for(int i = 0; i < N; i++)
{
for(int j = 0; j < N; j++)
{
for(int k = 0; k < N; k++)
{
printf("%f\n", out1[i*(N*N) + j*N + k]);
}
}
}
fftw_free(in);
free(out);
fftw_free(out1);
return 0;
}`
显然我的转换结果不一样。我不明白这是怎么回事?
如FFTW3手册所述:
These transforms are unnormalized, so an r2c followed by a c2r transform (or vice versa) will result in the original data scaled by the number of real data elements—that is, the product of the (logical) dimensions of the real data.
在您的情况下,实际维度的数量是 4^3
,将第一个结果 115474520512
除以该数字返回第一个输入 115474520512/(4^3) = 1804289383
。
您的 FFT 未归一化。你的输入和输出之间有一个常数因子。
看看here
These transforms are unnormalized, so an r2c followed by a c2r transform (or vice versa) will result in the original data scaled by the number of real data elements—that is, the product of the (logical) dimensions of the real data.
所以系数应该是N * N * N
。只需将您的数据除以该系数,您应该会得到与输入相同的数据。