Plotting power spectral density in Matplotlib

A Power Spectral Density (PSD) plot shows how the power of a signal is distributed across different frequencies. Matplotlib provides the psd() method to create these plots, which is useful for analyzing signal frequency content and noise characteristics.

Basic PSD Plot

Let's create a signal with noise and plot its power spectral density ?

import matplotlib.pyplot as plt
import numpy as np

# Set figure parameters
plt.rcParams["figure.figsize"] = [7.50, 3.50]
plt.rcParams["figure.autolayout"] = True

# Create time series data
dt = 0.01
t = np.arange(0, 10, dt)
nse = np.random.randn(len(t))
r = np.exp(-t / 0.05)
cnse = np.convolve(nse, r) * dt
cnse = cnse[:len(t)]
s = 0.1 * np.sin(2 * np.pi * t) + cnse

# Create subplots
fig, (ax0, ax1) = plt.subplots(2, 1)

# Plot original signal
ax0.plot(t, s)
ax0.set_title('Original Signal')
ax0.set_xlabel('Time (s)')
ax0.set_ylabel('Amplitude')

# Plot power spectral density
ax1.psd(s, 512, 1 / dt)
ax1.set_title('Power Spectral Density')

plt.tight_layout()
plt.show()

Understanding the PSD Parameters

The psd() method accepts several important parameters ?

import matplotlib.pyplot as plt
import numpy as np

# Generate a cleaner example signal
fs = 100  # Sampling frequency
t = np.linspace(0, 2, 2*fs)
signal = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*20*t) + 0.2*np.random.randn(len(t))

plt.figure(figsize=(10, 6))

# Plot with different NFFT values
plt.subplot(2, 2, 1)
plt.psd(signal, NFFT=64, Fs=fs)
plt.title('NFFT = 64')

plt.subplot(2, 2, 2)
plt.psd(signal, NFFT=256, Fs=fs)
plt.title('NFFT = 256')

plt.subplot(2, 2, 3)
plt.psd(signal, NFFT=512, Fs=fs)
plt.title('NFFT = 512')

plt.subplot(2, 2, 4)
plt.plot(t[:200], signal[:200])
plt.title('Original Signal (first 2 seconds)')
plt.xlabel('Time (s)')

plt.tight_layout()
plt.show()

Key Parameters

Practical Example with Multiple Frequencies

Here's a more detailed example showing how PSD helps identify frequency components ?

import matplotlib.pyplot as plt
import numpy as np

# Create a composite signal with known frequencies
fs = 1000  # Sampling frequency
t = np.linspace(0, 1, fs)

# Signal with 50Hz, 120Hz components and noise
freq1, freq2 = 50, 120
signal = (np.sin(2*np.pi*freq1*t) + 
          0.5*np.sin(2*np.pi*freq2*t) + 
          0.3*np.random.randn(len(t)))

plt.figure(figsize=(12, 8))

# Time domain plot
plt.subplot(3, 1, 1)
plt.plot(t[:200], signal[:200])
plt.title('Time Domain Signal (first 200 samples)')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

# PSD plot
plt.subplot(3, 1, 2)
plt.psd(signal, NFFT=1024, Fs=fs)
plt.title('Power Spectral Density')
plt.grid(True)

# Zoomed PSD
plt.subplot(3, 1, 3)
plt.psd(signal, NFFT=1024, Fs=fs)
plt.xlim(0, 200)
plt.title('PSD (Zoomed: 0-200 Hz)')
plt.grid(True)

plt.tight_layout()
plt.show()

Conclusion

Power Spectral Density plots are essential for frequency domain analysis of signals. Use higher NFFT values for better frequency resolution, and adjust the sampling frequency parameter to get correct frequency scales in your plots.