Reply by Steven G. Johnson●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
Reply by Rune Allnor●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 Marc2050●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 Rune Allnor●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 dbd●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"?
>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 dbd●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 steve●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.63
the first term is suppose to be just the sum, so matlab is correct,
fftw is not
Reply by Marc2050●May 18, 20112011-05-18
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