python - 使用 Python 进行不准确的傅里叶变换

Question

我的目标是对分布进行傅立叶变换。这是一个物理问题，我正在尝试将函数从位置空间转换为动量空间。然而，我发现当我尝试使用 scipys fft 进行傅立叶变换时，它会变成锯齿状，而预期会是平滑的形状。我认为这与采样有关，但我无法找出问题所在。

这是转换后的函数当前的样子：

这是它大致应该看起来的样子（它的宽度可能略有不同，但就平滑度而言，它应该看起来相似）：

这是用于生成蓝色图像的代码：

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

---

更新

---

如果有人可以提供一些进一步的建议，那就太好了，因为我仍然遇到麻烦。根据@ScottStensland 的评论，我试图找到正弦波的 FT 以查看是否发现任何问题，然后将示例改回为我最初的问题。

以下是 sin(x) 的 FT 的结果：

这是预期的（我认为）。但是当我通过初始示例将代码改回时，我得到以下信息（顶部图像是我的初始分布）：

sin(x) 示例的代码如下：

# 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()

对于我自己的尝试：

#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()

相关文件在这里： https ://github.com/georgedixon4321/NewDistribution.git

标签： pythonscipytransformfftcontinuous-fourier

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)

请注意，外推是危险的，它恰好在这个例子中起作用。在执行此操作之前，您始终要确保外推确实返回了您想要的曲线并且它不会弄乱任何东西。