convolution Archives - GaussianWaves

Cyclic Prefix in OFDM: hands-on demo in Matlab

Synopsis: Cyclic prefix in OFDM, tricks a natural channel to perform circular convolution. This simplifies equalizer design at the receiver. Hands-on demo in Matlab.

Cyclic Prefix-ed OFDM

A cyclic-prefixed OFDM (CP-OFDM) transceiver architecture is typically implemented using inverse discrete Fourier transform (IDFT) and discrete Fourier transform (DFT) blocks (refer Figure 13.3). In an OFDM transmitter, the modulated symbols are assigned to individual subcarriers and sent to an IDFT block. The output of the IDFT block, generally viewed as time domain samples, results in an OFDM symbol. Such OFDM symbols are then transmitted across a channel with certain channel impulse response (CIR). On the other hand, the receiver applies DFT over the received OFDM symbols for further demodulation of the information in the individual subcarriers.

This article is part of the book
● Wireless Communication Systems in Matlab (second edition), ISBN: 979-8648350779 available in ebook (PDF) format and Paperback (hardcopy) format.

In reality, a multipath channel or any channel in nature acts as a linear filter on the transmitted OFDM symbols. Mathematically, the transmitted OFDM symbol (denoted as gets linearly convolved with the CIR and gets corrupted with additive white gaussian noise – designated as . Denoting linear convolution as ‘‘, the received signal in discrete-time can be represented as

The idea behind using OFDM is to combat frequency selective fading, where different frequency components of a transmitted signal can undergo different levels of fading. The OFDM divides the available spectrum in to small chunks called sub-channels. Within the subchannels, the fading experienced by the individual modulated symbol can be considered flat. This opens-up the possibility of using a simple frequency domain equalizer to neutralize the channel effects in the individual subchannels.

Circular convolution and designing a simple frequency domain equalizer

From one of the DFT properties, we know that the circular convolution of two sequences in time domain is equivalent to the multiplication of their individual responses in frequency domain.

Let and are two sequences of length with their DFTs denoted as and respectively. Denoting circular convolution as ,

If we ignore channel noise in the OFDM transmission, the received signal is written as

We can note that the channel performs a linear convolution operation on the transmitted signal. Instead, if the channel performs circular convolution (which is not the case in nature) then the equation would have been

By applying the DFT property given in equation 2,

As a consequence, the channel effects can be neutralized at the receiver using a simple frequency domain equalizer (actually this is a zero-forcing equalizer when viewed in time domain) that just inverts the estimated channel response and multiplies it with the frequency response of the received signal to obtain the estimates of the OFDM symbols in the subcarriers as

Demonstrating the role of Cyclic Prefix

The simple frequency domain equalizer shown in equation 6 is possible only if the channel performs circular convolution. But in nature, all channels perform linear convolution. The linear convolution can be converted into circular convolution by adding Cyclic Prefix (CP) in the OFDM architecture. The addition of CP makes the linear convolution imparted by the channel appear as circular convolution to the DFT process at the receiver (Reference [1]).

Let’s understand this by demonstration.

To simply stuffs, we will create two randomized vectors for and . is of length , the channel impulse response is of length and we intend to use -point DFT/IDFT wherever applicable.

N=8; %period of DFT
s = randn(1,8);
h = randn(1,3);

>> s =   -0.0155    2.5770    1.9238   -0.0629   -0.8105    0.6727   -1.5924   -0.8007
>> h =   -0.4878   -1.5351    0.2355

Now, convolve the vectors and linearly and circularly. The outputs are plotted in Figure 1. We note that the linear convolution and the circular convolution does not yield identical results.

lin_s_h = conv(h,s) %linear convolution of h and s
cir_s_h = cconv(h,s,N) % circular convolution of h and s with period N

lin_s_h = 0.0076 -1.2332 -4.8981 -2.3158  0.9449  0.9013 -0.4468  2.9934  0.8542 -0.1885 
cir_s_h = 0.8618 -1.4217 -4.8981 -2.3158  0.9449  0.9013 -0.4468  2.9934

*Figure 1: Difference between linear convolution and circular convolution*

*Figure 1: Difference between linear convolution and circular convolution*

Let’s append a cyclic prefix to the signal by copying last symbols from and pasting it to its front. Since the channel delay is 3 symbols (CIR of length 3), we need to add atleast 2 CP symbols.

Ncp = 2; %number of symbols to copy and paste for CP
s_cp = [s(end-Ncp+1:end) s]; %copy last Ncp symbols from s and prefix it.

s_cp = -1.5924 -0.8007 -0.0155  2.5770  1.9238 -0.0629 -0.8105  0.6727 -1.5924 -0.8007

Lets assume that we send the cyclic-prefixed OFDM symbol through the channel which perform linear filtering

lin_scp_h = conv(h,s_cp) %linear convolution of CP-OFDM symbol s_cp and the CIR h

lin_scp_h = 0.7767  2.8350  0.8618 -1.4217 -4.8981 -2.3158  0.9449  0.9013 -0.4468  2.9934  0.8542 -0.1885

Compare the outputs due to cir_s_h and lin_scp_h . We can immediately recognize that that the middle section of lin_scp_h is exactly same as the cir_s_h vector. This is shown in the figure 2.

*Figure 2: Equivalence of circular convolution and linear convolution with cyclic prefixed signal*

cir_s_h   = 0.8618  -1.4217 -4.8981 -2.3158  0.9449  0.9013 -0.4468  2.9934
lin_scp_h = 0.7767   2.8350  0.8618 -1.4217 -4.8981 -2.3158  0.9449  0.9013 -0.4468 2.9934 0.8542 -0.1885

We have just tricked the channel to perform circular convolution by adding a cyclic extension to the OFDM symbol. At the receiver, we are only interested in the middle section of the lin_scp_h which is the channel output. The first block in the OFDM receiver removes the additional symbols from the front and back of the channel output. The resultant vector is the received symbol r after the removal of cyclic prefix in the receiver.

r = lin_scp_h(Ncp+1:N+Ncp)%cut from index Ncp+1 to N+Ncp

Verifying DFT property

The DFT property given in equation 5 can be re-written as

To verify this in Matlab, take N-point DFTs of the received signal and CIR. Then, we can see that the IDFT of product of DFTs of and will be equal to the N-point circular convolution of and

R = fft(r,N); %frequency response of received signal
H = fft(h,N); %frequency response of CIR
S = fft(s,N); %frequency response of OFDM signal (non CP)

r1 = ifft(S.*H); %IFFT of product of individual DFTs

display(['IFFT(DFT(H)*DFT(S)) : ',num2str(r1)])
display([cconv(s,h): ', numstr(r)])

IFFT(DFT(H)*DFT(S)) : 0.86175  -1.4217  -4.8981  -2.3158  0.94491  0.90128 -0.44682  2.9934
cconv(s,h):           0.86175  -1.4217  -4.8981  -2.3158  0.94491  0.90128 -0.44682  2.9934

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Reference:

[1] H. Sari, G. Karam, and I. Jeanclaude, Transmission Techniques for Digital Terrestrial TV Broadcasting, IEEE Commun. Magazine, Vol. 33, pp. 100-109, Feb. 1995.↗

Topics in this chapter

Orthogonal Frequency Division Multiplexing (OFDM)
● Introduction
● Understanding the role of cyclic prefix in a CP-OFDM system
□ Circular convolution and designing a simple frequency domain equalizer
□ Demonstrating the role of cyclic prefix
□ Verifying DFT property
● Discrete-time implementation of baseband CP-OFDM
● Performance of MPSK-CP-OFDM and MQAM-CP-OFDM on AWGN channel
● Performance of MPSK-CP-OFDM and MQAM-CP-OFDM on frequency selective Rayleigh channel

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) (172 votes, average: 3.66 out of 5) Checkout Added to cart	Digital Modulations using Python (PDF ebook) (127 votes, average: 3.58 out of 5) Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) (134 votes, average: 3.63 out of 5) Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Construct autocorrelation Matrix in Matlab & Python

Auto-correlation, also called series correlation, is the correlation of a given sequence with itself as a function of time lag. Cross-correlation is a more generic term, which gives the correlation between two different sequences as a function of time lag.

Given two sequences and , the cross-correlation at times separated by lag i is given by ( denotes complex conjugate operation)

Auto-correlation is a special case of cross-correlation, where x=y. One can use a brute force method (using for loops implementing the above equation) to compute the auto-correlation sequence. However, other alternatives are also at your disposal.

Method 1: Auto-correlation using xcorr function

Matlab

For a N-dimensional given vector x, the Matlab command xcorr(x,x) or simply xcorr(x) gives the auto-correlation sequence. For the input sequence x=[1,2,3,4], the command xcorr(x) gives the following result.

>> x=[1 2 3 4]
>> acf=xcorr(x)
acf= 4   11   20   30   20   11  4

Python

In Python, autocorrelation of 1-D sequence can be obtained using numpy.correlate function. Set the parameter mode=’full’ which is useful for calculating the autocorrelation as a function of lag.

import numpy as np
x = np.asarray([1,2,3,4])
np.correlate(x, x,mode='full')
# output: array([ 4, 11, 20, 30, 20, 11,  4])

Method 2: Auto-correlation using Convolution

Auto-correlation sequence can be computed as the convolution between the given sequence and the reversed (flipped) version of the conjugate of the sequence.The conjugate operation is not needed if the input sequence is real.

Matlab

>> x=[1 2 3 4]
>> acf=conv(x,fliplr(conj(x)))
acf= 4   11   20   30   20   11  4

Python

import numpy as np
x = np.asarray([1,2,3,4])
np.convolve(x,np.conj(x)[::-1])
# output: array([ 4, 11, 20, 30, 20, 11,  4])

Method 3: Autocorrelation using Toeplitz matrix

Autocorrelation sequence can be found using Toeplitz matrices. An example for using Toeplitz matrix structure for computing convolution is given here. The same technique is extended here, where one signal is set as input sequence and the other is just the flipped version of its conjugate. The conjugate operation is not needed if the input sequence is real.

Matlab

>> x=[1 2 3 4]
>> acf=toeplitz([x zeros(1,length(x)-1)],[x(1) zeros(1,length(conj(x))-1)])*fliplr(x).'
acf= 4   11   20   30   20   11  4

Python

import numpy as np
from scipy.linalg import toeplitz
x = np.asarray([1,2,3,4])
toeplitz(np.pad(x, (0,len(x)-1),mode='constant'),np.pad([x[0]], (0,len(x)-1),mode='constant'))@x[::-1]
# output: array([ 4, 11, 20, 30, 20, 11, 4])

Method 4: Auto-correlation using FFT/IFFT

Auto-correlation sequence can be found using FFT/IFFT pairs. An example for using FFT/IFFT for computing convolution is given here. The same technique is extended here, where one signal is set as input sequence and the other is just the flipped version of its conjugate.The conjugate operation is not needed if the input sequence is real.

Matlab

>> x=[1 2 3 4]
>> L = 2*length(x)-1
>> acf=ifft(fft(x,L).*fft(fliplr(conj(x)),L))
acf= 4   11   20   30   20   11  4

Python

import numpy as np
from scipy.fftpack import fft,ifft

x = np.asarray([1,2,3,4])

L = 2*len(x) - 1
ifft(fft(x,L)*fft(np.conj(x)[::-1],L))
#output: array([ 4+0j, 11+0j, 20+0j, 30+0j, 20+0j, 11+0j, 4+0j])

Note that in all the above cases, due to the symmetry property of auto-correlation function, the center element represents .

Construction the Auto-correlation Matrix

Auto-correlation matrix is a special form of matrix constructed from auto-correlation sequence. It takes the following form.

The auto-correlation matrix is easily constructed, once the auto-correlation sequence is known. The auto-correlation matrix is a Hermitian matrix as well as a Toeplitz matrix. This property is exploited in the following code for constructing the Auto-Correlation matrix.

Matlab

>> x=[1+j 2+j 3-j] %x is complex
>> acf=conv(x,fliplr(conj(x)))% %using Method 2 to compute Auto-correlation sequence
>>Rxx=acf(3:end); % R_xx(0) is the center element
>>Rx=toeplitz(Rxx,[Rxx(1) conj(Rxx(2:end))])

Python

import numpy as np
x = np.asarray([1+1j,2+1j,3-1j]) #x is complex
acf = np.convolve(x,np.conj(x)[::-1]) # using Method 2 to compute Auto-correlation sequence
Rxx=acf[2:]; # R_xx(0) is the center element
Rx = toeplitz(Rxx,np.hstack((Rxx[0], np.conj(Rxx[1:]))))

Result:

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

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”,Department of Electrical Engineering,Stanford University,Stanford 94305,USA.↗
[3] Matlab documentation help on Toeplitz command.↗

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

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 Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Convolution: understand the mathematics

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.

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.

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 are coefficients of the polynomial.

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

It is also represented by a vector of coefficients as .

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 and . The usual representation of such polynomials is given by

The product vector is represented as

Or equivalently as

Since the subscripts obey the equality , changing the subscript to gives

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

This operation is referred as convolution (linear convolution, to be precise), denoted as . 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 and terms, their multiplication produces terms.

The straightforward algorithm (shown below) that implements the above product expression will consume 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 and an input sequence , the output of the system is obtained by convolving the input sequence and impulse response.

where, the sequence is of length and is of length .

Assume that the sequence is of length 4 given by and the sequence is of length 3 given by . The convolution is given by

Computing each sample in the convolution,

Note the above result. See how the series and 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 ) is reversed. It is delayed to the extreme left where there are no overlaps between the two sequences. Now, the sample offset of is increased 1 step at a time. At each step, the overlapping portions of and are multiplied and summed. This process is repeated till the sequence is slid to the extreme right where no more overlaps between and are possible.

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

When the sequences and are represented as matrices, the convolution operation can be equivalently represented as

The matrix representing the incremental delays of 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 as ,

Now, the convolution of and is simply a matrix multiplication of Toeplitz matrix and the matrix representation of denoted as

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: Note: There is a rating embedded within this post, please visit this post to rate it.

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

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Cookie	Duration	Description
cookielawinfo-checbox-analytics	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Analytics".
cookielawinfo-checbox-analytics	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Analytics".
cookielawinfo-checbox-functional	11 months	The cookie is set by GDPR cookie consent to record the user consent for the cookies in the category "Functional".
cookielawinfo-checbox-functional	11 months	The cookie is set by GDPR cookie consent to record the user consent for the cookies in the category "Functional".
cookielawinfo-checbox-others	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Other.
cookielawinfo-checbox-others	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Other.
cookielawinfo-checkbox-necessary	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookies is used to store the user consent for the cookies in the category "Necessary".
cookielawinfo-checkbox-performance	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Performance".
viewed_cookie_policy	11 months	The cookie is set by the GDPR Cookie Consent plugin and is used to store whether or not user has consented to the use of cookies. It does not store any personal data.