# Doubt about a discrete Fourier Transform

1. Aug 25, 2017

### 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.

2. Relevant 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$

3. 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?

#### Attached Files:

• ###### fft.png
File size:
14.9 KB
Views:
22
2. Aug 25, 2017

### 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?

3. Aug 25, 2017

### scottdave

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

4. Aug 26, 2017

### Telemachus

Yes, what is wrong with that?

5. Aug 26, 2017

### vela

Staff Emeritus
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?

6. Aug 26, 2017

### 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:

Code (Fortran):

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 im 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

#### Attached Files:

File size:
14.4 KB
Views:
14
File size:
14.6 KB
Views:
15
File size:
12.9 KB
Views:
15
• ###### user2.png
File size:
17.4 KB
Views:
15
Last edited: Aug 26, 2017
7. Aug 26, 2017

### Telemachus

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?

8. Aug 26, 2017

### vela

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

9. Aug 26, 2017

### 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?

10. Aug 26, 2017

### vela

Staff Emeritus
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.

11. Aug 26, 2017

### Telemachus

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.

12. Aug 26, 2017

### vela

Staff Emeritus
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.

13. Aug 26, 2017

### 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

Code (Fortran):

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.

Last edited: Aug 26, 2017
14. Aug 26, 2017

### vela

Staff Emeritus
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.

15. Aug 26, 2017

### 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.

#### Attached Files:

File size:
13.5 KB
Views:
9
• ###### cuser2.png
File size:
11.9 KB
Views:
10
Last edited: Aug 26, 2017
16. Aug 27, 2017

### vela

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

17. Aug 27, 2017

### Telemachus

Yes, that was it. Thanks!

18. Aug 29, 2017

### 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.