Python Nitime: マルチテーパ法による PSD 推定

今回はマルチテーパ法(MTM)による PSD 推定についてメモ。

MTM は Matalb では関数 pmtm で使用することができます。 残念ながら SciPy には MTM の関数は含まれていないため、Nitime パッケージを使用します。 Nitime は pip でインストールでき、Python3 でも使用できます。

MTM の利点、欠点について簡単にまとめます。 周波数領域解析では、どのような窓処理を行なってもデータの境界付近の情報が失なわれ、必ずスペクトル漏れがあります。 PSD 推定において Welch 法 は窓処理したセグメントをオーバラップさせて平均化していますが、これもスペクトル漏れは避けられません。 MTM ではスペクトル漏れの許容値に対してそれを満たす複数の窓関数でデータを窓処理し、スペクトルを加重平均します。 窓関数群は こちら に示されているようなものになります。 これにより MTM ではスペクトル漏れの影響を抑えることができます。 また、ノイズによるスペクトル漏れを抑えるので、スペクトルのばらつきを低減することにもなります。 欠点としては平均化のためスペクトルの周波数分解能が下がってしまうこと等があります。

# coding: utf-8
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import nitime.algorithms as tsa

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)  # Periodogram
freq2, P2 = signal.welch(y, fs, nperseg=n/2)  # Welch`s method
freq3, P_mt, nu = tsa.multi_taper_psd(y, Fs=fs)  # MTM

plt.figure()
plt.plot(freq1, 10*np.log10(P1), "b", label="Periodogram")
plt.plot(freq2, 10*np.log10(P2), "r", linewidth=2, label="Welch(nseg=n/2)")
plt.plot(freq3, 10*np.log10(P_mt), "y", linewidth=2, label="MTM")
plt.ylim(-60, 0)
plt.legend(loc="upper right")
plt.xlabel("Frequency[Hz]")
plt.ylabel("Power/frequency[dB/Hz]")
plt.show()
15063001.png

コメント

Comments powered by Disqus
書籍更新情報
2017-02-18
Pythonによる科学技術計算 基礎編
1.4版への更新が可能になりました。
サポートページはこちら
電子書籍