How to perform spectral analysis in Python | DSP library
This page covers steps on how to perform spectral analysis in Python. It mentions python DSP library for performing spectral analysis on time domain signals.
Introduction: The parameters of any signal can be measured either in time domain or frequency domain. Signals can be in many shapes and sizes such as sinewave, square wave, rectangular wave, random binary sequence, pulse, random noise etc. Signals can be periodic or aperiodic based on their repetitions across time. Spectral analysis refers to measuring/plotting and analyzing signal spectrum in time as well as frequency domains.
Time domain analysis of signals usually can be done using Oscilloscope in which X-axis represents time and Y-axis represents or maps amplitude. Frequency domain analysis of signals can be done using Spectrum Analyzer in which X-axis represents frequency and Y-axis represents power or amplitude2.
Time domain and frequency domain analysis can be done in the computer using languages such as MATLAB, Labview, Python etc. IFFT (Inverse Fast Fourier Transform) function converts frequency domain signal vector into time domain vector. FFT (Fast Fourier Transform) or DFT (Discrete Fourier Transform) function converts time domain signal vector into frequency domain signal vector. Hence it is very easy to perform spectral analysis of signals in python or similar languages using built-in functions provided by the respective toolkits.
Spectral analysis python script
Following python code mentions time domain and frequency domain plots useful for spectral analysis of signals.
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fc=1000 # frequency of the carrier
N = 1e3
fs=10*fc # sampling frequency with oversampling factor=10
time = np.arange(N) / fs
x = np.sin(2*np.pi*fc*time) # time domain signal (real number)
plt.figure()
plt.plot(time, x)
plt.xlabel('Time')
plt.ylabel('Signal amplitude')
plt.title('Time domain signal amplitude versus time samples')
N1 = 256 # FFT size
X1 = fft(x, N1) # N-point complex DFT, output contains DC at index 0
# calculate frequency bins with FFT
# Power spectral density (PSD) using FFT and absolute squared
df=fs/N1 # frequency resolution
sampleIndex = np.arange(start = 0,stop = N1) # raw index for FFT plot
f=sampleIndex*df # x-axis index converted to frequencies
PSD = abs(X1)**2
n1 = len(f)//2
n2 = len(PSD)//2
plt.figure()
plt.plot(f[1:n1], PSD[1:n2])
plt.xlabel('Frequency indices')
plt.ylabel('Signal power')
plt.title('Frequency domain signal power vs frequency indices')
# Power spectral density (PSD) using signal.welch
f1, Pxx_spec = signal.welch(x, fs, 'flattop', 512, scaling='spectrum')
plt.figure()
plt.semilogy(f1, np.sqrt(Pxx_spec))
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.title('Power spectrum as per scipy.signal.welch function')
plt.show()
Spectrum analysis in python output plots
Following are the respective outputs at different points in above python code.
It mentions following plots.
1. Plot time domain waveform (Amplitude versus time)
2. Plot instantaneous power vs time waveform (Power vs time)
3. Plot frequency domain waveform (Power vs frequency)
• Using FFT
• Using signal.welch function
Amplitude vs time plot
Power vs frequency plot using FFT
Power vs frequency plot using Welch
There are other plots which provide more insights when signal contains complex baseband information with advanced modulation schemes such as OFDM, OFDMA, 16QAM, 64QAM, 256QAM, 1024QAM etc. These plots are EVM (Error Vector Magnitude) vs symbols, EVM vs subcarriers, Constellation IQ diagram, CCDF, channel frequency response, spectral flatness, jitter (measured in time domain), phase noise (measured in frequency domain) etc.