如何用FFT替代求和快速计算德拜公式模拟X射线衍射？

我试图通过Python编码来加深对X射线衍射的理解。对于位置为$R_i$的点集，德拜公式如下：
$I(\mathbf{q}) = \sum_{i,j} b_i b_j e^{-i \mathbf{q} \cdot (\mathbf{R}_i - \mathbf{R}_j)}$，其中指数中的$i$代表虚数单位，其余$i,j$为点的索引，为简化计算暂时令bᵢ = bⱼ = 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)

运行几秒后得到符合预期的结果，但900个点需计算810000个距离，手动求和本质效率偏低。

从求和形式看，FFT可大幅提速，但存在问题：离散傅里叶变换需将图像像素化，引入大量空白区域，效率不高；后续需移动点位置，规则网格采样无帮助，非均匀傅里叶变换仍需像素化。

能否从np.array格式的(x,y)坐标列表（可视为狄拉克δ函数）出发，用FFT或其他方法更快计算该求和？希望获取相关数学技术、Python函数或包的指引。

核心数学变形 #

先对德拜公式做关键简化：当$b_i = b_j = 1$时，原$O(N^2)$的求和可拆分为：
$$
\sum_{i,j} e^{-i\mathbf{q}\cdot(\mathbf{R}_i - \mathbf{R}j)} = \left| \sum{i} e^{-i\mathbf{q}\cdot\mathbf{R}_i} \right|^2
$$
这个变形直接把双重求和转化为单个求和的模平方，而后者可通过傅里叶变换高效计算，无需依赖规则网格。

两种高效实现路径 #

路径1：基于常规FFT的优化方案 #

若接受轻度像素化，可通过自适应网格减少空白区域：

1. 构建自适应密度场：
   • 计算点集的最小包围盒，根据点的最小间距设置网格分辨率
   • 用np.histogram2d将点映射到网格，生成仅包含点集区域的密度矩阵（避免全尺寸空白网格）
2. 执行FFT并计算强度：
   对密度矩阵做2D FFT，取结果的模平方即为德拜强度，再将FFT的频率轴映射到目标$\mathbf{q}$空间即可。

路径2：非均匀傅里叶变换（NUFFT） #

完全避免像素化，直接处理非均匀点集：

• 推荐工具包：
  • finufft：快速非均匀FFT的高效实现，支持2D/3D场景，直接输入点坐标和权重（此处全为1），即可输出指定$\mathbf{q}$点的结果
  • pyfftw：可配合自定义规划加速FFT/NUFFT计算

NUFFT代码示例

import numpy as np
import finufft

# 假设coords是(2, N)格式的输入坐标数组
x = coords[0]
y = coords[1]
N = x.size

# 定义目标q空间点
qs = 200
qx = np.linspace(-2, 8, qs)
qy = np.zeros(qs)

# 调用finufft计算非均匀傅里叶变换
result = finufft.finufft2d1(x, y, np.ones(N, dtype=np.complex128), qx, qy, eps=1e-6)

# 计算德拜强度（模平方）
diffrac_intensity = np.abs(result)**2

性能对比 #

• 900个点的场景下，显式求和为$O(N²⁾⁼⁸¹⁰4$次运算；FFT/NUFFT为$O(N \log N)$或线性复杂度，速度可提升1~2个数量级
• 当点数量增至10^4以上时，性能差距会更显著

内容的提问来源于stack exchange，提问作者andwerb