Trazado de la Densidad Espectral de Potencia con Matplotlib

Introducción

Este tutorial lo guiará a través del proceso de trazar la Densidad Espectral de Potencia (PSD, por sus siglas en inglés) utilizando la biblioteca Matplotlib en Python. La PSD es una representación comúnmente utilizada en el campo del procesamiento de señales. NumPy tiene muchas bibliotecas útiles para calcular una PSD, y mostraremos algunos ejemplos de cómo se puede realizar y visualizar con Matplotlib.

Consejos sobre la VM

Una vez que se haya completado la inicialización de la VM, haga clic en la esquina superior izquierda para cambiar a la pestaña Cuaderno y acceder a Jupyter Notebook para practicar.

A veces, es posible que tenga que esperar unos segundos a que Jupyter Notebook termine de cargarse. La validación de las operaciones no se puede automatizar debido a las limitaciones de Jupyter Notebook.

Si tiene problemas durante el aprendizaje, no dude en preguntar a Labby. Deje sus comentarios después de la sesión y lo resolveremos rápidamente para usted.

Trazar una PSD básica

Primero, trazaremos una PSD básica utilizando datos aleatorios. Crearemos una serie temporal, agregaremos ruido y luego trazaremos la PSD utilizando la función psd de la biblioteca matplotlib.mlab.

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab

np.random.seed(19680801)
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

fig, (ax0, ax1) = plt.subplots(2, 1)
ax0.plot(t, s)
ax1.psd(s, 512, 1 / dt)

plt.show()

Comparación con el código equivalente en MATLAB

Podemos comparar el código anterior con el código equivalente en MATLAB para lograr lo mismo:

dt = 0.01;
t = [0:dt:10];
nse = randn(size(t));
r = exp(-t/0.05);
cnse = conv(nse, r)*dt;
cnse = cnse(1:length(t));
s = 0.1*sin(2*pi*t) + cnse;

subplot(211)
plot(t, s)
subplot(212)
psd(s, 512, 1/dt)

Trazar la PSD con diferentes cantidades de relleno con ceros

A continuación, trazaremos la PSD con diferentes cantidades de relleno con ceros. Esto utiliza la serie temporal completa de una vez.

dt = np.pi / 100.
fs = 1. / dt
t = np.arange(0, 8, dt)
y = 10. * np.sin(2 * np.pi * 4 * t) + 5. * np.sin(2 * np.pi * 4.25 * t)
y = y + np.random.randn(*t.shape)

fig, axs = plt.subplot_mosaic([
    ['signal','signal','signal'],
    ['zero padding', 'block size', 'overlap'],
], layout='constrained')

axs['signal'].plot(t, y)
axs['signal'].set_xlabel('tiempo [s]')
axs['signal'].set_ylabel('señal')

axs['zero padding'].psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)
axs['zero padding'].psd(y, NFFT=len(t), pad_to=len(t) * 2, Fs=fs)
axs['zero padding'].psd(y, NFFT=len(t), pad_to=len(t) * 4, Fs=fs)

axs['block size'].psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)
axs['block size'].psd(y, NFFT=len(t) // 2, pad_to=len(t), Fs=fs)
axs['block size'].psd(y, NFFT=len(t) // 4, pad_to=len(t), Fs=fs)
axs['block size'].set_ylabel('')

axs['overlap'].psd(y, NFFT=len(t) // 2, pad_to=len(t), noverlap=0, Fs=fs)
axs['overlap'].psd(y, NFFT=len(t) // 2, pad_to=len(t),
                   noverlap=int(0.025 * len(t)), Fs=fs)
axs['overlap'].psd(y, NFFT=len(t) // 2, pad_to=len(t),
                   noverlap=int(0.1 * len(t)), Fs=fs)
axs['overlap'].set_ylabel('')
axs['overlap'].set_title('superposición')

for title, ax in axs.items():
    if title =='signal':
        continue

    ax.set_title(title)
    ax.sharex(axs['zero padding'])
    ax.sharey(axs['zero padding'])

plt.show()

Trazar la PSD de una señal compleja

Finalmente, trazaremos la PSD de una señal compleja para demostrar que las PSD complejas funcionan correctamente.

prng = np.random.RandomState(19680801)
fs = 1000
t = np.linspace(0, 0.3, 301)
A = np.array([2, 8]).reshape(-1, 1)
f = np.array([150, 140]).reshape(-1, 1)
xn = (A * np.exp(2j * np.pi * f * t)).sum(axis=0) + 5 * prng.randn(*t.shape)

fig, (ax0, ax1) = plt.subplots(ncols=2, layout='constrained')

yticks = np.arange(-50, 30, 10)
yrange = (yticks[0], yticks[-1])
xticks = np.arange(-500, 550, 200)

ax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,
        scale_by_freq=True)
ax0.set_title('Periodograma')
ax0.set_yticks(yticks)
ax0.set_xticks(xticks)
ax0.grid(True)
ax0.set_ylim(yrange)

ax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75,
        scale_by_freq=True)
ax1.set_title('Welch')
ax1.set_xticks(xticks)
ax1.set_yticks(yticks)
ax1.set_ylabel('')
ax1.grid(True)
ax1.set_ylim(yrange)

plt.show()

Resumen

En este tutorial, te hemos mostrado cómo trazar la Densidad Espectral de Potencia (PSD, por sus siglas en inglés) utilizando la biblioteca Matplotlib en Python. Comenzamos trazando una PSD básica, luego la comparamos con el código equivalente en MATLAB. A continuación, trazamos la PSD con diferentes cantidades de relleno con ceros, y después trazamos la PSD de una señal compleja.