Matplotlib PSD 플로팅

Beginner

This tutorial is from open-source community. Access the source code

소개

이 튜토리얼은 Python 의 Matplotlib 라이브러리를 사용하여 전력 스펙트럼 밀도 (PSD, Power Spectral Density) 를 플로팅하는 과정을 안내합니다. PSD 는 신호 처리 분야에서 일반적으로 사용되는 플롯입니다. NumPy 는 PSD 를 계산하는 데 유용한 많은 라이브러리를 가지고 있으며, Matplotlib 을 사용하여 이를 수행하고 시각화하는 몇 가지 예시를 보여드리겠습니다.

VM 팁

VM 시작이 완료되면, 왼쪽 상단을 클릭하여 Notebook 탭으로 전환하여 실습을 위해 Jupyter Notebook에 접속하십시오.

때로는 Jupyter Notebook 이 로딩을 완료하는 데 몇 초 정도 기다려야 할 수 있습니다. Jupyter Notebook 의 제한으로 인해 작업의 유효성 검사는 자동화될 수 없습니다.

학습 중 문제가 발생하면 언제든지 Labby 에게 문의하십시오. 세션 후 피드백을 제공해주시면 문제를 신속하게 해결해 드리겠습니다.

기본 PSD 플롯

먼저, 임의의 데이터를 사용하여 기본 PSD 를 플롯합니다. 시계열을 생성하고, 노이즈를 추가한 다음, matplotlib.mlab 라이브러리의 psd 함수를 사용하여 PSD 를 플롯합니다.

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()

동등한 MATLAB 코드와 비교

이전 코드를 동일한 작업을 수행하는 동등한 MATLAB 코드와 비교할 수 있습니다.

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)

서로 다른 양의 패딩으로 PSD 플롯

다음으로, 서로 다른 양의 제로 패딩 (zero padding) 으로 PSD 를 플롯합니다. 이 예제는 전체 시계열을 한 번에 사용합니다.

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('time [s]')
axs['signal'].set_ylabel('signal')

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('overlap')

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()

복소 신호의 PSD 플롯

마지막으로, 복소 PSD 가 제대로 작동하는지 보여주기 위해 복소 신호의 PSD 를 플롯합니다.

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('Periodogram')
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()

요약

이 튜토리얼에서는 Python 의 Matplotlib 라이브러리를 사용하여 전력 스펙트럼 밀도 (PSD, Power Spectral Density) 를 플롯하는 방법을 보여드렸습니다. 기본적인 PSD 를 플롯하는 것으로 시작하여, 이를 MATLAB 코드와 비교했습니다. 그런 다음 서로 다른 양의 제로 패딩 (zero padding) 으로 PSD 를 플롯하고, 복소 신호의 PSD 를 플롯했습니다.