Real-valued spreaders and interpolators

[mreineck]

I know that this functionality is already mentioned in the TODO file, but I wanted to extend the discussion a bit more to include the FFT part of the algorithm and see whether we can perhaps arrive at a fully "real-valued NUFFT".

It is obvious that the spreading step can be accelerated considerably if the user-provided strengths are real-valued instead of complex. The following FFT of the oversampled grid will then result in a complex array with certain symmetries that make half of its entries redundant. However, the array returned to the user is only a sub-array of this, and without special precautions this array will lose some of those properties.
For example, if we consider a 1D type 1 NUFFT with real-valued input and a requested output size that is even, we will most likely end up with an entry at the Nyquist frequency that has a non-vanishing imaginary component. And that would no longer correspond to the Fourier transform of a real-valued field, which is bad...
Is there a mathematically clean way to deal with this effect?

If so, we could perhaps go one step further and, instead of performing an FFT, perform a Fast Hartley transform (https://en.wikipedia.org/wiki/Discrete_Hartley_transform), and return a real-valued array of the requested size instead of a complex-valued one. (I admit this is a bit exotic, but has an undeniable neatness to it ...)
Overall, except for the coordinates, all arrays involved in the transform would have use half the memory of their complex-NUFFT counterparts, which is a nice additional bonus to te faster spreading.

Analogous things can be done for type 2 NUFFTs. It's probably not worth bothering with type 3, since that almost always requires multiplying the strengths with complex phases, which prohibits any performance gains.

As is already mentioned in the TODO file, all of this will lead to a lot of almost-code-duplication. So it's probably only worth thinking about if there are real-world use cases.