Fast evaluation of an integral convolution with an "expanding kernel"

165 Views Asked by At

Suppose I have a 1-D integral convolution transform like this: $$ g(x) = \int_{-\infty}^{+\infty} dy\, f(y)\, K(x-y). \qquad (1) $$ Say the kernel $K(x)$ is a known analytic function, and say we have an analytic form for its Fourier transform $$ \hat{K}(q) \equiv \int_{-\infty}^{+\infty}dx\, K(x)\, \exp(-iqx)\, . $$ Suppose furthermore that the function $f(x)$ is entirely supported on the interval $[-L/2, +L/2]$ of length $L$, and that instead of an analytical form for $f$ we have $N$ samples of $f$ at evenly spaced points $$ x_n \equiv -\frac{L}{2} + \frac{n L}{N}\, . $$

Given all this, and assuming $g(x)$ has the same support as $f(x)$ I might choose to approximate the integral (1) for $g(x)$ as a discrete sum: \begin{align} g_n \equiv g(x_n) &\approx \frac{L}{N} \sum_{m=0}^{N-1} f(x_m)\, K(x_n - x_m)\\ &= \frac{L}{N} \sum_{m=0}^{N-1} f_m\, \phi_{n - m}\, ,\qquad (2) \end{align} where $f_m \equiv f(x_m)$, $\phi_{n - m} \equiv K\left(\frac{L}{N}(n-m)\right) = K(x_n - x_m)$, and the $L/N$ term is to ensure the correct units.

We could evaluate the sum (2) directly by brute force for each of the $N$ points to get the $N$ samples ${\{g_n\}}$. This would be an $O(N^2)$ process. In practice, we would probably take advantage of the convolution theorem using Fast Fourier Transforms (FFTs): $$ g_n \approx \frac{L}{N} {\left\{{FFT}^{-1}\{\hat{f} \cdot \hat{\phi}\}\right\}}_n $$ Here $\hat{f}$ and $\hat{\phi}$ are the FFTs of $f$ and $\phi$, and we'd probably also take advantage of the fact that we have an analytic expression for $\hat{\phi}$ to save a little time. The overall calculation would now be $O(N\log N)$. Great.

Now, suppose instead that $K$ is an "expanding kernel", i.e. it spreads things out so that the support of $g(x)$ is now much larger. Say the support of $g$ is $[-L'/2, +L'/2]$, where $L' = \alpha L$ and $\alpha > 1$. To get the same number of samples of $g$ spread out over this wider space, we define another grid: $$ {x'}_n \equiv -\frac{L'}{2} + \frac{n L'}{N} $$ Armed with this, the discrete analog of the integral in Eq. (1) becomes: $$ g({x'}_n) \equiv g_n \approx \frac{L}{N}\sum_{m=0}^{N-1} f_m \, \varphi(\alpha\, n - m)\, ,\qquad (3) $$ where $\varphi(\alpha\, n - m) \equiv K({x'}_n - x_m)$. This is no longer a convolution. I could evaluate this via brute force in $O(N^2)$, but this is too slow if $N$ is large.

My question: Are there fast methods for evaluating a not-quite-a-convolution like Eq. (3)? Or are there other (fast) methods for evaluating Eq. (1), where $g$ has a much wider support than $f$, that I'm overlooking?

Thanks!