What are good/best practices to take the Fourier transform of an InterpolatingFunction?

Question

I have a function which I have obtained from numerical integration of a differential equation, and I would like to take its Fourier transform. What are good practices for doing this?

To make things more precise, suppose I have a function like the following:

solution =  u /. First@NDSolve[{
         u'[t] == -E^(-(t^2/2)) t Cos[5 t] - 5 E^(-(t^2/2)) Sin[5 t] t,
         u[0] == 1
      }, u, {t, 0, 50}]

(Note that, while this example has an analytical solution, my full problem does not.) I am using NDSolve since it appears to be the best way to obtain functions $F$ of the form $F(t)=\int_{t_0}^t f(\tau) \,\text d\tau$ where $f$ is known. I am interested in the Fourier transform of this solution.

I know of course that I can sample this function over a discrete set of points and then use the discrete Fourier transform Fourier to obtain a fair approximation:

ListPlot[
 Abs[
  Fourier[solution /@ Range[0, 50, 0.01]]
  ]
 , PlotRange -> {{0, 150}, {0, 1}}]

However, this completely ignores the fact that the InterpolatingFunction itself is already a discrete representation of the desired function. As it stands, the code above includes interpolation from a discrete data set into a new discrete data set, and this can only

• lose information, in the case that the sample grid is too sparse, or alternatively
• give the impression that it contains more information than was already in the original data, if the sample grid is too closely spaced.

Either way, this extra sampling step has no communication with the internals of the InterpolatingFunction and as such it cannot know how much information it contains or how good it is. I am looking for a method which uses this information as optimally as possible. The accuracy of the ODE solution, and therefore of its transform, are best controlled from the options of NDSolve, and the subsequent transformations of the data should take heed of this.

I am aware of NFourierTransform, which, as I understand, uses NIntegrate internally; I'm also aware that NIntegrate natively supports InterpolatingFunction objects, which is presumably still the case if they are multiplied by a further factor of $e^{i\omega t}$. (On the other hand, the behaviour of this factor can rather radically change the overall behaviour of the whole integrand from 'slow' into highly oscillatory, so that some attention to Method would be warranted.)

It could be, though, that the Fast nature of Fourier might be enough to outstrip such advantages, and still be faster than NFourierTransform even when over-blasted on the accuracy. Speed is definitely an issue, as I will be performing automated scans over many such spectra, so I want them to be accurate enough but not waste time overcalculating.

So: what are good ways of performing this transform?