Inaccurate Fourier Transform using Python

Question

I am aiming to take the fourier transform of a distribution. It is a physics problem and I am trying to transform the function from position space to momentum space. I am however finding that when I attempt to take the fourier transform using scipys fft, that it becomes jagged whereas a smooth shape is expected. I assume it is something to do with sampling, but I cannot work out what is wrong.

This is what the transformed function currently looks like:

This is what it is roughly supposed to look like (it may have a slightly different width, but in terms of smoothness it should look similar):

and here is the code used to generate the blue image:

from scipy.fft import fft, fftfreq, fftshift
import numpy as np
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
import scipy
from scipy import interpolate
from scipy import integrate
# number of signal points
x = np.load('xvalues.npy') #Previously generated x values
y=np.load('function_to_be_transformed.npy') #Previously generated function (with same number of values as x)
y = np.asarray(y).squeeze() 
f = interpolate.interp1d(x, y) #interpolating data to make accessible function

N = 80000
# sample spacing
T = 1.0 / 80000.0
x = np.linspace(-N*T, N*T, N)
y=f(x)
yf = fft(y)
xf = fftfreq(N, T)
xf = fftshift(xf)
yplot = fftshift(yf)
import matplotlib.pyplot as plt
plt.plot(x,np.abs(f(x))**2)
plt.xlabel('x')
plt.ylabel(r'$|\Psi(x)|^2$')
plt.savefig("firstPo.eps", format="eps")
plt.show()

plt.plot(xf, np.abs(1.0/N * np.abs(yplot))**2)
plt.xlim(right=100.0)  # adjust the right leaving left unchanged
plt.xlim(left=-100.0)  # adjust the left leaving right unchanged
#plt.grid()
plt.ylabel(r'$|\phi(p)|^2$')
plt.xlabel('p')
plt.savefig("firstMo.eps", format="eps")
plt.show()

Update

If anyone could offer some further advice, that'd be great because I am still having trouble. Following from @ScottStensland 's comment, I have attempted to find the FT of a sin wave to see if I find any problems and then retrofit the example back onto my initial problem.

Here are the results for the FT of sin(x):

This is as expected (I think). But when I retrofit the code back to by initial example I get the following (The top image is my initial distribution):

The code is as follows for the sin(x) example:

# sin wave
import numpy as np
from numpy import arange
from numpy.fft import rfft
from math import sin,pi
import matplotlib.pyplot as plt
def f(x):
    return sin(x)

N=1000
x=np.arange(0.0,1.0,1.0/N)
y=np.zeros(len(x))
for i in range(len(x)):
    y[i]=f(x[i])
#y=map(f,x)
#print(y)
c=rfft(y)       
plt.plot(abs(c))    
plt.xlim(0,100)     
plt.show()

and for the attempt at my own one:

#Interpolated Function
# sin wave
import numpy as np
from numpy import arange
from numpy.fft import rfft
from math import sin,pi
import matplotlib.pyplot as plt
x = np.linspace(-1.0,1.0,1001) #Previously generated x values
y=np.load('function_to_be_transformed.npy') #Previously generated function (with same number of values as x)
y = np.asarray(y).squeeze()
N=1001
x=np.arange(-1.0,1.0,2.0/N)

#y=map(f,x)
#print(y)
plt.plot(x,y)
plt.show()
c=rfft(y)       
plt.plot(abs(c))    
plt.show()

The relevant files are here: https://github.com/georgedixon4321/NewDistribution.git

when coming up to speed on Fourier Transforms its often very helpful to begin with a simple sin curve as the time domain input ... feed that into a fft call and plot this as the frequency domain data ... confirm you are seeing the known frequency which matches your known input ... then feed your freq domain data into an inverse Fourier Transform to synthesize a time domain signal which should match or be very close to your original sin curve input ... having handy this simple known input is very helpful to get the kinks out of the code and gain an appreciation as to the many options avail — Scott Stensland, Aug 06 '20 at 22:49
Try to plot your frequency spectrum without interpolation first. Also plot the interpolated signal, make sure it looks the way you think it does. I have no idea what your input x and y values look like, are you sure your input x goes from -1 to 1? If not, you might be cutting off part of your data, or even extrapolating. This whole interpolating data business smells, IMO. You shouldn't do that. — Cris Luengo, Aug 06 '20 at 23:07
@CrisLuengo Thanks for your response, I have checked and the interpolated plot is definitely identical to the array that I originally had plotted. The input definitely goes from 1 to -1, I only interpolated because I did not know how to transform the array directly. — George, Aug 06 '20 at 23:57
@ScottStensland Thanks for your response, I will have a play around with the sin wave and see if I can find a way to fix it — George, Aug 06 '20 at 23:58
@ScottStensland I have followed your advice and seem to have gotten a similar, jagged result (though it returns a smooth curve for the transform of the sin wave). I have updated the original post. — George, Aug 07 '20 at 12:54

Péter Leéh · Accepted Answer · 2020-08-08T09:19:29.963

The problem is that the resolution of the details you want to resolve is limited, no matter how big N is. You need to extend the limits of the original x, resampling with interpolation is not doing anything there. Here is a sample run: I created a similar dataset you have. Check out what happens if you set loc to 2, 50, 80 when leaving the limits of x.

from scipy.fftpack import fft, fftshift, fftfreq, ifft

loc = 2
x = np.linspace(-130, 130, 10000)

y1 = np.exp(-((x - loc) ** 2) / (2 ** 2))
y2 = np.exp(-((x + loc) ** 2) / (2 ** 2))
y = y1 + y2

plt.figure()
plt.plot(x, y)

xf = fftshift(fftfreq(len(x), np.diff(x)[0]))
yf = ifft(y)

plt.figure()
plt.plot(fftshift(xf), np.abs(yf))
plt.xlim(-0.5, .5)

As the spikes get further and further away from each other, you need to extend the limits of the domain to achieve the same resolution.

Applying this to your example:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.fftpack import fft, ifft, fftfreq, fftshift

x = np.load('xvalues.npy')
y = np.load('function_to_be_transformed.npy').ravel()

f = interp1d(x, y, fill_value="extrapolate")

N = 1000000

# I made a bigger domain
x = np.linspace(10*x[0], 10*x[-1], N)

y = f(x)


xf = fftshift(fftfreq(len(x), np.diff(x)[0]))
yf = ifft(y)

plt.figure()
plt.plot(fftshift(xf), np.abs(yf))
plt.xlim(-30, 30)

Note that extrapolating is dangerous, it just happened to work in this example. Before doing this, you always want to make sure the extrapolation does return the curve you want and it does not mess up anything.

Inaccurate Fourier Transform using Python

1 Answers1