Signal Processing Archives - Page 2 of 2

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

Sampling Theorem – Bandpass or Intermediate or Under Sampling

Prerequisite: Sampling theorem – baseband sampling

Intermediate Sampling or Under-Sampling

A signal is a bandpass signal if we can fit all its frequency content inside a bandwidth . Bandwidth is simply the difference between the lowest and the highest frequency present in the signal.

“In order for a faithful reproduction and reconstruction of a bandpass analog signal with bandwidth –

, the signal should be sampled at a Sampling frequency (

) that is greater than or equal to twice the maximum bandwidth of the signal:

Consider a bandpass signal extending from 150Hz to 200Hz. The bandwidth of this signal is . In order to faithfully represent the above signal in the digital domain the sampling frequency must be . Note that the sampling frequency 100Hz is far below the maximum content of the signal (which is 200Hz). That is why the bandpass sampling is also called “under-sampling”. As long as the sampling frequency is greater than or equal to twice the bandwidth of the signal, the reconstruction back to analog domain will be error free.

Going back to the aliasing zone figure, if the signal of interest is in the zone other than zone 1, it is called a bandpass signal and the sampling operation is called “Intermediate Sampling” or “Harmonic Sampling” or “Under Sampling” or “Bandpass Sampling”.

Note that zone 1 is a mirror image of zone 2 (with frequency reversal). Similarly zone 3 is a mirror image of zone 4 etc.., Also, any signal in zone 1 will be reflected in zone 2 with frequency reversal which inturn will be copied in zone 3 and so on.

Lets say the signal of interest lies in zone 2. This will be copied in all the other zones. Zone 1 also contains the sampled signal with frequency reversal which can be correct by reversing the order of FFT in digital domain.

No matter in which zone the signal of interest lies, zone 1 always contains the signal after sampling operation is performed. If the signal of interest lies in any of the even zones, zone 1 contains the sampled signal with frequency reversal. If the signal of interest lies in any of the odd zones, zone 1 contains the sampled signal without frequency reversal.

Example:

Consider an AM signal centered at carrier frequency 1MHz, with two components offset by 10KHz – 0.99 MHz and 1.01 MHz. So the AM signal contains three frequency components at 0.99 MHz, 1 MHz and 1.01 MHz.

Our desire is to sample the AM signal. As with the usual sampling theorem (baseband), we know that if we sample the signal at twice the maximum frequency i.e Fs>=2*1.01MHz=2.02 MHz there should be no problem in representing the analog signal in digital domain.

By the bandpass sampling theorem, we do not need to use a sampler running at Fs>=2.02 MHz. Faster sampler implies more cost. By applying the bandpass sampling theorem, we can use a slower sampler and reduce the cost of the system. The bandwidth of the signal is 1.01MHz-0.99 MHz = 20 KHz. So, just sampling at will convert the signal to digital domain properly and we can also avoid using an expensive high rate sampler (if used according to baseband sampling theorem).

Lets set the sampling frequency to be (which is 3 times higher than the minimum required sampling rate of 40KHz or oversampling rate =3).

Now we can easily find the position of the spectral components in the sampled output by using the aliasing zone figure as given above. Since , will be 60KHz. So the zone 1 will be from 0 to 60 KHz, zone 2 -> 60-120KHz and so on. The three spectral components at 0.99MHz, 1MHz and 1.01 MHz will fall at zone 17 ( how ? 0.99 MHz/60 KHz = 16.5 , 1MHz/60KHz = 16.67 and 1.01MHz/60KHz = 16.83 , all figures approximating to 17). By the aliasing zone figure, zone 16 contains a copy of zone 17, zone 15 contains a copy of zone 16, zone 14 contains a copy of zone 15 and so on… Finally zone 1 contains the copy of zone 2 (Frequency reversal also exist at even zones). In effect, zone 1 contains a copy of zone 17. Since the original spectral components are at zone 17, which is an odd zone, zone 1 contains the copy of spectral components at zone 17 without frequency reversal.

Since there is no frequency reversal, in zone 1 the three components will be at 30KHz, 40KHz and 50KHz (You can easily figure this out ).

This operation has downconverted our signal from zone 17 to zone 1 without distorting the signal components. The downconverted signal can be further processed by using a filter to select the baseband downconverted components. Following figure illustrates the concept of bandpass sampling.

Consider the same AM signal with three components at 0.99MHz, 1MHz and 1.01 MHz. Now we also have an “unwanted” fourth component at 1.2 MHz along with the incoming signal. If we sample the signal at 120 KHz, it will cause aliasing (because the bandwidth of the entire signal is 1.2-0.99 = 0.21 MHz = 210 KHz and the sampling frequency of 120KHz is below twice the bandwidth). In order to avoid anti-aliasing and to discard the unwanted component at 1.2 MHz, an anti-aliasing bandpass filter has to be used to select those desired component before performing the sampling operation at 120KHz. This is also called “pre-filtering”. The following figure illustrates this concept.

Sampling Theorem – Baseband Sampling

For Matlab demo of sampling see here.

“Nyquist-Shannon Sampling Theorem” is the fundamental base over which all the digital processing techniques are built. Processing a signal in digital domain gives several advantages (like immunity to temperature drift, accuracy, predictability, ease of design, ease of implementation etc..,) over analog domain processing.

Analog to Digital conversion:

In analog domain, the signal that is of concern is continuous in both time and amplitude. The process of discretization of the analog signal in both time domain and amplitude levels yields the equivalent digital signal. The conversion of analog to digital domain is a three step process 1) Discretization in time – Sampling 2) Discretization of amplitude levels – Quantization 3) Converting the discrete samples to digital samples – Coding/Encoding

The sampling operation samples (“chops”) the incoming signal at regular interval called “Sampling Rate” (denoted by ). Sampling Rate is determined by Sampling Frequency (denoted by ) as Lets consider the following logical questions: * Given a real world signal, how do we select the sampling rate in order to faithfully represent the signal in digital domain ? * Is there any criteria for selecting the sampling rate ? * Will there be any deviation if the signal is converted back to analog domain ? Answer : Consult the “Nyquist-Shannon Sampling Theorem” to select the sampling rate or sampling frequency.

Nyquist-Shannon Sampling Theorem:

The following sampling theorem is the exact reproduction of text from Shannon’s classic paper^[1],

“If a function f(t) contains no frequencies higher than W cps, it is completely determined by giving its ordinates at a series of points spaced 1/2W seconds apart.”

Sampling Theorem mainly falls into two categories : 1) Baseband Sampling – Applied for signals in the baseband (useful frequency components extending from 0Hz to some Fm Hz) 2) Bandpass Sampling – Applied for signals whose frequency components extent from some F1 Hz to F2Hz (where F2>F1) In simple terms, the Nyquist Shannon Sampling Theorem for Baseband can be explained as follows

Baseband Sampling:

If the signal is confined to a maximum frequency of Fm Hz, in other words, the signal is a baseband signal (extending from 0 Hz to maximum Fm Hz).

In order for a faithful reproduction and reconstruction of an analog signal that is confined to a maximum frequency Fm, the signal should be sampled at a Sampling frequency (Fs) that is greater than or equal to twice the maximum frequency of the signal:

Consider a 10Hz sine wave in analog domain. The maximum frequency present in this signal is Fm=10 Hz (obviously no doubt about it !!!). Now, to satisfy the sampling theorem that is stated above and to have a faithful representation of the signal in digital domain, the sampling frequency can be chosen as Fs >=20Hz. That is, we are free to choose any number above 20 Hz. Higher the sampling frequency higher is the accuracy of representation of the signal. Higher sampling frequency also implies more samples, which implies more storage space or more memory requirements. In time domain, the process of sampling can be viewed as multiplying the signal with a series of pulses (“pulse train) at regular interval – Ts. In frequency domain, the output of the sampling process gives the following components – Fm (original frequency content of the signal), Fs±Fm,2Fs±Fm,3Fs±Fm,4Fs±Fm and so on and on…

Now the sampled signal contains lots of unwanted frequency components (Fs±Fm,2Fs±Fm,…). If we want to convert the sampled signal back to analog domain, all we need to do is to filter out those unwanted frequency components by using a “reconstruction” filter (In this case it is a low pass filter) that is designed to select only those frequency components that are upto Fm Hz. The above process mentions only the sampling part which samples the incoming analog signal at regular intervals. Actually a quantizer will follow the sampler which will discretize (“quantize”) amplitude levels of the sampled signal. The quantized amplitude levels are sent to an encoder that converts the discrete amplitude levels to binary representation (binary data). So when converting the binary data back to analog domain, we need a Digital to Analog Converter (DAC) that converts the binary data to analog signal. Now the converted signal after the DAC contains the same unwanted frequencies as well as the wanted component. Thus a reconstruction filter with proper cut-off frequency has to placed after the DAC to filter out only the wanted components.

Aliasing and Anti-aliasing:

Consider a signal with two frequency components f1=10Hz – which is our desired signal and f2=20Hz – which is a noise. Let’s say we sample the signal at 30Hz. The first frequency component f1=10Hz will generate following frequency components at the output of the multiplier (sampler) – 10Hz,20Hz,40Hz,50Hz,70Hz and so on. The second frequency component f2=20Hz will generate the following frequency components at the output of the multiplier – 20Hz,10Hz,50Hz,40Hz,80Hz and so on…

Note the 10Hz component that is generated by f2=20Hz. This 10Hz component (which is a manifestation of noisy component f2=20Hz) will interfere with our original f1=10Hz component and are indistinguishable. This 10Hz component is called “alias” of the original component f2=20Hz (noise). Similarly the 20Hz component generated by f1=10Hz component is an “alias” of f1=10Hz component. This 20Hz alias of f1=10Hz will interfere with our original component f2=20Hz and are indistinguishable. We do not need to care about the interference that occurs at 20Hz since it is a noise and any way it has to be eliminated. But we need to do something about the aliasing component generated by the f2=20Hz. Since this is a noise component, the aliasing component generated by this noise will interfere with our original f1=10Hz component and will corrupt it. Aliasing depends on the sampling frequency and its relationship with the frequency components. If we sample a signal at Fs, all the frequency components from Fs/2 to Fs will be alias of frequency components from 0 to Fs/2 and vice versa. This frequency – Fs/2 is called “Folding frequency” since the frequency components from Fs/2 to Fs folds back itself and interferes with the components from 0Hz to Fs/2 Hz and vice versa. Actually the aliasing zones occur on the either sides of 0.5Fs, 1.5Fs, 2.5Fs,3.5Fs etc… All these frequencies are also called “Folding Frequencies” that causes frequency reversal. Similarly aliasing also occurs on either side of Fs,2Fs,3Fs,4Fs… without frequency reversals. The following figure illustrates the concept of aliasing zones.

In the above figure, zone 2 is just a mirror image of zone 1 with frequency reversal. Similarly zone 2 will create aliases in zone 3 (without frequency reversal), zone 3 creates mirror image in zone 4 with frequency reversal and so on… In the example above, the folding frequency was at Fs/2=15Hz, so all the components from 15Hz to 30Hz will be the alias of the components from 0Hz to 15Hz. Once the aliasing components enter our band of interest, it is impossible to distinguish between original components and aliased components and as a result, the original content of the signal will be lost. In order to prevent aliasing, it is necessary to remove those frequencies that are above Fs/2 before sampling the signal. This is achieved by using an “anti-aliasing” filter that precedes the analog to digital converter. An anti-aliasing filter is designed to restrict all the frequencies above the folding frequency Fs/2 and therefore avoids aliasing that may occur at the output of the multiplier otherwise.

Thus, a complete design of analog to digital conversion contains an anti-aliasing filter preceding the ADC and the complete design of digital to analog conversion contains a reconstruction filter succeeding the DAC.

Note: Remember that both the anti-aliasing and reconstruction filters are analog filters since they operate on analog signal. So it is imperative that the sampling rate has to be chosen carefully to relax the requirements for the anti-aliasing and reconstruction filters.

Effects of Sampling Rate:

Consider a sinusoidal signal of frequency Fm=2MHz. Lets say that we sample the signal at Fs=8MHz (Fs>=2*Fm). The factor Fs/Fm is called “over-sampling factor”. In this case we are over-sampling the signal by a factor of Fm=8MHz/2MHz = 4. Now the folding frequency will be at Fs/2 = 4MHz. Our anti-aliasing filter has to be designed to strictly cut off all the frequencies above 4MHz to prevent aliasing. In practice, ideal brick wall response for filters is not possible. Any filter will have a transition band between pass-band and stop-band. Sharper/faster roll off transition band (or narrow transition band) filters are always desired. But such filters are always of high orders. Since both the anti-aliasing and reconstruction filters are analog filters, high order filters that provide faster roll-off transition bands are expensive (Cost increases proportionally with filter order). The system also gets bulkier with increase in filter order.Therefore, to build a relatively cheaper system, the filter requirement in-terms of width of the transition band has to be relaxed. This can be done by increasing the sampling rate or equivalently the over-sampling factor. When the sampling rate (Fs) is increased, the distance between the maximum frequency content Fm and Fs/2 will increase. This increase in the gap between the maximum frequency content of the signal and Fs/2 will ease requirements on the transition bands of the anti-aliasing analog filter. Following figure illustrates this concept. If we use a sampling frequency of Fs=8MHz (over-sampling factor = 4), the transition band is narrower and it calls for a higher order anti-aliasing filter (which will be very expensive and bulkier). If we increase the sampling frequency to Fs=32MHz (over-sampling factor = 32MHz/2MHz=16), the distance between the desired component and Fs/2 has greatly increased that it facilitates a relatively inexpensive anti-aliasing filter with a wider transition band. Thus, increasing the sampling rate of the ADC facilitates simpler lower order anti-aliasing filter as well as reconstruction filter. However, increase in the sampling rate calls for a faster sampler which makes ADC expensive. It is necessary to compromise and to strike balance between the sampling rate and the requirement of the anti-aliasing/reconstruction filter.

Application : Up-conversion

In the above examples, the reconstruction filter was conceived as a low pass filter that is designed to pass only the baseband frequency components after reconstruction. Remember that any frequency component present in zone 1 will be reflected in all the zones (with frequency reversals in even zones and without frequency reversals in odd zones). So, if we design the reconstruction filter to be a bandpass filter selection the reflected frequencies in any of the zones expect zone 1, then we achieve up-conversion. In any communication system, the digital signal that comes out of a digital signal processor cannot be transmitted across as such. The processed signal (which is in the digital domain) has to be converted to analog signal and the analog signal has to be translated to appropriate frequency of operation that fits the medium of transmission. For example, in an RF system, a baseband signal is converted to higher frequency (up-conversion) using a multiplier and oscillator and then the high frequency signal is transmitted across the medium. If we have a band-pass reconstruction filter at the output of the DAC, we can directly achieve up-conversion (which saves us from using a multiplier and oscillator). The following figure illustrates this concept.

Reference:

[1] Claude E. Shannon, “Communication in the presence of noise” ,Proceedings of the IRE, Vol 37, no.1,pp.10-21,Jan 1949

Correlative Coding – Modified Duobinary Signaling

Modified Duobinary Signaling is an extension of duobinary signaling. It has the advantage of zero PSD at low frequencies (especially at DC ) that is suitable for channels with poor DC response. It correlates two symbols that are 2T time instants apart, whereas in duobinary signaling, symbols that are 1T apart are correlated.

The general condition to achieve zero ISI is given by

As discussed in a previous article, in correlative coding , the requirement of zero ISI condition is relaxed as a controlled amount of ISI is introduced in the transmitted signal and is counteracted in the receiver side

In the case of modified duobinary signaling, the above equation is modified as

which states that the ISI is limited to two alternate samples. Here a controlled or “deterministic” amount of ISI is introduced and hence its effect can be removed upon signal detection at the receiver.

Modified Duobinary Signaling:

The following figure shows the modified duobinary signaling scheme (click to enlarge).

Encoding Process:

1) a_n = binary input bit; a_n ∈ {0,1}.
2) b_n = NRZ polar output of Level converter in the precoder and is given by,

where a_k is the precoded output (before level converter).

3) y_n can be represented as

Note that the samples b_n are uncorrelated ( i.e either +d for “1” or -d for “0” input). On the other-hand,the samples y_n are correlated ( i.e. there are three possible values +2d,0,-2d depending on ak and a_k-2). Meaning that the modified duobinary encoding correlates present sample a_k and the previous input sample a_k-2.

4) From the diagram,impulse response of the modified duobinary encoder is computed as

Decoding Process:

5) The receiver consists of a modified duobinary decoder and a postcoder (inverse of precoder). The decoder implements the following equation (which can be deduced from the equation given under step 3 (see above))

This equation indicates that the decoding process is prone to error propagation as the estimate of present sample relies on the estimate of previous sample. This error propagation is avoided by using a precoder before modified-duobinary encoder at the transmitter and a postcoder after the modified-duobinary decoder. The precoder ties the present sample and the sample that precedes the previous sample ( correlates these two samples) and the postcoder does the reverse process.

6) The entire process of modified-duobinary decoding and the postcoding can be combined together as one algorithm. The following decision rule is used for detecting the original modified-duobinary signal samples {a_n} from {y_n}

Matlab code:

Check this book for full Matlab code.
Wireless Communication Systems in Matlab – by Mathuranathan Viswanathan

Simulation Results:

To know more on the simulation and results – visit this page – “Partial response signalling schemes – impulse and frequency responses”

Impulse response and frequency response of various Partial response (PR) signaling schemes

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

Recommended Books

Correlative coding – Duobinary Signaling

The condition for zero ISI (Inter Symbol Interference) is

which states that when sampling a particular symbol (at time instant nT=0), the effect of all other symbols on the current sampled symbol is zero.

As discussed in the previous article, one of the practical ways to mitigate ISI is to use partial response signaling technique ( otherwise called as “correlative coding”). In partial response signaling, the requirement of zero ISI condition is relaxed as a controlled amount of ISI is introduced in the transmitted signal and is counteracted in the receiver side.

By relaxing the zero ISI condition, the above equation can be modified as,

which states that the ISI is limited to two adjacent samples. Here we introduce a controlled or “deterministic” amount of ISI and hence its effect can be removed upon signal detection at the receiver.

Duobinary Signaling:

The following figure shows the duobinary signaling scheme.

Encoding Process:

1) a_n = binary input bit; a_n ∈ {0,1}.
2) b_n = NRZ polar output of Level converter in the precoder and is given by,

3) y_n can be represented as

Note that the samples b_n are uncorrelated ( i.e either +d for “1” or -d for “0” input). On the other-hand, the samples y_n are correlated ( i.e. there are three possible values +2d,0,-2d depending on a_n and a_n-1). Meaning that the duobinary encoding correlates present sample a_n and the previous input sample a_n-1_.

4) From the diagram,impulse response of the duobinary encoder is computed as

Decoding Process:

5) The receiver consists of a duobinary decoder and a postcoder (inverse of precoder).The duobinary decoder implements the following equation (which can be deduced from the equation given under step 3 (see above))

This equation indicates that the decoding process is prone to error propagation as the estimate of present sample relies on the estimate of previous sample. This error propagation is avoided by using a precoder before duobinary encoder at the transmitter and a postcoder after the duobinary decoder. The precoder ties the present sample and previous sample ( correlates these two samples) and the postcoder does the reverse process.

6) The entire process of duobinary decoding and the postcoding can be combined together as one algorithm. The following decision rule is used for detecting the original duobinary signal samples {a_n} from {y_n}

Matlab Code:

Check this book for full Matlab code and simulation results.
Wireless Communication Systems in Matlab – by Mathuranathan Viswanathan

Simulation and results

To know more on the simulation and results – visit this page – “Partial response signalling schemes – impulse and frequency responses”

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

External Resources:

[1] The care and feeding of digital, pulse-shaping filters – By Ken Gentile↗
[2] Inter Symbol Interference and Root Raised Cosine Filtering – Complex2real↗

Understand Moving Average Filter with Python & Matlab

The moving average filter is a simple Low Pass FIR (Finite Impulse Response) filter commonly used for smoothing an array of sampled data/signal. It takes samples of input at a time and takes the average of those -samples and produces a single output point. It is a very simple LPF (Low Pass Filter) structure that comes handy for scientists and engineers to filter unwanted noisy component from the intended data.

As the filter length increases (the parameter ) the smoothness of the output increases, whereas the sharp transitions in the data are made increasingly blunt. This implies that this filter has excellent time domain response but a poor frequency response.

The MA filter performs three important functions:

1) It takes input points, computes the average of those -points and produces a single output point
2) Due to the computation/calculations involved, the filter introduces a definite amount of delay
3) The filter acts as a Low Pass Filter (with poor frequency domain response and a good time domain response).

Implementation

The difference equation for a -point discrete-time moving average filter with input represented by the vector and the averaged output vector , is

\[y[n] = \displaystyle{\frac{1}{L} \sum_{k=0}^{L-1}x[n-k]} \quad\quad (1) \]

For example, a -point Moving Average FIR filter takes the current and previous four samples of input and calculates the average. This operation is represented as shown in the Figure 1 with the following difference equation for the input output relationship in discrete-time.

\[\begin{aligned} y[n] &= \displaystyle{\frac{1}{5} \left(x[n] + x[n-1] + x[n-2] + x[n-3] + x[n-4] \right) } \\ &= 0.2 \left(x[n] + x[n-1] + x[n-2] + x[n-3] + x[n-4] \right) \end{aligned} \quad\quad (2) \]

Figure 1: Discrete-time 5-point Moving Average FIR filter

The unit delay shown in Figure 1 is realized by either of the two options:

Representing the input samples as an array in the computer memory and processing them
Using D-Flip flop shift registers for digital hardware implementation. If each discrete value of the input is represented as a -bit signal line from ADC (analog to digital converter), then we would require 4 sets of 12-Flip flops to implement the -point moving average filter shown in Figure 1.

Z-Transform and Transfer function

In signal processing, delaying a signal by sample period (unit delay) is equivalent to multiplying the Z-transform by . By applying this idea, we can find the Z-transform of the -point moving average filter in equation (2) as

\[ \begin{aligned} y[n] &= 0.2 \left(x[n] + x[n-1] + x[n-2] + x[n-3] + x[n-4] \right) \\ Y[z] &= 0.2 \left(z^0 + z^{-1} + z^{-2} + z^{-3} + z^{-4} \right) X[z] \\ Y[z] &= 0.2 \left(1 + z^{-1} + z^{-2} + z^{-3} + z^{-4} \right) X[z] \quad\quad (3) \end{aligned}\]

Similarly, the Z-transform of the generic -sample Moving Average filter of equation (1) is

\[x Y(z) = \left(\displaystyle{\frac{1}{L} \sum_{k=0}^{L-1} z^{-k}} \right) X(z) \quad\quad (4) \]

The transfer function describes the input-output relationship of the system and for the -point Moving Average filter, the transfer function is given by

\[H(z) = \displaystyle{\frac{Y(z)}{X(z)}} =\displaystyle{\frac{1}{L} \sum_{k=0}^{L-1} z^{-k}} \quad\quad (5)\]

Simulating the filter in Matlab and Python

In Matlab, we can use the filter function or conv (convolution) to implement the moving average FIR filter.

In general, the Z-transform of a discrete-time filter’s output is related to the Z-transform of the input by

\[H(z) = \displaystyle{\frac{Y(z)}{X(z)}} = \displaystyle{\frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots + b_n z^{-n}}{a_0 + a_1 z^{-1} + a_2 z^{-2} + \cdots + a_m z^{-m}}} \quad\quad (6) \]

where, and are the filter coefficients and the order of the filter is the maximum of and

For implementing equation (6) using the filter function, the Matlab function is called as

B = [b_0, b_1, b_2, ..., b_n] %numerator coefficients
A = [a_0, a_1, a_2, ..., a_m] %denominator coefficients
y = filter(B,A,x) %filter input x and get result in y

We can note from the difference equation and transfer function of the -point moving average filter, that following values for the numerator coefficients and denominator coefficients .

\[\begin{aligned} \{b_i\} &= \{1/L, 1/L, …., 1/L\} &,i=0,1,\cdots,L-1 \\ \{a_i\} &= {1} &,i=0 \end{aligned} \quad\quad (7)\]

Therefore, the -point moving average filter can be coded as

B = [0.2, 0.2, 0.2, 0.2, 0.2] %numerator coefficients
A = [1] %denominator coefficients
y = filter(B,A,x) %filter input x and get result in y

The numerator coefficients for the moving average filter can be conveniently expressed in short notion as shown below

 L = 5
B = ones(1,L)/L %numerator coefficients
A = [1] %denominator coefficients
x = rand(1,10) %random samples for x
y = filter(B,A,x) %filter input x and get result in y

When using conv function to implement the moving average filter, the following code can be used

L = 5;
x = rand(1,10) %random samples for x;
y = conv(x,ones(1,L)/L)

There exists a difference between using conv function and filter function for implementing an FIR filter. The conv function gives the result of complete convolution and the length of the result is length(x)+ L -1. Whereas, the filter function gives the output that is of same length as that of the input .

In python, the filtering operation can be performed using the lfilter and convolve functions available in the scipy signal processing package. The equivalent python code is shown below.

import numpy as np
from scipy import signal
L=5 #L-point filter
b = (np.ones(L))/L #numerator co-effs of filter transfer function
a = np.ones(1)  #denominator co-effs of filter transfer function
x = np.random.randn(10) #10 random samples for x
y = signal.convolve(x,b) #filter output using convolution

y = signal.lfilter(b,a,x) #filter output using lfilter function

Pole-zero plot and frequency response

A pole-zero plot for a filter transfer function , displays the pole and zero locations in the z-plane. In the pole-zero plot, the zeros occur at locations (frequencies) where and the poles occur at locations (frequencies) where . Equivalently, zeros occurs at frequencies for which the numerator of the transfer function in equation 6 becomes zero and the poles occurs at frequencies for which the denominator of the transfer function becomes zero.

In a pole-zero plot, the locations of the poles are usually marked by cross () and the zeros are marked as circles (). The poles and zeros of a transfer function effectively define the system response and determines the stability and performance of the filtering system.

In Matlab, the pole-zero plot and the frequency response of the -point moving average can be obtained as follows.

L=11
zplane([ones(1,L)]/L,1); %pole-zero plot
w = -pi:(pi/100):pi; %to plot frequency response
freqz([ones(1,L)]/L,1,w); %plot magnitude and phase response

Figure 2: Pole-Zero plot for L=11-point Moving Average filter

Figure 3: Magnitude and phase response of L=11-point Moving Average filter

The magnitude and phase frequency responses can be coded in Python as follows

[updated] See my Google colab for full Python code

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

L=11 #L-point filter
b = (np.ones(L))/L #numerator co-effs of filter transfer function
a = np.ones(1)  #denominator co-effs of filter transfer function

w, h = signal.freqz(b,a)
plt.subplot(2, 1, 1)
plt.plot(w, 20 * np.log10(abs(h)))
plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [rad/sample]')

plt.subplot(2, 1, 2)
angles = np.unwrap(np.angle(h))
plt.plot(w, angles)
plt.ylabel('Angle (radians)')
plt.xlabel('Frequency [rad/sample]')
plt.show()

Figure 4: Impulse response, Pole-zero plot, magnitude and phase response of L=11 moving average filter

Case study:

Following figures depict the time domain & frequency domain responses of a -point Moving Average filter. A noisy square wave signal is driven through the filter and the time domain response is obtained.

*Figure 5: Processing a signal through Moving average filter of various lengths*

On the first plot, we have the noisy square wave signal that is going into the moving average filter. The input is noisy and our objective is to reduce the noise as much as possible. The next figure is the output response of a 3-point Moving Average filter. It can be deduced from the figure that the 3-point Moving Average filter has not done much in filtering out the noise. We increase the filter taps to 10-points and we can see that the noise in the output has reduced a lot, which is depicted in next figure.

With L=51 tap filter, though the noise is almost zero, the transitions are blunted out drastically (observe the slope on the either side of the signal and compare them with the ideal brick wall transitions in the input signal).

*Figure 6: Frequency response of moving average filters of various lengths*

From the frequency response of lower order filters (L=3, L=10), it can be asserted that the roll-off is very slow and the stop band attenuation is not good. Given this stop band attenuation, clearly, the moving average filter cannot separate one band of frequencies from another. As we increase the filter order to 51, the transitions are not preserved in time domain. Good performance in the time domain results in poor performance in the frequency domain, and vice versa. Compromise need for optimal filter design.

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

References:

[1] Filtering – OpenCourseWare – MIT.↗
[2] HC Chen et al.,”Moving Average Filter with its application to QRS detection”,IEEE Computers in Cardiology, 2003.↗

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

QPSK – Quadrature Phase Shift Keying

Quadrature Phase Shift Keying (QPSK) is a form of phase modulation technique, in which two information bits (combined as one symbol) are modulated at once, selecting one of the four possible carrier phase shift states.

*Figure 1: Waveform* simulation model for QPSK modulation

The QPSK signal within a symbol duration \(T_{sym}\) is defined as

\[s(t) = A \cdot cos \left[2 \pi f_c t + \theta_n \right], \quad \quad 0 \leq t \leq T_{sym},\; n=1,2,3,4 \quad \quad (1) \]

where the signal phase is given by

\[\theta_n = \left(2n – 1 \right) \frac{\pi}{4} \quad \quad (2)\]

Therefore, the four possible initial signal phases are \(\pi/4, 3 \pi/4, 5 \pi/4\) and \(7 \pi/4\) radians. Equation (1) can be re-written as

\[\begin{align} s(t) &= A \cdot cos \theta_n \cdot cos \left( 2 \pi f_c t\right) – A \cdot sin \theta_n \cdot sin \left( 2 \pi f_c t\right) \\ &= s_{ni} \phi_i(t) + s_{nq} \phi_q(t) \quad\quad \quad \quad \quad\quad \quad \quad \quad\quad \quad \quad \quad\quad \quad \quad (3) \end{align} \]

The above expression indicates the use of two orthonormal basis functions: \( \left\langle \phi_i(t),\phi_q(t)\right\rangle\) together with the inphase and quadrature signaling points: \( \left\langle s_{ni}, s_{nq}\right\rangle\). Therefore, on a two dimensional co-ordinate system with the axes set to \( \phi_i(t)\) and \(\phi_q(t)\), the QPSK signal is represented by four constellation points dictated by the vectors \(\left\langle s_{ni}, s_{nq}\right\rangle\) with \( n=1,2,3,4\).

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
All books available in ebook (PDF) and Paperback formats

The transmitter

The QPSK transmitter, shown in Figure 1, is implemented as a matlab function qpsk_mod. In this implementation, a splitter separates the odd and even bits from the generated information bits. Each stream of odd bits (quadrature arm) and even bits (in-phase arm) are converted to NRZ format in a parallel manner.

Refer Digital Modulations using Matlab : Build Simulation Models from Scratch for full Matlab code.
Refer Digital Modulations using Python for full Python code

File 1: qpsk_mod.m: QPSK modulator

function [s,t,I,Q] = qpsk_mod(a,fc,OF)
%Modulate an incoming binary stream using conventional QPSK
%a - input binary data stream (0's and 1's) to modulate
%fc - carrier frequency in Hertz
%OF - oversampling factor (multiples of fc) - at least 4 is better
%s - QPSK modulated signal with carrier
%t - time base for the carrier modulated signal
%I - baseband I channel waveform (no carrier)
%Q - baseband Q channel waveform (no carrier)
L = 2*OF;%samples in each symbol (QPSK has 2 bits in each symbol)
ak = 2*a-1; %NRZ encoding 0-> -1, 1->+1
I = ak(1:2:end);Q = ak(2:2:end);%even and odd bit streams
I=repmat(I,1,L).'; Q=repmat(Q,1,L).';%even/odd streams at 1/2Tb baud
I = I(:).'; Q = Q(:).'; %serialize
fs = OF*fc; %sampling frequency
t=0:1/fs:(length(I)-1)/fs; %time base
iChannel = I.*cos(2*pi*fc*t);qChannel = -Q.*sin(2*pi*fc*t);
s = iChannel + qChannel; %QPSK modulated baseband signal

The timing diagram for BPSK and QPSK modulation is shown in Figure 2. For BPSK modulation the symbol duration for each bit is same as bit duration, but for QPSK the symbol duration is twice the bit duration: \(T_{sym}=2T_b\). Therefore, if the QPSK symbols were transmitted at same rate as BPSK, it is clear that QPSK sends twice as much data as BPSK does. After oversampling and pulse shaping, it is intuitively clear that the signal on the I-arm and Q-arm are BPSK signals with symbol duration \(2T_b\). The signal on the in-phase arm is then multiplied by \(cos (2 \pi f_c t)\) and the signal on the quadrature arm is multiplied by \(-sin (2 \pi f_c t)\). QPSK modulated signal is obtained by adding the signal from both in-phase and quadrature arms.

Note: The oversampling rate for the simulation is chosen as \(L=2 f_s/f_c\), where \(f_c\) is the given carrier frequency and \(f_s\) is the sampling frequency satisfying Nyquist sampling theorem with respect to the carrier frequency (\(f_s \geq f_c\)). This configuration gives integral number of carrier cycles for one symbol duration.

*Figure 2: Timing diagram for BPSK and QPSK modulations*

The receiver

Due to its special relationship with BPSK, the QPSK receiver takes the simplest form as shown in Figure 3. In this implementation, the I-channel and Q-channel signals are individually demodulated in the same way as that of BPSK demodulation. After demodulation, the I-channel bits and Q-channel sequences are combined into a single sequence. The function qpsk_demod implements a QPSK demodulator as per Figure 3.

Read more about QPSK, implementation of their modulator and demodulator, performance simulation in these books:

*Figure 3: Waveform simulation model for QPSK demodulation*

Performance simulation over AWGN

The complete waveform simulation for the aforementioned QPSK modulation and demodulation is given next. The simulation involves, generating random message bits, modulating them using QPSK modulation, addition of AWGN channel noise corresponding to the given signal-to-noise ratio and demodulating the noisy signal using a coherent QPSK receiver. The waveforms at the various stages of the modulator are shown in the Figure 4.

*Figure 4: Simulated QPSK waveforms at the transmitter side*

The performance simulation for the QPSK transmitter-receiver combination was also coded in the code given above and the resulting bit-error rate performance curve will be same as that of conventional BPSK. A QPSK signal essentially combines two orthogonally modulated BPSK signals. Therefore, the resulting performance curves for QPSK – \(E_b/N_0\) Vs. bits-in-error – will be same as that of conventional BPSK.

QPSK variants

QPSK modulation has several variants, three such flavors among them are: Offset QPSK, π/4-QPSK and π/4-DQPSK.

Offset-QPSK

Offset-QPSK is essentially same as QPSK, except that the orthogonal carrier signals on the I-channel and the Q-channel are staggered (one of them is delayed in time). In OQPSK, the orthogonal components cannot change states at the same time, this is because the components change state only at the middle of the symbol periods (due to the half symbol offset in the Q-channel). This eliminates 180° phase shifts all together and the phase changes are limited to 0° or 90° every bit period.

Elimination of 180° phase shifts in OQPSK offers many advantages over QPSK. Unlike QPSK, the spectrum of OQPSK remains unchanged when band-limited [1]. Additionally, OQPSK performs better than QPSK when subjected to phase jitters [2]. Further improvements to OQPSK can be obtained if the phase transitions are avoided altogether – as evident from continuous modulation schemes like Minimum Shift Keying (MSK) technique.

π/4-QPSK and π/4-DQPSK

In π/4-QPSK, the signaling points of the modulated signals are chosen from two QPSK constellations that are just shifted π/4 radians (45°) with respect to each other. Switching between the two constellations every successive bit ensures that the phase changes are confined to odd multiples of 45°. Therefore, phase transitions of 90° and 180° are eliminated.

π/4-QPSK preserves the constant envelope property better than QPSK and OQPSK. Unlike QPSK and OQPSK schemes, π/4-QPSK can be differentially encoded, therefore enabling the use of both coherent and non-coherent demodulation techniques. Choice of non-coherent demodulation results in simpler receiver design. Differentially encoded π/4-QPSK is referred as π/4-DQPSK.

Read more about QPSK and its variants, implementation of their modulator and demodulator, performance simulation in these books:

Constellation diagram

The phase transition properties of the different variants of QPSK schemes, are easily investigated using constellation diagram. Refer this article that discusses the method to plot signal space constellations, for the various modulations used in the transmitter.

*Figure 5: Constellations plots for: (a) a = 0.3 RC-filtered QPSK, (b) α = 0.3 RC-filtered O-QPSK, (c) α = 0.3 RC-filtered π/4-DQPSK and (d) MSK*

Refer Digital Modulations using Matlab : Build Simulation Models from Scratch for full Matlab code.
Refer Digital Modulations using Python for full Python code

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

References

[1] S. A. Rhodes, “Effects of hardlimiting on bandlimited transmissions with conventional and offset QPSK modulation”, in Proc. Nat. TeIecommun. Conf., Houston, TX, 1972, PP. 20F/1-20F/7
[2] S. A. Rhodes, “Effect of noisy phase reference on coherent detection of offset QPSK signals”, IEEE Trans. Commun., vol. COM-22, PP. 1046-1055, Aug. 1974.↗

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

Topics in this chapter

Digital Modulators and Demodulators - Passband Simulation Models
● Introduction
● Binary Phase Shift Keying (BPSK)
□ BPSK transmitter
□ BPSK receiver
□ End-to-end simulation
● Coherent detection of Differentially Encoded BPSK (DEBPSK)
● Differential BPSK (D-BPSK)
□ Sub-optimum receiver for DBPSK
□ Optimum noncoherent receiver for DBPSK
● Quadrature Phase Shift Keying (QPSK)
□ QPSK transmitter
□ QPSK receiver
□ Performance simulation over AWGN
● Offset QPSK (O-QPSK)
● π/p=4-DQPSK
● Continuous Phase Modulation (CPM)
□ Motivation behind CPM
□ Continuous Phase Frequency Shift Keying (CPFSK) modulation
□ Minimum Shift Keying (MSK)
● Investigating phase transition properties
● Power Spectral Density (PSD) plots
● Gaussian Minimum Shift Keying (GMSK)
□ Pre-modulation Gaussian Low Pass Filter
□ Quadrature implementation of GMSK modulator
□ GMSK spectra
□ GMSK demodulator
□ Performance
● Frequency Shift Keying (FSK)
□ Binary-FSK (BFSK)
□ Orthogonality condition for non-coherent BFSK detection
□ Orthogonality condition for coherent BFSK
□ Modulator
□ Coherent Demodulator
□ Non-coherent Demodulator
□ Performance simulation
□ Power spectral density

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.

Method 1: Auto-correlation using xcorr function

Matlab

Python

Method 2: Auto-correlation using Convolution

Matlab

Python

Method 3: Autocorrelation using Toeplitz matrix

Matlab

Python

Method 4: Auto-correlation using FFT/IFFT

Matlab

Python

Construction the Auto-correlation Matrix

Matlab

Python

References:

Books by the author

Definition

Method 1 – Brute-Force Method:

Python

Matlab

Method 2 – Using Toeplitz Matrix:

Method 3: Using FFT:

Python

Matlab

Other Methods:

Standard Algorithms for Fast Convolution:

Matlab Code

Simulation output:

References:

Topics in this chapter

Books by the author

Polynomial functions:

Representing single variable polynomial functions:

Multiplication of polynomials and linear convolution:

Toeplitz Matrix and Convolution:

References:

Topics in this chapter

Books by the author

Prerequisite: Sampling theorem – baseband sampling

Intermediate Sampling or Under-Sampling

Example:

See also:

Analog to Digital conversion:

Nyquist-Shannon Sampling Theorem:

Baseband Sampling:

Aliasing and Anti-aliasing:

Effects of Sampling Rate:

Application : Up-conversion

Reference:

See also:

Modified Duobinary Signaling:

Encoding Process:

Decoding Process:

Matlab code:

Simulation Results:

See also :

Recommended Books

Duobinary Signaling:

Encoding Process:

Decoding Process:

Matlab Code:

Simulation and results

See also :

External Resources:

Implementation

Z-Transform and Transfer function

Simulating the filter in Matlab and Python

Pole-zero plot and frequency response

Case study:

References:

Similar topics

Books by the author

The transmitter

The receiver

Performance simulation over AWGN

QPSK variants

Offset-QPSK

π/4-QPSK and π/4-DQPSK

Constellation diagram