In the lectures and the problem sets, you learned how to evaluate convolutions analytically. In this notebook, I'll show you how to work them out numerically in python giving you a useful way to check homework answers. Then, I'll show some applications of the convolution theorem.

Finding Convolutions

To convolve two functions, you have to evaluate $$ (f \otimes g)(x) = \int_{-\infty}^\infty f(x) g(x-u) du.$$ As you've found in the homeworks, sometimes this can be tricky, but this is very easy numerically. In the example below, I convolve a Gaussin with a step function to give me a smoothed step function.

Now we'll take an example from the practice problems. We'll convolve \begin{equation} f(x) = \begin{cases} 1, \ |x| < 1 \\ 0, \ \text{otherwise} \end{cases} \end{equation} with \begin{equation} g(x) = \begin{cases} x, \ 0 < x < 1 \\ 0, \ \text{otherwise} . \end{cases} \end{equation} We know what to expect of this from the practice problem \begin{equation} (f\otimes g)(x) = \begin{cases} 0, \ x < 1 \\ \frac{1}{2}(x + 1)^2, \ -1 < x < 0 \\ \frac{1}{2} \ 0 < x < 1 \\ x - \frac{x^2}{2} \ 1 < x < 2 \\ 0 \ x > 2 \end{cases} \end{equation} but let's see what we get numerically. You'll notice that the normalisations are slightly different. When we work out convolutions numerically, the result is normalised to range from 0 to 1.

Application: Circulant Matrices

We also learned about the convolution theorem, which states that the Fourier transform of a convolution is just the product of the Fourier transforms of the two functions $$\mathcal{F}[f \otimes g] = \mathcal{F}[f] \mathcal{F}[g].$$ It turns out this can be used to speed up solving certain linear systems. There is a special type of matrix called a circulant matrix that has the structure \begin{equation} C = \begin{pmatrix} c_0 & c_{n-1} & \cdots & c_2 & c_1 \\ c_1 & c_0 & c_{n-1} & & c_2 \\ \vdots & c_1 & c_0 & \ldots & \vdots \\ c_{n-2} & & \ldots & \ldots & c_{n-1} \\ c_{n-1} & c_{n-2} & \cdots & c_1 & c_0 \end{pmatrix} . \end{equation} If you have the linear system $$\boldsymbol{C} \boldsymbol{x} = \boldsymbol{b},$$ where $\boldsymbol{C}$ is circulant we can note that the elements of $c_{ij}$ only depend upon the difference of their indices $c_{i-j}$ so that the linear system is $$\sum_j c_{i-j} x_j = b_i,$$ which is just a convolution! We can also realise that a circulant matrix can be generated from only the first column. This means that the linear system can be written as $$\boldsymbol{c} \otimes \boldsymbol{x} = \boldsymbol{b}$$ where $\boldsymbol{c}$ is the first column of the circulant matrix $C$. To this we can apply the colvolution theorem to find $\boldsymbol{x}$ \begin{align} \boldsymbol{c} \otimes \boldsymbol{x} &= \boldsymbol{b} \\ \mathcal{F}[\boldsymbol{c} \otimes \boldsymbol{x}] &= \mathcal{F}[\boldsymbol{b}] \\ \mathcal{F}[\boldsymbol{c}] \mathcal{F}[\boldsymbol{x}] &= \mathcal{F}[\boldsymbol{b}] \\ \boldsymbol{x} &= \mathcal{F}^{-1} \left[ \frac{\mathcal{F}[\boldsymbol{b}]}{\mathcal{F}[\boldsymbol{c}]} \right] . \end{align} This converts a matrix inversion into three Fourier transforms. Since numerically evaluating Fourier transforms has complexity $\mathcal{O}(N \log N)$, while solving a linear system using standard methods has complexity $\mathcal{O}(N^{3})$, this can provide a significant speed up when solving large systems providing the matrix is cirulant. This technique is, for example, to model materials made of several dipoles arranged in a regular array. In the code below we'll try this out and see where it really is any faster. Scipy's linear algebra library has a function for constructing circulant matrices.

We see that using two Fourier transforms rather than a direct solution speeds things up a lot!

Application: Image Smoothing

Another use of convolutions is to smooth images. We saw at the start of the notebook that convolving a Gaussian with a step function smoothed the step function; image smoothing uses exactly the same idea. Sharp features of the image are smeared by the convolution.

In the code below, we'll take a grainy image and smooth it out by convolving it with a Gaussain.