如何解释 cuFFT R2C 结果答案

【问题标题】：How to interpret cuFFT R2C result如何解释 cuFFT R2C 结果
【发布时间】：2022-01-11 22:07:42
【问题描述】：

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

我确定我只是遗漏了一些关于 cuFFT 中的 FFT 实现的明显内容，但我很难在 cuFFT 文档中找到它。

为了解决这个问题，我创建了 500 毫秒的数据，以 100 毫秒/秒的速度采样，并带有一些光谱分量。然后，我声明 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 的数据。

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

【问题讨论】：

您是否尝试过 rDev = gp.GPUArray(int(nPtsFFT/2+1), np.complex128 在您的原始（非交错）方法中？
@SleuthEye 天哪，我是个白痴，解决了它。如果你这样回答，我会接受的。

标签： fft cufft scikit-cuda

【解决方案1】：

由于 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 实现很好地对齐。

【讨论】：

【解决方案2】：

我还没有找到对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()

【讨论】：