Okay, what I'm trying to achieve is simple. Apply FFT on some random data and then apply the reverse algorithm on the output to get back the input. I'm using kissFFT
library for this.
Code:
const int fft_siz = 512;
const int inverse = 1;
kiss_fft_cpx* in = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * fft_siz);
kiss_fft_cpx* out = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * fft_siz);
kiss_fft_cpx* rec = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * fft_siz);
kiss_fft_cfg cfg = kiss_fft_alloc(fft_siz, !inverse, NULL, NULL);
kiss_fft_cfg icfg = kiss_fft_alloc(fft_siz, inverse, NULL, NULL);
srand((unsigned int)time(NULL));
for(int i = 0; i < fft_siz; i++)
{
in[i].r = rand() % 256;
in[i].i = rand() % 256;
}
kiss_fft(cfg, in, out);
// scaling
for(int i = 0; i < fft_siz; i++)
{
out[i].r /= fft_siz;
out[i].i /= fft_siz;
}
kiss_fft(icfg, out, rec);
unsigned int count = 0;
for(int i = 0; i < fft_siz; i++)
if(in[i].r != rec[i].r)
{
count++;
printf( "in[%3d].r does not match rec[%3d].r :: %3d :: %f\n",
i, i, count, in[i].r - rec[i].r);
}
else if(in[i].i != rec[i].i)
{
count++;
printf( "in[%3d].i does not match rec[%3d].i :: %3d :: %f\n",
i, i, count, in[i].i - rec[i].i);
}
free(in);
free(out);
free(rec);
free(cfg);
free(icfg);
kiss_fft_cleanup();
Output:
in[ 0]: 71.000000 85.000000 -- out[ 0]: 127.095703 124.541016
in[ 1]: 248.000000 27.000000 -- out[ 1]: -7.083314 0.072701
in[ 2]: 64.000000 18.000000 -- out[ 2]: -3.770610 2.682554
in[ 3]: 6.000000 96.000000 -- out[ 3]: -7.929140 -2.897723
in[ 4]: 98.000000 23.000000 -- out[ 4]: -0.719621 -5.854260
in[ 5]: 250.000000 188.000000 -- out[ 5]: 0.397226 -1.248124
in[ 6]: 231.000000 3.000000 -- out[ 6]: -7.934285 -2.367196
in[ 7]: 6.000000 105.000000 -- out[ 7]: -0.317480 -2.955601
in[ 8]: 172.000000 143.000000 -- out[ 8]: -4.236186 3.911616
in[ 9]: 16.000000 134.000000 -- out[ 9]: -0.162577 -5.353521
in[ 10]: 230.000000 112.000000 -- out[ 10]: -4.703711 7.791993
in[ 11]: 5.000000 26.000000 -- out[ 11]: -2.636305 0.188381
in[ 12]: 16.000000 127.000000 -- out[ 12]: 1.137413 4.576081
in[ 13]: 112.000000 86.000000 -- out[ 13]: 0.978051 -0.408992
in[ 14]: 40.000000 23.000000 -- out[ 14]: 5.231920 -2.347566
in[ 15]: 75.000000 26.000000 -- out[ 15]: 0.009981 -2.091559
note ::count::difference
--------------------------------------------------------
in[ 1].r does not match rec[ 1].r :: 1 :: -0.000031
in[ 3].r does not match rec[ 3].r :: 2 :: -0.000015
in[ 4].i does not match rec[ 4].i :: 3 :: -0.000004
in[ 6].i does not match rec[ 6].i :: 4 :: -0.000008
in[ 7].r does not match rec[ 7].r :: 5 :: -0.000002
in[ 9].r does not match rec[ 9].r :: 6 :: -0.000015
in[ 11].r does not match rec[ 11].r :: 7 :: -0.000015
in[ 12].r does not match rec[ 12].r :: 8 :: -0.000015
in[ 13].i does not match rec[ 13].i :: 9 :: -0.000008
in[ 14].i does not match rec[ 14].i :: 10 :: 0.000008
in[ 15].r does not match rec[ 15].r :: 11 :: -0.000015
Debug. If you go to the bottom, you'll see that there's 317 mismatches. I'm also outputting the difference between values i.e. (in[].r - rec[].r)
or (in[].i - rec[].i)
.
What I'm showing next is the input data where white dots represent the real part and red dots the imaginary part.
This is the output data of FFT represented in purple along with the reconstructed data in white and red.
Notice the small difference? I'm guessing this is related to floating point precision. How can I overcome this problem to get the exact same input data that I used FFT on ?