FFTと窓関数に関する考察 - TadaoYamaokaの開発日記

前回の日記で、ピークの周波数の周辺のパワースペクトルの和がエネルギー(全周波数のパワースペクトルの総和)に占める割合について書いたろころ、Chachayさんからスペクトル漏れについてご指摘を頂いた。

chachay.hatenablog.com

FFTについてのはじめの日記を書いたときに、窓関数について理解が不十分だったため、窓関数ありなしでも周波数分解能との整数倍とのずれによってFFTの値が影響を受けるという結論については同じ傾向であったため、以降の調査で窓関数なしのまま調べてしまった。

ここでは、窓関数を適用した場合に、解析対象の周波数が周波数分解能(サンプリング周波数をFFTのフレームサイズで割った値)の整数倍からずれていた場合にも、FFTのピークからの裾の広がりが抑えられることについて考察する。

周波数分解能の整数倍に近い440Hzと、整数倍とのずれが大きい435Hzの正弦波について、窓関数ありなしの場合の裾の広がり方は、以下のようになる。

440Hz

han = np.hanning(size)

f = 440
y = np.sin(2 * np.pi * f * t)
Y = np.fft.fft(y)
print(np.var(abs(Y[0:size/2])))
plt.plot(frq, abs(Y))
Y_han = np.fft.fft(y * han)
print(np.var(abs(Y_han[0:size/2])))
plt.plot(frq, abs(Y_han))
plt.axis([200, 700, 0, 2100])
plt.show()

f:id:TadaoYamaoka:20160903120019p:plain

435Hz

f = 435
y = np.sin(2 * np.pi * f * t)
Y = np.fft.fft(y)
print(np.var(abs(Y[0:size/2])))
plt.plot(frq, abs(Y))
Y_han = np.fft.fft(y * han)
print(np.var(abs(Y_han[0:size/2])))
plt.plot(frq, abs(Y_han))
plt.axis([200, 700, 0, 2100])
plt.show()

f:id:TadaoYamaoka:20160903120050p:plain

グラフの青が窓関数なし、緑が窓関数ありの場合である。窓関数はハニング窓とした。

窓関数なしの場合は、440Hzに比べて435Hzの方がピークから裾が広がっている。
対して、窓関数ありの場合は、裾の広がりがほぼ同じになっている。
また、ピークの高さが440Hzに比べて435Hzが低くなっているが、グラフを補完すると同じ形の曲線を描いているように見える。

スペクトルの分散は以下の通りであった。

	窓関数なし	窓関数あり
440Hz	2043.14	766.75
435Hz	2032.43	766.69

窓関数なしに比べて窓関数ありの場合、分散の差が抑えられている。（単純比較はできないですが、統計知識が不足しているので感覚的な解釈です。）

つまり、窓関数を使用することで、解析対象の周波数と周波数分解能の整数倍とのずれによらず、スペクトルは近い形になるといえる。
ただし、ピークの値は、周波数分解能の整数倍とずれている場合は、実際のスペクトルの曲線のピークの位置にはならないので注意が必要となる。

このピークのずれは、前回の日記で書いた周辺の周波数のパワースペクトルの和を取ることでほぼ抑えることができる。
前回の日記の検証を窓関数ありで行うと、以下の通りとなる。

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

f:id:TadaoYamaoka:20160903122050p:plain
※青：元のパワースペクトル、緑：周辺3つの和、赤：周辺5つの和