Problems using FFTW- SOLVED

General Usage
1

Hi!

I am using FFTW in Julia 1.1.

I am trying to do the Fourier transform of exp(-0.5*x^2), whose analytical solution is exactly itself but I obtain a complex array.

Does anyone know what I am doing wrong? You can find the code I use just below:

Thanks for your time!!ft_fftw

using FFTW
using Plots

gridpoints=256
function_of_t=zeros(gridpoints)

x=range(-5,stop=5,length=gridpoints)

for i=1:gridpoints
	function_of_t[i]=exp(-0.5*x[i]*x[i])
end

fourier_trans=fftshift(fft(function_of_t))

y=function_of_t

plot(x,y, title="Function of t", xlabel="Time", ylabel="f(t)")
savefig("analytic_function.png")

y=fourier_trans

plot(x,y, title="FFT by means of FFTW", xlabel="omega", ylabel="F(omega)")
savefig("ft_fftw.png")

SOLVED. EDIT.

Hi! Thanks so much for your help. I really appreciate it.

Let’s do the Fourier Transform of exp(-0.5x^2)/(sqrt(2pi)) whose analytical solution is exp(-2pipi*x^2). The array scaled is the analytical solution.

using FFTW
using Plots


gridpoints=256

x=range(-5,stop=5,length=gridpoints+1)[1:end-1]
function_of_t=exp.(-0.5*x.^2)/(sqrt(2*pi))

function_of_t=fftshift(function_of_t)

fourier_trans=fftshift(fft(function_of_t))

inverse_fft=fftshift(ifft(fourier_trans))

y=fourier_trans

plot(x,y, title="FFT by means of FFTW", xlabel="omega", ylabel="F(omega)")
savefig("ft_fftw.png")


plot(x,abs.(y), title="FFT by means of FFTW", xlabel="omega", ylabel="F(omega)")
savefig("ft_fftw_abs.png")

z=inverse_fft

plot(x,abs.(z), title="Inverse FFT by means of FFTW", xlabel="omega", ylabel="F(omega)")
savefig("ift_fftw_abs.png")


scaled=zeros(gridpoints)

for i=1:gridpoints
	scaled[i]=real(fourier_trans[i])*(x[length(x)]-x[length(1)])/gridpoints
end

y=scaled


plot(x,y, title="FFT by means of FFTW", xlabel="omega", ylabel="F(omega)")
savefig("ft_scaled_fftw_abs.png")


fft_of_analytics=exp.(-2*pi*pi*x.^2)
plot(x,y, title="FFT analytical", xlabel="omega", ylabel="F(omega)")
savefig("ft_scaled_fftw_analytical.png")
2

Shouldn’t the array you are passing to fft be shifted around 0 with fftshift? Also your analytical solution works for infinity, so there could be some effect from those missing tails. I’m only guessing, cannot run it, please quote code with “```”

3

Hi!

Thanks for your answer. When shift the array, I still obtain a complex array. I modify the post. Not sure if the code is more visible now. It is my first post :sweat_smile:

4

Welcome!

Your function is not symmetric around the origin. Try

x=range(-5,stop=5,length=gridpoints+1)[1:end-1]
function_of_t=exp.(-0.5*x.^2)
fft(function_of_t)

As for the code, you should quote it within three backticks, not a single one.

5

I don’t see a problem with:

using FFTW
using Plots

x = range(-5, stop=5, length=256)
function_of_t = @. exp(-0.5x^2)

fourier_trans = fftshift(fft(function_of_t))

y = function_of_t
plot(x, y, title = "Function of t", xlabel = "Time", ylabel = "f(t)")

y = fourier_trans
plot(x, abs.(y), title = "FFT by means of FFTW", xlabel = "omega", ylabel = "F(omega)")

Which gives:

1
11920×1160 57.3 KB

and

2
21920×1160 43.8 KB

6

Hi!

Thanks so much for your reply.

My point is that the Fourier transform of exp(-0.5t²) is exactly exp(-0.5t²). When FFT is done, I do not get the same result.

7

Hi!

When I do these changes, I still get the complex array:ft_fftw

using FFTW
using Plots


gridpoints=256

x=range(-5,stop=5,length=gridpoints+1)[1:end-1]
function_of_t=exp.(-0.5*x.^2)
fourier_trans=fft(function_of_t)






y=fourier_trans

plot(x,y, title="FFT by means of FFTW", xlabel="omega", ylabel="F(omega)")
savefig("ft_fftw.png")
8

You should plot the magnitude of y, plot(x, abs.(y), .. )

9

When I do that, I obtain the following plot that does not correspond with the analytical function.ft_fftw

10

The amplitude of the real part is of order 100, that of the imaginary part 1e-14: since a Float64 is 16 digits, the imaginary part is effectively zero. fft always compute a complex transform.

1 Like
11

But you’re not calculating the FT of exp(-0.5 *t²), but rather the DFT of a sampled version of that function, defined only for a finite amount of time.

Also, the numerical innacuracy of Floats often result in a small phase being returned by FFT, even when in theory the phase should be zero.

1 Like