ground0state · January 23, 2020 15:07 · Jan 23, 2020
diff --git a/power-spectral-density.py b/power-spectral-density.py
@@ -0,0 +1,61 @@
+# https://org-technology.com/posts/power-spectral-density.html
+"""
+適当な信号を作成して 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)
No results found