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
Using FFTW3 to calculate an Gaussian
Started by ●June 19, 2008
Reply by ●June 20, 20082008-06-20
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
Reply by ●June 20, 20082008-06-20
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
Reply by ●June 20, 20082008-06-20
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
Reply by ●June 21, 20082008-06-21
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
Reply by ●June 23, 20082008-06-23
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-39Shifting 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
Reply by ●June 24, 20082008-06-24
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
Reply by ●June 24, 20082008-06-24
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
Reply by ●June 24, 20082008-06-24
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
Reply by ●June 25, 20082008-06-25
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