Forums

matlab fft and fftw

Started by Marc2050 May 18, 2011
Hi.

I did a simple test (1d fft) of some data with fftw and matlab.
The results seem to differ. I wonder what am I missing?
I need to port my codes from matlab to a custom program.
Many thanks for any help.

fftw:
    fftw_complex *in, *out;
    fftw_plan p;
    int N = 16;
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    in[0][0] = -157;
    in[1][0] = -157;
    in[2][0] = -163;
    in[3][0] = -162;
    in[4][0] = -153;
    in[5][0] = -147;
    in[6][0] = -146;
    in[7][0] = -148;
    in[8][0] = -152;
    in[9][0] = -152;
    in[10][0] = -153;
    in[11][0] = -154;
    in[12][0] = -156;
    in[13][0] = -161;
    in[14][0] = -160;
    in[15][0] = -153;
    for (int i = 0; i < N; i++)
        in[i][1] = 0;
    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p);

This is what I got:
-2475
-40.6284
8.89949
27.4691
3
-5.52795
-10.8995
-5.31273
-7
-5.31273
-10.8995
-5.52795
3
27.4691
8.89949
-40.6284

compare this to the following I got from fft in matlab
-2474.00
-39.63
9.90
28.47
4.00
-4.53
-9.90
-4.31
-6.00
-4.31
-9.90
-4.53
4.00
28.47
9.90
-39.63




On May 18, 7:37&#2013266080;am, "Marc2050" <maarcc@n_o_s_p_a_m.gmail.com> wrote:
> Hi. > > I did a simple test (1d fft) of some data with fftw and matlab. > The results seem to differ. I wonder what am I missing? > I need to port my codes from matlab to a custom program. > Many thanks for any help. > > fftw: > &#2013266080; &#2013266080; fftw_complex *in, *out; > &#2013266080; &#2013266080; fftw_plan p; > &#2013266080; &#2013266080; int N = 16; > &#2013266080; &#2013266080; in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); > &#2013266080; &#2013266080; out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); > &#2013266080; &#2013266080; in[0][0] = -157; > &#2013266080; &#2013266080; in[1][0] = -157; > &#2013266080; &#2013266080; in[2][0] = -163; > &#2013266080; &#2013266080; in[3][0] = -162; > &#2013266080; &#2013266080; in[4][0] = -153; > &#2013266080; &#2013266080; in[5][0] = -147; > &#2013266080; &#2013266080; in[6][0] = -146; > &#2013266080; &#2013266080; in[7][0] = -148; > &#2013266080; &#2013266080; in[8][0] = -152; > &#2013266080; &#2013266080; in[9][0] = -152; > &#2013266080; &#2013266080; in[10][0] = -153; > &#2013266080; &#2013266080; in[11][0] = -154; > &#2013266080; &#2013266080; in[12][0] = -156; > &#2013266080; &#2013266080; in[13][0] = -161; > &#2013266080; &#2013266080; in[14][0] = -160; > &#2013266080; &#2013266080; in[15][0] = -153; > &#2013266080; &#2013266080; for (int i = 0; i < N; i++) > &#2013266080; &#2013266080; &#2013266080; &#2013266080; in[i][1] = 0; > &#2013266080; &#2013266080; p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); > &#2013266080; &#2013266080; fftw_execute(p); > > This is what I got: > -2475 > -40.6284 > 8.89949 > 27.4691 > 3 > -5.52795 > -10.8995 > -5.31273 > -7 > -5.31273 > -10.8995 > -5.52795 > 3 > 27.4691 > 8.89949 > -40.6284 > > compare this to the following I got from fft in matlab > -2474.00 > -39.63 > 9.90 > 28.47 > 4.00 > -4.53 > -9.90 > -4.31 > -6.00 > -4.31 > -9.90 > -4.53 > 4.00 > 28.47 > 9.90 > -39.63
the first term is suppose to be just the sum, so matlab is correct, fftw is not
On May 18, 4:37&#2013266080;am, "Marc2050" <maarcc@n_o_s_p_a_m.gmail.com> wrote:
> Hi. > > I did a simple test (1d fft) of some data with fftw and matlab. > The results seem to differ. I wonder what am I missing? > I need to port my codes from matlab to a custom program. > Many thanks for any help. > ...
The data sequence is not dft-even so it does not have a real only dft. The values you list for Matlab appear to be the real components rounded to 2 decimal places. The values you list for fftw appear wrong by -1.0 for the real components with different rounding. Dale B. Dalrymple
>On May 18, 4:37=A0am, "Marc2050" <maarcc@n_o_s_p_a_m.gmail.com> wrote: >> Hi. >> >> I did a simple test (1d fft) of some data with fftw and matlab. >> The results seem to differ. I wonder what am I missing? >> I need to port my codes from matlab to a custom program. >> Many thanks for any help. >> ... > >The data sequence is not dft-even so it does not have a real only dft. >The values you list for Matlab appear to be the real components >rounded to 2 decimal places. The values you list for fftw appear wrong >by -1.0 for the real components with different rounding. > >Dale B. Dalrymple > >
Sorry for my ignorant, how do you make the input data I showed above to be "dft-even"?
On May 18, 4:12&#2013266080;pm, "Marc2050" <maarcc@n_o_s_p_a_m.gmail.com> wrote:
> >On May 18, 4:37=A0am, "Marc2050" <maarcc@n_o_s_p_a_m.gmail.com> wrote: > >> Hi. > > >> I did a simple test (1d fft) of some data with fftw and matlab. > >> The results seem to differ. I wonder what am I missing? > >> I need to port my codes from matlab to a custom program. > >> Many thanks for any help. > >> ... > > >The data sequence is not dft-even so it does not have a real only dft. > >The values you list for Matlab appear to be the real components > >rounded to 2 decimal places. The values you list for fftw appear wrong > >by -1.0 for the real components with different rounding. > > >Dale B. Dalrymple > > Sorry for my ignorant, how do you make the input data I showed above to be > "dft-even"?
See the first column, second page of: http://web.mit.edu/xiphmont/Public/windows.pdf Dale B. Dalrymple
On May 18, 1:37&#2013266080;pm, "Marc2050" <maarcc@n_o_s_p_a_m.gmail.com> wrote:
> Hi. > > I did a simple test (1d fft) of some data with fftw and matlab. > The results seem to differ. I wonder what am I missing? > I need to port my codes from matlab to a custom program. > Many thanks for any help. > > fftw: > &#2013266080; &#2013266080; fftw_complex *in, *out; > &#2013266080; &#2013266080; fftw_plan p; > &#2013266080; &#2013266080; int N = 16; > &#2013266080; &#2013266080; in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); > &#2013266080; &#2013266080; out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); > &#2013266080; &#2013266080; in[0][0] = -157;
You only initialize the real omponent of the input array. Depending on how fftw_malloc works, you might have to explicitly initialize the imaginary component to 0, or you will do computations that include numerical 'trash' in the un-initialized arrays. Rune
>On May 18, 1:37=A0pm, "Marc2050" <maarcc@n_o_s_p_a_m.gmail.com> wrote: >> Hi. >> >> I did a simple test (1d fft) of some data with fftw and matlab. >> The results seem to differ. I wonder what am I missing? >> I need to port my codes from matlab to a custom program. >> Many thanks for any help. >> >> fftw: >> =A0 =A0 fftw_complex *in, *out; >> =A0 =A0 fftw_plan p; >> =A0 =A0 int N =3D 16; >> =A0 =A0 in =3D (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); >> =A0 =A0 out =3D (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); >> =A0 =A0 in[0][0] =3D -157; > >You only initialize the real omponent of the input array. >Depending on how fftw_malloc works, you might have to >explicitly initialize the imaginary component to 0, or >you will do computations that include numerical 'trash' >in the un-initialized arrays. > >Rune >
I did do a for loop that initialize the imag component as shown in the codes. Or did I understand wrong about "initialize the real component"?
On May 19, 2:09&#2013266080;pm, "Marc2050" <maarcc@n_o_s_p_a_m.gmail.com> wrote:
> >On May 18, 1:37=A0pm, "Marc2050" <maarcc@n_o_s_p_a_m.gmail.com> wrote: > >> Hi. > > >> I did a simple test (1d fft) of some data with fftw and matlab. > >> The results seem to differ. I wonder what am I missing? > >> I need to port my codes from matlab to a custom program. > >> Many thanks for any help. > > >> fftw: > >> =A0 =A0 fftw_complex *in, *out; > >> =A0 =A0 fftw_plan p; > >> =A0 =A0 int N =3D 16; > >> =A0 =A0 in =3D (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); > >> =A0 =A0 out =3D (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); > >> =A0 =A0 in[0][0] =3D -157; > > >You only initialize the real omponent of the input array. > >Depending on how fftw_malloc works, you might have to > >explicitly initialize the imaginary component to 0, or > >you will do computations that include numerical 'trash' > >in the un-initialized arrays. > > >Rune > > I did do a for loop that initialize the imag component as shown in the > codes. Or did I understand wrong about "initialize the real component"?- Hide quoted text - > > - Show quoted text -
Sorry, I didn't see that loop. Rune
I just tried your code snippet with FFTW and gcc on a GNU/Linux
system, and I can't reproduce your error.  I get the same results as
the Matlab results that you quote (indeed, Matlab uses FFTW
internally).  So, possibly you have a bug elsewhere in your program?
I have attached my complete program below (compile with gcc -std=c99
marc2050.c -lfftw3 -lm).  If running this program gives different
results for you, please send an email to fftw@fftw.org in case you
have encountered a compiler bug or something of that sort.

Regards,
Steven G. Johnson

----  test program, based on the code snippet you posted ----
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>
int main(void)
{
     fftw_complex *in, *out;
     fftw_plan p;
     int N = 16;
     in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
     out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
     in[0][0] = -157;
     in[1][0] = -157;
     in[2][0] = -163;
     in[3][0] = -162;
     in[4][0] = -153;
     in[5][0] = -147;
     in[6][0] = -146;
     in[7][0] = -148;
     in[8][0] = -152;
     in[9][0] = -152;
     in[10][0] = -153;
     in[11][0] = -154;
     in[12][0] = -156;
     in[13][0] = -161;
     in[14][0] = -160;
     in[15][0] = -153;
     for (int i = 0; i < N; i++)
	  in[i][1] = 0;
     p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
     fftw_execute(p);
     for (int i = 0; i < N; i++)
	  printf("out[%d] = %g%+gi\n", i, out[i][0], out[i][1]);
     fftw_destroy_plan(p);
     return 0;
}

Output:

out[0] = -2474+0i
out[1] = -39.6284-11.3717i
out[2] = 9.89949+21.3137i
out[3] = 28.4691+2.46767i
out[4] = 4+0i
out[5] = -4.52795+2.12453i
out[6] = -9.89949+1.31371i
out[7] = -4.31273+0.28515i
out[8] = -6+0i
out[9] = -4.31273-0.28515i
out[10] = -9.89949-1.31371i
out[11] = -4.52795-2.12453i
out[12] = 4+0i
out[13] = 28.4691-2.46767i
out[14] = 9.89949-21.3137i
out[15] = -39.6284+11.3717i