如何解释 cuFFT R2C 结果

Question

我正在使用 GPU 加速一些数据分析代码，目前正在 numpy.fft 库和 cuFFT（使用 skcuda.fft 包装器）之间进行一些分析和比较。

我确定我只是遗漏了一些关于 cuFFT 中 FFT 实现的明显信息，但我正在努力寻找 cuFFT 文档中的内容。

为了解决这个问题，我创建了 500 毫秒的数据，采样率为 100 MS/s，其中包含一些频谱分量。然后，我声明 GPU 阵列、袖带计划 (R2C) 和运行带有数据子集的 fft。最后，我与 numpy.fft.rfft:

进行了比较

time = np.linspace(0, 500E-3, int(50E6))
freq = 100E6

data = np.sin(2*np.pi*time*1E6)
data += np.sin(2*np.pi*time*2E6 + 0.5)
data += np.sin(2*np.pi*time*3E6 - 0.1)
data += np.sin(2*np.pi*time*4E6 - 0.9)
data += np.sin(2*np.pi*time*1E6 - 1.9)
data += np.sin(2*np.pi*time*15E6 - 2.1)
data += np.sin(2*np.pi*time*20E6 - 0.3)
data += np.sin(2*np.pi*time*25E6 - 0.3)

nPtsFFT = int(2**13)

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

rDev = gp.GPUArray(int(nPtsFFT/2+1), np.float64)

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

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()

我天真地假设它们大致相等，但我发现袖带峰都发生了偏移，其他所有点都低于预期。

这让我想起了 scipy.fftpack.rfft 的输出，所以检查那里的文档我发现了 Re 和 Im 部分的交错。所以，如果我将绘图修改为：

plt.loglog(freqs, np.abs(hfftRes), label='npfft')
plt.loglog(freqs[:-1:2]/2, np.abs(rHost[:-1:2] + 1j*rHost[1::2]), label='cufft')
plt.legend()
plt.show()

我现在得到了我所期望的，但最高只有 25 MHz，而在给定采样率的情况下，我应该能够获得最高 50 MHz 的数据。

有没有办法从此变换中提取高达奈奎斯特频率的数据？

Answer 1

我还没有找到 cufft 输出的解释，但我可以通过 cupyx.scipy.fft.rfft 得到我想要的行为，如果其他人发现同样的问题，这可能会有用。

deviceData = cp.array(data, dtype = np.float64)
fftPlan = cpfftPack.get_fft_plan(deviceData[:nPtsFFT], (nPtsFFT,),
                                 value_type = 'R2C')
dfft = cpfft.rfft(deviceData[:nPtsFFT], plan = fftPlan)
h_dfft = dfft.get()
hfft = npfft.rfft(data[:nPtsFFT])

plt.loglog(freqs, np.abs(hfft), label='npfft')
plt.loglog(freqs, np.abs(h_dfft), label='cupyx fft')
plt.legend()
plt.show()

Answer 2

由于 R2C 接口产生复杂的输出，您必须提供一个 np.complex128 类型的数组来获取整个 int(nPtsFFT/2+1) 复杂值，而不仅仅是 int(nPtsFFT/2+1) 浮点值（这将只对应了一半的数据量）。

这可以通过如下更改 rDev 定义（并保持其他所有内容相同）来完成：

rDev = gp.GPUArray(int(nPtsFFT/2+1), np.complex128)

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

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()

结果应该会如预期的那样一直上升到 50MHz 奈奎斯特频率，尖峰与参考 np.fft.rfft 实现很好地对齐。

如何解释 cuFFT R2C 结果

How to interpret cuFFT R2C result

fft

cufft

scikit-cuda