如何解读cuFFT R2C变换结果？为何其与numpy.fft.rfft输出存在差异

问题核心原因 #

你遇到的输出截断、需要手动拆分实虚部的问题，本质是cuFFT R2C变换的输出数组声明错误：

• cuFFT实序列到复数序列（R2C）变换的输出为复数类型，skcuda封装要求输出数组声明为对应复数精度（float32对应complex64，float64对应complex128）
• 你之前将输出数组声明为np.float64类型，且长度仅为nPtsFFT/2+1，仅能容纳前一半复数的实部、虚部交错值，因此最终仅能覆盖到半频段（25MHz）

解决方法 #

方案1：直接声明正确的复数输出数组（推荐） #

修改输出GPU数组的类型和长度声明，即可得到和numpy.fft.rfft完全格式对齐的结果，无需手动处理实虚部：

# 测试数据构造、输入GPU数组声明部分不变
nPtsFFT = int(2**13)

dDev = gp.GPUArray(nPtsFFT, np.float64)
dDev.set(data[:nPtsFFT])

# 修正：输出数组为complex128类型，长度和numpy.rfft输出一致
rDev = gp.GPUArray(nPtsFFT // 2 + 1, np.complex128)

plan = cufft.Plan(nPtsFFT, np.float64, np.complex128)
cufft.fft(dDev, rDev, plan)
rHost = rDev.get()

# 直接和numpy结果对比，频率轴覆盖到全奈奎斯特频率（50MHz）
freqs = np.fft.rfftfreq(nPtsFFT, 1/freq)
hfftRes = np.fft.rfft(data[:nPtsFFT])

plt.loglog(freqs, np.abs(hfftRes), label='npfft')
plt.loglog(freqs, np.abs(rHost), label='cufft')
plt.legend()
plt.show()

方案2：从float类型交错输出恢复全频段数据 #

如果有特殊场景必须使用float数组接收cuFFT输出，按如下步骤调整即可得到全频段结果：

1. 输出float数组长度设置为2*(nPtsFFT//2 +1)，容纳所有复数的实部、虚部交错值
2. 恢复复数时无需缩放频率轴，直接拼接实虚部即可
   示例代码：

# 输出数组长度改为2倍，类型保持float64
rDev = gp.GPUArray(2 * (nPtsFFT // 2 + 1), np.float64)
plan = cufft.Plan(nPtsFFT, np.float64, np.complex128)
cufft.fft(dDev, rDev, plan)
rHost = rDev.get()

# 拼接实虚部得到复数结果，频率轴直接使用原始生成值
cufft_res = rHost[::2] + 1j * rHost[1::2]
freqs = np.fft.rfftfreq(nPtsFFT, 1/freq)

plt.loglog(freqs, np.abs(hfftRes), label='npfft')
plt.loglog(freqs, np.abs(cufft_res), label='cufft')
plt.legend()
plt.show()

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