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
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
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
An oversampling factor of
Different representations of FFT:
Since FFT is just a numeric computation of
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
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()
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
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()
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
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()
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
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()
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
Plotting the PSD plot with y-axis on log scale, produces the most encountered type of PSD plot in signal processing.
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. |  Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. |  Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. |
| Hand-picked Best books on Communication Engineering Best books on Signal Processing |
||
Topics in this chapter
![\begin{aligned} G(f) &=F[g(t)]\\ &= \int_{-\infty }^{\infty } g(t)e^{-j2\pi ft}\, dt\\ &= \frac{1}{\sigma \sqrt{2 \pi } } \int_{-\infty }^{\infty } e^{- \frac{t^2}{2 \sigma^2}}e^{-j2\pi ft}\, dt\\ &=\frac{1}{\sigma \sqrt{2 \pi } } \int_{-\infty }^{\infty } e^{- \frac{1}{2 \sigma^2}\left[t^2+j4 \pi \sigma^2 ft \right]}\, dt\\ &=\frac{1}{\sigma \sqrt{2 \pi } } \int_{-\infty }^{\infty } e^{- \frac{1}{2 \sigma^2}\left[t^2+j4 \pi \sigma^2 ft + (j 2 \pi \sigma^2 f)^2 - (j 2 \pi \sigma^2 f)^2\right]}\, dt\\ &=e^{ \frac{1}{2 \sigma^2}(j 2 \pi \sigma^2 f)^2}\frac{1}{\sigma \sqrt{2 \pi } } \int_{-\infty }^{\infty } e^{- \frac{1}{2 \sigma^2}\left[t+j 2 \pi \sigma^2 f \right]^2}\, dt\\ &=e^{ \frac{1}{2 \sigma^2}(j 2 \pi \sigma^2 f)^2}=e^{ - \frac{1}{2}( 2 \pi \sigma f)^2}\end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+G%28f%29+%26%3DF%5Bg%28t%29%5D%5C%5C+%26%3D+%5Cint_%7B-%5Cinfty+%7D%5E%7B%5Cinfty+%7D+g%28t%29e%5E%7B-j2%5Cpi+ft%7D%5C%2C+dt%5C%5C+%26%3D+%5Cfrac%7B1%7D%7B%5Csigma+%5Csqrt%7B2+%5Cpi+%7D+%7D+%5Cint_%7B-%5Cinfty+%7D%5E%7B%5Cinfty+%7D+e%5E%7B-+%5Cfrac%7Bt%5E2%7D%7B2+%5Csigma%5E2%7D%7De%5E%7B-j2%5Cpi+ft%7D%5C%2C+dt%5C%5C+%26%3D%5Cfrac%7B1%7D%7B%5Csigma+%5Csqrt%7B2+%5Cpi+%7D+%7D+%5Cint_%7B-%5Cinfty+%7D%5E%7B%5Cinfty+%7D+e%5E%7B-+%5Cfrac%7B1%7D%7B2+%5Csigma%5E2%7D%5Cleft%5Bt%5E2%2Bj4+%5Cpi+%5Csigma%5E2+ft+%5Cright%5D%7D%5C%2C+dt%5C%5C+%26%3D%5Cfrac%7B1%7D%7B%5Csigma+%5Csqrt%7B2+%5Cpi+%7D+%7D+%5Cint_%7B-%5Cinfty+%7D%5E%7B%5Cinfty+%7D+e%5E%7B-+%5Cfrac%7B1%7D%7B2+%5Csigma%5E2%7D%5Cleft%5Bt%5E2%2Bj4+%5Cpi+%5Csigma%5E2+ft+%2B+%28j+2+%5Cpi+%5Csigma%5E2+f%29%5E2+-+%28j+2+%5Cpi+%5Csigma%5E2+f%29%5E2%5Cright%5D%7D%5C%2C+dt%5C%5C+%26%3De%5E%7B+%5Cfrac%7B1%7D%7B2+%5Csigma%5E2%7D%28j+2+%5Cpi+%5Csigma%5E2+f%29%5E2%7D%5Cfrac%7B1%7D%7B%5Csigma+%5Csqrt%7B2+%5Cpi+%7D+%7D+%5Cint_%7B-%5Cinfty+%7D%5E%7B%5Cinfty+%7D+e%5E%7B-+%5Cfrac%7B1%7D%7B2+%5Csigma%5E2%7D%5Cleft%5Bt%2Bj+2+%5Cpi+%5Csigma%5E2+f+%5Cright%5D%5E2%7D%5C%2C+dt%5C%5C+%26%3De%5E%7B+%5Cfrac%7B1%7D%7B2+%5Csigma%5E2%7D%28j+2+%5Cpi+%5Csigma%5E2+f%29%5E2%7D%3De%5E%7B+-+%5Cfrac%7B1%7D%7B2%7D%28+2+%5Cpi+%5Csigma+f%29%5E2%7D%5Cend%7Baligned%7D+&bg=ffffff&fg=000&s=2&c=20201002)
![\displaystyle{\frac{1}{\sigma \sqrt{2 \pi} }\int_{-\infty }^{\infty }e^{- \frac{1}{2 \sigma^2}\left[t+j 2 \pi \sigma^2 f \right]^2}\, dt =\frac{1}{\sigma \sqrt{2 \pi } }\int_{-\infty }^{\infty }e^{- \frac{1}{2 \sigma^2} u^2}\, du =1}](https://s0.wp.com/latex.php?latex=%5Cdisplaystyle%7B%5Cfrac%7B1%7D%7B%5Csigma+%5Csqrt%7B2+%5Cpi%7D+%7D%5Cint_%7B-%5Cinfty+%7D%5E%7B%5Cinfty+%7De%5E%7B-+%5Cfrac%7B1%7D%7B2+%5Csigma%5E2%7D%5Cleft%5Bt%2Bj+2+%5Cpi+%5Csigma%5E2+f+%5Cright%5D%5E2%7D%5C%2C+dt+%3D%5Cfrac%7B1%7D%7B%5Csigma+%5Csqrt%7B2+%5Cpi+%7D+%7D%5Cint_%7B-%5Cinfty+%7D%5E%7B%5Cinfty+%7De%5E%7B-+%5Cfrac%7B1%7D%7B2+%5Csigma%5E2%7D+u%5E2%7D%5C%2C+du+%3D1%7D+&bg=ffffff&fg=000&s=2&c=20201002)
