TadaoYamaokaの日記

山岡忠夫 Home で公開しているプログラムの開発ネタを中心に書いていきます。

FFTと窓関数に関する考察

前回の日記で、ピークの周波数の周辺のパワースペクトルの和がエネルギー(全周波数のパワースペクトルの総和)に占める割合について書いたろころ、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つの和

周辺3の和を取ることで、周波数分解能の整数倍とずれの影響はほとんど受けなくなっている。