Methods to compute linear convolution

Mathematical details of convolution, its relationship to polynomial multiplication and the application of Toeplitz matrices in computing linear convolution are discussed in the previous article. A short survey of different techniques to compute discrete linear convolution (with Matlab code) is given here.

Definition

Given an 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] \ast x[n] = \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\).

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

Method 1 – Brute-Force Method:

The above equation can simply be implemented using nested for-loops, which consume the most computational time of all the methods given here. Typically, the computational complexity is \(O(n^2)\) time.

Python

def conv_brute_force(x,h):
	"""
	Brute force method to compute convolution
	Parameters:
	x, h : numpy vectors
	Returns:
		y : convolution of x and h
	"""
	N=len(x)
	M=len(h)
	y = np.zeros(N+M-1) #array filled with zeros
	
	for i in np.arange(0,N):
		for j in np.arange(0,M):
			y[i+j] = y[i+j] + x[i] * h[j]
return y

Matlab

for i = 0:n+m,
   yi = 0;
for i = 0 :n,
   for k = 0 to m,
      y[i+k] = y[i+k] + x[i] · h[k];

Method 2 – Using Toeplitz Matrix:

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

\[y=h \ast x=x \ast h=Toeplitz(h) X\]

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] \ast x[n] = \displaystyle{\sum_{i=-\infty}^{\infty}x[i]h[k-i]}, \quad\quad k=0,1,\cdots,5\]

Equivalent representation of the above convolution can be written as

\[\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}\]

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

For python code: refer the book – Digital modulations using Python

Matlab has inbuilt function to compute Toeplitz matrix from given vector. Toeplitz matrix of sequence \(h\) is given as

\[Toeplitz (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}\]

In Matlab, it is constructed by calling Toeplitz function[3] with two arguments. The first argument is the first column of the above matrix -\( [h_0,h_1,h_2,h_3,0,0]\) and the second argument is the first row of the above matrix – \([h_0,0,0]\).

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

In the generalized case, where \(x\) and \(h\) are of any arbitrary lengths (say \(N\) and \(M\)), the code can be modified as follows [here it is assumed that \( length (x) \geq length(h)\) ]

k=toeplitz([h zeros(1,length(x)-1)],[h(1) zeros(1,length(x)-1])*x.'

Method 3: Using FFT:

Computation of convolution using FFT (Fast Fourier Transform) has the advantage of reduced computational complexity when the length of inputs are large. To compute convolution, take FFT of the two sequences \(x\) and \(h\) with FFT length set to convolution output length \(length (x)+length(h)-1\), multiply the results and convert back to time-domain using IFFT (Inverse Fast Fourier Transform). Note that FFT is a direct implementation of circular convolution in time domain. Here we are attempting to compute linear convolution using circular convolution (or FFT) with zero-padding either one of the input sequence. This causes inefficiency when compared to circular convolution. Nevertheless, this method still provides \(O\left(\frac{N}{log_2N}\right)\) savings over brute-force method.

\[y(n) = IFFT[FFT_L (X) \ast FFT_L (H) ], \quad \; L=2^{nextpower2(N+M-1)}\]

Usually, the following algorithm is suffice which ignores additional zeros in the output terms.

\[y(n) = IFFT [ FFT_L (X) \ast FFT_L (H)], \quad \; L=N+M-1\]

Python

from scipy.fftpack import fft,ifft
y=ifft(fft(x,L)*(fft(h,L))) #Convolution using FFT and IFFT

Matlab

y=ifft(fft(x,L).*(fft(h,L)))

Other Methods:

If the input sequence is of infinite length or very very large as in many real time applications, block processing methods like Overlap-Add and Overlap-Save methods can be used to compute convolution in a faster and efficient way.

Standard Algorithms for Fast Convolution:

Standard algorithms for computing convolution that greatly reduce the computational complexity are

1) Cook-Toom Algorithm
2) Modified Cook-Toom Algorithm
3) Winograd Algorithm
4) Modified Winograd Algorithm
5) Iterated Convolution

Refer [4] -Fast Algorithms for Signal Processing by Richard E. Blahut for more details on the above algorithms (see reference section below).

Matlab Code

Implementing method 2 (convolution using Toeplitz matrix transformation) and method 3 (convolution using FFT) and comparing against Matlab’s standard conv function:

%Create random vectors for test
x=rand(1,7)+1i*rand(1,7);
h=rand(1,3)+1i*rand(1,3);
L=length(x)+length(h)-1;

%Convolution Using Toeplitz matrix
y1=toeplitz([h zeros(1,length(x)-1)],[h(1) zeros(1,length(x)-1)])*x.'
%Convolution FFT
y2=ifft(fft(x,L).*(fft(h,L))).'
%Matlab's standard function
y3=conv(h,x).'

Simulation output:

Comparing method 2 and method 3 against Matlab’s standard convolution function returns identical results

\[y1 =\begin{bmatrix} 0.2428 + 0.4536i\\ 0.1512 + 0.6668i\\ 0.8663 + 1.1188i\\ 0.5332 + 1.0162i\\ 0.8589 + 0.9822i\\ 0.4633 + 1.1006i\\ 0.4024 + 1.4706i\\ 0.3225 + 0.9790i\\ 0.2324 + 0.8227i\end{bmatrix}; y2 = \begin{bmatrix}0.2428 + 0.4536i\\ 0.1512 + 0.6668i\\ 0.8663 + 1.1188i\\ 0.5332 + 1.0162i\\ 0.8589 + 0.9822i\\ 0.4633 + 1.1006i\\ 0.4024 + 1.4706i\\ 0.3225 + 0.9790i\\ 0.2324 + 0.8227i\\ \end{bmatrix};y3 = \begin{bmatrix}0.2428 + 0.4536i\\ 0.1512 + 0.6668i\\ 0.8663 + 1.1188i\\ 0.5332 + 1.0162i\\ 0.8589 + 0.9822i\\ 0.4633 + 1.1006i\\ 0.4024 + 1.4706i\\ 0.3225 + 0.9790i\\ 0.2324 + 0.8227i\end{bmatrix}\]

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.↗
[4] Richard E. Blahut,”Fast Algorithms for Signal Processing”,ISBN-13: 978-0521190497,Cambridge University Press,August 16, 2010.↗

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 (162 votes, average: 3.78 out of 5)

Digital modulations using Python
Digital Modulations using Python
(PDF ebook)

PoorBelow averageAverageGoodExcellent (123 votes, average: 3.60 out of 5)

digital_modulations_using_matlab_book_cover
Digital Modulations using Matlab
(PDF ebook)

PoorBelow averageAverageGoodExcellent (126 votes, average: 3.70 out of 5)

Hand-picked Best books on Communication Engineering
Best books on Signal Processing

2 thoughts on “Methods to compute linear convolution”

Post your valuable comments !!!