はじめに
このチュートリアルでは、Python の Matplotlib ライブラリを使ってパワー スペクトル密度 (PSD) をプロットする方法を説明します。PSD は信号処理の分野で一般的に使われるプロットです。NumPy には PSD を計算するための便利なライブラリがたくさんあり、Matplotlib を使ってこれをどのように実行し、可視化するかのいくつかの例を示します。
VM のヒント
VM の起動が完了したら、左上隅をクリックして ノートブック タブに切り替え、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 をプロットする
次に、異なる量のゼロ パディングで 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) をプロットする方法を紹介しました。まず、基本的な PSD をプロットし、その後、同等の MATLAB コードと比較しました。次に、異なる量のゼロ パディングで PSD をプロットし、その後、複素信号の PSD をプロットしました。