SciPy でパワースペクトル密度(PSD) を推定する方法についてメモ。

関数について

SciPy では表 1 の関数で PSD を推定することができます。 SciPy にはノンパラメトリック法で PSD を推定する関数しかありません。個人的には計測データから PSD を求められれば十分なので困ったことはありません。

表1: PSD 推定の関数
関数	説明
periodogram	Periodogram 法(修正 Periodogram 法)
welch	Welch 法(平均 Periodogram 法)
lombscargle	Lomb-Scargle 法

periodogram は最も簡単な方法で、単純に測定データの離散時間フーリエ変換を求め、その結果の大きさを二乗するものです。計算は高速ですが、原理的に推定誤差が大きい方法です。オプション引数で FFT の前処理で使用する窓関数を指定できます。(修正 Periodogram 法)

測定データが欠けのない等間隔なサンプリングデータであれば、大抵は welch を使います。 Welch 法はデータを(何% か重ね合せるように)いくつかのセグメントに分割し、各セグメントについて修正 Periodogram 法で PSD を推定し、最後にその結果から平均を計算するものです。そのため、welch のオプション引数でセグメント長をデータ長と同じに設定し、同じ窓関数を使用すれば、periodogram の結果と同じになります。

lombscargle は欠損があったりして不等間隔でサンプリングされたデータから PSD を推定するのに使用します。

PSD の推定方法の種類と特徴についてはここがまとまっていると思いました。

Welch 法について最初はオーバーラップの概念がわかりにくですが、ここの説明がわかりやすかったです。そこで紹介されていた論文が窓関数について詳しく解説しています。窓関数を選ぶ際に大変参考になります。

[2015-06-30 火] マルチテーパ法 (MTM)による PSD 推定について記事を書きました。

例

適当な信号を作成して periodogram と welch で PSD を推定してみます。下の例ではセグメント長 nseg はデフォルトで 256 なので、ちょうど n/4 の長さです。いくつかセグメント長を変えて推定した結果をプロットしています。なおオーバーラップ noverlap はデフォルトのままで nseg/2、つまり 50% オーバーラップして計算しています。 periodogram では窓関数はデフォルトではなし(Boxcar と同じ)、welch では Hanning です。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

n = 1024
dt = 0.001
fs = 1/dt
f1 = 120
f2 = 150
t = np.linspace(1, n, n)*dt-dt
y = np.sin(2*np.pi*f1*t)+2*np.sin(2*np.pi*f2*t)+0.1*np.random.randn(t.size)

freq1, P1 = signal.periodogram(y, fs)
freq2, P2 = signal.welch(y, fs)
freq3, P3 = signal.welch(y, fs, nperseg=n/2)
freq4, P4 = signal.welch(y, fs, nperseg=n/8)

plt.figure()
plt.plot(freq1, 10*np.log10(P1), "b", label="periodogram")
plt.plot(freq2, 10*np.log10(P2), "r", linewidth=2, label="nseg=n/4")
plt.plot(freq3, 10*np.log10(P3), "c", linewidth=2, label="nseg=n/2")
plt.plot(freq4, 10*np.log10(P4), "y", linewidth=2, label="nseg=n/8")
plt.ylim(-60, 0)
plt.legend(loc="upper right")
plt.xlabel("Frequency[Hz]")
plt.ylabel("Power/frequency[dB/Hz]")
# plt.show()

オプション引数で様々な窓関数を指定できます。

freq5, P5 = signal.welch(y, fs, window="boxcar")
freq6, P6 = signal.welch(y, fs, window="flattop")

plt.figure()
plt.plot(freq2, 10*np.log10(P2), "r", linewidth=2, label="hanning")
plt.plot(freq5, 10*np.log10(P5), "g", linewidth=2, label="boxcar")
plt.plot(freq6, 10*np.log10(P6), "m", linewidth=2, label="flattop")
plt.ylim(-60, 0)
plt.legend(loc="upper right")
plt.xlabel("Frequency[Hz]")
plt.ylabel("Power/frequency[dB/Hz]")
# plt.show()

以下のようにして推定した PSD から実効値も簡単に求められます。

df = 1/dt/n
np.sqrt(np.sum(P2)*df)

0.79160410834650941