Confused by FFTW output

From: Tom (flurboglarf_at_mailinator.com)
Date: 03/15/05 

Date: 15 Mar 2005 08:00:52 -0800

Hi,
I am trying to get acquainted with the latest FFTW and wrote a small
test program to ensure correct usage by comparing a computed spectrum
with an analytical solution. For a start, I chose a Gaussian,
f(t)=exp(-a*t^2), whose spectrum is a purely real function, namely
F(w)=sqrt(pi/(4a))*exp(-w^2/(4a)). However, when computing the spectrum
with FFTW, I get a nonvanishing imaginary part in the transform, the
real part is definitely something different (even when allowing for
rounding errors), and it is the absolute value which reproduces the
closed-form solution. For further checking, I compared the derivative
of f(t) with the result of the inverse FT of iw*F(w) (i being the
imaginary unit), and the result is correct, so it seems to me that I do
call the routine in the right way. But I don't understand why the
numerical spectrum does not match the analytical result.
Does someone have an idea what I get wrong? My test code is appended
below; I tried different values for n to ensure that aliasing is not an
issue.
Thanks,
Tom

program testfft
implicit none
include "fftw3.f"
integer, parameter :: n=32
integer :: i
integer*8 :: plan
real, parameter :: pi=3.1415926535898
real :: t,dt,f,df,ranlsp,a,ttot
real, dimension(n) :: sig
complex, dimension(n/2+1) :: spec
a=0.2
ttot=10.
dt=ttot/float(n-1)
df=1./(2.*dt*float(n/2))
! signal
do i=1,n
   t=(i-1)*dt-ttot/2
   sig(i)=exp(-a*t*t)
   write(10,*) i,t,sig(i)
end do
! forward FT
call sfftw_plan_dft_r2c_1d(plan,n,sig,spec,FFTW_ESTIMATE)
call sfftw_execute(plan)
call sfftw_destroy_plan(plan)
spec=spec/sqrt(float(n)) ! normalize
! numerical spectrum and analytic solution (ranlsp)
do i=1,n/2+1
   f=(i-1)*df
   ranlsp=0.5*sqrt(pi/a)*exp(-(2.*pi*f)**2/(4.*a))
   write(20,*) i,f,real(spec(i)),aimag(spec(i)),abs(spec(i)),ranlsp
   spec(i)=2.*pi*f*(0.,1.)*spec(i) ! derivative
end do
! inverse FT
call sfftw_plan_dft_c2r_1d(plan,n,spec,sig,FFTW_ESTIMATE)
call sfftw_execute(plan)
call sfftw_destroy_plan(plan)
sig=sig/sqrt(float(n)) ! normalize
! derivative of signal (numerical and analytic)
do i=1,n
   t=(i-1)*dt-ttot/2
   write(30,*) i,t,sig(i),-2.*a*t*exp(-a*t*t)
end do
end program testfft