FFT Archives - GaussianWaves

Window function – figure of merits

Key focus: Window function smooths the observed signal over the edges. Analysis of some important parameters to help select the window for an application.

Spectral leakage

As we know, the DFT operation can be viewed as processing a signal through a set of filter banks with bandwidth Δf centered on the bin (frequency) of interest (Figure 1). To accurately resolve the frequency component in each bin, we desire the resolution bandwidth Δf to be as small as possible. In effect, each filter filter bank is viewed as filter having a brick-wall response with narrow bandwidth. If the signal under study contains several spectral components and broadband noise, the computed amplitude of the DFT output in each bin (filter bank) is significantly affected by the noise contained in the corresponding bandwidth of the filter bank.

*Figure 1: DFT (Discrete Fourier Transform) viewed as processing through a set of filter banks*

*Figure 1: DFT (Discrete Fourier Transform) viewed as processing through a set of filter banks*

In reality, signals are of time-limited nature and nothing can be known about the signal beyond the measured interval. When the time-limited signal slice is subjected to analysis using FFT algorithm (used for computing DFT), the FFT implicitly assumes that the signal essentially repeats itself after the observed interval. This may lead to discontinuities at the edges of each slice, that causes the energy contained in each frequency bin to spill into other bins. This phenomenon is called spectral leakage.

Equivalent Noise Bandwidth:^[2]

To suppress the spectral leakage, the signals are multiplied with a window function so as to smooth the discontinuity at the edges of the FFT slices. As a result, the choice of window function affects the amount of signal and noise that goes inside each filter bank. Hence the amount of noise that gets accumulated into each filter bank is affected by the bandwidth of the chosen window function. In effect, the amplitude estimate in each frequency bin is influenced by the amount of accumulated noise contributed by the window function of choice. In other words, the noise floor displayed at the FFT output, varies according to the window of choice.

Equivalent noise bandwidth (ENBW), often specified in terms of FFT bins, is the bandwidth of a fictitious brick-wall filter such that the amount of noise power accumulated inside the brick-wall filter is same as the noise power accumulated when processing white noise through the selected window.

*Figure 2: Illustrating equivalent noise bandwidth (ENBW) (figure not to scale*)

In discrete domain, the equivalent noise bandwidth B_enbw can be computed using time domain samples as

\[B_{enbw} = L \frac{\sum_{n}|w[n]|^2}{ \left| \sum_{n} w[n] \right|^2}\]

where L is the length of the window.

ENBW is a commonly used metric to characterize the variation in the displayed noise floor if a non rectangular window is used. If noise floor measurements are desired, the ENBW correction factor need to be applied to account for the variation in the noise floor [3]. Table 1 below, lists the ENBW (bins) and ENBW correction factor for some well known window functions.

Following figure illustrates the equivalent noise bandwidth illustrated for Hann window.

*Figure 3: Equivalent noise bandwidth (ENBW) illustrated for Hann window*

Coherent power gain:

In time-domain, windowing a signal reduces the amplitude of the signal at its boundaries (Figure 4). This causes a reduction in the amplitude of the spectral component when FFT used to visualize the signal in frequency domain. This reduction in amplitude in the spectral components is characterized as coherent power gain (which actually is a loss). So, when FFT is computed on a windowed signal, in order to compensate for the loss in the amplitude, we simply add the coherent power gain to the FFT output.

*Figure 4: Windowing a 10 Hz sinusoidal signal with Hann window of length L = 151 (sampling frequency F_s=160 Hz)*

Coherent power gain of a window w[n] of length L is given by

\[\text{Coherent power gain (dB)} = 20 \; log_{10} \left( \frac{\sum_n w[n]}{L} \right)\]

Following plot depicts the coherent power gain (i.e, the reduction in the FFT magnitude when the input signal is processed with a window) of a windowed sinusoidal signal of frequency 10 Hz. The length of the window is L = 151 points and the simulation assumes an oversampling factor of 16 (i.e, F_s=160 Hz). The FFT length is N_FFT =2048.

Table 1 below, lists the coherent power gain for some well known window functions.

*Figure 5: Coherent power gain – illustrated for various window functions*

Scalloping loss

As discussed above, the FFT (equivalently the DFT) operation can be seen as processing a signal through a set of filter banks with bandwidth Δf centered on the successive bin frequencies. The frequency resolution (Δf) of FFT depends on the size of FFT (N_FFT) and the sampling frequency F_s.

\[\Delta f = \frac{F_s}{N_{FFT}}\]

Therefore, the FFT process can measure amplitude of tones that are multiples of the frequency resolution Δf. That is, if the frequency of the observed signal matches a particular bin frequency, the signal amplitude will be maximum at that bin.

However, if the tone of the signal falls between two bins, the signal power is spread between the adjacent bins and hence the displayed signal amplitude/power is reduced.

Scalloping loss quantifies the worst case reduction in the signal level, when the tone of the observed signal falls exactly mid way between two bin frequencies. The worst case scalloping loss is calculated as

\[\text{Scalloping loss} = \frac{\sum_n w[n] e^{\left( -j \frac{\pi}{N_{FFT}}n \right)}}{\sum_n w[n]}\]

In reality, scalloping loss depends on the type of window and the relationship between the signal tones to the bin frequencies. Table 1 below, lists the scalloping loss for some well known window functions.

*Figure 6: Scalloping loss illustrated for processing a signal with Boxcar and Bartlett window*

Figure of merits for window functions

*Table 1: Figure of merits for various window functions*

List of window functions

Boxcar (Rectangular)

\[w[n] = \begin{cases} 1, \quad 0 \leq n \leq L-1 \\ 0, \quad \text{otherwise} \end{cases}\]

Figure 1: Boxcarr (rectangular) window *(L=51)*

Barthann

\[ w[n] = 0.62 – 0.48 \left| \frac{n}{L – 1} – 0.5 \right| + 0.38 cos \left(2 \pi \left| \frac{n}{L – 1} – 0.5 \right| \right), \quad 0 \leq n \leq L-1\]

Bartlett

\[w[n] = \frac{2}{L-1} \left(\frac{L-1}{2} – \left|n – \frac{L-1}{2}\right| \right),\quad 0 \leq n \leq L-1\]

Blackman

\[w[n] = 0.42 – 0.5 \cos(2\pi n/L) + 0.08 \cos(4\pi n/L), \quad 0 \leq n \leq L-1\]

Blackman-Harris

\[w[n]=a_0 – a_1 \cos \left ( \frac{2 \pi n}{L} \right)+ a_2 \cos \left ( \frac{4 \pi n}{L} \right)- a_3 \cos \left ( \frac{6 \pi n}{L} \right), \quad 0 \leq n \leq L-1\] \[a_0=0.35875;\quad a_1=0.48829;\quad a_2=0.14128;\quad a_3=0.01168\]

Bohman

\[w[n] = \left( 1 – |x|\right) cos \left( \pi |x| \right) + \frac{1}{\pi} sin \left( \pi |x| \right) , \quad -1 \leq x \leq 1 \]

Cosine

\[w[n] = sin \left( \frac{\pi}{ L } \left(n + 0.5 \right) \right),\quad 0 \leq n \leq L-1\]

Flattop

\[w[n]=a_0 – a_1 \cos \left ( \frac{2 \pi n}{L} \right)+ a_2 \cos \left ( \frac{4 \pi n}{L} \right)- a_3 \cos \left ( \frac{6 \pi n}{L} \right), \quad 0 \leq n \leq L-1 \\ a_0=0.21557895;\; a_1=0.41663158 \\ a_2=0.277263158;\; a_3=0.006947368 \]

Hamming

\[w[n] = 0.54 – 0.46 \cos\left(\frac{2\pi{n}}{L-1}\right), \quad 0 \leq n \leq L-1\]

Hann

\[w[n] = 0.5 – 0.5 \cos\left(\frac{2\pi{n}}{L-1}\right), \quad 0 \leq n \leq L-1\]

Nuttall

\[ w[n]=a_0 – a_1 \cos \left ( \frac{2 \pi n}{L} \right)+ a_2 \cos \left ( \frac{4 \pi n}{L} \right)- a_3 \cos \left ( \frac{6 \pi n}{L} \right), \quad 0 \leq n \leq L-1 \\ a_0=0.3635819;\; a_1=0.4891775;\\ a_2=0.1365995;\; a_3=0.0106411\]

Parzen

\[w[n] = \begin{cases} 1- 6 \left( \frac{|n|}{L/2} \right)^2 + 6 \left( \frac{|n|}{L/2} \right)^3, & 0 \leq |n| \leq (L-1)/4 \\ 2 \left( 1 – \frac{|n|}{L/2} \right)^3, & (L-1)/4 < |n| \leq (L-1)/2 \end{cases}\]

Triangular

\[w[n] = 1 – \frac{|n|}{L/2}, \quad n = -L/2, \cdots, -1, 0, 1, \cdots, L/2\]

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

References:

[1] Harris, “On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform”, proceedings of IEEE, vol. 66, no. 1, January 1978 .↗
[2] Stefan Scholl, “Exact Signal Measurements using FFT Analysis”,Microelectronic Systems Design Research Group, TU Kaiserslautern, Germany.↗

Additional Resources:

[1] If you are particularly interested in calculating the Equivalent Noise Bandwidth of Analog Filters refer – Thomas J. Karras,”Equivalent Noise Bandwidth Analysis from Transfer Functions”,NASA Technical Note,NASA TN D-2842.↗
[2] Strategies for Choosing Windows : Michael Cerna and Audrey F. Harvey,”The Fundamentals of FFT-Based Signal Analysis and Measurement”,National Instruments Application Notes 041.↗

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

Equivalent noise bandwidth (ENBW) of window functions

Key focus: Equivalent noise bandwidth (ENBW), is the bandwidth of a fictitious brick-wall filter that allows same amount of noise as a window function. Learn how to calculate ENBW in applications involving window functions and FFT operation.

FFT and spectral leakage

Here is an article with illustrations on – FFT slices and the phenomenon of spectral leakage,

Hence, to suppress the spectral leakage, the signals are multiplied with a window function so as to smooth the discontinuity at the edges of the FFT slices. As a result, the choice of window function affects the amount of signal and noise that goes inside each filter bank. Hence the amount of noise that gets accumulated into each filter bank is affected by the bandwidth of the chosen window function. In effect, the amplitude estimate in each frequency bin is influenced by the amount of accumulated noise contributed by the window function of choice. In other words, the noise floor displayed at the FFT output, varies according to the window of choice.

Equivalent noise bandwidth

In other words, in frequency domain, the power accumulated in the chosen window is given by

\[P_{window} = \int_{-f}^{f} |W(f)|^2 df\]

Then we wish to find the equivalent noise bandwidth B_enbw, such that

\[P_{brick-wall} = B_{enbw} |W(f_0)|^2 = P_{window}\]

where f₀ is the frequency at which the power peaks and P_brick-wall is the power accumulated in the fictitious rectangular window of bandwidth B_enbw.

Therefore, the equivalent noise bandwidth B_enbw is given by

\[B_{enbw} = \frac{\int_{-f}^{f} |W(f)|^2 df}{|W(f_0)|^2}\]

Translating to discrete domain, the equivalent noise bandwidth can be computed using DFT samples as

\[B_{enbw} =\frac{\sum_{k=0}^{N-1}|W[k]|^2}{|W[k_0]|^2}\]

where, k₀ is the index at which the magnitude of FFT output is maximum and N is the window length. Applying Parseval’s theorem and with a window of length L, B_enbw can also be computed using time domain samples as

\[B_{enbw} = L \frac{\sum_{n}|w[n]|^2}{ \left| \sum_{n} w[n] \right|^2}\]

Following figure illustrates the equivalent noise bandwidth illustrated for Hann window.

Python script

ENBW can be calculated from the time domain samples in a straightforward manner. Following Python 3.0 script calculates the ENBW for some well-known window functions provided as part of Scipy module↗.

#Author : Mathuranathan Viswanathan for gaussianwaves.com
#Date: 13 Sept 2020
#Script to calculate equivalent noise bandwidth (ENBW) for some well known window functions

import numpy as np
import pandas as pd
from scipy import signal

def equivalent_noise_bandwidth(window):
    #Returns the Equivalent Noise BandWidth (ENBW)
    return len(window) * np.sum(window**2) / np.sum(window)**2

def get_enbw_windows():
    #Return ENBW for all the following windows as a dataframe
    window_names = ['boxcar','barthann','bartlett','blackman','blackmanharris','bohman','cosine','exponential','flattop','hamming','hann','nuttall','parzen','triang']
    
    df = pd.DataFrame(columns=['Window','ENBW (bins)','ENBW correction (dB)'])
    for window_name in window_names:
        method_name = window_name
        func_to_run = getattr(signal, method_name) #map window names to window functions in scipy package
        L = 16384 #Number of points in the output window
        window = func_to_run(L) #call the functions
        
        enbw = equivalent_noise_bandwidth(window) #compute ENBW
        
        df = df.append({'Window': window_name.title(),'ENBW (bins)':round(enbw,3),'ENBW correction (dB)': round(10*np.log10(enbw),3)},ignore_index=True)
                       
    return df

df = get_enbw_windows() #call the function
df.head(14) #display dataframe

The resulting output is displayed in the following table (dataframe)

*Table 1: Table of equivalent noise bandwidth for some well known window functions*

As an example, the following plot depicts the difference in the noise floor of FFT output of a noise added sine wave that is processed through Boxcar and Flattop window.

*Figure 4: Noise floor and ENBW for flattop & boxcar window (FFT output) for noise added 10 Hz sinewave (oversampling factor = 16, Fs = 160 Hz, window length L =2048, SNR = 30 dB)*

Continue reading on Window function and figure of merits…

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

Reference

[1] Stefan Scholl, “Exact Signal Measurements using FFT Analysis”,Microelectronic Systems Design Research Group, TU Kaiserslautern, Germany.↗

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

Plot FFT using Python – FFT of sine wave & cosine wave

Key focus: Learn how to plot FFT of sine wave and cosine wave using Python. Understand FFTshift. Plot one-sided, double-sided and normalized spectrum using FFT.

Introduction

Numerous texts are available to explain the basics of Discrete Fourier Transform and its very efficient implementation – Fast Fourier Transform (FFT). Often we are confronted with the need to generate simple, standard signals (sine, cosine, Gaussian pulse, squarewave, isolated rectangular pulse, exponential decay, chirp signal) for simulation purpose. I intend to show (in a series of articles) how these basic signals can be generated in Python and how to represent them in frequency domain using FFT. If you are inclined towards Matlab programming, visit here.

This article is part of the book Digital Modulations using Python, ISBN: 978-1712321638 available in ebook (PDF) and Paperback (hardcopy) formats

Sine Wave

In order to generate a sine wave, the first step is to fix the frequency f of the sine wave. For example, we wish to generate a sine wave whose minimum and maximum amplitudes are -1V and +1V respectively. Given the frequency of the sinewave, the next step is to determine the sampling rate.

For baseband signals, the sampling is straight forward. By Nyquist Shannon sampling theorem, for faithful reproduction of a continuous signal in discrete domain, one has to sample the signal at a rate higher than at-least twice the maximum frequency contained in the signal (actually, it is twice the one-sided bandwidth occupied by a real signal. For a baseband signal bandwidth ( to ) and maximum frequency in a given band are equivalent).

For Python implementation, let us write a function to generate a sinusoidal signal using the Python’s Numpy library. Numpy is a fundamental library for scientific computations in Python. In order to use the numpy package, it needs to be imported. Here, we are importing the numpy package and renaming it as a shorter alias np.

import numpy as np

Next, we define a function for generating a sine wave signal with the required parameters.

def sine_wave(f,overSampRate,phase,nCyl):
	"""
	Generate sine wave signal with the following parameters
	Parameters:
		f : frequency of sine wave in Hertz
		overSampRate : oversampling rate (integer)
		phase : desired phase shift in radians
		nCyl : number of cycles of sine wave to generate
	Returns:
		(t,g) : time base (t) and the signal g(t) as tuple
	Example:
		f=10; overSampRate=30;
		phase = 1/3*np.pi;nCyl = 5;
		(t,g) = sine_wave(f,overSampRate,phase,nCyl)
	"""
	fs = overSampRate*f # sampling frequency
	t = np.arange(0,nCyl*1/f-1/fs,1/fs) # time base
	g = np.sin(2*np.pi*f*t+phase) # replace with cos if a cosine wave is desired
	return (t,g) # return time base and signal g(t) as tuple

We note that the function sine wave is defined inside a file named signalgen.py. We will add more such similar functions in the same file. The intent is to hold all the related signal generation functions, in a single file. This approach can be extended to object oriented programming. Now that we have defined the sine wave function in signalgen.py, all we need to do is call it with required parameters and plot the output.

"""
Simulate a sinusoidal signal with given sampling rate
"""
import numpy as np
import matplotlib.pyplot as plt # library for plotting
from signalgen import sine_wave # import the function

f = 10 #frequency = 10 Hz
overSampRate = 30 #oversammpling rate
fs = f*overSampRate #sampling frequency
phase = 1/3*np.pi #phase shift in radians
nCyl = 5 # desired number of cycles of the sine wave

(t,x) = sine_wave(f,overSampRate,phase,nCyl) #function call

plt.plot(t,x) # plot using pyplot library from matplotlib package
plt.title('Sine wave f='+str(f)+' Hz') # plot title
plt.xlabel('Time (s)') # x-axis label
plt.ylabel('Amplitude') # y-axis label
plt.show() # display the figure

Python is an interpreter based software language that processes everything in digital. In order to obtain a smooth sine wave, the sampling rate must be far higher than the prescribed minimum required sampling rate, that is at least twice the frequency – as per Nyquist-Shannon theorem. Hence, we need to sample the input signal at a rate significantly higher than what the Nyquist criterion dictates. Higher oversampling rate requires more memory for signal storage. It is advisable to keep the oversampling factor to an acceptable value.

An oversampling factor of is chosen in the previous function. This is to plot a smooth continuous like sine wave. Thus, the sampling rate becomes . If a phase shift is desired for the sine wave, specify it too.

*Figure 1: A 10Hz sinusoidal wave with 5 cycles and phase shift 1/3π radians*

Different representations of FFT:

Since FFT is just a numeric computation of -point DFT, there are many ways to plot the result. The FFT, implemented in Scipy.fftpack package, is an algorithm published in 1965 by J.W.Cooley and
J.W.Tuckey for efficiently calculating the DFT.

The SciPy functions that implement the FFT and IFFT can be invoked as follows

from scipy.fftpack import fft, ifft
X = fft(x,N) #compute X[k]
x = ifft(X,N) #compute x[n]

1. Plotting raw values of DFT:

The x-axis runs from to – representing sample values. Since the DFT values are complex, the magnitude of the DFT is plotted on the y-axis. From this plot we cannot identify the frequency of the sinusoid that was generated.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

NFFT=1024 #NFFT-point DFT      
X=fft(x,NFFT) #compute DFT using FFT    

fig1, ax = plt.subplots(nrows=1, ncols=1) #create figure handle
nVals = np.arange(start = 0,stop = NFFT) # raw index for FFT plot
ax.plot(nVals,np.abs(X))      
ax.set_title('Double Sided FFT - without FFTShift')
ax.set_xlabel('Sample points (N-point DFT)')        
ax.set_ylabel('DFT Values')
fig1.show()

*Figure 2: Double sided FFT – without FFTShift*

2. FFT plot – plotting raw values against normalized frequency axis:

In the next version of plot, the frequency axis (x-axis) is normalized to unity. Just divide the sample index on the x-axis by the length of the FFT. This normalizes the x-axis with respect to the sampling rate . Still, we cannot figure out the frequency of the sinusoid from the plot.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

NFFT=1024 #NFFT-point DFT  
X=fft(x,NFFT) #compute DFT using FFT     

fig2, ax = plt.subplots(nrows=1, ncols=1) #create figure handle
   
nVals=np.arange(start = 0,stop = NFFT)/NFFT #Normalized DFT Sample points         
ax.plot(nVals,np.abs(X))     
ax.set_title('Double Sided FFT - without FFTShift')        
ax.set_xlabel('Normalized Frequency')
ax.set_ylabel('DFT Values')
fig2.show()

*Figure 3: Double sided FFT with normalized x-axis (0 to 1)*

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

As you know, in the frequency domain, the values take up both positive and negative frequency axis. In order to plot the DFT values on a frequency axis with both positive and negative values, the DFT value at sample index has to be centered at the middle of the array. This is done by using FFTshift function in Scipy Python. The x-axis runs from to where the end points are the normalized ‘folding frequencies’ with respect to the sampling rate .

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,fftshift

NFFT=1024 #NFFT-point DFT      
X=fftshift(fft(x,NFFT)) #compute DFT using FFT  

fig3, ax = plt.subplots(nrows=1, ncols=1) #create figure handle
    
fVals=np.arange(start = -NFFT/2,stop = NFFT/2)/NFFT #DFT Sample points        
ax.plot(fVals,np.abs(X))
ax.set_title('Double Sided FFT - with FFTShift')
ax.set_xlabel('Normalized Frequency')
ax.set_ylabel('DFT Values');
ax.autoscale(enable=True, axis='x', tight=True)
ax.set_xticks(np.arange(-0.5, 0.5+0.1,0.1))
fig.show()

*Figure 4: Double sided FFT with normalized x-axis (-0.5 to 0.5)*

4. FFT plot – Absolute frequency on the x-axis vs. magnitude on y-axis:

Here, the normalized frequency axis is just multiplied by the sampling rate. From the plot below we can ascertain that the absolute value of FFT peaks at and . Thus the frequency of the generated sinusoid is . The small side-lobes next to the peak values at and are due to spectral leakage.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,fftshift

NFFT=1024     
X=fftshift(fft(x,NFFT))

fig4, ax = plt.subplots(nrows=1, ncols=1) #create figure handle

fVals=np.arange(start = -NFFT/2,stop = NFFT/2)*fs/NFFT
ax.plot(fVals,np.abs(X),'b')
ax.set_title('Double Sided FFT - with FFTShift')
ax.set_xlabel('Frequency (Hz)')         
ax.set_ylabel('|DFT Values|')
ax.set_xlim(-50,50)
ax.set_xticks(np.arange(-50, 50+10,10))
fig4.show()

*Figure 5: Double sided FFT – Absolute frequency on the x-axis vs. magnitude on y-axis*

5. Power Spectrum – Absolute frequency on the x-axis vs. power on y-axis:

The following is the most important representation of FFT. It plots the power of each frequency component on the y-axis and the frequency on the x-axis. The power can be plotted in linear scale or in log scale. The power of each frequency component is calculated as

Where is the frequency domain representation of the signal . In Python, the power has to be calculated with proper scaling terms.

*Figure 6: Power spectral density using FFT*

Plotting the PSD plot with y-axis on log scale, produces the most encountered type of PSD plot in signal processing.

*Figure 7: Power spectral density (y-axis on log scale) using FFT*

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

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

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

Phase demodulation via Hilbert transform: Hands-on

Key focus: Demodulation of phase modulated signal by extracting instantaneous phase can be done using Hilbert transform. Hands-on demo in Python & Matlab.

Phase modulated signal:

The concept of instantaneous amplitude/phase/frequency are fundamental to information communication and appears in many signal processing application. We know that a monochromatic signal of form x(t) = a cos(ω t + ɸ) cannot carry any information. To carry information, the signal need to be modulated. Different types of modulations can be performed – amplitude modulation, phase modulation / frequency modulation.

In amplitude modulation, the information is encoded as variations in the amplitude of a carrier signal. Demodulation of an amplitude modulated signal, involves extraction of envelope of the modulated signal. This was discussed and demonstrated 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

In phase modulation, the information is encoded as variations in the phase of the carrier signal. In its generic form, a phase modulated signal is expressed as an information-bearing sinusoidal signal modulating another sinusoidal carrier signal

\[x(t) = A cos \left[ 2 \pi f_c t + \beta + \alpha sin \left( 2 \pi f_m t + \theta \right) \right] \quad \quad \quad (1)\]

where, m(t) = α sin (2 π f_m t + θ ) represents the information-bearing modulating signal, with the following parameters

α – amplitude of the modulating sinusoidal signal
f_m – frequency of the modulating sinusoidal signal
θ – phase offset of the modulating sinusoidal signal

The carrier signal has the following parameters

A – amplitude of the carrier
f_c – frequency of the carrier and f_c>>f_m
β – phase offset of the carrier

Demodulating a phase modulated signal:

The phase modulated signal shown in equation (1), can be simply expressed as

\[x(t) = A cos \left[ \phi(t) \right] \quad\quad\quad (2) \]

Here, ɸ(t) is the instantaneous phase that varies according to the information signal m(t).

A phase modulated signal of form x(t) can be demodulated by forming an analytic signal by applying Hilbert transform and then extracting the instantaneous phase. This method is explained here.

We note that the instantaneous phase is ɸ(t) = 2 π f_c t + β + α sin (2 π f_m t + θ) is linear in time, that is proportional to 2 π f_c t . This linear offset needs to be subtracted from the instantaneous phase to obtained the information bearing modulated signal. If the carrier frequency is known at the receiver, this can be done easily. If not, the carrier frequency term 2 π f_c t needs to be estimated using a linear fit of the unwrapped instantaneous phase. The following Matlab and Python codes demonstrate all these methods.

Matlab code

%Demonstrate simple Phase Demodulation using Hilbert transform
clearvars; clc;
fc = 240; %carrier frequency
fm = 10; %frequency of modulating signal
alpha = 1; %amplitude of modulating signal
theta = pi/4; %phase offset of modulating signal
beta = pi/5; %constant carrier phase offset 
receiverKnowsCarrier= 'False'; %If receiver knows the carrier frequency & phase offset

fs = 8*fc; %sampling frequency
duration = 0.5; %duration of the signal
t = 0:1/fs:1-1/fs; %time base

%Phase Modulation
m_t = alpha*sin(2*pi*fm*t + theta); %modulating signal
x = cos(2*pi*fc*t + beta + m_t ); %modulated signal

figure(); subplot(2,1,1)
plot(t,m_t) %plot modulating signal
title('Modulating signal'); xlabel('t'); ylabel('m(t)')

subplot(2,1,2)
plot(t,x) %plot modulated signal
title('Modulated signal'); xlabel('t');ylabel('x(t)')

%Add AWGN noise to the transmitted signal
nMean = 0; %noise mean
nSigma = 0.1; %noise sigma
n = nMean + nSigma*randn(size(t)); %awgn noise
r = x + n;  %noisy received signal

%Demodulation of the noisy Phase Modulated signal
z= hilbert(r); %form the analytical signal from the received vector
inst_phase = unwrap(angle(z)); %instaneous phase

%If receiver don't know the carrier, estimate the subtraction term
if strcmpi(receiverKnowsCarrier,'True')
    offsetTerm = 2*pi*fc*t+beta; %if carrier frequency & phase offset is known
else
    p = polyfit(t,inst_phase,1); %linearly fit the instaneous phase
    estimated = polyval(p,t); %re-evaluate the offset term using the fitted values
    offsetTerm = estimated;
end
    
demodulated = inst_phase - offsetTerm;

figure()
plot(t,demodulated); %demodulated signal
title('Demodulated signal'); xlabel('n'); ylabel('\hat{m(t)}');

Python code

import numpy as np
from scipy.signal import hilbert
import matplotlib.pyplot as plt
PI = np.pi

fc = 240 #carrier frequency
fm = 10 #frequency of modulating signal
alpha = 1 #amplitude of modulating signal
theta = PI/4 #phase offset of modulating signal
beta = PI/5 #constant carrier phase offset 
receiverKnowsCarrier= False; #If receiver knows the carrier frequency & phase offset

fs = 8*fc #sampling frequency
duration = 0.5 #duration of the signal
t = np.arange(int(fs*duration)) / fs #time base

#Phase Modulation
m_t = alpha*np.sin(2*PI*fm*t + theta) #modulating signal
x = np.cos(2*PI*fc*t + beta + m_t ) #modulated signal

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,m_t) #plot modulating signal
plt.title('Modulating signal')
plt.xlabel('t')
plt.ylabel('m(t)')
plt.subplot(2,1,2)
plt.plot(t,x) #plot modulated signal
plt.title('Modulated signal')
plt.xlabel('t')
plt.ylabel('x(t)')

#Add AWGN noise to the transmitted signal
nMean = 0 #noise mean
nSigma = 0.1 #noise sigma
n = np.random.normal(nMean, nSigma, len(t))
r = x + n  #noisy received signal

#Demodulation of the noisy Phase Modulated signal
z= hilbert(r) #form the analytical signal from the received vector
inst_phase = np.unwrap(np.angle(z))#instaneous phase

#If receiver don't know the carrier, estimate the subtraction term
if receiverKnowsCarrier:
    offsetTerm = 2*PI*fc*t+beta; #if carrier frequency & phase offset is known
else:
    p = np.poly1d(np.polyfit(t,inst_phase,1)) #linearly fit the instaneous phase
    estimated = p(t) #re-evaluate the offset term using the fitted values
    offsetTerm = estimated
                           
demodulated = inst_phase - offsetTerm 

plt.figure()
plt.plot(t,demodulated) #demodulated signal
plt.title('Demodulated signal')
plt.xlabel('n')
plt.ylabel('\hat{m(t)}')

Results

*Figure 1: Phase modulation – modulating signal and modulated (transmitted) signal*

*Figure 2: Demodulated signal from the noisy received signal*

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

For further reading

[1] V. Cizek, “Discrete Hilbert transform”, IEEE Transactions on Audio and Electroacoustics, Volume: 18 , Issue: 4 , December 1970.↗

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

Extract envelope, phase using Hilbert transform: Demo

Key focus: Learn how to use Hilbert transform to extract envelope, instantaneous phase and frequency from a modulated signal. Hands-on demo using Python & Matlab.

If you would like to brush-up the basics on analytic signal and how it related to Hilbert transform, you may visit article: Understanding Analytic Signal and Hilbert Transform

Introduction

The concept of instantaneous amplitude/phase/frequency are fundamental to information communication and appears in many signal processing application. We know that a monochromatic signal of form cannot carry any information. To carry information, the signal need to be modulated. Take for example the case of amplitude modulation, in which a positive real-valued signal modulates a carrier . That is, the amplitude modulation is effected by multiplying the information bearing signal with the carrier signal .

Here, is the angular frequency of the signal measured in radians/sec and is related to the temporal frequency as . The term is also called instantaneous amplitude.

Similarly, in the case of phase or frequency modulations, the concept of instantaneous phase or instantaneous frequency is required for describing the modulated signal.

Here, is the constant amplitude factor (no change in the envelope of the signal) and is the instantaneous phase which varies according to the information. The instantaneous angular frequency is expressed as the derivative of instantaneous phase.

Definition

Generalizing these concepts, if a signal is expressed as

The instantaneous amplitude or the envelope of the signal is given by
The instantaneous phase is given by
The instantaneous angular frequency is derived as
The instantaneous temporal frequency is derived as

Problem statement

An amplitude modulated signal is formed by multiplying a sinusoidal information and a linear frequency chirp. The information content is expressed as and the linear frequency chirp is made to vary from to . Given the modulated signal, extract the instantaneous amplitude (envelope), instantaneous phase and the instantaneous frequency.

Solution

We note that the modulated signal is a real-valued signal. We also take note of the fact that amplitude/phase and frequency can be easily computed if the signal is expressed in complex form. Which transform should we use such that the we can convert a real signal to the complex plane without altering the required properties ?? Answer: Apply Hilbert transform and form the analytic signal on the complex plane. Figure 1 illustrates this concept.

*Figure 1: Converting a real-valued signal to complex plane using Hilbert Transform*

If we express the real-valued modulated signal as an analytic signal, it is expressed in complex plane as

where, (HT[.]) represents the Hilbert Transform operation. Now, the required parameters are very easy to obtain.

The instantaneous amplitude (envelope extraction) is computed in the complex plane as

The instantaneous phase is computed in the complex plane as

The instantaneous temporal frequency is computed in the complex plane as

Once we know the instantaneous phase, the carrier can be regenerated as . The regenerated carrier is often referred as Temporal Fine Structure (TFS) in Acoustic signal processing [1].

Matlab

fs = 600; %sampling frequency in Hz
t = 0:1/fs:1-1/fs; %time base
a_t = 1.0 + 0.7 * sin(2.0*pi*3.0*t) ; %information signal
c_t = chirp(t,20,t(end),80); %chirp carrier
x = a_t .* c_t; %modulated signal

subplot(2,1,1); plot(x);hold on; %plot the modulated signal

z = hilbert(x); %form the analytical signal
inst_amplitude = abs(z); %envelope extraction
inst_phase = unwrap(angle(z));%inst phase
inst_freq = diff(inst_phase)/(2*pi)*fs;%inst frequency

%Regenerate the carrier from the instantaneous phase
regenerated_carrier = cos(inst_phase);

plot(inst_amplitude,'r'); %overlay the extracted envelope
title('Modulated signal and extracted envelope'); xlabel('n'); ylabel('x(t) and |z(t)|');
subplot(2,1,2); plot(cos(inst_phase));
title('Extracted carrier or TFS'); xlabel('n'); ylabel('cos[\omega(t)]');

Python

import numpy as np
from scipy.signal import hilbert, chirp
import matplotlib.pyplot as plt

fs = 600.0 #sampling frequency
duration = 1.0 #duration of the signal
t = np.arange(int(fs*duration)) / fs #time base

a_t =  1.0 + 0.7 * np.sin(2.0*np.pi*3.0*t)#information signal
c_t = chirp(t, 20.0, t[-1], 80) #chirp carrier
x = a_t * c_t #modulated signal

plt.subplot(2,1,1)
plt.plot(x) #plot the modulated signal

z= hilbert(x) #form the analytical signal
inst_amplitude = np.abs(z) #envelope extraction
inst_phase = np.unwrap(np.angle(z))#inst phase
inst_freq = np.diff(inst_phase)/(2*np.pi)*fs #inst frequency

#Regenerate the carrier from the instantaneous phase
regenerated_carrier = np.cos(inst_phase)

plt.plot(inst_amplitude,'r'); #overlay the extracted envelope
plt.title('Modulated signal and extracted envelope')
plt.xlabel('n')
plt.ylabel('x(t) and |z(t)|')
plt.subplot(2,1,2)
plt.plot(regenerated_carrier)
plt.title('Extracted carrier or TFS')
plt.xlabel('n')
plt.ylabel('cos[\omega(t)]')

Results

*Figure 2: Amplitude modulation using a chirp signal and extraction of envelope and TFS*

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

Reference

[1] Moon, Il Joon, and Sung Hwa Hong. “What Is Temporal Fine Structure and Why Is It Important?” Korean Journal of Audiology 18.1 (2014): 1–7. PMC. Web. 24 Apr. 2017.↗

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

Understanding Analytic Signal and Hilbert Transform

Key focus of this article: Understand the relationship between analytic signal, Hilbert transform and FFT. Hands-on demonstration using Python and Matlab.

Introduction

Fourier Transform of a real-valued signal is complex-symmetric. It implies that the content at negative frequencies are redundant with respect to the positive frequencies. In their works, Gabor [1] and Ville [2], aimed to create an analytic signal by removing redundant negative frequency content resulting from the Fourier transform. The analytic signal is complex-valued but its spectrum will be one-sided (only positive frequencies) that preserved the spectral content of the original real-valued signal. Using an analytic signal instead of the original real-valued signal, has proven to be useful in many signal processing applications. For example, in spectral analysis, use of analytic signal in-lieu of the original real-valued signal mitigates estimation biases and eliminates cross-term artifacts due to negative and positive frequency components [3].

Continuous-time analytic signal

Let be a real-valued non-bandlimited finite energy signal, for which we wish to construct a corresponding analytic signal . The Continuous Time Fourier Transform of is given by

Lets say the magnitude spectrum of is as shown in Figure 1(a). We note that the signal is a real-valued and its magnitude spectrum is symmetric and extends infinitely in the frequency domain.

*Figure 1: (a) Spectrum of continuous signal x(t) and (b) spectrum of analytic signal z(t)*

As mentioned in the introduction, an analytic signal can be formed by suppressing the negative frequency contents of the Fourier Transform of the real-valued signal. That is, in frequency domain, the spectral content of the analytic signal is given by

The corresponding spectrum of the resulting analytic signal is shown in Figure 1(b).

Since the spectrum of the analytic signal is one-sided, the analytic signal will be complex valued in the time domain, hence the analytic signal can be represented in terms of real and imaginary components as . Since the spectral content is preserved in an analytic signal, it turns out that the real part of the analytic signal in time domain is essentially the original real-valued signal itself . Then, what takes place of the imaginary part ? Who is the companion to x(t) that occupies the imaginary part in the resulting analytic signal ? Summarizing as equation,

It is interesting to note that Hilbert transform [4] can be used to find a companion function (imaginary part in the equation above) to a real-valued signal such that the real signal can be analytically extended from the real axis to the upper half of the complex plane . Denoting Hilbert transform as , the analytic signal is given by

From these discussion, we can see that an analytic signal for a real-valued signal , can be constructed using two approaches.

● Frequency domain approach: The one-sided spectrum of is formed from the two-sided spectrum of the real-valued signal by applying equation (2)
● Time domain approach: Using Hilbert transform approach given in equation (4)

One of the important property of an analytic signal is that its real and imaginary components are orthogonal

Discrete-time analytic signal

Since we are in digital era, we are more interested in discrete-time signal processing. Consider a continuous real-valued signal gets sampled at interval seconds and results in real-valued discrete samples , i.e, . The spectrum of the continuous signal is shown in Figure 2(a). The spectrum of that results from the process of periodic sampling is given in Figure 2(b) (Refer here more details on the process of sampling). The spectrum of discrete-time signal can be obtained by Discrete-Time Fourier Transform (DTFT).

*Figure 2: (a) CTFT of continuous signal x(t), (b) Spectrum of x[n] resulted due to periodic sampling and (c) Periodic one-sided spectrum of analytical signal z[n]*

At this point, we would like to construct a discrete-time analytic signal from the real-valued sampled signal . We wish the analytic signal is complex valued and should satisfy the following two desired properties

● The real part of the analytic signal should be same as the original real-valued signal.
● The real and imaginary part of the analytic signal should satisfy the following property of orthogonality

In Frequency domain approach for the continuous-time case, we saw that an analytic signal is constructed by suppressing the negative frequency components from the spectrum of the real signal. We cannot do this for our periodically sampled signal . Periodic mirroring nature of the spectrum prevents one from suppressing the negative components. If we do so, it will vanish the entire spectrum. One solution to this problem is to set the negative half of each spectral period to zero. The resulting spectrum of the analytic signal is shown in Figure 2(c).

Given a record of samples of even length , the procedure to construct the analytic signal is as follows. This method satisfies both the desired properties listed above.

● Compute the -point DTFT of using FFT
● N-point periodic one-sided analytic signal is computed by the following transform

● Finally, the analytic signal (z[n]) is obtained by taking the inverse DTFT of

Matlab

The given procedure can be coded in Matlab using the FFT function. Given a record of real-valued samples , the corresponding analytic signal can be constructed as given next. Note that the Matlab has an inbuilt function to compute the analytic signal. The in-built function is called hilbert.

function z = analytic_signal(x)
%x is a real-valued record of length N, where N is even %returns the analytic signal z[n]
x = x(:); %serialize
N = length(x);
X = fft(x,N);
z = ifft([X(1); 2*X(2:N/2); X(N/2+1); zeros(N/2-1,1)],N);
end

To test this function, we create a 5 seconds record of a real-valued sine signal. The analytic signal is constructed and the orthogonal components are plotted in Figure 3. From the plot, we can see that the real part of the analytic signal is exactly same as the original signal (which is the cosine signal) and the imaginary part of the analytic signal is phase shifted version of the original signal. We note that the imaginary part of the analytic signal is a cosine function with amplitude scaled by which is none other than the Hilbert transform of sine function.

t=0:0.001:0.5-0.001;
x = sin(2*pi*10*t); %real-valued f = 10 Hz
subplot(2,1,1); plot(t,x);%plot the original signal
title('x[n] - original signal'); xlabel('n'); ylabel('x[n]');

z = analytic_signal(x); %construct analytic signal
subplot(2,1,2); plot(t, real(z), 'k'); hold on;
plot(t, imag(z), 'r');
title('Components of Analytic signal'); 
xlabel('n'); ylabel('z_r[n] and z_i[n]');
legend('Real(z[n])','Imag(z[n])');

Python

Equivalent code in Python is given below (tested with Python 3.6.0)

import numpy as np
def main():
    t = np.arange(start=0,stop=0.5,step=0.001)
    x = np.sin(2*np.pi*10*t)
    
    import matplotlib.pyplot as plt
    plt.subplot(2,1,1)
    plt.plot(t,x)
    plt.title('x[n] - original signal')
    plt.xlabel('n')
    plt.ylabel('x[n]')
    
    z = analytic_signal(x)
    
    plt.subplot(2,1,2)
    plt.plot(t,z.real,'k',label='Real(z[n])')
    plt.plot(t,z.imag,'r',label='Imag(z[n])')
    plt.title('Components of Analytic signal')
    plt.xlabel('n')
    plt.ylabel('z_r[n] and z_i[n]')
    plt.legend()

def analytic_signal(x):
    from scipy.fftpack import fft,ifft
    N = len(x)
    X = fft(x,N)
    h = np.zeros(N)
    h[0] = 1
    h[1:N//2] = 2*np.ones(N//2-1)
    h[N//2] = 1
    Z = X*h
    z = ifft(Z,N)
    return z

if __name__ == '__main__':
    main()

*Figure 3: Components of analytic signal for a real-valued sine function*

Hilbert Transform using FFT

We should note that the hilbert function in Matlab returns the analytic signal $latex z[n]$ not the hilbert transform of the signal . To get the hilbert transform, we should simply get the imaginary part of the analytic signal. Since we have written our own function to compute the analytic signal, getting the hilbert transform of a real-valued signal goes like this.

x_hilbert = imag(analytic_signal(x))

In the coming posts, we will some of the applications of constructing an analytic signal. For example: Find the instantaneous amplitude and phase of a signal, envelope detector for an amplitude modulated signal, detecting phase changes in a sine wave.

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

References:

[1] D. Gabor, “Theory of communications”, Journal of the Inst. Electr. Eng., vol. 93, pt. 111, pp. 42-57, 1946. See definition of complex signal on p. 432.↗
[2] J. A. Ville, “Theorie et application de la notion du signal analytique”, Cables el Transmission, vol. 2, pp. 61-74, 1948.↗
[3] S. M. Kay, “Maximum entropy spectral estimation using the analytical signal”, IEEE transactions on Acoustics, Speech, and Signal Processing, vol. 26, pp. 467-469, October 1978.↗
[4] Frank R. Kschischang, “The Hilbert Transform”, University of Toronto, October 22, 2006.↗
[5] S. L. Marple, “Computing the discrete-time ‘analytic’ signal via FFT,” Conference Record of the Thirty-First Asilomar Conference on Signals, Systems and Computers , Pacific Grove, CA, USA, 1997, pp. 1322-1325 vol.2.↗

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

Interpret FFT results – obtaining magnitude and phase information

In the previous post, Interpretation of frequency bins, frequency axis arrangement (fftshift/ifftshift) for complex DFT were discussed. In this post, I intend to show you how to interpret FFT results and obtain magnitude and phase information.

Outline

For the discussion here, lets take an arbitrary cosine function of the form $x(t)= A cos \left(2 \pi f_c t + \phi \right)$ and proceed step by step as follows

● Represent the signal $x(t)$ in computer (discrete-time) and plot the signal (time domain)
● Represent the signal in frequency domain using FFT ($ X[k]$)
● Extract amplitude and phase information from the FFT result
● Reconstruct the time domain signal from the frequency domain samples

This article is part of the book Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885 available in ebook (PDF) format (click here) and Paperback (hardcopy) format (click here)
Wireless Communication Systems in Matlab, ISBN: 978-1720114352 available in ebook (PDF) format (click here) and Paperback (hardcopy) format (click here).

Discrete-time domain representation

Consider a cosine signal of amplitude $A=0.5$, frequency $f_c=10 Hz$ and phase $phi= \pi/6$ radians (or $30^{\circ}$ )

\[x(t) = 0.5 cos \left( 2 \pi 10 t + \pi/6 \right)\]

In order to represent the continuous time signal $x(t)$ in computer memory, we need to sample the signal at sufficiently high rate (according to Nyquist sampling theorem). I have chosen a oversampling factor of $32$ so that the sampling frequency will be $f_s = 32 \times f_c $, and that gives $640$ samples in a $2$ seconds duration of the waveform record.

A = 0.5; %amplitude of the cosine wave
fc=10;%frequency of the cosine wave
phase=30; %desired phase shift of the cosine in degrees
fs=32*fc;%sampling frequency with oversampling factor 32
t=0:1/fs:2-1/fs;%2 seconds duration

phi = phase*pi/180; %convert phase shift in degrees in radians
x=A*cos(2*pi*fc*t+phi);%time domain signal with phase shift

figure; plot(t,x); %plot the signal

Represent the signal in frequency domain using FFT

Lets represent the signal in frequency domain using the FFT function. The FFT function computes $N$-point complex DFT. The length of the transformation $N$ should cover the signal of interest otherwise we will some loose valuable information in the conversion process to frequency domain. However, we can choose a reasonable length if we know about the nature of the signal.

For example, the cosine signal of our interest is periodic in nature and is of length $640$ samples (for 2 seconds duration signal). We can simply use a lower number $N=256$ for computing the FFT. In this case, only the first $256$ time domain samples will be considered for taking FFT. No need to worry about loss of information in this case, as the $256$ samples will have sufficient number of cycles using which we can calculate the frequency information.

N=256; %FFT size
X = 1/N*fftshift(fft(x,N));%N-point complex DFT

In the code above, $fftshift$ is used only for obtaining a nice double-sided frequency spectrum that delineates negative frequencies and positive frequencies in order. This transformation is not necessary. A scaling factor $1/N$ was used to account for the difference between the FFT implementation in Matlab and the text definition of complex DFT.

3a. Extract amplitude of frequency components (amplitude spectrum)

The FFT function computes the complex DFT and the hence the results in a sequence of complex numbers of form $X_{re} + j X_{im}$. The amplitude spectrum is obtained

\[|X[k]| = \sqrt{X_{re}^2 + X_{im}^2 } \]

For obtaining a double-sided plot, the ordered frequency axis (result of fftshift) is computed based on the sampling frequency and the amplitude spectrum is plotted.

df=fs/N; %frequency resolution
sampleIndex = -N/2:N/2-1; %ordered index for FFT plot
f=sampleIndex*df; %x-axis index converted to ordered frequencies
stem(f,abs(X)); %magnitudes vs frequencies
xlabel('f (Hz)'); ylabel('|X(k)|');

3b. Extract phase of frequency components (phase spectrum)

Extracting the correct phase spectrum is a tricky business. I will show you why it is so. The phase of the spectral components are computed as

\[\angle X[k] = tan^{-1} \left( \frac{X_{im}}{X_{re}} \right)\]

That equation looks naive, but one should be careful when computing the inverse tangents using computers. The obvious choice for implementation seems to be the $atan$ function in Matlab. However, usage of $atan$ function will prove disastrous unless additional precautions are taken. The $atan$ function computes the inverse tangent over two quadrants only, i.e, it will return values only in the $[-\pi/2 , \pi/2]$ interval. Therefore, the phase need to be unwrapped properly. We can simply fix this issue by computing the inverse tangent over all the four quadrants using the $atan2(X_{img},X_{re})$ function.

Lets compute and plot the phase information using $atan2$ function and see how the phase spectrum looks

phase=atan2(imag(X),real(X))*180/pi; %phase information
plot(f,phase); %phase vs frequencies

The phase spectrum is completely noisy. Unexpected !!!. The phase spectrum is noisy due to fact that the inverse tangents are computed from the $ratio$ of imaginary part to real part of the FFT result. Even a small floating rounding off error will amplify the result and manifest incorrectly as useful phase information (read how a computer program approximates very small numbers).

To understand, print the first few samples from the FFT result and observe that they are not absolute zeros (they are very small numbers in the order $10^{-16}$. Computing inverse tangent will result in incorrect results.

>> X(1:5)
ans =
   1.0e-16 *
  -0.7286            -0.3637 - 0.2501i  -0.4809 - 0.1579i  -0.3602 - 0.5579i   0.0261 - 0.4950i
>> atan2(imag(X(1:5)),real(X(1:5)))
ans =
    3.1416   -2.5391   -2.8244   -2.1441   -1.5181

The solution is to define a tolerance threshold and ignore all the computed phase values that are below the threshold.

X2=X;%store the FFT results in another array
%detect noise (very small numbers (eps)) and ignore them
threshold = max(abs(X))/10000; %tolerance threshold
X2(abs(X)<threshold) = 0; %maskout values that are below the threshold
phase=atan2(imag(X2),real(X2))*180/pi; %phase information
plot(f,phase); %phase vs frequencies

The recomputed phase spectrum is plotted below. The phase spectrum has correctly registered the $30^{\circ}$ phase shift at the frequency $f=10 Hz$. The phase spectrum is anti-symmetric ($\phi=-30^{\circ}$ at $f=-10 Hz$ ), which is expected for real-valued signals.

Reconstruct the time domain signal from the frequency domain samples

Reconstruction of the time domain signal from the frequency domain sample is pretty straightforward

x_recon = N*ifft(ifftshift(X),N); %reconstructed signal
t = [0:1:length(x_recon)-1]/fs; %recompute time index 
plot(t,x_recon);%reconstructed signal

The reconstructed signal has preserved the same initial phase shift and the frequency of the original signal. Note: The length of the reconstructed signal is only $256$ sample long ($\approx 0.8$ seconds duration), this is because the size of FFT is considered as $N=256$. Since the signal is periodic it is not a concern. For more complicated signals, appropriate FFT length (better to use a value that is larger than the length of the signal) need to be used.

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

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

Interpret FFT, complex DFT, frequency bins & FFTShift

Key focus: Interpret FFT results, complex DFT, frequency bins, fftshift and ifftshift. Know how to use them in analysis using Matlab and Python.

Four types of Fourier Transforms:

Often, one is confronted with the problem of converting a time domain signal to frequency domain and vice-versa. Fourier Transform is an excellent tool to achieve this conversion and is ubiquitously used in many applications. In signal processing, a time domain signal can be continuous or discrete and it can be aperiodic or periodic. This gives rise to four types of Fourier transforms.

From Table 1, we note note that when the signal is discrete in one domain, it will be periodic in other domain. Similarly, if the signal is continuous in one domain, it will be aperiodic (non-periodic) in another domain. For simplicity, let’s not venture into the specific equations for each of the transforms above. We will limit our discussion to DFT, that is widely available as part of software packages like Matlab, Scipy(python) etc.., however we can approximate other transforms using DFT.

Real and Complex versions of transforms

For each of the listed transforms above, there exist a real version and complex version. The real version of the transform, takes in a real numbers and gives two sets of real frequency domain points – one set representing coefficients over $cosine$ basis function and the other set representing the co-efficient over $sine$ basis function. The complex version of the transforms represent positive and negative frequencies in a single array. The complex versions are flexible that it can process both complex valued signals and real valued signals. The following figure captures the difference between real DFT and complex DFT.

Real DFT:

Consider the case of N-point $real$ DFT , it takes in N samples of (real-valued) time domain waveform $x[n]$ and gives two arrays of length $N/2+1$ each set projected on cosine and sine functions respectively.

\[\begin{align} X_{re}[k] &= \phantom{-}\frac{2}{N} \displaystyle{ \sum_{n=0}^{N-1} x[n] \cdot cos\left( \frac{2 \pi k n}{N} \right)} \\ X_{im}[k] &= -\frac{2}{N} \sum_{n=0}^{N-1} \displaystyle{ x[n] \cdot sin\left( \frac{2 \pi k n}{N} \right)} \end{align}\]

Here, the time domain index $n$ runs from $0 \rightarrow N$, the frequency domain index $k$ runs from $0 \rightarrow N/2$

The real-valued time domain signal $x [n]$ can be synthesized from the real DFT pairs as

\[x[n] = \displaystyle{ \sum_{k=0}^{N/2} X_{re}[K] \cdot cos\left( \frac{2 \pi k n}{N} \right) – X_{im}[K] \cdot sin\left( \frac{2 \pi k n}{N} \right)}\]

Caveat: When using the synthesis equation, the values $X_{re}[0]$ and $X_{re}[N/2] $ must be divided by two. This problem is due to the fact that we restrict the analysis to real-values only. These type of problems can be avoided by using complex version of DFT.

Complex DFT:

Consider the case of N-point complex DFT, it takes in N samples of complex-valued time domain waveform $x[n]$ and produces an array $X[k]$ of length N.

\[ X[k]= \displaystyle{\sum_{n=0}^{N-1} x[n] e^{-j2 \pi k n/N}}\]

The arrays values are interpreted as follows

● $X[0]$ represents DC frequency component
● Next $N/2$ terms are positive frequency components with $X[N/2]$ being the Nyquist frequency (which is equal to half of sampling frequency)
● Next $N/2-1$ terms are negative frequency components (note: negative frequency components are the phasors rotating in opposite direction, they can be optionally omitted depending on the application)

The corresponding synthesis equation (reconstruct $x[n]$ from frequency domain samples $X[k]$) is

\[x[n]= \displaystyle{ \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j2 \pi k n/N}} \]

From these equations we can see that the real DFT is computed by projecting the signal on cosine and sine basis functions. However, the complex DFT projects the input signal on exponential basis functions (Euler’s formula connects these two concepts).

When the input signal in the time domain is real valued, the complex DFT zero-fills the imaginary part during computation (That’s its flexibility and avoids the caveat needed for real DFT). The following figure shows how to interpret the raw FFT results in Matlab that computes complex DFT. The specifics will be discussed next with an example.

Figure 2: Interpretation of frequencies in complex DFT output

Fast Fourier Transform (FFT)

The FFT function in Matlab is an algorithm published in 1965 by J.W.Cooley and J.W.Tuckey for efficiently calculating the DFT. It exploits the special structure of DFT when the signal length is a power of 2, when this happens, the computation complexity is significantly reduced. FFT length is generally considered as power of 2 – this is called $radix-2$ FFT which exploits the twiddle factors. The FFT length can be odd as used in this particular FFT implementation – Prime-factor FFT algorithm where the FFT length factors into two co-primes.

FFT is widely available in software packages like Matlab, Scipy etc.., FFT in Matlab/Scipy implements the complex version of DFT. Matlab’s FFT implementation computes the complex DFT that is very similar to above equations except for the scaling factor. For comparison, the Matlab’s FFT implementation computes the complex DFT and its inverse as

\[ \begin{align} X[k] &= \displaystyle{ \phantom {\frac{1}{N}}\sum_{n=0}^{N-1} x[n] e^{-j2 \pi k n/N}} \\ x[n] &= \displaystyle{ \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j2 \pi k n/N}} \end{align}\]

The Matlab commands that implement the above equations are FFT and IFFT) respectively. The corresponding syntax is as follows

X = fft(x,N) %compute X[k]
x = ifft(X,N) %compute x[n]

Interpreting the FFT results

Lets assume that the $x[n]$ is the time domain cosine signal of frequency $f_c=10Hz$ that is sampled at a frequency $f_s=32*fc$ for representing it in the computer memory.

fc=10;%frequency of the carrier
fs=32*fc;%sampling frequency with oversampling factor=32
t=0:1/fs:2-1/fs;%2 seconds duration
x=cos(2*pi*fc*t);%time domain signal (real number)
subplot(3,1,1); plot(t,x);hold on; %plot the signal
title('x[n]=cos(2 \pi 10 t)'); xlabel('n'); ylabel('x[n]');

Figure 3: A 2 seconds record of 10 Hz cosine wave

Lets consider taking a $N=256$ point FFT, which is the $8^{th}$ power of $2$.

N=256; %FFT size
X = fft(x,N);%N-point complex DFT
%output contains DC at index 1, Nyquist frequency at N/2+1 th index
%positive frequencies from index 2 to N/2
%negative frequencies from index N/2+1 to N

Note: The FFT length should be sufficient to cover the entire length of the input signal. If $N$ is less than the length of the input signal, the input signal will be truncated when computing the FFT. In our case, the cosine wave is of 2 seconds duration and it will have 640 points (a $10Hz$ frequency wave sampled at 32 times oversampling factor will have $2 \times 32 \times 10 = 640$ samples in 2 seconds of the record). Since our input signal is periodic, we can safely use $N=256$ point FFT, anyways the FFT will extend the signal when computing the FFT (see additional topic on spectral leakage that explains this extension concept).

Due to Matlab’s index starting at 1, the DC component of the FFT decomposition is present at index 1.

>>X(1)
 1.1762e-14   (approximately equal to zero)

That’s pretty easy. Note that the index for the raw FFT are integers from $1 \rightarrow N$. We need to process it to convert these integers to $frequencies$. That is where the $sampling$ frequency counts.

Each point/bin in the FFT output array is spaced by the frequency resolution $\Delta f$ that is calculated as

\[ \Delta f = \frac{f_s}{N} \]

where, $f_s$ is the sampling frequency and $N$ is the FFT size that is considered. Thus, for our example, each point in the array is spaced by the frequency resolution

\[\Delta f = \frac{f_s}{N} = \frac{32*f_c}{256} = \frac{320}{256} = 1.25 Hz \]

Now, the $10 Hz$ cosine signal will leave a spike at the 8th sample ($10/1.25=8$), which is located at index 9 (See next figure).

>> abs(X(8:10)) %display samples 7 to 9
ans = 0.0000  128.0000    0.0000

Therefore, from the frequency resolution, the entire frequency axis can be computed as

%calculate frequency bins with FFT
df=fs/N %frequency resolution
sampleIndex = 0:N-1; %raw index for FFT plot
f=sampleIndex*df; %x-axis index converted to frequencies

Now we can plot the absolute value of the FFT against frequencies as

subplot(3,1,2); stem(sampleIndex,abs(X)); %sample values on x-axis
title('X[k]'); xlabel('k'); ylabel('|X(k)|');
subplot(3,1,3); plot(f,abs(X)); %x-axis represent frequencies
title('X[k]'); xlabel('frequencies (f)'); ylabel('|X(k)|');

The following plot shows the frequency axis and the sample index as it is for the complex FFT output.

Figure 4: Magnitude response from FFT output plotted against sample index (top) and computed frequencies (bottom)

After the frequency axis is properly transformed with respect to the sampling frequency, we note that the cosine signal has registered a spike at $10 Hz$. In addition to that, it has also registered a spike at $256-8=248^{th}$ sample that belongs to negative frequency portion. Since we know the nature of the signal, we can optionally ignore the negative frequencies. The sample at the Nyquist frequency ($f_s/2$) mark the boundary between the positive and negative frequencies.

>> nyquistIndex=N/2+1;
>> X(nyquistIndex-2:nyquistIndex+2).'
ans =
1.0e-13 *
  -0.2428 + 0.0404i
  -0.1897 + 0.0999i
  -0.3784          
  -0.1897 - 0.0999i
  -0.2428 - 0.0404i

Note that the complex numbers surrounding the Nyquist index are complex conjugates and they represent positive and negative frequencies respectively.

FFTShift

From the plot we see that the frequency axis starts with DC, followed by positive frequency terms which is in turn followed by the negative frequency terms. To introduce proper order in the x-axis, one can use FFTshift function Matlab, which arranges the frequencies in order: negative frequencies $\rightarrow$ DC $\rightarrow$ positive frequencies. The fftshift function need to be carefully used when $N$ is odd.

For even N, the original order returned by FFT is as follows (note: all indices below corresponds to Matlab’s index)

● $X[1]$ represents DC frequency component
● $X[2]$ to $X[N/2]$ terms are positive frequency components
● $X[N/2+1]$ is the Nyquist frequency ($F_s/2$) that is common to both positive and negative frequencies. We will consider it as part of negative frequencies to have the same equivalence with the fftshift function.
● $X[N/2+1]$ to $X[N]$ terms are considered as negative frequency components

FFTshift shifts the DC component to the center of the spectrum. It is important to remember that the Nyquist frequency at the (N/2+1)th Matlab index is common to both positive and negative frequency sides. FFTshift command puts the Nyquist frequency in the negative frequency side. This is captured in the following illustration.

Figure 5: Role of FFTShift in ordering the frequencies

Therefore, when $N$ is even, ordered frequency axis is set as

\[f = \Delta f \left[ -\frac{N}{2}:1:\frac{N}{2}-1 \right] = \frac{f_s}{N} \left[ -\frac{N}{2}:1:\frac{N}{2}-1 \right]\]

When (N) is odd, the ordered frequency axis should be set as

\[ f = \Delta f \left[ -\frac{N-1}{2}:1:\frac{N+1}{2}-1 \right] = \frac{f_s}{N} \left[ -\frac{N-1}{2}:1:\frac{N-1}{2}-1 \right]\]

The following code snippet, computes the fftshift using both the manual method and using the Matlab’s in-build command. The results are plotted by superimposing them on each other. The plot shows that both the manual method and fftshift method are in good agreement.

%two-sided FFT with negative frequencies on left and positive frequencies on right
%following works only if N is even, for odd N see equation above
X1 = [(X(N/2+1:N)) X(1:N/2)]; %order frequencies without using fftShift
X2 = fftshift(X);%order frequencies by using fftshift

df=fs/N; %frequency resolution
sampleIndex = -N/2:N/2-1; %raw index for FFT plot
f=sampleIndex*df; %x-axis index converted to frequencies
%plot ordered spectrum using the two methods
figure;subplot(2,1,1);stem(sampleIndex,abs(X1));hold on; stem(sampleIndex,abs(X2),'r') %sample index on x-axis
title('Frequency Domain'); xlabel('k'); ylabel('|X(k)|');%x-axis represent sample index
subplot(2,1,2);stem(f,abs(X1)); stem(f,abs(X2),'r') %x-axis represent frequencies
xlabel('frequencies (f)'); ylabel('|X(k)|');

Figure 6: Magnitude response of FFT result after applying FFTShift : plotted against sample index (top) and against computed frequencies (bottom)

Comparing the bottom figures in the Figure 4 and Figure 6, we see that the ordered frequency axis is more meaningful to interpret.

IFFTShift

One can undo the effect of fftshift by employing ifftshift function. The ifftshift function restores the raw frequency order. If the FFT output is ordered using fftshift function, then one must restore the frequency components back to original order BEFORE taking IFFT. Following statements are equivalent.

X = fft(x,N) %compute X[k]
x = ifft(X,N) %compute x[n]

X = fftshift(fft(x,N)); %take FFT and rearrange frequency order (this is mainly done for interpretation)
x = ifft(ifftshift(X),N)% restore raw frequency order and then take IFFT

Some observations on FFTShift and IFFTShift

When $N$ is odd and for an arbitrary sequence, the fftshift and ifftshift functions will produce different results. However, when they are used in tandem, it restores the original sequence.

>>  x=[0,1,2,3,4,5,6,7,8]
0     1     2     3     4     5     6     7     8
>>  fftshift(x)
5     6     7     8     0     1     2     3     4
>>  ifftshift(x)
4     5     6     7     8     0     1     2     3
>>  ifftshift(fftshift(x))
0     1     2     3     4     5     6     7     8
>>  fftshift(ifftshift(x))
0     1     2     3     4     5     6     7     8

When $N$ is even and for an arbitrary sequence, the fftshift and ifftshift functions will produce the same result. When they are used in tandem, it restores the original sequence.

>>  x=[0,1,2,3,4,5,6,7]
0     1     2     3     4     5     6     7
>>  fftshift(x)
4     5     6     7     0     1     2     3
>>  ifftshift(x)
4     5     6     7     0     1     2     3
>>  ifftshift(fftshift(x))
0     1     2     3     4     5     6     7
>>  fftshift(ifftshift(x))
0     1     2     3     4     5     6     7

For Python code, please check the following book: Digital Modulations using Python – by Mathuranathan Viswanathan

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

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

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

Chirp Signal – FFT & PSD in Matlab & Python

Key focus: Know how to generate a Chirp signal, compute its Fourier Transform using FFT and power spectral density (PSD) in Matlab & Python.

Introduction

All the signals discussed so far do not change in frequency over time. Obtaining a signal with time-varying frequency is of main focus here. A signal that varies in frequency over time is called “chirp”. The frequency of the chirp signal can vary from low to high frequency (up-chirp) or from high to low frequency (low-chirp).

Observation

Chirp signals/signatures are encountered in many applications ranging from radar, sonar, spread spectrum, optical communication, image processing, doppler effect, motion of a pendulum, as gravitation waves, manifestation as Frequency Modulation (FM), echo location ^[1] etc.

Mathematical Description:

A linear chirp signal sweeps the frequency from low to high frequency (or vice-versa) linearly. One approach to generate a chirp signal is to concatenate a series of segments of sine waves each with increasing(or decreasing) frequency in order. This method introduces discontinuities in the chirp signal due to the mismatch in the phases of each such segments. Modifying the equation of a sinusoid to generate a chirp signal is a better approach.

The equation for generating a sinusoidal (cosine here) signal with amplitude A, angular frequency and initial phase is

This can be written as a function of instantaneous phase

where is the instantaneous phase of the sinusoid and it is linear in time. The time derivative of instantaneous phase is equal to the angular frequency of the sinusoid – which in case is a constant in the above equation.

Instead of having the phase linear in time, let’s change the phase to quadratic form and thus non-linear.

for some constant .

Therefore, the equation for chirp signal takes the following form,

The first derivative of the phase, which is the instantaneous angular frequency becomes a function of time, which is given by

The time-varying frequency in Hertz is given by

In the above equation, the frequency is no longer a constant, rather it is of time-varying nature with initial frequency given by . Thus, from the above equation, given a time duration T, the rate of change of frequency is given by

where, is the starting frequency of the sweep, is the frequency at the end of the duration T.

Substituting (7) & (8) in (6)

From (6) and (8)

where is a constant which will act as the initial phase of the sweep.

Thus the modified equation for generating a chirp signal (from equations (5) and (10)) is given by

where the time-varying frequency function is given by

Generation of Chirp signal, computing its Fourier Transform using FFT and power spectral density (PSD) in Matlab is shown as example, for Python code, please refer the book Digital Modulations using Python.

Generating a chirp signal without using in-built “chirp” Function in Matlab:

Implement a function that describes the chirp using equation (11) and (12). The starting frequency of the sweep is and the frequency at time is . The initial phase forms the final part of the argument in the following function

function x=mychirp(t,f0,t1,f1,phase)
%Y = mychirp(t,f0,t1,f1) generates samples of a linear swept-frequency
%   signal at the time instances defined in timebase array t.  The instantaneous
%   frequency at time 0 is f0 Hertz.  The instantaneous frequency f1
%   is achieved at time t1.
%   The argument 'phase' is optional. It defines the initial phase of the
%   signal degined in radians. By default phase=0 radian
    
if nargin==4
    phase=0;
end
    t0=t(1);
    T=t1-t0;
    k=(f1-f0)/T;
    x=cos(2*pi*(k/2*t+f0).*t+phase);
end

The following wrapper script utilizes the above function and generates a chirp with starting frequency at the start of the time base and at which is the end of the time base. From the PSD plot, it can be ascertained that the signal energy is concentrated only upto 25 Hz

fs=500; %sampling frequency
t=0:1/fs:1; %time base - upto 1 second

f0=1;% starting frequency of the chirp
f1=fs/20; %frequency of the chirp at t1=1 second
x = mychirp(t,f0,1,f1); 
subplot(2,2,1)
plot(t,x,'k');
title(['Chirp Signal']);
xlabel('Time(s)');
ylabel('Amplitude');

FFT and power spectral density

As with other signals, describes in the previous posts, let’s plot the FFT of the generated chirp signal and its power spectral density (PSD).

L=length(x);
NFFT = 1024;
X = fftshift(fft(x,NFFT));
Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

subplot(2,2,2)
plot(f,abs(X)/(L),'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
xlim([-50 50])


Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
subplot(2,2,3)
plot(f,10*log10(Pxx),'r');
title('Double Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
xlim([-100 100])

X = fft(x,NFFT);
X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
Pxx=X.*conj(X)/(NFFT*NFFT);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
subplot(2,2,4)
plot(f,10*log10(Pxx),'r');
title('Single Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');

For Python code, please refer the book Digital Modulations using Python

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

References:

[1] Patrick Flandrin,“Chirps everywhere”,CNRS — Ecole Normale Supérieure de Lyon

Topics in this chapter

Books by the author

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

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

Spectral leakage

Equivalent Noise Bandwidth:[2]

Coherent power gain:

Scalloping loss

Figure of merits for window functions

List of window functions

Boxcar (Rectangular)

Barthann

Bartlett

Blackman

Blackman-Harris

Bohman

Cosine

Flattop

Hamming

Hann

Nuttall

Parzen

Triangular

References:

Similar articles

Additional Resources:

Books by the author

FFT and spectral leakage

Equivalent noise bandwidth

Python script

Reference

Similar articles

Books by the author

Introduction

Sine Wave

Different representations of FFT:

1. Plotting raw values of DFT:

2. FFT plot – plotting raw values against normalized frequency axis:

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

4. FFT plot – Absolute frequency on the x-axis vs. magnitude on y-axis:

5. Power Spectrum – Absolute frequency on the x-axis vs. power on y-axis:

Books by the author

Topics in this chapter

Phase modulated signal:

Demodulating a phase modulated signal:

Matlab code

Python code

Results

For further reading

Topics in this chapter

Books by the author

Introduction

Definition

Problem statement

Solution

Matlab

Python

Results

Reference

Topics in this chapter

Books by the author

Introduction

Continuous-time analytic signal

Discrete-time analytic signal

Matlab

Python

Hilbert Transform using FFT

References:

Topics in this chapter

Books by the author

Outline

Discrete-time domain representation

Represent the signal in frequency domain using FFT

3a. Extract amplitude of frequency components (amplitude spectrum)

3b. Extract phase of frequency components (phase spectrum)

Reconstruct the time domain signal from the frequency domain samples

Topics in this chapter

Books by the author

Four types of Fourier Transforms:

Real and Complex versions of transforms

Real DFT:

Complex DFT:

Fast Fourier Transform (FFT)

Interpreting the FFT results

Equivalent Noise Bandwidth:^[2]