Thanks for the help.
But why when I increase N, the results curve then get "sharper"? If my N
is larger, that means I have more terms, shouldn't the results look more
like a Gaussian? why it only make the curve narrower than expected? Is it
mean something still wrong with my code?
Reply by glen herrmannsfeldt●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
Reply by dbd●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 jv_ju●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 jv_ju●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 glen herrmannsfeldt●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.
(snip of zeros)
Notice the nice shape if you change the sign on
every other value of the output!
-- glen
Reply by Steven G. Johnson●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 jv_ju●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 jv_ju●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 glen herrmannsfeldt●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