Fast computation of convolution integral of a gaussian function

Question

Given a convolution integral $$ 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 gaussian function
• $f:x\in[a,b] \to \Bbb R$ a known function

I'm seeking a fast algorithm that allows to accurately approximate the function $g(y)$ over an interval $y\in [c,d]$.

I think this problem may be known in image processing or signal processing. Perhaps there exists a clever method using

• the property that $\varphi$ is a Gaussian function and/or
• $g(y)$ is a convolution of two function $\phi$ and $f\cdot \mathbb{I}_{[a,b]}$.

My attempt:

Denote $\Delta y= \frac{d-c}{N_y} $ , $\Delta x= \frac{b-a}{N_x} $, $(y_i,x_i) = (c+i\cdot \Delta y, a+ i\cdot \Delta x)$

Approximate the integal g(y) by $(g(y_i))_{i=1,..,N}$, with $$g(y_i)=\Delta x \cdot\sum_{j=1}^{N_x} \varphi(y_i -x_j) f(x_j) \text{for } i=1,...,N_y$$

The complexity of this method is $\mathcal{O}(N_x N_y)$, which is slow.

I tried to set $\Delta x = \Delta y$ for arranging the terms $\varphi(y_i-x_j)_{i,j}$ but it's difficult because there are cases where $(d-c) \gg (b-a)$ (or $(d-c) \ll (b-a)$) so if I fixed $N_y$ for example, $N_x$ becomes too many or too little.

Answer 1

Convolution with a Gaussian kernel of an $n$-point function has $n^2$ complexity, while Fourier transformation (FFT), multiplication, and inverse Fourier transformation is only of complexity $n\log n$. Here is a Python code for the two-dimensional case.