# Convolution

Let $y[n]$ represent the discrete time signal that results when $x[n]$ is convolved with $h[n]$. Its value at time $n$ is computed as:

\begin{equation} \begin{aligned} y[n] = x[n] * h[n] &= \sum_k x[k]h[n-k] = \sum_k h[k]x[n-k] = \\ &= \dots + hx[n-2] + hx[n-1] + hx[n] + h[-1]x[n+1] + h[-2]x[n+2] + \dots \end{aligned} \end{equation}

For all the examples below we will consider that the signal $x[n]$ is finite of length $N$ and the filter $h[n]$ is finite of length $M$. We will use the convolve() function from the scipy.signal library which features three different modes that we will further discuss in the table below. Note that the blue-shaded boxes indicate the alignment point between x and the flipped filter h for the computation of the first (i.e. index 0) and last samples of y.

Mode Visual Representation Explanation
full It is the default option of the function and corresponds to the mathematical definition of convolution.

The resulting array has length $N+M-1$.
valid The output consists only of those elements that do not rely on the zero-padding.

The length of the resulting array is decreased by $M-1$ with respect to the original signal.
same The filter has its initial position centered around the first element of the signal.

Hence, the output is the same size $N$ as the input signal.

In order to illustrate the concept of convolution we will take the discrete-time sinusoid of frequency F0 = 0.005,

\begin{equation} x[n] = 10 + 2\cos(2\pi 0.005 n) \end{equation}

and we will add white gaussian noise of $\sigma^2 = 1$ to the pure signal.

import numpy as np
import matplotlib.pyplot as plt

F0 = 0.005
L = 1000
n = np.arange(L)
pure = 10 + 2*np.cos(2*np.pi*F0*n)

noise = np.random.normal(0, 1, pure.shape)
noisy = pure + noise

plt.title(r'Discrete-time sinusoid of $F_0$=0.005')
plt.plot(noisy, '-b', linewidth=0.5, label='Noisy signal')
plt.plot(pure, '--r', linewidth=4, label='Original signal')
plt.xlabel('n')
plt.legend(loc='upper right')
plt.axis('tight') To smooth out these fluctuations and highlight the trend of the signal, we are interested in performing a moving average which is the unweighted mean of $M$ samples. Hence, we will consider the filter

\begin{equation} h[n] = \bigg[\underbrace{\frac{1}{M}, \frac{1}{M}, \dots, \frac{1}{M}}_{M}\bigg] \end{equation}

to convolve it with the noisy signal, taking $M=67$:

from scipy.signal import convolve

M = 67
h = np.repeat(1/M, M)
y = convolve(noisy, h, mode='full')

plt.figure()
plt.plot(y, '-b', linewidth=1, label='Convolved signal')
plt.plot(pure, '--r', linewidth=1, label='Original signal')
plt.xlabel('n')
plt.legend(loc='lower left')
plt.axis('tight') The length of the output array is $1066=1000+67-1$. Since we have used mode='full', the output at a given position corresponds to the average of the input at the same time instant with the $M-1$ previous samples. Hence, y has been computed using exactly $66$ values that are part of the zero padding, and

\begin{equation} \begin{aligned} y &= \frac{x}{M},\\ y &= \frac{1}{M}\bigg(x + x\bigg) \\ &\ \ \vdots \end{aligned} \end{equation}

Comparing the shapes of both input and ouput, we observe that the output is delayed with respect to the input (a phase shift has ocurred).

If we now change the mode to 'same',

y = convolve(noisy, h, mode='same')


the length of the output array is equal to the length of the input array: The first thing we notice is that the original signal and the output of the convolution look aligned. However, both the first and last values are largely underestimated. Since we have used mode='same', the output for a given position corresponds to the average of the previous $33$ samples, itself and the next $33$. Hence, y has been computed using exactly $33$ values that are part of the zero padding.

We can also use mode='valid'. In such a case, the output array will have length $934 = 1000 - 67 + 1$. The first computed sample in this case, i.e. y, will be equal to the 66th sample when mode='full' is used, and to the 33rd sample when mode='same' is used.

Finally, if instead of taking the averaging filter we consider h = [-1, 1], we actually remove the expected value of the signal and end up having white noise: This figure has been obtained with mode='valid' and the length of the output array is $999 = 1000 - 2 + 1$.

## Convolution in 2 dimensions

Analogous to one-dimensional convolutions, we have the convolve2d() function for the two-dimensional case. If we take for example the image of this circuit and aspire to detect the horizontal contours of the image, we can use the following filter,

$$h[m, n] = \begin{pmatrix} 1 & 1 & 1 \\ 0 & 0 & 0 \\ -1 & -1 & -1 \end{pmatrix}$$

also known as the Prewitt filter. If we want the ouput image to have the same dimensions as the input image we will use the convolution with mode='same'. In such a case, a dark border may appear if we use zero padding (although for images other types of padding are also common). If we use mode='valid' the resulting image will have smaller dimensions.

from scipy.signal import convolve2d

hor = np.array([[1,1,1],
[0,0,0],
[-1,-1,-1]])
circuit_hor = convolve2d(circuit, hor, mode='valid')

plt.figure()
plt.imshow(circuit_hor, cmap='gray')
plt.show() Given that we have used mode='valid', we can indeed check that the shape of the circuit has been reduced by $M-1 = 2$ in each dimension:

>>> circuit.shape
(449, 960)
>>> circuit_hor.shape
(447, 958)


Just like before, if we compute the horizontal gradient of the image we will detect vertical contours. The effect of this filter

$$h[m, n] = \begin{pmatrix} 1 & 0 & -1 \\ 1 & 0 & -1 \\ 1 & 0 & -1 \end{pmatrix}$$

on the circuit is displayed below: Finally, the averaging filter

$$h[m, n] = \begin{pmatrix} \frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\ \frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\ \frac{1}{9} & \frac{1}{9} & \frac{1}{9} \end{pmatrix},$$

which can be expressed as np.ones((3, 3))/9, has the effect of blurring the image.  Lliçons.jutge.org