verito - 1 year ago 63

C Question

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 Source

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.