Fast computation of a convolution integral with Gaussian kernel

Question

Given a convolution integral with gaussian kernel $$ g(y) =\int_a^b\varphi(y-x)f(x)dx=\int_{-\infty}^{+\infty}\varphi(y-x)f(x)\mathbb{I}_{[a,b]}(x)dx $$ where

• $\varphi(x)= \frac{1}{\sqrt{2\pi}}\exp{\left(-\frac{x^2}{2}\right)}$ a unidimensional gaussian function
• $f:x\in[a,b] \to \Bbb R$ a known function

Remark: After searching on the Internet, it seems to me that this problem is a simpler version of Gaussian blur, but I'm not sure.

I'm seeking a fast algorithm (for instance, $\mathcal{N \log (N)}$ of complexity) that allows to accurately approximate the function $g(y)$ over an interval $y\in [c,d]$ (because later, I need to compute a lot of function $g(y)$ with different parameters $a,b,c,d$).

For example, I need an output $(g(y_k))_{k=0,1,...,N}$ with $c = y_0 < y_1 <...< y_N = d$. I don't need values of $g(y)$ with $y \not \in [c,d]$.

Thank you in advance.

---

Currently, I have two attempts

---

First attempt: Fast Fourier Transform method

In theory, this method is simple. We notice that $g = \varphi * (f\cdot \mathbb{I}_{[a,b]})$ Hence $$ \begin{align} \mathcal{F}(g)(\omega) &= \mathcal{F}(\varphi) (\omega) \cdot\mathcal{F} (\mathbb{I}_{[a,b]})(\omega) \\ &= \varphi(\omega) \cdot\mathcal{F} (\mathbb{I}_{[a,b]})(\omega) \end{align} $$ $$\implies g(y) = \mathcal{F}^{-1}\mathcal{F}(g)(y) = \mathcal{F}^{-1}(\varphi(\omega) \cdot\mathcal{F} (\mathbb{I}_{[a,b]})(\omega) ) (y)$$

It suffices to use 2 FFT for the Fourier Transform $\mathcal{F} (\mathbb{I}_{[a,b]})$ and an the Inverse Fourier Transform $\mathcal{F}^{-1} (\mathcal{F}(g))$.

But in reality, it's quite hard for me. It's tricky to determine the right support of the function $f \mathbb{I}_{[a,b]}$ in Python (for example) (please look my first question here).

I have also doubt that FFT is not an accurate method for this kind of problem (please look my second question here). But if it is true, why can people use FFT for Gaussian blur problem?

---

Second attempt: Convolution method .

The integral can be discretized on a equidistance grid $(x_j)_{j=1,...,M}$ as follows $$g(y) \approx \Delta \cdot \sum_{j=1}^M \varphi(y-x_j) f(x_j)$$ with $x_j = a + k\Delta$ for $j = 1,...,M$

If we choose $y_i$ on this grid, for example $c = y_0 = x_0 + L \Delta$ (with $L \in \Bbb Z$) and $y_i = y_0 + i\Delta$ for $i = 1,... \left[\frac{d-c}{\Delta} \right]$.

With this method, it suffices to compute $M$ values of $f(x_k)$, $M+N-1$ values of $\varphi(x_k)$ (instead of $MN$ values of $f(x_k)$ and $\varphi(x_k)$) for $(g(y_i))_{i=0,...,N}$.

However, we can't control $M$, $N$ and $\Delta$, so, the accuracy. Indeed, $M,N$ and $Delta$ depends on the ratio $\rho = \frac{d-c}{b-a}$ and so $$\frac{N}{M} \approx \rho$$

If $\rho$ is too large and if we set $N = N_{\text{max}}$ then $M$ becomes very small and estimation of $g(x_k)$ is not sufficiently accurate (as we compute the integral with a small $M$).

If $\rho$ is too small and if we set $M = M_{\text{max}}$ then $N$ becomes very small and we have only a small $N$ number of $g(y_k)$ which is not sufficient to reconstruct the function $g(y)$ with $y$ in $[c,d]$.

I don't know whether it's a recurring problem in Signal processing or not. And how people deal with it.

Answer 1

"Overlap add" or "overlap save" should work just fine here. See https://en.wikipedia.org/wiki/Overlap%E2%80%93add_method

You'll have to truncate the Gaussian both in time and in frequency but since a Gaussian decays really fast, it's easy to find a length that's "good enough" for your application.

You will also have to pick a sample rate. Since you are doing something numerically in a computer, you are working with discrete signals, so the sampling theorem applies. The sample rate should be at least twice the bandwidth of your signal or the Gaussian, which ever is larger. In practice you want it a little bit higher so it's easier to manage the "transition band" between your highest frequency of interest and the Nyquist frequency.

1. Sample the Gaussian at the chosen sample rate and truncate to desired accuracy. It's best to pick an FFT-friendly length, e.g. a power of two or 3 times a power of two. Let's call this $N$
2. Zero pad to $2N$ and perform and FFT. That's your filter transfer function
3. Sample your input chop it into frames of length N
4. For each frame, zero pad to $2N$ and FFT, multiply with the transfer function, do inverse FFT.
5. The inverse FFT produces $2N$ samples. The output is the first $N$ samples of the current frame plus the last $N$ samples of the previous frame.
6. Continue until all your input is processed

This is 100% accurate in a sense that it produces the exact same result as a discrete convolution sum.

Probably the trickiest part is to get the sampling parameters dialed in to make sure there is no excessive aliasing in either domain, but you need to solve this problem any numerical algorithm.