Forums

Using FFTW3 to calculate an Gaussian

Started by jv_ju June 19, 2008
Hello,

I am a new user of FFTW, and I am using Fortran. My problem is to try to
get the Fourier Transform of a Gaussian function: f(x)=exp(-x**2). The
results should be another Gaussian: 1/sqrt(4*pi)*exp(-omega**2/4). I think
the out put files should be all real numbers, no imagine parts.

the parts doing FFTW in my code:

program test

  implicit none

  include "fftw3.f"

  integer N
  parameter(N=64)

  real*8, dimension(N) :: in
  real*8, dimension(N) :: out

  integer i
  real*8 :: x,xmax,xmin,dx
  integer ::  plan

  OPEN(UNIT=7,FILE='fftin')
  OPEN(UNIT=8,FILE='fftout')

 ! write(*,*) 'Input array:'

  xmax=10.
  xmin=-10.
  dx=(xmax-xmin)/N

  write(6,*) 'XMAX',xmax,xmin,dx,N
  do i = 1,N
     x=xmin+dx*i
     in(i)=exp(-x**2)
     write(7,*) x, in(i)
  enddo
  close(7)


   call dfftw_plan_r2r_1d (plan, N, in, out, fftw_r2hc, FFTW_ESTIMATE)
  !call dfftw_plan_dft_r2c_1d(plan,N,in,out,FFTW_ESTIMATE)
   !call dfftw_plan_dft_1d (plan, N, in, out, FFTW_FORWARD, FFTW_ESTIMATE
)

  call dfftw_execute(plan)

  !write(*,*) 'Output array:'

  do i = 1,N
     write(8,*) i,2*3.1415/10*i,out(i)
  enddo

  STOP
end
 I tried different plans( r2r,r2c ,etc.), but every time I can not get the
results as another Gaussian. I just wonder what's wrong?

Shuangyi



jv_ju wrote:

> I am a new user of FFTW, and I am using Fortran. My problem is to try to > get the Fourier Transform of a Gaussian function: f(x)=exp(-x**2). The > results should be another Gaussian: 1/sqrt(4*pi)*exp(-omega**2/4). I think > the out put files should be all real numbers, no imagine parts.
Remember, the FFT is periodic. Also, the zero term is the first element of the array, element 1 in your example. For a real transform, the input must be symmetric, which means in(i+1) must equal in(N+1-i) for arrays with origin 1. do i = 1,N/2+1 x=xmin+dx*(i-1) in(i)=exp(-x**2) if(i.gt.1) in(N+2-i)=in(i) write(7,*) x, in(i) enddo (For the discrete case it is an infinite train of gaussians transforms to an infinite train of gaussians.) -- glen
Thanks for the help.

Does it mean if I use dfftw_plan_dft_1d (plan, N, in, out, FFTW_FORWARD,
FFTW_ESTIMATE ), instead of dfftw_plan_r2r_1d, I still need to make the
input data symmetric?

If I increase the N like make N=1024, will I get a result with smaller
imagine parts?

Shuangyi

ain of gaussians
>transforms to an infinite train of gaussians.) > >-- glen > >
Hi, I tried to make the input data symmetric, but still get some negative output number. The following are my input and output: i input 1 3.720075976020800E-44 2 6.754396503887189E-39 3 5.614728092387960E-34 4 2.136871186570251E-29 5 3.723363121750522E-25 6 2.970297015029749E-21 7 1.084855264042939E-17 8 1.814057958631677E-14 9 1.388794386496399E-11 10 4.867793902108204E-09 11 7.811489408304497E-07 12 5.739088873946876E-05 13 1.930454136227708E-03 14 2.972921638615876E-02 15 0.2096113871510978 16 0.6766338461617289 17 1.000000000000000 18 0.6766338461617289 19 0.2096113871510978 20 2.972921638615876E-02 21 1.930454136227708E-03 22 5.739088873946876E-05 23 7.811489408304497E-07 24 4.867793902108204E-09 25 1.388794386496399E-11 26 1.814057958631677E-14 27 1.084855264042939E-17 28 2.970297015029749E-21 29 3.723363121750522E-25 30 2.136871186570251E-29 31 5.614728092387960E-34 32 6.754396503887189E-39 i output 1 2.835926161509187 2 -2.766808697386488 3 2.569400307801379 4 -2.271186709394578 5 1.910919236656792 6 -1.530383451170318 7 1.166612575289569 8 -0.8464896583457231 9 0.5846365717001540 10 -0.3843466132844360 11 0.2405191010038918 12 -0.1433000663466105 13 8.135894685384915E-02 14 -4.420947540452302E-02 15 2.346801579405655E-02 16 -1.327532866732173E-02 17 1.024432829142996E-02 18 0.000000000000000E+00 19 0.000000000000000E+00 20 0.000000000000000E+00 21 0.000000000000000E+00 22 0.000000000000000E+00 23 0.000000000000000E+00 24 0.000000000000000E+00 25 0.000000000000000E+00 26 0.000000000000000E+00 27 0.000000000000000E+00 28 0.000000000000000E+00 29 0.000000000000000E+00 30 0.000000000000000E+00 31 0.000000000000000E+00 32 0.000000000000000E+00
On Jun 20, 3:15 pm, "jv_ju" <shuangy...@hotmail.com> wrote:
> Hi, I tried to make the input data symmetric, but still get some negative > output number.
You are making a common mistake. For a DFT (what an FFT computes), the origin is not in the middle of the array, but at the beginning of the array (the 0th / DC element). Again, remember that the boundary conditions are periodic. So, this means that for symmetric data y[n], where n = 0 .. N-1, it should satisfy y[n] = y[N-n]. For example, here is an N=12 set of inputs containing a symmetric Gaussian: 1.0000e+00 3.6788e-01 1.8316e-02 1.2341e-04 1.1254e-07 1.3888e-11 2.3195e-16 1.3888e-11 1.1254e-07 1.2341e-04 1.8316e-02 3.6788e-01 The result of FFT'ing this data will be a Gaussian like you expect. At least, it will be approximately a Gaussian. The Fourier transform of a Gaussian is a Gaussian, but this is not true of the discrete Fourier transform, which ony approximates a Fourier integral. So, The result will only be approximately Gaussian. When you center your Gaussian at the middle of the array, by the shift theorem this multiplies the FFT by a linear phase, and in fact the phase that gets multiplied is an oscillating sign pattern (work it out). Regards, Steven G. Johnson
jv_ju wrote:

> Hi, I tried to make the input data symmetric, but still get some negative > output number.
Normally, it would be symmetric around the first input value.
> The following are my input and output:
> i input
> 1 3.720075976020800E-44 > 2 6.754396503887189E-39 > 3 5.614728092387960E-34 > 4 2.136871186570251E-29 > 5 3.723363121750522E-25 > 6 2.970297015029749E-21 > 7 1.084855264042939E-17 > 8 1.814057958631677E-14 > 9 1.388794386496399E-11 > 10 4.867793902108204E-09 > 11 7.811489408304497E-07 > 12 5.739088873946876E-05 > 13 1.930454136227708E-03 > 14 2.972921638615876E-02 > 15 0.2096113871510978 > 16 0.6766338461617289 > 17 1.000000000000000 > 18 0.6766338461617289 > 19 0.2096113871510978 > 20 2.972921638615876E-02 > 21 1.930454136227708E-03 > 22 5.739088873946876E-05 > 23 7.811489408304497E-07 > 24 4.867793902108204E-09 > 25 1.388794386496399E-11 > 26 1.814057958631677E-14 > 27 1.084855264042939E-17 > 28 2.970297015029749E-21 > 29 3.723363121750522E-25 > 30 2.136871186570251E-29 > 31 5.614728092387960E-34 > 32 6.754396503887189E-39
Shifting like this adds a pi term, such that successive transform output values change sign.
> i output > 1 2.835926161509187 > 2 -2.766808697386488 > 3 2.569400307801379 > 4 -2.271186709394578 > 5 1.910919236656792 > 6 -1.530383451170318 > 7 1.166612575289569 > 8 -0.8464896583457231 > 9 0.5846365717001540 > 10 -0.3843466132844360 > 11 0.2405191010038918 > 12 -0.1433000663466105 > 13 8.135894685384915E-02 > 14 -4.420947540452302E-02 > 15 2.346801579405655E-02 > 16 -1.327532866732173E-02 > 17 1.024432829142996E-02
(snip of zeros) Notice the nice shape if you change the sign on every other value of the output! -- glen
Hello,

Thanks for all the help, it seems the code works better, but the output
results seems to "sharp" and do not has the right amplitude, shall I divide
the output by N?

I changed the input for f(x)=exp(-x**2) as:
 i in(i)
 1 1.000000000000000 0.000000000000000E+00
 2 0.8406237433345053 0.000000000000000E+00
 3 0.4993517885992761 0.000000000000000E+00
 4 0.2096113871510978 0.000000000000000E+00
 5 6.217652402211627E-02 0.000000000000000E+00
 6 1.303290744850934E-02 0.000000000000000E+00
 7 1.930454136227708E-03 0.000000000000000E+00
 8 1.303290744850934E-02 0.000000000000000E+00
 9 6.217652402211627E-02 0.000000000000000E+00
 10 0.2096113871510978 0.000000000000000E+00
 11 0.4993517885992761 0.000000000000000E+00
 12 0.8406237433345053 0.000000000000000E+00

and here is my output:

 1 4.251523155247238 0.000000000000000E+00
 2 2.868674186073873 0.000000000000000E+00
 3 0.8748360179956545 0.000000000000000E+00
 4 0.1237190167094526 0.000000000000000E+00
 5 5.968265034016262E-03 0.000000000000000E+00
 6 1.815434807990846E-03 0.000000000000000E+00
 7 -1.548996489212140E-03 0.000000000000000E+00
 8 1.815434807990846E-03 0.000000000000000E+00
 9 5.968265034016262E-03 0.000000000000000E+00
 10 0.1237190167094526 0.000000000000000E+00
 11 0.8748360179956545 0.000000000000000E+00
 12 2.868674186073873 0.000000000000000E+00

Thanks!

Shuangyi


I also found a strange thing, if I change the N, the output results will
change. The following is the modified code:

  program test

  implicit none

  include "fftw3.f"

  integer N
  parameter(N=12)

  double complex in
  double complex out
  dimension out(N)
  dimension in(N)

  integer i
  real*8 :: x,xmax,xmin,dx
  integer*8  plan

  OPEN(UNIT=7,FILE='fftin')
  OPEN(UNIT=8,FILE='fftout')

 ! write(*,*) 'Input array:'

  xmax=10.
  xmin=0.
  dx=(xmax-xmin)/N

  write(6,*) 'XMAX',xmax,xmin,dx,N

   do i = 1,N/2+1
      x=xmin+dx*(i-1)
      in(i)=exp(-x**2)*(-1)**(i-1)
      if(i.gt.1) in(N+2-i)=in(i)
   enddo
   do i=1,N
      x=xmin+dx*(i-1)

      write(7,*) i,x,real(in(i)), imag(in(i))
   enddo

  close(7)



   call dfftw_plan_dft_1d (plan, N, in, out, FFTW_FORWARD, FFTW_ESTIMATE)

  call dfftw_execute(plan)

  !write(*,*) 'Output array:'

  do i = 1,N
 
     write(8,*) i, real(out(i)),imag(out(i))
  enddo
  close(8)
  STOP
end


On Jun 24, 5:12 pm, "jv_ju" <shuangy...@hotmail.com> wrote:

> I also found a strange thing, if I change the N, the output results will > change. > ...
You have already been given the reason: quote: ... The result of FFT'ing this data will be a Gaussian like you expect. At least, it will be approximately a Gaussian. The Fourier transform of a Gaussian is a Gaussian, but this is not true of the discrete Fourier transform, which ony approximates a Fourier integral. So, The result will only be approximately Gaussian. ... Regards, Steven G. Johnson end quote With different size transforms, you truncate the input Gaussian differently and get different results. For DFTs there are functions that transform more closely to themselves than the Gaussian. Dale B. Dalrymple
jv_ju wrote:

> I also found a strange thing, if I change the N,
> the output results will change. Remember, DFT (and so FFT) is periodic. You are not really transforming gaussian to gaussian, but a train of gaussians to a train of gaussians. The tails will be a little different, as the tail of the next one crosses into it. The effect will be more noticeable for smaller N. -- glen