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
matlab fft and fftw
Started by ●May 18, 2011
Reply by ●May 18, 20112011-05-18
On May 18, 7:37�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: > � � 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.63the first term is suppose to be just the sum, so matlab is correct, fftw is not
Reply by ●May 18, 20112011-05-18
On May 18, 4:37�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
Reply by ●May 18, 20112011-05-18
>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"?
Reply by ●May 18, 20112011-05-18
On May 18, 4:12�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
Reply by ●May 19, 20112011-05-19
On May 18, 1:37�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: > � � 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;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
Reply by ●May 19, 20112011-05-19
>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"?
Reply by ●May 19, 20112011-05-19
On May 19, 2:09�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
Reply by ●June 21, 20112011-06-21
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