How do I efficiently evaluate a convolution integral between different-sized domains? 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&gt;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 &gt; 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
&amp;= \int_{-a/2}^{+a/2}dx'\, g(x')\, K(X_n - x') \\
&amp;\approx \sum_{m = 0}^{N-1} g_m\, K(X_n - x_m) \qquad \text{(discretize)}\\
&amp;= \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 &gt; 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.
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.