Неточное преобразование Фурье с использованием Python

Question

Я стремлюсь использовать преобразование Фурье распределения. Это физическая проблема, и я пытаюсь преобразовать функцию из пространства позиций в пространство импульсов. Однако я обнаружил, что когда я пытаюсь выполнить преобразование Фурье с помощью scipys fft, оно становится неровным, тогда как ожидается гладкая форма. Я предполагаю, что это как-то связано с выборкой, но я не могу понять, что не так.

Вот как сейчас выглядит преобразованная функция:

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

Answer 1

Проблема в том, что разрешение деталей, которые вы хотите разрешить, ограничено, независимо от размера N. Вам нужно расширить пределы исходного x, передискретизация с интерполяцией там ничего не делает. Вот пример выполнения: я создал аналогичный набор данных, который есть у вас. Посмотрите, что произойдет, если вы установите loc на 2, 50, 80 при выходе за пределы 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)

По мере того, как шипы удаляются все дальше и дальше друг от друга, вам нужно увеличивать ограничения домена для достижения того же разрешения.

Применяя это к вашему примеру:

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)

Обратите внимание, что экстраполяция опасна, она просто сработала в этом примере. Перед тем как это сделать, вы всегда должны быть уверены, что экстраполяция вернет нужную кривую и ничего не испортит.