DFT of sine(x) using FFTW in Fortran muddled output

Question

The following is the code I have written to find the DFT of sine(x) over a period.

program fftw_test

implicit none

INTEGER FFTW_MEASURE
    PARAMETER (FFTW_MEASURE=0)
INTEGER FFTW_ESTIMATE
    PARAMETER (FFTW_ESTIMATE=64)
INTEGER FFTW_FORWARD
    PARAMETER (FFTW_FORWARD=-1)

integer, parameter :: n = 8
integer :: i
double complex, dimension(0:n-1) :: input, output
double precision, parameter :: pi = 3.141592653, h = 2.0d0*pi/(n)
integer*8 :: plan

call dfftw_plan_dft_1d(plan, n, input, output, fftw_forward, fftw_measure)

do i = 0, n-1
    input(i) = cmplx(sin(h*i), 0)
end do

call dfftw_execute_dft(plan, input, output)

output = output/n 

output(0) = cmplx(0,0) ! setting oddball wavenumber to be 0

call dfftw_destroy_plan(plan)

do i = -n/2, n/2-1, 1
    write(*, *) i, output(i+(n/2))
end do

end program

I am aware of the r2c (real to complex) function in the FFTW library. But I was advised to use the normal c2c function. So I defined the input function as a complex number with real part = sine(x) and complex part 0.

The DFT of sine(x) is supposed to be fk(-1) = cmplx(0, -0.5) and fk(1) = cmplx(0, 0.5) where fk(k) means the fourier coefficient of the k wavenumber

The output I received is as follows.

      -4 (  0.0000000000000000     ,  0.0000000000000000     )
      -3 (  3.2001271327131496E-008,-0.49999998518472011     )
      -2 ( -1.0927847071684482E-008,  1.4901161193847656E-008)
      -1 ( -1.0145577183762535E-008,  1.4815279864022202E-008)
       0 ( -1.0927847071684482E-008,  0.0000000000000000     )
       1 ( -1.0145577183762535E-008, -1.4815279864022202E-008)
       2 ( -1.0927847071684482E-008, -1.4901161193847656E-008)
       3 (  3.2001271327131496E-008, 0.49999998518472011     )

I am getting fk(-3) = cmplx(~0, -0.5) and fk(3) = cmplx(~0, 0.5). If I increase the grid size to 16, 32 or so I get -n/2 -1 and n/2 -1 wavenumbers with the required values instead of the -1 and 1 wavenumbers.

Does this have something to do with the way FFTW stores the output in the output array ? Or am I going wrong anywhere else ?

Also, I don't seem to be getting 'proper 0' where I should be. It is instead numbers of the order of 10^(-8) which I believe is the smallest my datatype double can hold. Is that something I should be worried about ?

10^(-8) is certainly not the smallest value a double can hold. You must however take into account how FFT works and the relative difference of these values to your maximum value. — Vladimir F
The guy who advised you to use the c2c transforms had some reasons for that? You must also take into account, that the negative frequencies are stored in the upper half of the array. — Vladimir F
If your input data are purely real, it's almost certainly more efficient to use r2c or r2hc (real to half-complex). It seems strange that someone would recommend c2c in this case. — Yossarian

Stefan Stefan · Accepted Answer · 2015-02-12T10:26:19

Like @VladimirF already said, the ordering of the values is a bit different, than you might expect. The first half of the array holds the positive frequencies, the second half holds the negative frequencies in reverse order (see this link). And you might have to check the sign convention used by FFTW.

The problem with accuracy stems from your single precision value for pi and the use of cmplx which produces single precision complex numbers (use the keyword argument kind). In this case you could simply assign your real value to the complex variables. Applying these two changes yields a precision of ~1e-10. This can be improved by supplying a better approximation for pi (i.e. more than 10 digits).

E.g. the value pi = 3.141592653589793d0 yields results with accuracy of 1e-16.

DFT of sine(x) using FFTW in Fortran muddled output

1 Answers