verito verito - 3 months ago 18
C Question

Computing 3D FFT and Inverse FFT in C

I want to compute FFT and back transform to check if it works the same. I have an application of large 3D matrix in my code I tried to test it with 4*4*4 matrix and here is my code

`

#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;

}`


Apparently my transformed results are not the same. I don't understand what is going wrong?

Answer

Your FFT is not normalized. There is a constant factor between your input and output.

Take a look 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.

So the factor should be N * N * N. Just divide your data by this factor and you should get back the same data as your input.