用FFT计算指数复数和以模拟衍射，而非求和？
背景
我想通过在Python中编写代码来更好地理解X射线衍射。对于一组位置为R_i的点，德拜公式是这样的：
这里的指数中的i是指复数，其他的i是索引。为了简单起见，暂时设定b_i = b_j = 1。
现在我尝试对一组我有坐标的点进行显式计算：
import numpy as np
# set up grid
dims = 2
side = 30
points = np.power(side, dims)
coords = np.zeros((dims, points))
xc, yc = np.meshgrid(np.arange(side), np.arange(side))
coords[0, :] = xc.reshape((points))
coords[1, :] = yc.reshape((points))
# calculate diffraction
xdist = np.subtract.outer(coords[0], coords[0])
ydist = np.subtract.outer(coords[1], coords[1])
rdist = np.stack((xdist, ydist))
rdist = rdist.reshape(2, rdist.shape[1]*rdist.shape[2])
qs = 200
qspace = np.stack((np.linspace(-2, 8, qs), np.zeros(qs)))
diffrac = np.sum(np.exp(-1j * np.tensordot(qspace.T, rdist, axes=1)), axis=1)
经过几秒钟，我得到了以下结果：
结果看起来符合预期（周期性为2π，因为点之间的间距为1）。这也很合理，因为对于900个点，需要计算810000个距离。我没有使用循环，所以我认为代码在效率上还不错，但仅仅是手动计算这个总和似乎本身就很慢。
想法
现在看起来，如果我能使用离散快速傅里叶变换（FFT），事情会大大加快，因为这个总和的形状。然而：
- 对于离散傅里叶变换，我仍然需要将图像“像素化”（根据我的理解），以在信号中的点之间包含很多空白区域。就像我想要转换我分享的第一张图片的像素。这似乎也不太高效（例如，因为采样的问题）。
- 我希望之后能移动这些点，所以第一张图片是一个网格并且采样规律这一点并没有特别帮助。看起来非均匀傅里叶变换可能对我有帮助，但这仍然需要我“像素化”图像并将一些值设置为0。
问题
有没有办法使用FFT（或其他方法）更快地计算这个总和，从一组np.array坐标（x,y）开始？（如果你愿意，可以看作是狄拉克δ函数...）。
特别希望能指点一些相关的数学技巧/Python函数/Python库。我对实际应用中的傅里叶变换不太熟悉，但我在网上找到的大部分资料似乎都不相关。所以我可能在错误的方向上，或者我的理解有些欠缺。任何帮助都非常感谢！
（第一张图片是来自https://www.ill.eu/fileadmin/user_upload/ILL/6_Careers/1_All_our_vacancies/PhD_recruitment/Student_Seminars/2017/19-2017-05-09_Fischer_Cookies.pdf的截图，因为似乎在SO上没有数学符号，或者我没有找到）
1 个回答
这个回答提供了一种方法，可以让代码运行得更高效，充分利用你的CPU计算能力，从而大大加快速度。
在运行过程中，有超过90%的时间都花在了np.exp上，因为计算复杂数的指数是非常耗费资源的。
为了加快这个过程，可以使用多线程（因为Numpy本身不支持多线程）。此外，我们还可以使用更快的np.exp实现（通常是利用CPU的SIMD单元）。这两者都可以通过Numexpr轻松实现。
接着，我们可以通过矩阵乘法qspace.T @ rdist来加速np.tensordot操作，因为Numpy的实现效率不高。
import numexpr as ne
# Equivalent of the last line of the code:
tmp1 = qspace.T @ rdist
tmp2 = ne.evaluate('exp(-1j * tmp1)')
diffrac = np.sum(tmp2, axis=1)
性能评估
以下是我在i5-9600KF CPU（6核）上的性能测试结果：
Initial code: 9.3 s
New proposed code: 1.1 s
因此，新的实现速度是原来的8.5倍。大部分时间仍然花在计算复杂数的指数上（超过60%）。