-1

I'm trying to use fftw3.3.5 for DFT in C. but i currently have trouble while applying DFT to real data, it results in something strange, like extremly large numbers, as well as something that isn't symmetric.

here is the test code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>

//#define DEBUG

int main()
{
    fftw_complex *out;
    fftw_plan p;
    double *in;
    int N = 8;
    int i;

    in = (double*) fftw_malloc(sizeof(double) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

    for (i = 0; i < N; i++)
    {
        in[i] = 1.0;
    }

    p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
    fftw_execute(p);

    for (i = 0; i < N; i++)
    {
        printf("%6.2f+j%6.2f\n", out[i][0],out[i][1]);
    }
    printf("\n");
    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
    return 0;
}

with the compile argument:

gcc -o test.exe test.c -L../lib/fftw-3.3.5-dll32 -lfftw3-3 -lm

it supposed to give the output as:

 8+j0
 0+j0
 0+j0
 0+j0
 0+j0
 0+j0
 0+j0
 0+j0

but what i get is:

  8.00+j  0.00
  0.00+j  0.00
  0.00+j  0.00
  0.00+j  0.00
  0.00+j  0.00
403786757349850460000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.00+j944948400857723240000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.00
191040283242158910000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.00+j89520024970316631000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.00
11715519029813030000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.00+j49940585465016918000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.00

also if i try data with 128 samples like:

0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   -3,05175781250000e-05   0   -3,05175781250000e-05   0   -3,05175781250000e-05   0   -3,05175781250000e-05   0   -3,05175781250000e-05   0   -3,05175781250000e-05   0   -3,05175781250000e-05   0   0   -3,05175781250000e-05

i get (use sqrt(real^2+imag^2)/len*2) as amplitude:

3.81481e-006 3.71065e-006 3.40881e-006 2.93997e-006 2.35125e-006 1.70090e-006 1.05247e-006 4.74860e-007 1.86058e-007 4.41965e-007 6.28417e-007 6.74806e-007 5.98803e-007 4.41757e-007 2.77456e-007 2.51045e-007 3.64966e-007 4.63350e-007 4.86514e-007 4.21271e-007 2.76131e-007 7.43513e-008 1.51306e-007 3.64232e-007 5.29849e-007 6.21102e-007 6.23034e-007 5.36511e-007 3.84268e-007 2.43994e-007 3.02517e-007 4.96143e-007 6.74370e-007 7.82784e-007 7.98526e-007 7.19781e-007 5.70181e-007 4.21114e-007 4.20667e-007 5.92154e-007 7.92975e-007 9.36630e-007 9.83523e-007 9.20893e-007 7.62605e-007 5.64734e-007 4.74754e-007 6.25706e-007 8.81107e-007 1.10263e-006 1.22349e-006 1.21059e-006 1.05585e-006 7.81027e-007 4.81996e-007 5.15201e-007 9.35378e-007 1.42703e-006 1.88275e-006 2.25805e-006 2.53386e-006 2.71137e-006 2.80810e-006 2.85021e-006 2.86111e-006 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000

as my input is real, such a non-symmetry output should never happen...

Can anyone point out what i do wrong and how to fix it?

tyker1
  • 123
  • 8
  • Did you check if the out-vector is entirely written to, for example by zeroing it after allocation? – Ctx Jul 18 '18 at 10:25
  • i checked it, its not zero, as malloc dont initialize the memory..is that the problem? and why will that cause problem? – tyker1 Jul 18 '18 at 10:27
  • No, it is not really a problem if everything works as intended. But you could check, if fftw really fills the out vector completely by setting it to some sentinel value before calling the fftw functions. The last three values look quite random. – Ctx Jul 18 '18 at 10:29
  • i checked it too, it seems that it simply left the symmetric part unfilled.. – tyker1 Jul 18 '18 at 10:35

1 Answers1

2

From the manual:

fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out, unsigned flags);

Here, n is the “logical” size of the DFT, not necessarily the physical size of the array. In particular, the real (double) array has n elements, while the complex (fftw_complex) array has n/2+1 elements (where the division is rounded down)

So, the symmetric part is not written (it is redundant anyway).

In your case, this matches exactly the 5 elements you read correctly, the other 3 are left uninitialized.

valgrind can help you detect the usage of uninitialized values in the future.

Ctx
  • 18,090
  • 24
  • 36
  • 51