FFTのピークでは基本周波数を測れないわけ

囲碁ばかりだったので久しぶり音ネタについて書きます。

AndroidとiOSでリリースしている「Vocal Pitch Monitor」では、基本周波数(ピッチ)の測定に自己相関関数を使用している。

基本周波数はFFTのピークで測れるのではないかと思われがちだが、FFTのピークでは正確に測れない。
理由はいくつかあるが、その一つである周波数分解能とピークの関係について解説する。

まず、単純な場合として、音声データが正弦波(sin波)の場合を考える。

440Hzの正弦波をFFTすると以下のようなグラフになる。

import numpy as np
import matplotlib.pyplot as plt
fs = 44100.0
size = 4096
t = np.arange(0, size) / fs
frq = np.fft.fftfreq(size, 1/fs)

f = 440
y = np.sin(2 * np.pi * f * t)
Y = np.fft.fft(y)
plt.plot(frq, abs(Y))
plt.axis([0, 2000, 0, 2500])
plt.show()

f:id:TadaoYamaoka:20160827115044p:plain

別の周波数435HzをFFTすると以下のようなグラフになる。

f = 435
y = np.sin(2 * np.pi * f * t)
Y = np.fft.fft(y)
plt.plot(frq, abs(Y))
plt.axis([0, 2000, 0, 2500])
plt.show()

f:id:TadaoYamaoka:20160827115408p:plain

それぞれのFFTピークの値は、以下の通りである。

440Hz	1999.5
435Hz	1527.2

同じ正弦波でも、1.3倍ほどの差がある。
また、グラフの形も435Hzの方がピークが低くなった分、裾が広がった形になっている。

FFTの値は、基底関数に掛ける係数であるため、ピークの基底関数が実際のピークとずれている場合、ピークからずれた基底関数を足し合わせて実際のピークの関数を表現することになる。
ずれが大きいほど多くの関数の足し合わせが必要になり、係数の値がピークから裾広がりになる。

なお、面積(係数の値の合計値)が同じなのかというと、積分値はそれぞれ以下の通りであり一致しない。

440Hz	5728.1
435Hz	8966.5

音声データは、基本周波数とその倍音から構成されている。
FFTのピークの値は基本周波数と倍音がどの周波数かによって変わってくるため、実際のピークの周波数でFFTの値もピークになるとは限らない。
試してみると頻繁に倍音を誤検知してしまうことがわかる。
また、特定の音階で誤検知しやすいという現象があらわれる。

では、周波数によって、FFTのピークがどのように変化するだろうか。

正弦波の周波数を0.1Hzずつ変化させたときのFFTピークをグラフにすると以下のようになる。

peaks = []
fset = np.linspace(400, 500, 1001)
for f in fset:
	y = np.sin(2 * np.pi * f * t)
	Y = np.fft.fft(y)
	peaks.append(max(abs(Y)))
plt.plot(fset, peaks)
plt.show()

f:id:TadaoYamaoka:20160827121852p:plain

約10Hz置きに周期的にFFTのピークの値が変化していることがわかる。
最小値と最大値では、1.57倍の差がある。

最小値	1303.9
最大値	2048.0	最小値の1.57倍

このFFTのピークの周期は、周波数分解能と関係している。
上の方で、基底関数と実際のピークがずれている場合にピークの値が小さくなることを説明したが、FFTの基底関数(sinとcos)は最小周波数(サンプリング周波数をFFTのフレームサイズで割った値)の関数の整数倍の周波数の関数であるため、ピークとのずれもFFTの最小周波数と同じ周期となる。
FFTの基底関数は最小周波数は、周波数分解能と呼ばれ、上記例ではサンプリング周波数44.1kHz、フレームサイズ4096であるため、44.1kHz / 4096 = 10.76Hzである。

最小周波数の整数倍とのずれを先ほどのグラフに重ねると以下のようになる。

fig, ax1 = plt.subplots()
ax1.plot(fset, peaks)
ax2 = ax1.twinx()
ax2.plot(fset, fset/(fs/size) - (fset/(fs/size)).astype(np.int), color='r', linewidth=0.2)
plt.show()

f:id:TadaoYamaoka:20160827122754p:plain