使用 Python 的傅里叶变换不准确

Question

我的目标是对分布进行傅立叶变换。这是一个物理问题，我正在尝试将函数从位置 space 转换为动量 space。然而，我发现当我尝试使用 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

Answer 1

问题是，无论N多大，你要解决的细节分辨率都是有限的。您需要扩展原始 x 的限制，使用插值重新采样不会在那里做任何事情。这是一个示例运行：我创建了一个与您相似的数据集。查看如果在离开 x.

的限制时将 loc 设置为 2、50、80 会发生什么

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)

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

使用 Python 的傅里叶变换不准确

Inaccurate Fourier Transform using Python

python

scipy

transform

fft

continuous-fourier