Accurate way to calculate convolution integrals in arbitrary grids

126 Views Asked by At

The convolution integral is defined by: \begin{eqnarray} (f*g)(y)=\int_{-\infty}^{+\infty}dx f(x)g(y-x). \end{eqnarray} Moreover, assume both $\lim_{x\to \pm \infty} f(x),g(x)=0$. Now, suppose we have our functions $f,g$ discretized on a set of $N$ points $\{x_{i}\}_{i=0,...,N-1}$, so we have $f_{i}\equiv f[x_{i}]$. Suppose the set $\{x_{i}\}$ ranges values of $\mathbb{R}$ from scales $x_{\text{max}}~ 10^{5}$ to $x_{min}\sim 10^{-5}$, and that we want to compute $(f*g)[y_{i}]$, with $y_{i}\equiv x_{i}$. How can we perform such operation in an efficient way? A FastFourier Transform method is only of use if the distance between points $\Delta x_{i}=x_{i+1}-x_{i}>0$ is constant; but if the points $\{x_{i}\}$ are unevenly spaced, the method can't be applied. How to compute the integral in that case?

My approach: I have mapped the real axis to the interval $[-1,1]$ by using: \begin{eqnarray} u =\tanh(x/s), \end{eqnarray} where $s$ is a scale factor. Then, I use linear interpolation to calculate the points on a dense set in $u$, and convolve the functions using FFT. To retrieve the result, I interpolate back to $x$ with the inverse mapping. The problem is, the method is not giving satisfactory results. I would like to know how this can be overcomed.