自己相関関数を使用した基本周波数の測定

昨日の日記で、FFTのピークは、周波数分解能によって決まる基底関数の最小周波数の整数倍とのずれによって値が変わることを説明した。

ここでは、自己相関関数を使用すると周波数分解能によらず安定して基本周波数を測定できることを示す。

まず、FFTで倍音を含む音声データでピークを測定した場合に、倍音が誤検知されることがある例を示す。

基本周波数を440Hzとして、振幅が基本周波数の95%の2倍音と、90%の3倍音を含む音声データをFFTすると以下のようになる。

import numpy as np
import matplotlib.pyplot as plt
fs = 44100
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) + np.sin(2 * np.pi * 2 * f * t) * 0.95 + np.sin(2 * np.pi * 3 * f * t) * 0.9
Y = np.fft.fft(y)
plt.plot(frq, abs(Y))
plt.axis([0, 2000, 0, 2500])
plt.show()

440Hzでは基本周波数がピークとなっている。

次に、基本周波数を435Hzとした場合は以下のようになる。

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

435Hzでは2倍音がピークになっている。次に高いピークは3倍音で、基本周波数のピークが一番低い値になっている。
これが、昨日の日記で説明した周波数分解能に起因する問題である。
一番高いピークではなく、一番の左側のピークを検出すればよさそうであるが、実際の音声データでは、ノイズなどで基本周波数とは無関係の位置に山の形が現れることがあり、どれをピークと見做すかは試行錯誤が必要であり、実際にやってみると精度を上げるのは難しいことがわかる。

では、同じ倍音を含む435Hzの音声データを自己相関関数を使用して基本周波数を測定するとどうなるだろうか。

自己相関関数は、2つの同じ音声データを1サンプルずつずらしながら全ての位置の値どうしを掛け算した和を計算する関数である。
周期性のあるデータの場合は、自己相関関数も周期的な関数となる。
パワースペクトルを逆フーリエ変換することで高速に計算できる。(ウィーナー＝ヒンチンの定理 - Wikipedia)

f = 435
y = np.sin(2 * np.pi * f * t) + np.sin(2 * np.pi * 2 * f * t) * 0.95 + np.sin(2 * np.pi * 3 * f * t) * 0.9
Y = np.fft.fft(y)
acf = np.fft.ifft(abs(Y)**2)
plt.plot(acf[0:size/2])
plt.axis([0, size/2, min(acf), max(acf)])
plt.show()

グラフの横軸は、ずらしたサンプル数であり、横軸の値をnとすると、fs(サンプリング周波数) / n が、周波数となる。

音声データには、周期的なデータが複数回現れるので、自己相関関数にもピークが複数回現れる。
グラフが右側に行くほど減衰する。（減衰する理由は、こちらの日記を参照）
よって、1番目のピークが基本周波数となる。(0番の値を除く)

ピークは、以下のようなピークピッキングによって検出する。

#ピークピッキング
def pick_peak(data):
	peaks_val = []
	peaks_index = []
	for i in range(2, data.size):
		if data[i-1] - data[i-2] >= 0 and data[i] - data[i-1] < 0:
			peaks_val.append(data[i-1])
			peaks_index.append(i-1)
	max_index = peaks_val.index(max(peaks_val))
	return peaks_index[max_index]

先ほどの435Hzのデータの自己相関関数のピークは以下のようになる。

n = pick_peak(acf[0:size/2])
f0 = fs/n

結果は、n = 101、f0 = 436.6 となり、基本周波数435Hzに近い値が求まっている。

FFTのピーク検出では倍音を誤検知することがあるが、自己相関関数のピークでは安定して基本周波数を求めることができる。
その理由は、自己相関関数のピークを計算する式に、FFTのフレームサイズ(size)が現れないためである。

なお、自己相関関数のピークが必ず基本周波数になるとは限らず、元から倍音が強い音声データでは、倍音が検出されてしまうため、基本周波数を検出するには工夫が必要になる。
また、自己相関関数で求まる周波数は、横軸nの逆数であるため高音になるほど(原点に近いほど)分解能が低くなる。
上記の例でも正しい基本周波数が435Hzのところ検出された周波数は436.6Hzであり、1.6Hzの誤差がある。
より正確な周波数を求めたい場合には工夫が必要となる。