# Fast computation of a convolution integral with Gaussian kernel

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]$$.

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.

• Convolving with a Gaussian is Gaussian blur. There are many different methods, but direct computation of the convolution sum is the fastest option if the sigma of the Gaussian is small w.r.t. the sample rate (eg sigma < 10 samples). There is an O(N) algorithm (independent of the sigma), it yields a pretty good approximation. The FFT approach is only useful IMO for very small sigma, where sampling the Gaussian is just not precise enough (eg sigma < 1 sample). Jul 18, 2021 at 0:48
• @CrisLuengo Thank you for the comment. There exists an $\mathcal{O}(N)$ algorithm? Could you please tell me the name of this algorithm?
– NN2
Jul 18, 2021 at 1:19
• See the literature for this function. Jul 18, 2021 at 1:24
• – uhoh
Aug 16, 2021 at 15:28
• @uhoh It's really helpful. I'll read it. Thank you very much!
– NN2
Aug 18, 2021 at 18:41

"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.

• Thank you very much for the answer. I'm trying to use the FFT. But I have some questions. It seems to me that FFT is only efficient for long signals, isn't it? Indeed, the complexity is $\mathcal{O}(N \log N)$ but for $N$ large and the real algorithm complexity is $C \cdot N \log N$ with a large constant $C$. Is my understanding is correct? And so for the case where the number of data points for each signal is less than 1000, do you think it's worthwise to use FFT? And are overlap methods efficient for small data?
– NN2
Jul 20, 2021 at 13:32
• If it's less than a 1000 points, why do you care about efficiency? Any processor these days will solve this in a fraction of a second no matter what algorithm you use . Jul 21, 2021 at 13:07
• In fact, I should apply this many times for a multilayer problem. Precisely, I need to compute consecutively $$f_{n+1}(x)=\int_{a_n}^{b_n}\varphi(x-t)f_n(t)dt$$ for $n=1,…,N$ ($N$ around 20). And calculate the multilayer problem for 10 millions times for each point in a grid of 1 thousand points. That’s reason why I need something fast enough.
– NN2
Jul 21, 2021 at 15:24
• I see, that makes sense. In this case the "optimum" solution will depend a lot on the details: what's the spectral content or you data, what's the cutoff frequency of your gaussian, how much numerical error can you accept, etc. Jul 22, 2021 at 15:09
• Another thought: can you do the entire problem in the frequency domain? Jul 22, 2021 at 15:17