Compute signal power in Matlab

Calculating the energy and power of a signal was discussed in one of the previous posts. Here, we will verify the calculation of signal power using Discrete Fourier Transform (DFT) in Matlab. Check here to know more on the concept of power and energy.

The total power of a signal can be computed using the following equation

For other forms of equations: refer here.

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

Case Study:

Let x(t) be a sine wave of amplitude A and frequency fc represented by the following equation.

When represented in frequency domain, it will look like the one on the right side plot in the following figure. This is evident from the fact that the sinewave can be mathematically represented by applying Euler’s formula.

Taking the Fourier transform of x(t) to represent it in frequency domain,

When considering the amplitude part,the above decomposition gives two spikes of amplitude A/2 on either side of the frequency domain at fc and -fc.

Figure 1: Time domain representation and frequency spectrum of a pure sinewave

Squaring the amplitudes gives the magnitude of power of the individual spikes/frequency components. The power spectrum is plotted below.

Figure 2: Power spectrum of a pure sinewave

Thus if the pure sinewave is of amplitude A=1V and frequency=100Hz, the power spectrum will have two spikes of value at 100 Hz and -100 Hz frequencies. The total power will be

Let’s verify this through simulation.

Simulation and Verification

A sine wave of 100 Hz frequency and amplitude 1V is taken for the experiment.

A=1; %Amplitude of sine wave
Fc=100; %Frequency of sine wave
Fs=1000; %Sampling rate - oversampled by the rate of 10
Ts=1/Fs; %Sampling period
nCycles=200; %Number of cycles of the sinewave

subplot(2,1,1);
t=0:Ts:nCycles/Fc-Ts; %Time base
x=A*sin(2*pi*Fc*t); %Sinusoidal function
stem(t,x); %Plot command

A sinusoidal wave of 10 cycles is plotted here

Figure 3: Sine wave generated in Matlab

Matlab’s Norm function:

Matlab’s basic installation comes with “norm” function. The p-norm in Matlab is computed as

By default, the single argument norm function computed 2-norm given as

To compute the total power of the signal x[n] (as in equation (1) above), all we have to do is – compute norm(x), square it and divide by the length of the signal.

L=length(x);
P=(norm(x)^2)/L;
sprintf('Power of the Signal from Time domain %f',P);

The above given piece of code will result in the following output

Power of the Signal from Time domain 0.500000

Verifying the total Power by DFT : Frequency Domain

Here, the total power is verified by applying DFT on the sinusoidal sequence. The sinusoidal sequence x[n] is represented in frequency domain X[f] using Matlab’s FFT function. The power associated with each frequency point is computed as

Finally, the total power is calculated as the sum of all the points in the frequency domain representation.

X = fft(x);
Px=sum(X.*conj(X))/(L^2); %Compute power with proper scaling.
subplot(2,1,2)
% Plot single-sided amplitude spectrum.
stem(Px);
sprintf('Total Power of the Signal from DFT %f',P);

Result:

Total Power of the Signal from DFT 0.500000
Figure 4: Power spectrum of a pure sinewave simulated in Matlab

A word on Matlab’s FFT: Matlab’s FFT is optimized for faster performance if the transform length is a power of 2. The following snippet of code simply calls “fft” without the transform length. In this case, the window length and the transform length are the same. This is to simplify the calculation of power. You can re-write the call to the FFT routine with transform length set to next power of two which is greater than or equal to the window length (sequence length). Then the step to compute total power will be differing slightly in the denominator. But that will not improve resolution (Remember : zero padding to compute FFT will not improve resolution).

Also note that in the above simulation we are using a pure sinusoid. The entire sequence of sinusoid defined all the cycles completely. There are no discontinuities in the sequence. If you call FFT with the transform length set to next power of 2 (as given in Matlab manuals), it will pad additional zeros to the sequence and creates discontinuities in the FFT computation. This will lead to spectral leakage. FFT and spectral leakage is discussed here.

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

References:

[1] Sanjay Lall,”Norm and Vector spaces”,Information Systems Laboratory,Stanford University.↗

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)

(175 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_book_cover
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

Power and Energy of a Signal : Demystified

Key focus: Clearly understand the terms: power and energy of a signal, their mathematical definition, physical significance and computation in signal processing context.

Energy of a signal:

Defining the term “size”:

In signal processing, a signal is viewed as a function of time. The term “size of a signal” is used to represent “strength of the signal”. It is crucial to know the “size” of a signal used in a certain application. For example, we may be interested to know the amount of electricity needed to power a LCD monitor as opposed to a CRT monitor. Both of these applications are different and have different tolerances. Thus the amount of electricity driving these devices will also be different.

A given signal’s size can be measured in many ways. Given a mathematical function (or a signal equivalently), it seems that the area under the curve, described by the mathematical function, is a good measure of describing the size of a signal. A signal can have both positive and negative values. This may render areas that are negative. Due to this effect, it is possible that the computed values cancel each other totally or partially, rendering incorrect result. Thus the metric function of “area under the curve” is not suitable for defining the “size” of a signal. Now, we are left with two options : either 1) computation of the area under the absolute value of the function or 2) computation of the area under the square of the function. The second choice is favored due to its mathematical tractability and its similarity to Euclidean Norm which is used in signal detection techniques (Note: Euclidean norm – otherwise called L2 norm or 2-norm[1] – is often considered in signal detection techniques – on the assumption that it provides a reasonable measure of distance between two points on signal space. It is computed as Euclidean distance in detection theory).

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

Going by the second choice of viewing the “size” as the computation of the area under the square of the function, the energy of a continuous-time complex signal \(x(t)\) is defined as

\[E_x = \int_{-\infty}^{\infty} |x(t)|^2 dt\]

If the signal \(x(t)\) is real, the modulus operator in the above equation does not matter.

This is called “Energy” in signal processing terms. This is also a measure of signal strength. This definition can be applied to any signal (or a vector) irrespective of whether it possesses actual energy (a basic quantitative property as described by physics) or not. If the signal is associated with some physical energy, then the above definition gives the energy content in the signal. If the signal is an electrical signal, then the above definition gives the total energy of the signal (in Joules) dissipated over a 1 Ohm resistor.

Actual Energy – physical quantity:

To know the actual energy of the signal \(E\), one has to know the value of load \(Z\) the signal is driving and also the nature the electrical signal (voltage or current). For a voltage signal, the above equation has to be scaled by a factor of \(1/Z\).

\[E = \frac{E_x}{Z} = \frac{1}{Z} \int_{-\infty}^{\infty} |x(t)|^2 dt \]

For current signal, it has to be scaled by \(Z\).

\[E = ZE_x = Z \int_{-\infty}^{\infty} |x(t)|^2 dt\]

Here, \(Z\) is the impedance driven by the signal \(x(t)\) , \(E_x\) is the signal energy (signal processing term) and \(E\) is the Energy of the signal (physical quantity) driving the load \( Z\)

Energy in discrete domain:

In the discrete domain, the energy of the signal is given by

\[E_x = \displaystyle{ \sum_{n=-\infty}^{\infty} |x(n)|^2}\]

The energy is finite only if the above sum converges to a finite value. This implies that the signal is “squarely-summable“. Such a signal is called finite energy signal.

What if the given signal does not decay with respect to time (as in a continuous sine wave repeating its cycle infinitely) ? The energy will be infinite  and such a signal is “not squarely-summable” in other words. We need another measurable quantity to circumvent this problem.This leads us to the notion of “Power”

Power:

Power is defined as the amount of energy consumed per unit time. This quantity is useful if the energy of the signal goes to infinity or the signal is “not-squarely-summable”. For “non-squarely-summable” signals, the power calculated by taking the snapshot of the signal over a specific interval of time as follows

1) Take a snapshot of the signal over some finite time duration
2) Compute the energy of the signal \(E_x\)
3) Divide the energy by number of samples taken for computation \(N\)
4) Extend the limit of number of samples to infinity \(N\rightarrow \infty \). This gives the total power of the signal.

In discrete domain, the total power of the signal is given by

\[P_x = \lim_{N \rightarrow \infty } \frac{1}{2N+1} \displaystyle{\sum_{n=-N}^{n=+N} |x(n)|^2} \]

Following equations are different forms of the same computation found in many text books. The only difference is the number of samples taken for computation. The denominator changes according to the number of samples taken for computation.

\[\begin{align} P_x &= \lim_{N \rightarrow \infty } \frac{1}{2N} \displaystyle{\sum_{n=-N}^{n=N-1} |x(n)|^2} \\ P_x &= \lim_{N \rightarrow \infty } \frac{1}{N} \displaystyle{\sum_{n=0}^{n=N-1} |x(n)|^2} \\ P_x &= \lim_{N \rightarrow \infty } \frac{1}{N_1 – N_0 + 1} \displaystyle{\sum_{n=N_0}^{n=N_1} |x(n)|^2} \end{align} \]

Classification of Signals:

A signal can be classified based on its power or energy content. Signals having finite energy are energy signals. Power signals have finite and non-zero power.

Energy Signal :

A finite energy signal will have zero TOTAL power. Let’s investigate this statement in detail. When the energy is finite, the total power will be zero. Check out the denominator in the equation for calculating the total power. When the limit \(N\rightarrow \infty\), the energy dilutes to zero over the infinite duration and hence the total power becomes zero.

Power Signal:

Signals whose total power is finite and non-zero. The energy of the power signal will be infinite. Example: Periodic sequences like sinusoid. A sinusoidal signal has finite, non-zero power but infinite energy.

A signal cannot be both an energy signal and a power signal.

Neither an Energy signal nor a Power signal:

Signals can also be a cat on the wall – neither an energy signal nor a power signal. Consider a signal of increasing amplitude defined by,

\[x(n) = n\]

For such a signal, both the energy and power will be infinite. Thus, it cannot be classified either as an energy signal or as a power signal.

Calculation of power and verifying it through Matlab is discussed next…

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

References:

[1] Sanjay Lall,”Norm and Vector spaces”,Information Systems Laboratory,Stanford University.↗

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

Similar Topics

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

Breaking the Bandwidth Barrier – Li-Fi turns LED bulbs into High Speed Internet Hotspots

Note: There is a rating embedded within this post, please visit this post to rate it.
Researchers, working on Ultra-Parallel Visible Light Communications (UP-VLC) project, have demonstrated 3.5Gpbs free space data transmission via three separate micro-LEDs that emit red, blue and green colours.  In effect, combining the three colours (that makes up for white colour), researchers say that achieving data rates over 10Gbps is feasible.

The technology enabler here is called – “Li-Fi”, a term coined by Prof. Harold Haas of University of Edinburgh, uses LED lights as medium of data transmission.  Li-Fi, academically referred as “Visible Light Communication”, aims to use existing micro- LED light bulbs for both illumination and communication.

With the advent of this technology, soon you will see all the illuminated spots in offices, houses or any other place turning into internet hotspots streaming high quality data. It has become a hot topic and soon many big corporations will jump on to the Li-Fi bandwagon.

It all started with Prof. Hass’s research group demonstrating the proof-of-concept results that exploited the higher peak-to-average ratio (PAR) property of Orthogonal Frequency Division Multiplexing (OFDM) systems.  PAR, otherwise called “Crest Factor”, is a disadvantage in RF communications. The research group have used this property to turn commercially available LED bulbs to ultra-high speed wireless hotspots. Later, this work was featured in TIME magazine’s “50 Best inventions of the year 2011” and Nobel Laureate Prof. Hänsch listed this work in his book titled “100 ground-breaking ideas” that could shape the next century.

Using OFDM, micro-LED light bulbs are enabled to handing millions of light-intensity-changes per second. Essentially, the micro-LED bulbs are switched on and off depending on the input bits –‘0’ or ‘1’. The switching is done at such a high speed that it is undetectable by human eye. Thus both un-flickering-illumination and high-speed-communication are achievable under this condition.

Li-Fi carries the promise of breaking the bandwidth-crunch suffered by existing wireless systems. The demand for more bandwidth is skyrocketing day by day with the proliferation of wireless enabled devices. The crunch faced by Radio-Frequency (RF) systems, is also compounded by complex government regulations, needs of a growing base of customers and risking of the performance issues due to these challenges.

Li-Fi is essentially a wireless system. It uses the visible light region of the electromagnetic spectrum which is 10,000 times larger than the microwave region of the electromagnetic spectrum. This opens up the utilizable frequencies to the order of terahertz level. It will not interfere with existing devices and it can be used in areas where there is extensive RF noise or in places where RF frequencies are restricted (like in airplanes).

Watch the TED talk – “Wireless data from every light bulb”- where Prof Haas demonstrates the capability of Li-Fi technology system prototype (using a desk-lamp) that steams a HD movie in real-time.

Symbol Timing Recovery for QPSK (digital modulations)

The goal of timing recovery is to estimate and correct the sampling instants and phase at the receiver, such that it allows the receiver to decode the transmitted symbols reliably.

What is Symbol timing Recovery :

When transmitting data across a communication system, three things are important: frequency of transmission, phase information and the symbol rate.

In coherent detection/demodulation, both the transmitter and receiver posses the knowledge of exact symbol sampling timing and symbol phase (and/or symbol frequency). While everything is set at the transmitter, the receiver is at the mercy of recovery algorithms to regenerate these information from the incoming signal itself. If the transmission is a passband transmission, the carrier recovery algorithm also recovers the carrier frequency. For phase sensitive systems like BPSK, QPSK etc.., the carrier recovery algorithm recovers the symbol phase so that it is synchronous with the transmitted symbol.

The first part in such a receiver architecture of a MPSK transmitting system is multiplying the incoming signal with sine and cosine components of the carrier wave.

The sine and cosine components are generated using a carrier recovery block (Phase Lock Loop (PLL) or setting a local oscillator to track the variations).

Once the in-phase and quadrature signals are separated out properly, the next task is to match each symbol with the transmitted pulse shape such that the overall SNR of the system improves.

Implementing this in digital domain, the architecture described so far would look like this (Note: the subscript of the incoming signal has changed from analog domain to digital domain – i.e. to )

In the digital architecture above, the matched filter is implemented as a simple finite impulse response (FIR) filter whose impulse response is matched to that of the transmitter pulse shape. It helps the receiver in timing recovery and also it improves the overall SNR of the system by suppressing some amount of noise. The incoming signal up to the point before the matched filter, may have fluctuations in the amplitude. The matched filter also behaves like an averaging filter that smooths out the variations in the signal.

Note that in this digital version, the incoming signal is already a sampled signal. It has already passed through an analog to digital converter that sampled the signal at some sampling rate. From the symbol perspective, the symbols have to be sampled at optimum sampling instant to extract its content properly.

This requires a re-sampler, which resamples the averaged signal at the optimum sampling instant. If the original sampling instant is before or after the optimum sampling point, the timing recovery signal will help to re-sample (re-adjust sampling times) accordingly.

Let’s take a simple BPSK transmitter for illustration. This would be equivalent to any of the single arms (in-phase and quadrature phase arms) of a QPSK transmitter or receiver.

An alternate data pattern (symbols) – [+1,-1,+1,+1,\cdots,] is transmitted across the channel. Assume that each symbol occupies Tsym=8 sample time.

clear all; clc;
n=10; %Number of data symbols
Tsym=8; %Symbol time interms of sample time or oversampling rate equivalently
%data=2*(rand(n,1)<0.5)-1;
data=[1 -1 1 -1 1 -1 1 -1 1 -1]'; %BPSK data
bpsk=reshape(repmat(data,1,Tsym)',n*Tsym,1); %BPSK signal

figure('Color',[1 1 1]);
subplot(3,1,1);
plot(bpsk);
title('Ideal BPSK symbols');
xlabel('Sample index [n]');
ylabel('Amplitude')
set(gca,'XTick',0:8:80);
axis([1 80 -2 2]); grid on;

Lets add white gaussian noise (awgn). A random noise of standard deviation 0.25 is generated and added with the generated BPSK symbols.

noise=0.25*randn(size(bpsk)); %Adding some amount of noise
received=bpsk+noise; %Received signal with noise

subplot(3,1,2);
plot(received);
title('Transmitted BPSK symbols (with noise)');
xlabel('Sample index [n]');
ylabel('Amplitude')
set(gca,'XTick',0:8:80);
axis([1 80 -2 2]); grid on;

From the first plot, we see that the transmitted pulse is a rectangular pulse that spans Tsym samples. In the illustration, Tsym=8. The best averaging filter (matched filter) for this case is a rectangular filter (though they are not preferred in practice, I am just using it for simplifying the illustration) that spans 8 samples. Such a rectangular pulse can be mathematically represented in terms of unit step function as

(Another type of averaging filter – “Moving Average Filter” is implemented here)

The resulting rectangular pulse will have a value of 0.5 at the edges of the sampling instants (index 0 and 7) and a value of ‘1’ at the remaining indices in between the edges. Such a rectangular function is indicated below.

The incoming signal is convolved with the averaging filter and the resultant output is given below

impRes=[0.5 ones(1,6) 0.5]; %Averaging Filter -> u[n]-u[n-Tsamp]
yy=conv(received,impRes,'full');
subplot(3,1,3);
plot(yy);
title('Matched Filter (Averaging Filter) output');
xlabel('Sample index [n]');
ylabel('Amplitude');

set(gca,'XTick',0:8:80);
axis([1 80 -10 10]); grid on;

We can note that the averaged output peaks at the locations where the symbol transition occurs. Thus, when the signal is sampled at those ideal locations, the BPSK symbols [+1,-1,+1, …] can be recovered perfectly.

In practice, a square root raised cosine (SRRC) filter is used both at the transmitter and the receiver (as a matched filter) to mitigate inter-symbol interference. An implementation of SRRC filter in Matlab is given here

But the problem here is: “How does the receiver know the ideal sampling instants?”. The solution is “someone has to supply those ideal sampling instants”. A symbol time recovery circuit is used for this purpose.

Coming back to the receiver architecture, lets add a symbol time recovery circuit that supplies the recovered timing instants. The signal will be re-sampled at those instants supplied by the recovery circuit.

The Algorithm behind Symbol Timing Recovery:

Different algorithms exist for symbol timing recovery and synchronization. An “Early/Late Symbol Recovery algorithm” is illustrated here.

The algorithm starts by selecting an arbitrary sample at some time (denoted by T). It captures the two adjacent samples (on either side of the sampling instant T) that are separated by δ seconds. The sample at the index T-δ is called Early Sample and the sample at the index T+δ is called Late Sample. The timing error is generated by comparing the amplitudes of the early and late samples. The next symbol sampling time instant is either advanced or delayed based on the sign of difference between the early and late sample.

1) If the Early Sample = Late Sample : The peak occurs at the on-time sampling instant T. No adjustment in the timing is needed.
2) If |Early Sample| > |Late Sample| : Late timing, the sampling time is offset so that the next symbol is sampled T-δ/2 seconds after the current sampling time.
3) If |Early Sample| < |Late Sample| : Early timing, the sampling time is offset so that the next symbol is sampled T+δ/2 seconds after the current sampling time.

These three situations are shown next.

There exist many variations to the above mentioned algorithm. The Early/Late synchronization technique given here is the simplest one taken for illustration.

Let’s complete the architecture with a signal quantization and constellation de-mapping block which gives out the estimated demodulated symbols.

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

For Further Reading:

[1] Technique for implementing an Early-Late Gate Synchronization structure for DPSK.↗
[2] Ying Li et al,”Hardware Implementation of Symbol Synchronization for Underwater FSK”, IEEE International Conference on Sensor Networks, Ubiquitous, and Trustworthy Computing (SUTC), p 82 – 88, 2010.↗
[3] Heinrich Meyr & Gerd Ascheid,”Digital Communication Receivers: Synchronization in Digital Communication Volume I, Phase-, Frequency-Locked Loops, and Amplitude Control (Wiley and Signal Processing)”,John Wiley & Sons; Volume 1 edition (March 1990),ISBN-13: 978-0471501930.↗
[4] Umberto Mengali,”Synchronization Techniques for Digital Receivers (Applications of Communications Theory)”,Springer; 1997 edition (October 31, 1997),ISBN-13: 978-0306457258.↗

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

Understanding Fourier Series

Understand Fourier Series, Fourier Cosine Series, Fourier Sine Series, partial sums, even odd symmetry. Hands-on simulation with Matlab code given.

Fourier analysis and Fourier Synthesis:

Fourier analysis – a term named after the French mathematician Joseph Fourier, is the process of breaking down a complex function and expressing it as a combination of simpler functions. The reverse process of combining simpler functions to reconstruct the complex function is termed as Fourier Synthesis.

Mostly, the simpler functions are chosen to be sine and cosine functions. Thus, the term “Fourier analysis” expresses a complex function in terms of sine and cosine terms and the term “Fourier Analysis” reconstructs the complex function from the sine and cosine terms.

Frequency is the measure of number of repetitive occurrences of a particular event. By definition, a sine wave is a smooth curve that repeats at a certain frequency. Thus, the term “frequency” and sine are almost synonymous. A cosine wave is also a sine wave but with 90* phase shift. Therefore, when you talk about sine and cosine functions, you are taking in terms of “frequencies”. That is why in signal processing, the Fourier analysis is applied in frequency (or spectrum) analysis.

Fourier series, Continuous Fourier Transform, Discrete Fourier Transform, and Discrete Time Fourier Transform are some of the variants of Fourier analysis.

Fourier series:

Applied on functions that are periodic. A periodic function is broken down and expressed in terms of sine and cosine terms. In mathematics, the term “series” represents a sum of sequence of numbers. For example we can make a series with a sequence of numbers that follows Geometric Progression (common ratio between the numbers)

Common ratio =3 : 1+ 3 + 9 + 27 + …

An infinite series is a series that has infinite number of terms. If the elements of the infinite series has a common ratio less than 1, then there is a possibility of the sum converging at a particular value. Fourier series falls under the category of trigonometric infinite series, where the individual elements of the series are expressed trigonometrically. The construct of the Fourier series is given by

Here f(x) is the complex periodic function we wish to break down in terms of sine and cosine basis functions. The coefficients a0, a1,… and b1, b2,… can be found by

Functions and Symmetry:

It is necessary to classify the functions according to its symmetry properties. Doing so will save computation time and effort. Functions either fall into odd symmetry or even symmetry or no symmetry category. Symmetry can be ascertained by plotting the function in a graph paper and folding it along the y axis. Symmetry of a function is always with respect to y axis.

Even Symmetry:

Figure 1: A function exhibiting even-symmetry

Mathematically depicted as f(x) = f(-x). The value of the given function f(x) at a given positive value x is same at corresponding negative value –x.  If plotted on a graph paper and folded along the y-axis, the left half and the right half of the function matches with each other (mirror image).

For even symmetry functions, only the cosine terms exist in Fourier Series expansion. The bn coefficients vanishes all-together (i.e, no sine basis). This leads to what is called Fourier Cosine Series.

Odd Symmetry:

Mathematically depicted as f(x) = -f(-x). The value of the given function f(x) at a given positive value x is same but with a sign change at corresponding negative value –x.  If plotted on a graph paper and folded along the y-axis, the left half of the graph will look like inverted (upside down) mirror image of the right half.

Figure 2: A function exhibiting odd-symmetry


For odd symmetry functions, only the sine terms exists in Fourier Series expansion. The an coefficients vanishes all-together (no cosine basis). This leads to what is called Fourier Sine Series.

Thus, knowing the symmetry could save us a lots of computation time and effort, as we do not have to calculate half the number of coefficients if symmetry exists.

Table: Fourier Series and Function Symmetry

Partial Sum and Convergence of Fourier Series:

Fourier Series is a class of infinite series, meaning that there are infinite terms in the expansion.We cannot go on calculating the terms indefinitely. To decompose a complex function using Fourier Series expansion, one has to limit the number of terms we wish to obtain and this process affects convergence. Convergence is based on certain criteria. There exists a separate branch of mathematics called Classical Harmonic Analysis that deals with this subject. Convergence is usually calculated over a partial sum – the sum of all terms upto which we have calculated the coefficients.

Example:

Consider the following periodic function :

Investigation of the function plot reveals that this function exhibits anti-symmetry (odd symmetry). So it is enough if we compute only the sine terms in the Fourier expansion.

Figure 3: Example for odd-symmetry

Computing the Fourier Sine Series:

Thus the complete Fourier expansion of the given function f(x) is given by

Note that the sine term vanishes when n is even (n=0,2,4,…,). Thus the above expansion can be simplified to

Matlab Simulation:

The following Matlab simulation computes the Fourier series expansion of the above mentioned function. The partial sum is plotted till an error criteria is satisfied.

The function f(x) either stays at +1 or at -1. The partial sum is calculated for each iteration and compared with either +1 or -1 and till the error reaches a small value of 0.01.

%Author Mathuranathan Viswanathan for https://gaussianwaves.com
%Creative Commons CC-BY-NC-SA
%If you use this piece of code you must attribute the author

clearvars;clc;
time=linspace(-pi,pi,1000);
partial_sum=0;

%Complex Function represented in terms of time and amplitude value
t=[-pi,-pi,0,0,pi,pi];
value=[0,-1,-1,1,1,-1];
handle1=line(t,value,'color','r','linewidth',2);
grid on;hold on; 
axis([-pi pi -1.5 1.5])
%Since the given complex function exhibits odd periodic extension
%only Bn term is valid with n=1,3,5,...

for n=1:2:200 %Odd terms to consider for partial sums
    %Plot 1 period of the given function    
    partial_sum=partial_sum+(4/(n*pi))*sin(n*time); %Fourier Series Expansion using Sine terms
    error=mean((abs(partial_sum)-1).^2); %Error Criteria
    handle2=plot(time,partial_sum,'k','linewidth',2);
    title(['Square Wave Partial Sum:  n = ',num2str(n),'  Error = ',num2str(error)])
    pause
    set(handle2,'Visible','off');
    if error<0.01
        break
    end
end

The plots below shows that the error gets minimized as more and more terms are included in the series expansion.

Figure 4: Simulated plots illustrating the role played by partial sums in Fourier Series expansion

Understanding the Plots:

In the first plot, the original square wave (red color) is decomposed into first three terms (n=3) of the Fourier Series. The plot in black color shows how the reconstructed (Fourier Synthesis) signal will look like if the three terms are combined together. As you progress further by increasing the number of terms ( n= 7, 15, 41, …) the plot in the black color increasingly resembles the original square wave.

Note the ringing effect at the corners of the black plot as the number of decomposed terms (n) is increased. This phenomenon is called Gibbs Phenomenon. Remember that the Fourier Series is an infinite series with indefinite number of terms. Since we cannot calculate all the infinite number of terms we have to stop at some point. This truncation of the number of decomposed terms leads to Gibbs Phenomenon. Read more on Gibss Phenomenon and its simulation in Matlab here.

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

For further reading

[1] Arthur L. Schoenstadt, “An Introduction to Fourier Analysis : Fourier Series, Partial Differential Equations and Fourier Transforms”, Department of Applied Mathematics, Naval Postgraduate School, Monterey, California, August 2005.↗

Similar articles

[1] Understanding Fourier Series
[2] Introduction to digital filter design
[3] Design FIR filter to reject unwanted frequencies
[4] FIR or IIR ? Understand the design perspective
[5] How to Interpret FFT results – complex DFT, frequency bins and FFTShift
[6] How to interpret FFT results – obtaining magnitude and phase information
[7] Analytic signal, Hilbert Transform and FFT
[8] FFT and spectral leakage
[9] Moving average filter in Python and Matlab

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

Folding Frequencies and Aliasing Zones

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.

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.

Bandpass sampling with pre-filtering

See also:

[1] Oversampling, ADC – DAC Conversion,pulse shaping and Matched Filter
[2] Sampling Theorem Basics and Baseband Sampling

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

Components of ADC

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…

Baseband Sampling

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…

Aliasing and Anti-aliasing

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.

Folding Frequencies and 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.

A complete design of ADC and DAC

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.

ADC and DAC chain

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.

Effects of Sampling Rate

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.

Application : Upconversion

Reference:

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

See also:

[1] Oversampling, ADC – DAC Conversion,pulse shaping and Matched Filter [2] Bandpass Sampling

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

Modified DuoBinary Signaling

Encoding Process:

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

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

3) yn can be represented as

Note that the samples bn are uncorrelated ( i.e either +d for “1” or -d for “0” input). On the other-hand,the samples yn are correlated ( i.e. there are three possible values +2d,0,-2d depending on ak and ak-2). Meaning that the modified duobinary encoding correlates present sample ak and the previous input sample ak-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 {an} from {yn}

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.

See also :

[1] Correlative Coding – Duobinary Signaling
[2] Introduction to Inter Symbol Interference

Recommended Books

FFT and Spectral Leakage

Spectral leakage due to FFT is caused by: mismatch between desired tone and chosen frequency resolution, time limiting an observation. Understand the concept using hands-on examples.

Limits of frequency domain studies

Frequency Transform is used to study a signal’s frequency domain characteristics. When using FFT to study the frequency domain characteristics of a signal, there are two limits : 1) The detect-ability of a small signal in the presence of a larger one; 2) frequency resolution – which distinguishes two different frequencies.

In practice, the measured signals are limited in time and the FFT calculates the frequency transform over a certain number of discrete frequencies called bins.

Spectral Leakage:

In reality, signals are of time-limited nature and nothing can be known about the signal beyond the measured interval. For example, if the measurement of a never ending continuous train of sinusoidal wave is of interest, at some point of time we need to terminate our observation to do further analysis. The limit on the time is also posed by limitations of the measurement system itself (example: buffer size), besides other factors.

It is often said that the FFT implicitly assumes that the signal essentially repeats itself after the measured interval and hence the FFT assumes the signal to be continuous (conceptually, juxtapose the measured signal repetitively). This lead to glitches in the assumed signal (see Figure 1). When the measurement time is purposefully made to be a non-integral multiple of the actual signal rate, these sharp discontinuities will spread out in the frequency domain leading to spectral leakage. This explanation for spectral leakage need to be carefully investigated.

Figure 1:Impact of observation interval on FFT

Experiment 1: Effect of FFT length and frequency resolution

Consider a pure sinusoidal signal of frequency \(f_x = 10 \;Hz\) and to represent in computer memory, the signal is observed for 1 second and sampled at frequency \(F_s=100 \; Hz\). Now, there will be 100 samples in the buffer and the buffer will contain integral number of waveform cycles (10 cycles in this case). The signal samples are analyzed using N-point DFT. Two cases are considered here for investigation : 1) The FFT size \(N\) is same as the length of the signal samples, i.e, N=100 and 2) FFT size set to next power of 2 that fits the signal samples i.e, N=128. The result are plotted next.

Fx=10; %Frequency of the sinusoid
Fs=100; %Sampling Frequency
observationTime = 1; %observation time in seconds
t=0:1/Fs:observationTime-1/Fs; %time base
x=sin(2*pi*Fx*t);%sampled sine wave

N=100; %DFT length same as signal length
X1 = 1/N*fftshift(fft(x,N));%N-point complex DFT of x
f1=(-N/2:1:N/2-1)*Fs/N; %frequencies on x-axis, Fs/N is the frequency resolution

N=128; %DFT length
X2 = 1/N*fftshift(fft(x,N));%N-point complex DFT of x
f2=(-N/2:1:N/2-1)*Fs/N; %frequencies on x-axis, Fs/N is the frequency resolution

figure;
subplot(3,1,1);stem(x,'r')
title('Time domain');xlabel('Sample index (n)');ylabel('x[n]');
subplot(3,1,2);stem(f1,abs(X1)); %magnitudes vs frequencies
xlim([-16,16]);title(['FFT, N=100, \Delta f=',num2str(Fs/100)]);xlabel('f (Hz)'); ylabel('|X(k)|');
subplot(3,1,3);stem(f2,abs(X2)); %magnitudes vs frequencies
xlim([-16,16]);title(['FFT, N=128, \Delta f=',num2str(Fs/128)]);xlabel('f (Hz)'); ylabel('|X(k)|');
Figure 2: Effect of FFT length and frequency resolution

One might wonder, even though the buffered samples contain integral number of waveform cycles, why the frequency spectrum registered a distinct spike at \(10 \; Hz\) when \(N=100\) and not when \(N=128\). This is due to the different frequency resolution, the measure of ability to resolve two different adjacent frequencies

For case 1, the frequency resolution is \(\Delta f = f_s/N = 100/10 = 1 Hz\). This means that the frequency bins are spaced 1 Hz apart and that is why it is able to hit the bull’s eye at 10 Hz peak.

For case 2, the frequency resolution is \(\Delta f = f_s/128 = 100/128 = 0.7813 Hz\). With this frequency resolution, the x-axis of the frequency plot cannot have exact value of 10 Hz. Instead, the nearest adjacent frequency bins are 9.375 Hz and 10.1563 Hz respectively. Therefore, the frequency spectrum cannot represent 10 Hz and the energy of the signal gets leaked to adjacent bins, leading to spectral leakage.

Experiment 2: Effect of time-limited observation

In the previous experiment, the signal wave observed for 1 second duration and that fetched whole 10 cycles in the signal buffer. Now, reduce the observation time to 0.91 second and re-run the same code, results below. In this case, the signal buffer will have 9.1 cycles of the sinewave, which is not a whole number. For case 1, the frequency resolution is 1 Hz and the FFT plot has registered a distinct peak at 10 Hz. Careful investigation of the plot will reveal very low spectral leakage even in case 1 (observe the non-zero amplitude values for the rest of the bins). This is primarily due to the change in the observation interval leading to non-integral number of cycles within the observed window. The spectral leakage in case 2, when N=128, is predominantly due to mismatch in the frequency resolution.

Figure 3: Effect of limited time observation on spectral leakage in time domain

From these two experiments, we can say that
1) The mismatch between the tone of the signal and the chosen frequency resolution (result of sampling frequency and the FFT length) leads to spectral leakage (experiment 1).
2) Time-limiting an observation (at inappropriate times), may lead to spectral leakage (experiment 2).
3) Hence, the spectral leakage from a larger signal component, if present, may significantly overshadow other smaller signals making them difficult to identify or detect.

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

\[\]

Reference

[1] Fredric J Harris, “On the use of windows for Harmonic Analysis with the Discrete Wavelet Transform”, Proceedings of IEEE, Vol 66,No 1, January 1978.↗

Similar articles

[1] Understanding Fourier Series
[2] Introduction to digital filter design
[3] Design FIR filter to reject unwanted frequencies
[4] FIR or IIR ? Understand the design perspective
[5] How to Interpret FFT results – complex DFT, frequency bins and FFTShift
[6] How to interpret FFT results – obtaining magnitude and phase information
[7] Analytic signal, Hilbert Transform and FFT
[8] FFT and spectral leakage
[9] Moving average filter in Python and Matlab
[10] Window function – figure of merits
[11] Equivalent noise bandwidth of window functions

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

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

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

Similar topics

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