Matplotlibでパワースペクトル密度（PSD）をプロットする方法

Matplotlibを使ってパワースペクトル密度（Power Spectral Density: PSD）をプロットするには、以下の手順に従います。

手順

• 図のサイズを設定し、サブプロット間および周囲の余白（パディング）を調整します。
• 時間刻み幅を表す変数 dt を初期化します。
• NumPyを使用して、t（時間）、nse（正規乱数ノイズ）、r（指数関数的に減衰する応答）、cnse（畳み込み後のノイズ）、s（信号）などのデータ点を作成します。
• 図と一連のサブプロットを作成します。
• plot() メソッドを使って、t と s のデータをプロットします。
• psd() メソッドを使ってパワースペクトル密度をプロットします。
• show() メソッドで図を表示します。

コード例

import matplotlib.pyplot as plt
import numpy as np

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

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

コードの解説

このコードでは、まず dt = 0.01 としてサンプリング間隔を定義し、np.arange() で0から10秒までの時間軸を作成しています。次に np.random.randn() でホワイトノイズを生成し、それを指数関数的に減衰するインパルス応答と畳み込むことで、色付きノイズ（相関のあるノイズ）を作り出しています。

最終的な信号 s は、振幅0.1の正弦波（1Hz）とこのノイズ成分を合成したものです。上段のサブプロットには時間領域の波形を表示し、下段のサブプロットでは psd(s, 512, 1 / dt) を呼び出して周波数領域におけるパワースペクトル密度を可視化します。第2引数の512はNFFT（FFTポイント数）、第3引数の 1 / dt はサンプリング周波数を表します。

実行結果

実行すると、上段に時間領域の信号波形が、下段に対応するパワースペクトル密度が表示されます。PSDのグラフからは、1Hz付近に正弦波成分によるピークが現れ、低周波域ほどノイズ成分のパワーが大きくなる傾向が確認できます。