Statement of the problem:
The convolution of two functions is defined by: \begin{eqnarray} (f*g)(y)=\int_{-\infty}^{+\infty}dx f(x)g(y-x), \end{eqnarray} for functions $f(x),g(x)$ defined over the complex plane $\mathbb{C}$; but note the integral is performed over the real axis only. Now, we assume that our functions are just defined over a specific grid of points $\{x_{j}\}, j=0,...,N-1$, and we approximate the integral by its Riemman sum, \begin{eqnarray} (f*g)[k]\approx \sum_{j=1}^{N-1}\Delta x_{j}f[x_{j}]g[y_{k}-x_{j}]. \end{eqnarray} Moreover, assume that $\Delta x_{j}$ are arbitrary, in the sense that the points are not equidistantly spaced. In addition, assume that the points $\{x_{j}\}$ define a vast variation of the scales of the real axis, meaning that generally $|\Delta x_{0}|>> |\Delta x_{N/2}|$, i.e. the differences between grid points tend to decrease towards the middle point of the grid.
If all points were equidistantly spaced, then one can use the DFT to perform a FastFourier Transform to calculate the convolution integral. In our case, since points are not equidistant, we use the NUFFT: \begin{eqnarray} f[x_{j}]=\frac{1}{\sqrt{N}}\sum_{k}e^{-i2\pi \omega_{k}x_{j}}f[\omega_{k}]. \end{eqnarray} Then we may write: \begin{eqnarray} (f*g)(y_{k})\approx \frac{1}{N}\sum_{k_{1},k_{2}}f[\omega_{k_{1}}]g[\omega_{k_{2}}]e^{-i2\pi y_{k}\omega_{k_{2}}}\left(\sum_{j}\Delta x_{j} e^{-i2\pi (\omega_{k_{1}}-\omega_{k_{2}})x_{j}}\right). \end{eqnarray}
Question: The problem comes on finding a transformation of the original grid that allows to compute the integral in a way similar to the FFT, whose complexity is $\mathcal{O}(N\log N)$. It is clear that such lower bound cannot be reached in this case, but I am looking for a way to perform this efficiently. Are there any good methods on how to solve this problem? Or at least some indications on how to proceed.