Doubt about a discrete Fourier Transform

[Telemachus]

Hi. I was checking the library for the discrete Fourier transform, fftw. So, I was using a functition $f(x)=sin(kx)$, which when transformed must give a delta function in k. When I transform, and then transform back, I effectively recover the function, so I think I am doing something right. But when I try to plot the function in reciprocal space, the delta function doesn't appear where it should be.

Homework Equations

The discrete Fourier transform of a function $f(x)$ is $\hat f_k=h \sum_{j=1}^N e^{-ikx_j}f(x_j)$ $k=0,1,2...,N-1$.
The inverse transform is:
$f_j=\frac{1}{2\pi}\sum_{k=0}^N e^{ikx_j}\hat f_k$, $j=1,2...,N$

The Attempt at a Solution

For example, here I have used the case with k=2:

I am using for the discrete case: $x=i*\Delta x$, with i from 0 to N, and similarly $k=i\Delta x$. Is my choice for the discrete k wrong? and how should I fix this?

 

[Telemachus]

I think I've fixed it by taking $\Delta k=\frac{2\pi}{L}$, where L is the length of the interval. Now, I was trying to represent a function using this Fourier transform, and the difficulty I encounter now is that I have a real function expressed as a set of real numbers and complex exponential. Is there anyway to visualize this as a real function?

 

[scottdave]

It looks like you are using k as the frequency of sine, and also using k as the index of the summation?

 

[Telemachus]

Yes, what is wrong with that?

 

[vela]

You see nothing wrong with $\hat{f}_k = \sum_{j=0}^{N-1} e^{-ikx}\sin(kx)$ in the context of the problem you were talking about in the original post?

 

[Telemachus]

Sorry, I don't see what you mean. Perhaps I'm committing the same mistake in this problem.

I have tried to expand a function in a Fourier series: I have $f(x)=e^{\sin(x)}(\sin(x)^2+\sin(x))$.

I have computed the Fourier coefficients $f_k$ using the fftw in fortran. Then I approximate $\displaystyle f(x_j) \sim \frac{1}{N} \sum_{k=0}^{N-1}f_k e^{ikx_j}$ where $x_j=j\frac{2\pi}{N}$. Where N is the number of points in the grid. This is the result I obtain for N=64:

If I increase N, I obtain a better representation.

However, due to the fact that the function is smooth, I hoped to have much better convergence than what I do.

I also have an analytic solution of a differential equation:
$u(x)-u''(x)=f(x)$
With periodic boundary conditions: $u(x+2\pi)=u(x)$

Using a Fourier expansion $\displaystyle f(x)=e^{\sin(x)}(\sin(x)^2+\sin(x))=\sum_{k=-\infty}^{\infty}\hat f_k e^{ikx}$, and similarly for $\displaystyle u(x)=\sum_{k=-\infty}^{\infty}\hat u_k e^{ikx}$

I arrived to: $\displaystyle u(x)=\sum_{k=-\infty}^{\infty}\frac{1}{1+k^2}\hat f_k e^{ikx}$.

In the discrete approximation I've used:
$\displaystyle u(x_j)\sim \frac{1}{N}\sum_{k=0}^{N-1}\frac{1}{1+k^2}\hat f_k e^{ikx_j}, \, x_j=j\frac{2\pi}{N}$ (1)

Now there is some problem, and I don't know what it is. The analytic solution is: $u(x)=e^{sin(x)}$, here is a plot of both, the analytic, and the one obtained by the Fourier series (1):

the problem is possibly in the summation, as the solution is purely real, but if I plot the imaginary part for the resultant of the summation (1) I get:

I have also tried to run with k from $-N/2$ to $N/2$, but I get zero.

This is the loop I use for the last problem:

     do j=1,N
       x = j*dx
       kv=zero
       do i=1,N
        kn=i-1
        ct=one/(N*(one+kn**2))
        kv=ct*zexp(im*kn*x)*zout(i)+kv
       enddo
      user = dreal(kv)
      cuser = dimag(kv)
      write(16,*) x,' ',user
      write(18,*) x,' ',cuser
     enddo

zout(i) are the Fourier coefficients obtained by fftw, and I am is the imaginary unit.

I don't know if this is related to the problem you've both mentioned before, or if this has something to do with what is explained in here with respect to the ordering of the Fourier coefficients: http://www.fftw.org/fftw3_doc/The-1d-Discrete-Fourier-Transform_0028DFT_0029.html

Note also that we use the standard “in-order” output ordering—the k-th output corresponds to the frequency k/n (or k/T, where T is your total sampling period). For those who like to think in terms of positive and negative frequencies, this means that the positive frequencies are stored in the first half of the output and the negative frequencies are stored in backwards order in the second half of the output. (The frequency -k/n is the same as the frequency (n-k)/n.)

 

[Telemachus]

You see nothing wrong with $\hat{f}_k = \sum_{j=0}^{N-1} e^{-ikx}\sin(kx)$ in the context of the problem you were talking about in the original post?

Is this related to the fact that the function is $2\pi$ periodic? that gives a bound in k, and actually set it to be a discrete set even for continuous x, is that it?

 

[vela]

No, it's that the $k$ in $\sin kx$ isn't the same $k$ as the $k$ in $e^{-ikx}$.

 

[Telemachus]

I wasn't aware of that. From Euler Formula I have:

$\sin(kx)=\frac{e^{ikx}-e^{-ikx}}{2i}$, So, when I do the projection I see that I will have: $\hat f_k=\sum_{j=0}^{N-1}e^{-ik}\frac{e^{ikx}-e^{-ikx}}{2i}$
so the k in the exponential is twice the k in the sine function?

 

[vela]

In the original post, you were looking at a function $f(x) = \sin kx$ with the specific case where $k=2$. In your expression for $\hat{f}_k$, you said $k=0,1,2,\dots,N-1$. And then in your expression for $f_j$, $k$ is the summation index. Using the same symbol to mean different things is confusing at best. More likely, it's going to lead to a mistake.

 

[Telemachus]

In the original post, you were looking at a function $f(x) = \sin kx$ with the specific case where $k=2$. In your expression for $\hat{f}_k$, you said $k=0,1,2,\dots,N-1$. And then in your expression for $f_j$, $k$ is the summation index. Using the same symbol to mean different things is confusing at best. More likely, it's going to lead to a mistake.

Yes, you are right, I was actually seeking for the wrong k in the transformed function. I didn't thought of it, thank you. Should I create a new thread for the question in post #6? I didn't do it because it is about discrete Fourier transforms too.

 

[vela]

No, I don't think that's necessary. Post 6 provides some context for what you were asking originally.

For what it's worth, I tried the calculation you did (in post 6) in Mathematica and got the results you were expecting. Because of rounding errors, the values of $u(x_j)$ had tiny imaginary parts, but ignoring those and plotting just the real part appears to have reproduced the analytical solution. I used values of $k$ running from $-N/2$ to $N/2-1$.

It's probably a good idea to make sure your transformations are working for simpler cases, as you seem to be trying now, to pinpoint where your problem lies.

 

[Telemachus]

Yes, that gives me zero. I'll post the whole code in case that anyone wants to try. I'm using a unix system, and compile with gfortran.

To compile I use:
gfortran program.for -L/usr/local/lib -lfftw3 -lm -o program.x

You have to install fftw first: http://www.fftw.org/download.html

    program fourser
c...FFTW transformation
       double complex zin, zout
       parameter (nmax=2**10)
       dimension zin(nmax), zout(nmax)
       integer*8 plan
   double complex zero,zone,eye,kv
   real*8 rzero,one
   real*8 rmin,rmax,dx
   integer N,kn
   real*8 x,ct,kr,user,u
   data rzero,one,two/0.0d0,1.0d0,2.0d0/
   include "/usr/local/include/fftw3.f"
   zero = dcmplx(rzero,rzero)
   eye = dcmplx(rzero,one)
   zone = dcmplx(one,rzero)
       Pi = two*dasin(one)
   rmin = rzero
   rmax = two*Pi
   N = 1000
   dx = (rmax-rmin)/N
c... output files
   open(unit=15,file='fx2.dat',status='unknown')
   open(unit=16,file='user2.dat',status='unknown')
   open(unit=17,file='uan2.dat',status='unknown')
   open(unit=18,file='cuser2.dat',status='unknown')
c   open(unit=25,file='Invfk.dat',status='unknown')
   do i=1,N
     x = i*dx
     zin(i) = dexp(dsin(x))*(dsin(x)**2+dsin(x))
     rz = dreal(zin(i))
     ri = dimag(zin(i))
     ra = cdabs(zin(i))
         u = dexp(dsin(x))
     write(15,*) x,' ',rz
     write(17,*) x,' ',u
   enddo
c... Direct Fast Fourier Transform
   call dfftw_plan_dft_1d(plan,N,zin,zout,
     +                        FFTW_FORWARD,FFTW_ESTIMATE)
   call dfftw_execute(plan)
   call dfftw_destroy_plan(plan)
c... Solve problem
      do j=1,N
       x = j*dx
       kv=zero
       do i=1,N
        kn=-N/2+1-i
        ct=one/(N*(one+kn**2))
        kv=ct*zexp(eye*kn*x)*zout(i)+kv
       enddo
     user = dreal(kv)
     cuser = dimag(kv)
    write(16,*) x,' ',user
    write(18,*) x,' ',cuser
      enddo
   close(15)
   close(16)
   close(17)
   close(18)
   stop
   end

I think that possible the problem is in how the coefficients are ordered, but I think I have tried all the possibilities.

 

[vela]

Rather than doing the summation for the inverse DFT yourself, why not let FFTW do it? You just need one loop to calculate $\hat{u}_k$ and the send the resulting array to FFTW.

 

[Telemachus]

Oh, I'll try that, its a good idea.

Yes, that was a good idea. Now I have some constant messing around, but the convergence seems good.

Sorry, after a check, I still having the complex component for the solution.

 

[vela]

When you calculate $\hat{u}_k$, you're probably using the wrong value for $k$ in $1/(1+k^2)$.

 

[Telemachus]

Yes, that was it. Thanks!

 

[Telemachus]

Sorry, I saw the stupid thing I did with the k the other day and I laughed. Anyway, of course in the program was well set.