Convolution: understand the mathematics

by

Convolution operation is ubiquitous in signal processing applications. The mathematics of convolution is strongly rooted in operation on polynomials. The intent of this text is to enhance the understanding on mathematical details of convolution.

Polynomial functions:

Polynomial functions are expressions consisting of sum of terms, where each term includes one or more variables raised to a non-negative power and each term may be scaled by a coefficient. Addition, Subtraction and multiplication of polynomials are possible.

This article is part of the following books
● Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
● Digital Modulations using Python ISBN: 978-1712321638
● Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

Polynomial functions can involve one or more variables. For example, following polynomial expression is a function of variable x. It involves sum of 3 terms where each term is scaled by a coefficient. Polynomial expression involving two variables ‘x’ and ‘y’ is given next.

f(x)=x^{2}+2x+1

f(x,y)=2x^{4}-4x^{2}y+3xy^{2}+8y^{2}+5

Representing single variable polynomial functions:

Polynomial functions involving single variable is of specific interest here. In general a single variable (say ‘x’) polynomial is expressed in the following sum of terms form, where a_0,a_1,a_2,...,a_{n-1} are coefficients of the polynomial.

f(x)=a_0+a_{1}x+ a_{2}x^{2}+ a_{3}x^{3}+...+ a_{n-1}x^{n-1}

The degree or order of the polynomial function is the highest power of x with a non-zero coefficient. The above equation can be written as

 f(x)=\displaystyle{\sum_{i=0}^{n-1}a_i x^i}

It is also represented by a vector of coefficients as a=[a_0,a_1,a_2,...,a_{n-1}].

Polynomials can also be represented using their roots which is a product of linear terms form, as explained later.

Multiplication of polynomials and linear convolution:

As indicated earlier, mathematical operations like additions, subtractions and multiplications can be performed on polynomial functions. Addition or subtraction of polynomials is straight forward. Multiplication of polynomials is of specific interest in the context of subject discussed here.

Computing the product of two polynomials represented by the coefficient vectors a=[a_0,a_1,a_2,...,a_n] and b=[b_0,b_1,b_2,...,b_n]. The usual representation of such polynomials is given by

p(x)= a_0+a_{1}x+ a_{2}x^{2}+ a_{3}x^{3}+\cdots+ a_{n}x^{n}

q(x)=b_0+b_{1}x+ b_{2}x^{2}+ a_{3}x^{3}+\cdots+ b_{m}x^{m}

The product vector is represented as

p(x).q(x) = a_0b_0+(a_1b_0+b_1a_0)x+\cdots+a_{n}b_{m}x^{n+m}

Or equivalently as

p(x) . q(x)=(p\ast q)(x)=a.b=[c_0,c_1,c_2,\cdots,c_{n+m}] \; , \;where

c_k=\displaystyle{\sum_{i,j:i+j=k}a_i b_j}   \quad\quad k=0,1,\cdots,n+m

Since the subscripts obey the equality i+j=k, changing the subscript j to k-i gives

c_k=\displaystyle{\sum_{i=-\infty}^{\infty}a_ib_{k-i}} \quad\quad k=0,1,\cdots,n+m

Which, when written in terms of indices, provides the most widely used form seen in signal processing text books.

c[k]=\displaystyle{\sum_{i=-\infty}^{\infty}a[i]b[k-i]} \quad\quad k=0,1,\cdots,n+m

This operation is referred as convolution (linear convolution, to be precise), denoted as \ast . It is very closely related to other operations on vectors – cross-correlation, auto-correlation and moving average computation.
Thus when we are computing convolution, we are actually multiplying two polynomials.

Note, that if the polynomials have N and M terms, their multiplication produces N+M-1 terms.

The straightforward algorithm (shown below) that implements the above product expression will consume O(n^{2}) time. A better algorithm is a necessity. It turns out that viewing the polynomials as product of linear factors of complex roots will save much of the computations. This forms the basis for Fast Fourier Transform, where nth root of unity is utilized recursively to do the computation in frequency domain (in signal processing perspective). More about this later.

for i = 0:n+m,
   ci = 0;
for i = 0 :n,
   for k = 0 to m,
      c[i+k] = c[i+k] + a[i] · b[k];

Toeplitz Matrix and Convolution:

Convolution operation of two sequences can be viewed as multiplying two matrices as explained next. Given a LTI (Linear Time Invariant) system with impulse response h[n] and an input sequence x[n], the output of the system y[n] is obtained by convolving the input sequence and impulse response.

y[k]=h[n]*x[n] = \displaystyle{\sum_{i=-\infty}^{\infty}x[i]h[k-i]}

where, the sequence x[n] is of length N and h[n] is of length M.

Assume that the sequence h[n] is of length 4 given by h[n]=[h_0,h_1,h_2,h_3] and the sequence x[n] is of length 3 given by x[n]=[x_0,x_1,x_2]. The convolution h[n] \ast x[n] is given by

y[k]=h[n]*x[n] = \displaystyle{\sum_{i=-\infty}^{\infty}x[i]h[k-i]} \quad k=0,1,\cdots,5

Computing each sample in the convolution,

mathematics of convolution

Note the above result. See how the series x[n] and h[n] multiply with each other from the opposite directions. This gives the reason on why we have to reverse one of the sequences and shift one step a time to do the convolution operation (This is often depicted graphically on standard signal processing texts.. You might have wondered why we have to do this… Now, rejoice that you have understood the reason behind it)

Thus, graphically, in convolution, one of the sequences (say h[n]) is reversed. It is delayed to the extreme left where there are no overlaps between the two sequences. Now, the sample offset of h[n] is increased 1 step at a time. At each step, the overlapping portions of h[n] and x[n] are multiplied and summed. This process is repeated till the sequence h[n] is slid to the extreme right where no more overlaps between h[n] and x[n] are possible.

Writing the final result in the previous equation in matrix form,

\begin{bmatrix} y[0]\\y[1]\\y[2]\\y[3]\\y[4]\\y[5] \end{bmatrix}=\begin{bmatrix}h[0] &0 & 0 \\ h[1] & h[0] & 0 \\ h[2] & h[1] & h[0] \\ h[3] & h[2] & h[1] \\ 0 & h[3] & h[2] \\ 0 & 0 & h[3] \end{bmatrix} \begin{bmatrix}x[0]\\x[1]\\x[2] \end{bmatrix}

When the sequences h[n] and x[n] are represented as matrices, the convolution operation can be equivalently represented as

y=h*x=x*h=\begin{bmatrix}h_0 &0 & 0 \\ h_1 & h_0 & 0 \\ h_2 & h_1 & h_0 \\ h_3 & h_2 & h_1 \\ 0 & h_3 & h_2 \\ 0 & 0 & h_3 \end{bmatrix}\begin{bmatrix}x_0\\x_1\\x_2\end{bmatrix}

The matrix representing the incremental delays of h[n] used in the above equation is a special form of matrix called Toeplitz matrix. Toeplitz matrix have constant entries along their diagonals. Toeplitz matrices are used to model systems that posses shift invariant properties. The property of shift invariance is evident from the matrix structure itself. Since we are modelling a Linear Time Invariant system[1], Toeplitz matrices are our natural choice. On a side note, a special form of Toeplitz matrix called “circulant matrix” is used in applications involving circular convolution and Discrete Fourier Transform (DFT)[2].

Representing the operation used to construct the Toeplitz matrix out of the sequence h as T,

T(h)=\begin{bmatrix} h_0 &0 & 0 \\ h_1 & h_0 & 0 \\ h_2 & h_1 & h_0 \\ h_3 & h_2 & h_1 \\ 0 & h_3 & h_2 \\ 0 & 0 & h_3 \\ \end{bmatrix}

Now, the convolution of h and x is simply a matrix multiplication of Toeplitz matrix T(h) and the matrix representation of x denoted as X

y=h*x=x*h=\begin{bmatrix} h_0 &0 & 0 \\ h_1 & h_0 & 0 \\ h_2 & h_1 & h_0 \\ h_3 & h_2 & h_1 \\ 0 & h_3 & h_2 \\ 0 & 0 & h_3 \\ \end{bmatrix}\begin{bmatrix}x_0\\x_1\\x_2\end{bmatrix}=T(h).X

One can quickly vectorize the convolution operation in matlab by using Toeplize matrices as shown below.

y=toeplitz([h0 h1 h2 h3 0 0],[h0 0 0])*x.';

Continue reading on “methods to compute linear convolution“…

For understanding the usage of toeplitz command in Matlab refer [3]

Rate this article: PoorBelow averageAverageGoodExcellent (18 votes, average: 4.44 out of 5)

References:

[1] Reddi.S.S,”Eigen Vector properties of Toeplitz matrices and their application to spectral analysis of time series”, Signal Processing, Vol 7,North-Holland, 1984,pp 46-56.↗
[2] Robert M. Gray,”Toeplitz and circulant matrices – an overview”,Deptartment of Electrical Engineering,Stanford University,Stanford 94305,USA.↗
[3] Matlab documentation help on Toeplitz command.↗

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
● Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
● Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
● Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author

Wireless Communication Systems in Matlab
Wireless Communication Systems in Matlab
Second Edition(PDF)
PoorBelow averageAverageGoodExcellent (168 votes, average: 3.71 out of 5)

Digital modulations using Python
Digital Modulations using Python
(PDF ebook)
PoorBelow averageAverageGoodExcellent (124 votes, average: 3.60 out of 5)

digital_modulations_using_matlab_book_cover
Digital Modulations using Matlab
(PDF ebook)
PoorBelow averageAverageGoodExcellent (129 votes, average: 3.71 out of 5)
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Post your valuable comments !!!