Juliaで基音周波数の推定(続き)

先日、Juliaを使って基音周波数を推定する記事を書きました。

marui.hatenablog.com

その中で、基音よりも倍音成分のほうが振幅が大きいときのために自己相関関数を使用して対処するという説明をしましたが、それではうまくいかない音があったので、別の方法を使うことにしたよというのが今日の記事です。

手元にあったエレキベースの4弦の解放音(E1の音、約41.20 Hz)です。時間波形を見ると30秒近くまでサスティンが伸びているのがわかります。

f:id:amarui:20210924125255p:plain

これを先日と同じようにFFTをかけて周波数スペクトルを見てみます。

f:id:amarui:20210924125325p:plain

すると、第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)

f:id:amarui:20210924125607p:plain

相関値が最も高いのは基音ではなく第2倍音となってしまいました。ここでfindmax()でピーク周波数を求めると83 Hzとなり、1オクターブずれた結果になってしまいます。

自己相関関数だとうまくいかない場合もあるということで、やりかたを変えてみます。ネットで少し調べたところHarmonic Product Spectrumというテクニックがあることがわかりました。下の図はリンク先からお借りしたものですが、スペクトルを周波数方向に1/2、1/3……と縮めながら積をとることで基音周波数の成分だけが突出するようになるという概念図になっています。

f:id:amarui:20210924130444p:plain

これをそのまま実装したのが次のコードです。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が以下のグラフです。 f:id:amarui:20210924131511p:plain

一か所だけが突出しているのがわかります。ピーク部分付近だけを拡大すると…… f:id:amarui:20210924131540p:plain

ここから基音周波数は41.52 Hzだと推定されました。正確なE1からは13~14セントくらいずれていますので、チューニングをしっかりやらないとですね。