Fft() function does not return the correct analytic result

Specific Domains Signal and Image Processing
,
1

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:

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
image963×84 11.4 KB
,
with a=2 being the width of the rectangle.

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

image
image730×491 16 KB

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

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

image
image725×484 11.1 KB

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:

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

image

image
image724×486 12.3 KB

Does anyone knows why this happens?

Best regards,
Tiago

2 
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.
2 Likes
3

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
A value far from machine precision. Is this a matter of some kind of normalization that I’m not taking into account?

Cheers,
Tiago

4

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

1 Like
5

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

6

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

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.

3 Likes