FFTとエネルギーについての考察(続き2)

前回の続きです。

ピークの周波数の周辺のパワースペクトルの和がエネルギー(全周波数のパワースペクトルの総和)に占める割合を、400Hzから500Hzの範囲で0.1Hzおきに調べてみた。

fset = np.linspace(400, 500, 1001)
range = [1, 3, 5]
for r in range:
	peaks = []
	for f in fset:
		y = np.sin(2 * np.pi * f * t)
		Y = np.fft.fft(y)
		P = abs(Y[0:size/2])**2
		P_around = np.convolve(P, np.ones(r), 'valid')
		peaks.append(max(P_around) / sum(P))
	plt.plot(fset, peaks)
plt.axis([400, 500, 0, 1.1])
plt.show()

青が元のパワースペクトル、緑が周辺3つの和、赤が周辺5つの和のときのピークの値である。

それぞれの最大値、最小値は以下の通り。

	最大値	最小値
元のパワースペクトル	0.99	0.40
周辺3つの和	0.99	0.85
周辺5つの和	0.99	0.91

周波数分解能の整数倍とのずれが最大のとき、元のパワースペクトルではエネルギーの40%しか測れないのに対して、周辺3つの和をとることで85%のエネルギーが測れている。周辺5つでは91%となっている。

なお、例として400Hzから500Hzの範囲を調べているが、他の周波数帯でもグラフが平行移動するだけで結果は同じとなる。

それにしても、python+numpy+matplotlibはこういう調査するには便利ですね。
Jupyter QtConsoleを使って、対話的にいろいろ調べることができます。
一時期OctaveとかR言語も試してみましたが、pythonの方が覚えやすいし使いやすい。
それ以前はC言語でFFTのプログラム書いて調べていましたけど、めんどくさくて途中で投げ出したくなってました。