如何用Python/Numpy从坐标与数值创建大型3D数组并计算傅里叶变换

核心结论：直接构建(6308,6308,6308)数组完全不可行 #

别浪费时间尝试这个操作——就算用单精度float32存储，这个数组的内存需求约为1TB，普通设备根本无法承载。你的6308个有效点只是全网格的极小部分，必须换思路处理。

可行方案1：稀疏数组+稀疏傅里叶变换 #

利用数据的稀疏性（仅6308个有效点），用坐标形式存储数据，避免全网格填充：

• 先将原始(x,y,z)物理坐标映射为连续整数网格索引（如果坐标本身是离散整数则直接使用），配合对应数值保存为坐标-值对。
• 使用专门的稀疏FFT工具处理，无需构建全尺寸网格，比如基于scipy.sparse的扩展工具，或直接用非均匀FFT库。

可行方案2：重采样到小网格（最易实现） #

如果可以接受降低分辨率，这是最快落地的方法：

1. 统计x、y、z坐标的最大/最小值，选择设备可承载的网格分辨率（比如256或512，256³的float32数组仅占64MB）。
2. 通过插值将离散点填充到小网格，再执行傅里叶变换：

import numpy as np
from scipy.interpolate import griddata

# 假设data是加载后的N行4列数组，格式为[x,y,z,value]
points = data[:, :3]
values = data[:, 3]

# 定义目标网格分辨率
res = 256
x_grid = np.linspace(points[:,0].min(), points[:,0].max(), res)
y_grid = np.linspace(points[:,1].min(), points[:,1].max(), res)
z_grid = np.linspace(points[:,2].min(), points[:,2].max(), res)
X, Y, Z = np.meshgrid(x_grid, y_grid, z_grid, indexing='ij')

# 插值填充网格，可选linear/nearest/cubic方法
grid_values = griddata(points, values, (X, Y, Z), method='linear')
# 执行3D傅里叶变换
fft_result = np.fft.fftn(grid_values)

可行方案3：非均匀FFT（NUFFT）处理离散点 #

如果不能降低分辨率，直接用NUFFT处理非均匀采样点，无需构建网格：

• 使用finufft这类轻量库，直接输入离散坐标和数值即可输出傅里叶变换结果：

import finufft

# 将坐标归一化到NUFFT要求的[-π, π]范围
x_range = points[:,0].max() - points[:,0].min()
y_range = points[:,1].max() - points[:,1].min()
z_range = points[:,2].max() - points[:,2].min()

x_norm = (points[:,0] - points[:,0].min())/x_range * 2*np.pi - np.pi
y_norm = (points[:,1] - points[:,1].min())/y_range * 2*np.pi - np.pi
z_norm = (points[:,2] - points[:,2].min())/z_range * 2*np.pi - np.pi

# 设置傅里叶空间的分辨率，示例为1024×1024×1024
n = 1024
fft_result = finufft.nufft3d1(x_norm, y_norm, z_norm, values, n, n, n)

关于xarray的补充用法 #

xarray确实依赖numpy数组作为基础，但你可以将上述生成的插值网格数组或稀疏数组传入xarray.DataArray，示例：

import xarray as xr

da = xr.DataArray(grid_values, coords=[x_grid, y_grid, z_grid], dims=['x', 'y', 'z'])
fft_da = xr.DataArray(np.fft.fftn(da), coords=da.coords, dims=da.dims)

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