Cómo interpretar el resultado R2C de cuFFT

Estoy en el proceso de acelerar algún código de análisis de datos con GPU y actualmente estoy haciendo algunos perfiles y comparaciones entre la biblioteca numpy.fft y cuFFT (utilizando el envoltorio skcuda.fft).

Estoy seguro de que solo estoy pasando por alto algo obvio sobre la implementación de FFT en cuFFT, pero estoy luchando por encontrar qué es en la documentación de cuFFT.

Para configurar el problema, creo 500 ms de datos muestreados a 100 MS/s con algunos componentes espectrales. Luego, declaro los arrays de GPU, el plan cufft (R2C) y ejecuto la fft con un subconjunto de los datos. Finalmente, hago una comparación con numpy.fft.rfft:

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

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

nPtsFFT = int(2**13)

dDev = gp.GPUArray(nPtsFFT, np.float64)
dDev.set(datos[: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/frecuencia)
hfftRes = np.fft.rfft(datos[:nPtsFFT])

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

Naive comparison of cuFFT and npfft

Esto me recordó a la salida de scipy.fftpack.rfft, así que al verificar la documentación allí, encontré la intercalación de las partes Re e Im. Por lo tanto, si modifico la trama a:

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

Comparison compensating for interleave Re and Im output

Ahora obtengo lo que espero, pero solo hasta 25 MHz, cuando debería poder obtener datos de hasta 50 MHz dada la frecuencia de muestreo.

¿Hay alguna manera de extraer los datos hasta la frecuencia de Nyquist de esta transformación?