How do I efficiently evaluate a convolution integral between different-sized domains?

Question

This is a question that I've previously asked over on math.stackexchange, and I have yet to receive a useful answer. It was suggested that I post this here.

The problem itself originally comes from simulation of wave propagation, and involves large two-dimensional arrays of complex data. I'll explain it here as a simplified one-dimensional problem, however.

Suppose I have the following integral relationship involving a convolution: $$ f(x) = \int_{-\infty}^{+\infty}dx'\, g(x')\, K(x - x') $$ Assume that the functional form of the kernel $K$ is known, and suppose that the function $g$ has compact support on $[-a/2, +a/2]$ for some width $a>0$. Assume further that I have data about the function $g$ in the form of $N$ equally-spaced samples $g_n$: $$ g_n = g(x_n)\, , $$ where $$ x_n = -\frac{a}{2} + n\frac{a}{N}\, , \qquad n = 0,..., N-1 $$ My goal is to generate some $N$ evenly-spaced samples of the function $f$.

If we suppose that we want to evaluate $f$ on the same domain on which $g$ lives, and for which we have samples of $g$, then it's easy. I discretize everything in the usual way, and the whole thing becomes a discrete convolution: $$ f_n = f(x_n) = \sum_{m = 0}^{N-1} g_m K_{n - m}\, , $$ where $$ K_{n - m} = K(x_n - x_m) = K\left(\frac{a}{N}(n-m)\right)\, . $$ Everything is simple since I can use the discrete convolution theorem: I FFT the $\{g_n\}$ and $\{K_n\}$ sequences, multiply them, do an inverse FFT, and I'm done. I have my $\{f_n\}$ data. The whole thing will be an $O(N\log N)$ process.

Suppose instead that I have to evaluate $f$ on a different-sized domain. For example, perhaps the $g$ data is narrow, but the kernel $K$ is very wide. Let's say I have to generate $N$ samples of $f$ on the domain $[-L/2, +L/2]$, for some domain width $L > a$. $L$ could even be much larger than a.

If I try to discretize everything in a similar way, I'd have to first define a sequence of points $\{X_n\}$ in the wider domain: $$ X_n = -\frac{L}{2} + n\frac{L}{N}\, , \qquad n = 0,..., N-1 $$ Then the data I seek is $\{f_n\}$ where $f_n = f(X_n)$. The original integral is still over the smaller domain, so we have: \begin{align} f_n &= \int_{-a/2}^{+a/2}dx'\, g(x')\, K(X_n - x') \\ &\approx \sum_{m = 0}^{N-1} g_m\, K(X_n - x_m) \qquad \text{(discretize)}\\ &= \sum_{m = 0}^{N-1} g_m\, K\left(\frac{L-a}{2} + \frac{a}{N}(\alpha\, n - m)\right)\, , \end{align} where $\alpha = L/a > 1$ is the ratio of the two domain sizes.

This is no longer a discrete convolution, since the second term is no longer a function of $n-m$, but instead of $\alpha\, n - m$. So I can no longer use the convolution theorem. I can imagine evaluating the sum above by brute force, but that would be an $O(N^2)$ process, which is prohibitive for my problem.

Is there any way of efficiently calculating the discrete version of this convolution integral between two different domain sizes? Thanks.

Answer 1

Essentially, to my understanding this problems boils down to a fast-convolution problem:

You have a function $g(x)$, which has a small domain and you have a function $K(\tau)$, which has a big domain and you want to calculate

$$f(x) = \int_{-\infty}^{\infty}g(x')K(x-x')dx'$$

in a discretized fashion, where you want to calculate as efficient as possible. Especially, you are interested in $f(-L/2+nL/N)$ for $n=0,\dots,N-1$.

When considering the discretization you need to take care of the bandwidth of $g$ and $K$ and you need to sample with the according Nyquist frequency at least. Let's assume that $g$ can be sufficiently sampled with $g(-\frac{a}{2}+n\frac{a}{N})$ (i.e. the sampling is fast enough regarding Nyquist frequency). Let's also assume, that $K(\tau)$ has the same or smaller bandwith, i.e. the sampling rate for $g$ is also sufficient for $K$.

Then, let's further assume $L=qa, q\in\mathbb{N}$, i.e. L is a multiple of $a$. Define the following functions:

$$g[n]=g(-\frac{L}{2}+n\frac{L}{qN}), \quad n=0,\dots,qN-1$$ and $$K[n]=K(-\frac{L}{2}+n\frac{L}{qN}), \quad n=0,\dots,qN-1$$

Note, that $g[n]$ is zero for most $n$, but that is fine here.

Now, you have

$$f[n]=\sum_n'g[n']*K[n-n'], \quad n=0,\dots,qN-1$$

Which is a higher-resolution version of your original $f_n$ (It has qN instead of only N samples within the interval $\pm\frac{L}{2}$). So, what you actually want to have is a downsampled version of this $f[n]$ given by

$f_n=f[qn]$ (where $f_n$ is the $f_n$ from your post, i.e. with $n=0,\dots,N-1$.

So, this is also a discrete convolution problem, and it can be solved via FFT (note that the FFT actually considers circular convolution, if you don't zero-pad). However, since your $g[n]$ has a much smaller domain that $K[n]$, you can resort to fast-convolution algorithms (Overlap-save or Overlap-add).

To summarize: The trick is actually just to sample both signals with the same sampling frequency (accepting that $g[n]$ contains a lot of zeros) and then understanding this problem as a simple discrete convolution.