先日、Juliaを使って基音周波数（≒ピッチ）を推定する記事を書きました。
その中で、基音よりも倍音成分のほうが振幅が大きいときのために自己相関関数を使用して対処するという説明をしましたが、それでもうまくいかない音があったので、別の方法を使うことにしたよというのが今日の記事です。
手元にあったエレキベースの4弦の解放音（E1の音、約41.20 Hz）です。時間波形を見ると30秒近くまでサスティンが伸びているのがわかります。
これを先日と同じようにFFTをかけて周波数スペクトルを見てみます。
すると、第2倍音・第4倍音が大きく出ているのがわかります。以下のように前回と同じコードを使って自己相関をとって20 Hz以下を無視したグラフを描いたところ……
yc = xcorr(y, y); yc = yc[length(yc)÷2+1 : end]; yc = yc ./ maximum(yc); lowfreq = 20; yc2 = yc[f .>= lowfreq]; fc2 = f[f .>= lowfreq]; plot(fc2[1:50000], yc2[1:50000], xlabel="Frequency (Hz)", ylabel="Correlation", legend=false)
相関値が最も高いのは基音ではなく第2倍音となってしまいました。ここでfindmax()でピーク周波数を求めると83 Hzとなり、1オクターブずれた結果になってしまいます。
自己相関関数だとうまくいかない場合もあるということで、やりかたを変えてみます。ネットで少し調べたところHarmonic Product Spectrumというテクニックがあることがわかりました。（最初の報告は A. Michael Noll. “Pitch Determination of Human Speech by the Harmonic Product Spectrum, the Harmonic Sum Spectrum and a Maximum Likelihood Estimate,” In Proceedings of the Symposium on Computer Processing in Communications, Vol. XIX, Polytechnic Press: Brooklyn, New York, (1970), pp. 779–797. らしいです。シンポジウムじたいは1969年に開催されたようです。）
下の図はリンク先からお借りしたものですが、スペクトルを周波数方向に1/2、1/3……と縮めながら積をとることで基音周波数の成分だけが突出するようになるという概念図になっています。
これをそのまま実装したのが次のコードです。DSP.resample()関数を使ってスペクトルを徐々に縮めて、ycに周波数スペクトルの積を入れていきます。第5倍音くらいまででも十分なようですが、以下では第10倍音までの積を取るようにしています。
yc = copy(y); for n = 2:10 yr = resample(y, 1/n); yr[isnan.(yr)] .= 0; yc[1:length(yr)] .*= yr; end
これで計算したHarmonic Product Spectrumが以下のグラフです。
一か所だけが突出しているのがわかります。ピーク部分付近だけを拡大すると……
ここから基音周波数は41.52 Hzだと推定されました。周波数から音名を計算したところ、正確なE1からは13～14セントくらいずれていますので、チューニングをしっかりやらないとですね。
