# Fft() function does not return the correct analytic result

**URL:** <https://discourse.julialang.org/t/fft-function-does-not-return-the-correct-analytic-result/50644>\
**Category:** Signal and Image Processing\
**Tags:** fftw, dsp\
**Created:** [November 23, 2020, 6:40pm UTC](https://discourse.julialang.org/t/fft-function-does-not-return-the-correct-analytic-result/50644 "2020-11-23T18:40:00Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [November 23, 2020, 6:40pm UTC](https://discourse.julialang.org/t/fft-function-does-not-return-the-correct-analytic-result/50644/1 "2020-11-23T18:40:01Z")

</div>

Hello, currently I’m having issues understanding how to perform with success a Fourier Transformation using the `fft()` function of the `fftw.jl` package.

Let’s pick the step function:

```julia
using FFTW, Plots
x = range(-5,stop=5,length=1000)
f = [abs(i) <= 1 ? 1 : 0 for i in x]

```

When applying a Fourier Transformation to this function we have the following analytic function:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/d/4d8bf71530f887b0d3a6d0e1ded84d4759363206.png) ,  
with `a=2` being the width of the rectangle.

```julia
Fanalyt = (1/pi)*sinc.(x)
plot(x,Fanalyt)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/2/223c468310d30373742a3e6f05b455dccafd6e11.png)

If I, naively, perform a fourier transformation to `f`, I’ll get this wrong solution:

```julia
Fnumeric = fftshift(fft(f))
plot(x,real(Fnumeric))

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/3/a3e92e62259886bb39665a3b6f1f1c5ee462bee1.png)

Neither the scale is correct, neither the ‘form’ of the solution.  
It’ s easy to properly scale the output of `fft()`, if I multiply `Fnumeric` by a factor of `dx/2pi`, which is the factor that appears in the first order approximation of the F.T., this scales the things down correctly. But i’m still having in incorrect behaviour, as you can see:

```julia
dx = x[2]-x[1]
Fnum_scaled = Fnumeric*(dx/(2pi))
maximum(Fanalyt)
maximum(Fnum_scaled)
plot(x,Fnum_scaled)

```

![image](https://global.discourse-cdn.com/julialang/original/3X/e/a/eafceef9ac583d8e1e7f1e68ff5da49e1fb31791.png)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/b/bb7f5b6abab7b2c78251006e37858451d689d77f.png)

Does anyone knows why this happens?

Best regards,  
Tiago

---

<div class="post-metadata">

**Author:** ![fedoroff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedoroff/32/53209_2.png) [@fedoroff](https://discourse.julialang.org/u/fedoroff)\
**Post date:** [November 23, 2020, 8:43pm UTC](https://discourse.julialang.org/t/fft-function-does-not-return-the-correct-analytic-result/50644/2 "2020-11-23T20:43:25Z")

</div>

```julia
import FFTW
using PyPlot

Nx = 1024
x = range(-10, 10, length=Nx)
dx = x[2] - x[1]
k = 2 * pi * FFTW.fftfreq(Nx, 1/dx)
k = FFTW.ifftshift(k)

a = 1.0
f = [abs(i) <= a / 2 ? 1 : 0 for i in x]

sa = @. a * sinc(k * a / 2 / pi) # sinc(x) = sin(pi * x) / (pi * x)

s = FFTW.fftshift(f) # compensate osciallations
s = FFTW.fft(s) * dx
s = FFTW.ifftshift(s)

plot(k, sa)
plot(k, real.(s))
show()

```

1. Note that sinc function in Julia is defined as normalized sinc function: `sinc(x) = sin(pi * x) / (pi * x)`.
2. Spectrum has extra oscillations due to the way how fft algorithm stores arrays. You should apply extra `fftshift` before calling `fft` in order to avoid these oscillations (shift in space = frequency modulation in spectrum).
3. By default FFTW does not normalize the integrals by 2pi.

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [November 25, 2020, 10:30am UTC](https://discourse.julialang.org/t/fft-function-does-not-return-the-correct-analytic-result/50644/3 "2020-11-25T10:30:09Z")

</div>

> [@fedoroff](#):
>
> ```julia
> plot(k, sa)
> plot(k, real.(s))
> 
> ```

Thank you for the explanation!  
I did not know that I had to write a variable this variable: `k = 2 * pi * FFTW.fftfreq(Nx, 1/dx)`, in order to plot the output of `fft()`.

But i have another question, why the imaginary part of `s` is so high? Theoretically it should be 0, so I was expecting to have values close to 0. But it’s not what I’m getting.  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/c/bc85ff4ae230623eaeaee3c4ce311e404ece5098.png)  
A value far from machine precision. Is this a matter of some kind of normalization that I’m not taking into account?

Cheers,  
Tiago

---

<div class="post-metadata">

**Author:** ![fedoroff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedoroff/32/53209_2.png) [@fedoroff](https://discourse.julialang.org/u/fedoroff)\
**Post date:** [November 25, 2020, 10:42am UTC](https://discourse.julialang.org/t/fft-function-does-not-return-the-correct-analytic-result/50644/4 "2020-11-25T10:42:03Z")

</div>

The imaginary part is nonzero because the sinc function decays very slow and touches the boundaries of the spectral grid (restricted by the Nyquist frequency). In order to improve the situation you need to increase the size of the spectral grid by decreasing the step size in the real domain (try to increase Nx twice).

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [November 25, 2020, 1:21pm UTC](https://discourse.julialang.org/t/fft-function-does-not-return-the-correct-analytic-result/50644/5 "2020-11-25T13:21:49Z")

</div>

Ohh i got it! thank you very much, this helped me a lot!!

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [November 26, 2020, 7:21pm UTC](https://discourse.julialang.org/t/fft-function-does-not-return-the-correct-analytic-result/50644/6 "2020-11-26T19:21:58Z")

</div>

One small remark. The even-symmetry of the input signal to the DFT could be further improved by defining `x` more carefully:

```julia
Nx = 1024
x = range(-10,10,length = Nx + 1)
x = x[1:Nx]

```

Then the code provided in the solution by @fedoroff, will show the imaginary part of the FFT close to zero and thus an even better fit to the analytic solution will result after taking the real part.  
PS: the fact that the FT of an even real function is real is exploited.